Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2020 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionRegister.h"
23 : #include "core/ActionShortcut.h"
24 : #include "multicolvar/MultiColvarShortcuts.h"
25 : #include "tools/IFile.h"
26 : #include "core/ActionSetup.h"
27 :
28 : //+PLUMEDOC MCOLVARF PAMM
29 : /*
30 : Probabilistic analysis of molecular motifs.
31 :
32 : Probabilistic analysis of molecular motifs (PAMM) was introduced in the papers in the bibliography.
33 : The essence of this approach involves calculating some large set of collective variables
34 : for a set of atoms in a short trajectory and fitting this data using a Gaussian Mixture Model.
35 : The idea is that modes in these distributions can be used to identify features such as hydrogen bonds or
36 : secondary structure types.
37 :
38 : The assumption within this implementation is that the fitting of the Gaussian mixture model has been
39 : done elsewhere by a separate code. You thus provide an input file to this action which contains the
40 : means, covariance matrices and weights for a set of Gaussian kernels, $\{\phi\}$. The values and
41 : derivatives for the following set of quantities is then computed:
42 :
43 : $$
44 : s_k = \frac{ \phi_k}{ \sum_i \phi_i }
45 : $$
46 :
47 : Each of the $\phi_k$ is a Gaussian function that acts on a set in quantities calculated that might be calculated
48 : using a [TORSION](TORSION.md), [DISTANCE](DISTANCE.md) or [ANGLE](ANGLE.md) action for example.
49 : These quantities are then inserted into the set of $n$ kernels that are in the input file. This will be done for multiple sets of values
50 : for the input quantities and a final quantity will be calculated by summing the above $s_k$ values or
51 : some transformation of the above. This sounds less complicated than it is and is best understood by
52 : looking through the example given below, which can be expanded to show the full set of operations that PLUMED is performing.
53 :
54 : \warning Mixing input variables that are periodic with variables that are not periodic has not been tested
55 :
56 : ## Examples
57 :
58 : In this example I will explain in detail what the following input is computing:
59 :
60 : ```plumed
61 : #SETTINGS MOLFILE=regtest/pamm/rt-pamm-periodic/M1d.pdb INPUTFILES=regtest/pamm/rt-pamm-periodic/2D-testc-0.75.pammp
62 : MOLINFO MOLTYPE=protein STRUCTURE=regtest/pamm/rt-pamm-periodic/M1d.pdb
63 : psi: TORSION ATOMS1=@psi-2 ATOMS2=@psi-3 ATOMS3=@psi-4
64 : phi: TORSION ATOMS1=@phi-2 ATOMS2=@phi-3 ATOMS3=@phi-4
65 : p: PAMM ...
66 : ARG=phi,psi MEAN
67 : CLUSTERS=regtest/pamm/rt-pamm-periodic/2D-testc-0.75.pammp
68 : ...
69 : PRINT ARG=p-1_mean,p-2_mean FILE=colvar
70 : ```
71 :
72 : The best place to start our explanation is to look at the contents of the `2D-testc-0.75.pammp` file, which you can do
73 : by clicking on the links in the annotated input above. This files contains the parameters of two two-dimensional Gaussian functions.
74 : Each of these Gaussian kernels has a weight, $w_k$, a vector that specifies the position of its center, $\mathbf{c}_k$, and a covariance matrix, $\Sigma_k$.
75 : The $\phi_k$ functions that we use to calculate our PAMM components are thus:
76 :
77 : $$
78 : \phi_k = \frac{w_k}{N_k} \exp\left( -(\mathbf{s} - \mathbf{c}_k)^T \Sigma^{-1}_k (\mathbf{s} - \mathbf{c}_k) \right)
79 : $$
80 :
81 : In the above $N_k$ is a normalization factor that is calculated based on $\Sigma$. The vector $\mathbf{s}$ is a vector of quantities
82 : that are calculated by the input [TORSION](TORSION.md) actions. This vector must be two dimensional and in this case each component is the value of a
83 : torsion angle. If we look at the two TORSION actions in the above we are calculating the $\phi$ and $\psi$ backbone torsional
84 : angles in a protein (Note the use of [MOLINFO](MOLINFO.md) to make specification of atoms straightforward). We thus calculate the values of our
85 : 2 $ \{ \phi \} $ kernels 3 times. The first time we use the $\phi$ and $\psi$ angles in the second residue of the protein,
86 : the second time it is the $\phi$ and $\psi$ angles of the third residue of the protein and the third time it is the $\phi$ and $\psi$ angles
87 : of the fourth residue in the protein. The final two quantities that are output by the print command, p.mean-1 and p.mean-2, are the averages
88 : over these three residues for the quantities:
89 :
90 : $$
91 : s_1 = \frac{ \phi_1}{ \phi_1 + \phi_2 }
92 : $$
93 :
94 : and
95 :
96 : $$
97 : s_2 = \frac{ \phi_2}{ \phi_1 + \phi_2 }
98 : $$
99 :
100 : There is a great deal of flexibility in this input. We can work with, and examine, any number of components, we can use any set of collective variables
101 : and compute these PAMM variables and we can transform the PAMM variables themselves in a large number of different ways when computing these sums. Furthermore,
102 : by expanding the shortcuts in the example above we can obtain insight into how the PAMM method operates.
103 : */
104 : //+ENDPLUMEDOC
105 :
106 : namespace PLMD {
107 : namespace pamm {
108 :
109 : class PAMM : public ActionShortcut {
110 : public:
111 : static void registerKeywords( Keywords& keys );
112 : explicit PAMM(const ActionOptions&);
113 : };
114 :
115 : PLUMED_REGISTER_ACTION(PAMM,"PAMM")
116 :
117 6 : void PAMM::registerKeywords( Keywords& keys ) {
118 6 : ActionShortcut::registerKeywords( keys );
119 6 : keys.add("compulsory","ARG","the vectors from which the pamm coordinates are calculated");
120 6 : keys.add("compulsory","CLUSTERS","the name of the file that contains the definitions of all the clusters");
121 6 : keys.add("compulsory","REGULARISE","0.001","don't allow the denominator to be smaller then this value");
122 6 : keys.add("compulsory","KERNELS","all","which kernels are we computing the PAMM values for");
123 6 : multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
124 6 : keys.needsAction("KERNEL");
125 6 : keys.needsAction("COMBINE");
126 6 : keys.addDOI("10.1063/1.4900655");
127 6 : keys.addDOI("10.1021/acs.jctc.7b00993");
128 6 : }
129 :
130 2 : PAMM::PAMM(const ActionOptions& ao) :
131 : Action(ao),
132 2 : ActionShortcut(ao) {
133 : // Must get list of input value names
134 2 : std::vector<std::string> valnames;
135 4 : parseVector("ARG",valnames);
136 : // Create input values
137 2 : std::string argstr=" ARG=" + valnames[0];
138 3 : for(unsigned j=1; j<valnames.size(); ++j) {
139 2 : argstr += "," + valnames[j];
140 : }
141 :
142 : // Create actions to calculate all pamm kernels
143 : unsigned nkernels = 0;
144 : double h;
145 : std::string fname;
146 2 : parse("CLUSTERS",fname);
147 2 : IFile ifile;
148 2 : ifile.open(fname);
149 2 : ifile.allowIgnoredFields();
150 : for(unsigned k=0;; ++k) {
151 22 : if( !ifile.scanField("height",h) ) {
152 : break;
153 : }
154 : // Create a kernel for this cluster
155 : std::string num, wstr, ktype;
156 9 : Tools::convert( k+1, num );
157 9 : Tools::convert(h,wstr);
158 9 : ifile.scanField("kerneltype",ktype);
159 54 : readInputLine( getShortcutLabel() + "_kernel-" + num + ": KERNEL NORMALIZED" + argstr + " NUMBER=" + num + " REFERENCE=" + fname + " WEIGHT=" + wstr + " TYPE=" + ktype );
160 9 : nkernels++;
161 9 : ifile.scanField();
162 9 : }
163 2 : ifile.close();
164 :
165 : // And add on the regularization
166 : std::string regparam;
167 4 : parse("REGULARISE",regparam);
168 : // Now combine all the PAMM objects with the regparam
169 2 : std::string paramstr, cinput = getShortcutLabel() + "_ksum: COMBINE PERIODIC=NO";
170 11 : for(unsigned k=0; k<nkernels; ++k) {
171 : std::string num;
172 9 : Tools::convert( k+1, num );
173 9 : if( k==0 ) {
174 : cinput += " ARG=";
175 4 : paramstr=" PARAMETERS=-" + regparam;
176 : } else {
177 : cinput += ",";
178 : paramstr += ",0";
179 : }
180 18 : cinput += getShortcutLabel() + "_kernel-" + num;
181 : }
182 4 : readInputLine( cinput + paramstr );
183 :
184 : // And now compute all the pamm kernels
185 : std::string kchoice;
186 4 : parse("KERNELS",kchoice);
187 : std::map<std::string,std::string> keymap;
188 2 : multicolvar::MultiColvarShortcuts::readShortcutKeywords( keymap, this );
189 2 : if( kchoice=="all" ) {
190 11 : for(unsigned k=0; k<nkernels; ++k) {
191 : std::string num;
192 9 : Tools::convert( k+1, num );
193 45 : readInputLine( getShortcutLabel() + "-" + num + ": CUSTOM ARG=" + getShortcutLabel() + "_kernel-" + num + "," + getShortcutLabel() + "_ksum FUNC=x/y PERIODIC=NO");
194 27 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel() + "-" + num, getShortcutLabel() + "-" + num, "", keymap, this );
195 : }
196 : } else {
197 0 : std::vector<std::string> awords=Tools::getWords(kchoice,"\t\n ,");
198 0 : Tools::interpretRanges( awords );
199 0 : for(unsigned k=0; k<awords.size(); ++k) {
200 0 : readInputLine( getShortcutLabel() + "-" + awords[k] + ": CUSTOM ARG=" + getShortcutLabel() + "_kernel-" + awords[k] + "," + getShortcutLabel() + "_ksum FUNC=x/y PERIODIC=NO");
201 0 : multicolvar::MultiColvarShortcuts::expandFunctions( getShortcutLabel() + "-" + awords[k], getShortcutLabel() + "-" + awords[k], "", keymap, this );
202 : }
203 0 : }
204 6 : }
205 :
206 : }
207 : }
|
