Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2018 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/ActionShortcut.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/ActionWithArguments.h"
26 : #include "tools/PDB.h"
27 : #include "Path.h"
28 :
29 : namespace PLMD {
30 : namespace mapping {
31 :
32 : //+PLUMEDOC COLVAR PCAVARS
33 : /*
34 : Projection on principal component eigenvectors or other high dimensional linear subspace
35 :
36 : As discussed in the documentation for [RMSD](RMSD.md) and [the refdist module](module_refdist.md) there are various different ways
37 : of calculating the distance between the instantaneous structure adopted by the system and some high-dimensional, reference configuration. The
38 : problem with all these methods, as one gets further and further from the reference configuration, the
39 : distance from it becomes a progressively poorer and poorer collective variable. This happens because
40 : the ``number" of structures at a distance $d$ from a reference configuration is proportional to $d^N$ in
41 : an $N$ dimensional space. Consequently, when $d$ is small the distance from the reference configuration
42 : may well be a good collective variable. However, when $d$ is large it is unlikely that the distance from the reference
43 : structure is a good CV. When the distance is large there will almost certainly be markedly different
44 : configurations that have the same CV value and hence barriers in transverse degrees of
45 : freedom.
46 :
47 : For these reasons dimensionality reduction is often employed so a projection $\mathbf{s}$ of a high-dimensional configuration
48 : $\mathbf{X}$ in a lower dimensionality space is found using a function:
49 :
50 : $$
51 : \mathbf{s} = F(\mathbf{X}-\mathbf{X}^{ref})
52 : $$
53 :
54 : where here we have introduced some high-dimensional reference configuration $\mathbf{X}^{ref}$. By far the simplest way to
55 : do this is to use some linear operator for $F$. That is to say we find a low-dimensional projection
56 : by rotating the basis vectors using some linear algebra:
57 :
58 : $$
59 : \mathbf{s}_i = \sum_k A_{ik} ( X_{k} - X_{k}^{ref} )
60 : $$
61 :
62 : Here $\mathbf{X}-\mathbf{X}^{ref}$ is the displacement from the refernece configuration in the high dimenisonal space, which you can
63 : calculate using the [RMSD](RMSD.md) action or by simply calculating difference between the instantaneous values for a collection of variables
64 : and the values of the variables for a particular reference configuration. $A$ is then a $d$ by $D$ matrix where $D$ is the dimensionality
65 : of the high dimensional space and $d$ is the dimensionality of the lower dimensional subspace. This matrix, $A$,
66 : can be found by various means including principal component analysis ([PCA](PCA.md)) and normal mode analysis. In both these methods the
67 : rows of $A$ would be the principle eigenvectors of a square matrix. For PCA the covariance while for normal modes the Hessian.
68 :
69 : ## Examples
70 :
71 : The following input calculates a projection on a linear subspace where the displacements
72 : from the reference configuration are calculated using [RMSD](RMSD.md) and TYPE=OPTIMAL. Consequently,
73 : both translation of the center of mass of the atoms and rotation of the reference
74 : frame are removed from the displacements that appear in the equatiosn above. The matrix $A$ and the reference
75 : configuration $R^{ref}$ are specified in the pdb input file reference.pdb and the
76 : value of all projections (and the residual) are output to a file called colvar2.
77 :
78 : ```plumed
79 : #SETTINGS INPUTFILES=regtest/mapping/rt-pca/reference.pdb
80 : pca2: PCAVARS REFERENCE=regtest/mapping/rt-pca/reference.pdb TYPE=OPTIMAL
81 : PRINT ARG=pca2.* FILE=colvar2
82 : ```
83 :
84 : As you can see from the example above, the reference configurations are specified using a pdb file. The first configuration that you provide is the reference configuration,
85 : which is referred to in the above expressions as $X^{ref}$. Subsequent configurations give the directions of row vectors that are contained in
86 : the matrix $A$ above. These directions are specified by giving a second configuration that describes the components of $A$ explicitly.
87 :
88 : When you use an input like the one in the example above PLUMED assumes that the atoms for which RMSD displacements are being computed together form a molecule.
89 : A procedure akin to that in [WHOLEMOLECULES](WHOLEMOLECULES.md) is thus performed to ensure that any bonds that are broken by the periodic boundary conditions are
90 : reformed. If you would like to turn this feature off for any reason you add the NOPBC flag to the input line as shown below:
91 :
92 : ```plumed
93 : #SETTINGS INPUTFILES=regtest/mapping/rt-pca/reference.pdb
94 : pca2: PCAVARS REFERENCE=regtest/mapping/rt-pca/reference.pdb NOPBC TYPE=OPTIMAL
95 : PRINT ARG=pca2.* FILE=colvar2
96 : ```
97 :
98 : Notice that the PCAVARS command is a shortcut. If you look at how the shortcut in the above input is expanded you should be able to see how the command works
99 : by calculating the RMSD displacement between the instantaneous and reference configuration and how those displacements are then projected on the eigenvector that
100 : was specified in the second frame of the pdb input above. Understanding the expanded version of this shortcut command allows you to recognise that you can project
101 : the displacement vector on any arbitrary vector. For example in the input below two reference structures are provided in the pdb file. PLUMED calculates the RMSD
102 : distance between these two reference configurations during setup and sets up a unit vector called eig that points along the director connecting the two RMSD structure.
103 : During the calculation the vector connecting the instantaneous configuration and the first of the two reference configurations in computed. This vector is then projected
104 : on the unit vector connecting the two initial structures:
105 :
106 : ```plumed
107 : #SETTINGS INPUTFILES=regtest/mapping/rt-pca-two-frames/two-frames.pdb
108 :
109 : # Read in two reference configuratioms from PDB file
110 : ref1: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pca-two-frames/two-frames.pdb NUMBER=1
111 : ref1T: TRANSPOSE ARG=ref1
112 : ref2: PDB2CONSTANT REFERENCE=regtest/mapping/rt-pca-two-frames/two-frames.pdb NUMBER=2
113 : # Calculate the displacement vector that takes you from ref1 to ref2
114 : eigu: RMSD ARG=ref1T,ref2 DISPLACEMENT TYPE=OPTIMAL
115 : # Normalise the reference vector
116 : eigu2: CUSTOM ARG=eigu.disp FUNC=x*x PERIODIC=NO
117 : eign2: SUM ARG=eigu2 PERIODIC=NO
118 : eig: CUSTOM ARG=eigu.disp,eign2 FUNC=x/sqrt(y) PERIODIC=NO
119 : eigT: TRANSPOSE ARG=eig
120 : # Everything prior to this point is only done in setup. From here onwards we have the commands that are done during the main calculate loop.
121 :
122 : # Calculate the RMSD displacement between the instantaneous structure and the first reference structure
123 : rmsd: RMSD REFERENCE=regtest/mapping/rt-pca-two-frames/two-frames.pdb NUMBER=1 TYPE=OPTIMAL DISPLACEMENT SQUARED
124 : # Project the displacement computed above on the director of the vector that connects reference structure ref1 to refeference structure ref2
125 : pca: MATRIX_VECTOR_PRODUCT ARG=eigT,rmsd.disp
126 :
127 : # Print the final CV to a file
128 : PRINT ARG=pca FILE=colvar
129 : ```
130 :
131 : You can also project vectors of differences of arguments on reference vectors. For example, the input below can be used to look at the projection
132 : of the vector connecting the instantanous configuraiton to a reference point in CV on a reference vector that has been specified in the PDB file.
133 :
134 : ```plumed
135 : #SETTINGS INPUTFILES=regtest/mapping/rt-pca-args/epath.pdb
136 :
137 : d1: DISTANCE ATOMS=1,2
138 : d2: DISTANCE ATOMS=3,4
139 : d3: DISTANCE ATOMS=5,6
140 : pca: PCAVARS ARG=d1,d2,d3 REFERENCE=regtest/mapping/rt-pca-args/epath.pdb
141 : PRINT ARG=pca_eig-1,pca_residual FILE=colvar
142 : ```
143 :
144 : The first set of argument values in this input file are the reference values for the arguments. The second any subsquent sets of arguments give the
145 : coefficients that should be used when constructing linear combinations.
146 :
147 : Notice, lastly, that you can also use a combination of argument values and atomic positions when specifying the reference configuration and the reference
148 : directions. If you are doing something this complicated, however, you are perhaps better working with the PLUMED input directly rather than this shortcut
149 : as you will need to take special measures to ensure that all your CVs are in the same units.
150 :
151 : */
152 : //+ENDPLUMEDOC
153 :
154 : class PCAVars : public ActionShortcut {
155 : public:
156 : static void registerKeywords(Keywords& keys);
157 : explicit PCAVars(const ActionOptions&);
158 : };
159 :
160 :
161 : PLUMED_REGISTER_ACTION(PCAVars,"PCAVARS")
162 :
163 11 : void PCAVars::registerKeywords( Keywords& keys ) {
164 11 : ActionShortcut::registerKeywords( keys );
165 11 : keys.add("compulsory","REFERENCE","a pdb file containing the set of reference configurations");
166 11 : keys.add("compulsory","TYPE","OPTIMAL-FAST","the manner in which distances are calculated. More information on the different "
167 : "metrics that are available in PLUMED can be found in the documentation for RMSD");
168 22 : keys.addInputKeyword("optional","ARG","scalar/vector","if there are arguments to be used specify them here");
169 11 : keys.addFlag("NOPBC",false,"do not use periodic boundary conditions when computing this quantity");
170 22 : keys.addOutputComponent("eig","default","vector","the projections on the eigenvalues");
171 22 : keys.addOutputComponent("residual","default","scalar","the residual distance that is not projected on any of the eigenvalues");
172 11 : keys.needsAction("RMSD");
173 11 : keys.needsAction("PDB2CONSTANT");
174 11 : keys.needsAction("TRANSPOSE");
175 11 : keys.needsAction("EUCLIDEAN_DISTANCE");
176 11 : keys.needsAction("CONCATENATE");
177 11 : keys.needsAction("COMBINE");
178 11 : keys.needsAction("CONSTANT");
179 11 : keys.needsAction("COMBINE");
180 11 : keys.needsAction("MATRIX_VECTOR_PRODUCT");
181 11 : keys.needsAction("CUSTOM");
182 11 : keys.needsAction("SUM");
183 11 : keys.needsAction("SELECT_COMPONENTS");
184 11 : }
185 :
186 5 : PCAVars::PCAVars( const ActionOptions& ao ):
187 : Action(ao),
188 5 : ActionShortcut(ao) {
189 : std::string reference;
190 10 : parse("REFERENCE",reference);
191 : // Create the object that holds the atomic positions by reading the first frame
192 5 : FILE* fp=std::fopen(reference.c_str(),"r");
193 5 : PDB pdb;
194 5 : if(!fp) {
195 0 : error("could not open reference file " + reference );
196 : }
197 5 : bool do_read=pdb.readFromFilepointer(fp,false,0.1);
198 5 : if( !do_read ) {
199 0 : plumed_merror("missing file " + reference );
200 : }
201 : std::string mtype;
202 5 : parse("TYPE",mtype);
203 :
204 5 : if( pdb.getPositions().size()>0 ) {
205 : // And now create the rmsd object
206 3 : std::string rmsd_line = getShortcutLabel() + "_at: RMSD DISPLACEMENT SQUARED NUMBER=1 REFERENCE=" + reference;
207 : bool nopbc;
208 3 : parseFlag("NOPBC",nopbc);
209 3 : if(nopbc) {
210 : rmsd_line += " NOPBC";
211 : }
212 : // Now create the RMSD object
213 6 : readInputLine( rmsd_line + " TYPE=" + mtype );
214 : }
215 5 : std::vector<std::string> argnames;
216 10 : parseVector("ARG",argnames);
217 : unsigned nargs=0;
218 : std::string instargs, refargs;
219 5 : std::vector<Value*> theargs;
220 5 : if( argnames.size()>0 ) {
221 2 : ActionWithArguments::interpretArgumentList( argnames, plumed.getActionSet(), this, theargs );
222 : }
223 9 : for(unsigned i=0; i<theargs.size(); ++i) {
224 4 : std::string iargn = Path::fixArgumentName( theargs[i]->getName() );
225 4 : nargs += theargs[i]->getNumberOfValues();
226 4 : if( theargs[i]->getNumberOfValues()>1 ) {
227 3 : readInputLine( getShortcutLabel() + "_ref_" + iargn + "T: PDB2CONSTANT NUMBER=1 REFERENCE=" + reference + " ARG=" + theargs[i]->getName() );
228 5 : readInputLine( getShortcutLabel() + "_ref_" + iargn + ": TRANSPOSE ARG=" + getShortcutLabel() + "_ref_" + iargn + "T");
229 : } else {
230 12 : readInputLine( getShortcutLabel() + "_ref_" + iargn + ": PDB2CONSTANT NUMBER=1 REFERENCE=" + reference + " ARG=" + theargs[i]->getName() );
231 : }
232 4 : if( i==0 ) {
233 2 : instargs=" ARG1=" + theargs[i]->getName();
234 6 : refargs=" ARG2=" + getShortcutLabel() + "_ref_" + iargn;
235 : } else {
236 2 : instargs +="," + theargs[i]->getName();
237 6 : refargs +="," + getShortcutLabel() + "_ref_" + iargn;
238 : }
239 : }
240 5 : if( theargs.size()>0 ) {
241 4 : readInputLine( getShortcutLabel() + "_argdist: EUCLIDEAN_DISTANCE SQUARED" + instargs + refargs );
242 : }
243 5 : if( pdb.getPositions().size()>0 && theargs.size()>0 ) {
244 0 : readInputLine( getShortcutLabel() + ": CONCATENATE ARG=" + getShortcutLabel() + "_at.disp," + getShortcutLabel() + "_argdist_diffT");
245 0 : readInputLine( getShortcutLabel() + "_dist: COMBINE ARG=" + getShortcutLabel() + "_at.dist," + getShortcutLabel() + "_argdist PERIODIC=NO");
246 : }
247 :
248 : // Get the displace stuff
249 5 : std::vector<double> displace( pdb.getBeta() );
250 : double dtot = 0;
251 26 : for(unsigned i=0; i<displace.size(); ++i) {
252 21 : dtot += displace[i];
253 : }
254 26 : for(unsigned i=0; i<displace.size(); ++i) {
255 21 : displace[i] = displace[i] / dtot;
256 : }
257 :
258 : // Now read in the directions and create matheval objects to compute the pca components
259 5 : unsigned nfram=0, ncomp=0;
260 : std::string pvec;
261 13 : while( do_read ) {
262 13 : std::vector<double> argdir(nargs);
263 13 : PDB mypdb;
264 13 : do_read=mypdb.readFromFilepointer(fp,plumed.usingNaturalUnits(),0.1/plumed.getUnits().getLength());
265 13 : if( do_read ) {
266 8 : nfram++;
267 : // Normalize the eigenvector in the input
268 : double norm=0;
269 50 : for(unsigned i=0; i<mypdb.getPositions().size(); ++i) {
270 42 : norm += mypdb.getPositions()[i][0]*mypdb.getPositions()[i][0];
271 42 : norm += mypdb.getPositions()[i][1]*mypdb.getPositions()[i][1];
272 42 : norm += mypdb.getPositions()[i][2]*mypdb.getPositions()[i][2];
273 : }
274 : unsigned k=0;
275 12 : for(unsigned i=0; i<theargs.size(); ++i) {
276 4 : std::vector<double> argval( theargs[i]->getNumberOfValues() );
277 4 : if( !mypdb.getArgumentValue(theargs[i]->getName(), argval) ) {
278 0 : error("argument " + theargs[i]->getName() + " was not set in pdb input");
279 : }
280 10 : for(unsigned j=0; j<argval.size(); ++j) {
281 6 : argdir[k] = argval[j];
282 6 : norm += argdir[k]*argdir[k];
283 6 : k++;
284 : }
285 4 : }
286 8 : norm = sqrt( norm );
287 8 : std::vector<double> normed_coeffs( 3*mypdb.getPositions().size() );
288 50 : for(unsigned i=0; i<mypdb.getPositions().size(); ++i) {
289 42 : if( mtype=="SIMPLE" ) {
290 28 : normed_coeffs[0*mypdb.getPositions().size()+i] = mypdb.getPositions()[i][0] / norm;
291 28 : normed_coeffs[1*mypdb.getPositions().size()+i] = mypdb.getPositions()[i][1] / norm;
292 28 : normed_coeffs[2*mypdb.getPositions().size()+i] = mypdb.getPositions()[i][2] / norm;
293 : } else {
294 14 : normed_coeffs[0*mypdb.getPositions().size()+i] = sqrt(displace[i])*mypdb.getPositions()[i][0] / norm;
295 14 : normed_coeffs[1*mypdb.getPositions().size()+i] = sqrt(displace[i])*mypdb.getPositions()[i][1] / norm;
296 14 : normed_coeffs[2*mypdb.getPositions().size()+i] = sqrt(displace[i])*mypdb.getPositions()[i][2] / norm;
297 : }
298 : }
299 : std::string coeff1;
300 8 : if( mypdb.getPositions().size()>0 ) {
301 6 : if( nfram==1 ) {
302 3 : Tools::convert( normed_coeffs[0], pvec );
303 : } else {
304 3 : Tools::convert( normed_coeffs[0], coeff1 );
305 3 : pvec += "," + coeff1;
306 : }
307 126 : for(unsigned i=1; i<normed_coeffs.size(); ++i) {
308 120 : Tools::convert( normed_coeffs[i], coeff1 );
309 120 : pvec += "," + coeff1;
310 : }
311 6 : for(unsigned i=0; i<argdir.size(); ++i) {
312 0 : Tools::convert( argdir[i] / norm, coeff1 );
313 0 : pvec += "," + coeff1;
314 : }
315 2 : } else if( theargs.size()>0 ) {
316 2 : if( nfram==1 ) {
317 2 : Tools::convert( argdir[0] / norm, pvec );
318 : } else {
319 0 : Tools::convert( argdir[0] / norm, coeff1 );
320 0 : pvec += "," + coeff1;
321 : }
322 6 : for(unsigned i=1; i<argdir.size(); ++i) {
323 4 : Tools::convert( argdir[i] / norm, coeff1 );
324 4 : pvec += "," + coeff1;
325 : }
326 : }
327 8 : ncomp = 3*mypdb.getPositions().size() + nargs;
328 8 : } else {
329 : break;
330 : }
331 13 : }
332 5 : std::fclose(fp);
333 : std::string neig, ncols;
334 5 : Tools::convert( nfram, neig );
335 5 : Tools::convert( ncomp, ncols );
336 15 : readInputLine( getShortcutLabel() + "_peig: CONSTANT VALUES=" + pvec + " NROWS=" + neig + " NCOLS=" + ncols );
337 5 : if( pdb.getPositions().size()>0 && theargs.size()>0 ) {
338 0 : readInputLine( getShortcutLabel() + "_eig: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_peig," + getShortcutLabel() );
339 5 : } else if( pdb.getPositions().size()>0 ) {
340 12 : readInputLine( getShortcutLabel() + "_eig: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_peig," + getShortcutLabel() + "_at.disp");
341 2 : } else if( theargs.size()>0 ) {
342 8 : readInputLine( getShortcutLabel() + "_eig: MATRIX_VECTOR_PRODUCT ARG=" + getShortcutLabel() + "_peig," + getShortcutLabel() + "_argdist_diffT");
343 : }
344 13 : for(unsigned i=0; i<nfram; ++i) {
345 : std::string num;
346 8 : Tools::convert( i+1, num );
347 32 : readInputLine( getShortcutLabel() + "_eig-" + num + ": SELECT_COMPONENTS ARG=" + getShortcutLabel() + "_eig COMPONENTS=" + num );
348 : }
349 10 : readInputLine( getShortcutLabel() + "_eig2: CUSTOM ARG=" + getShortcutLabel() + "_eig FUNC=x*x PERIODIC=NO");
350 10 : readInputLine( getShortcutLabel() + "_eigsum2: SUM ARG=" + getShortcutLabel() + "_eig2 PERIODIC=NO");
351 5 : if( pdb.getPositions().size()>0 && theargs.size()>0 ) {
352 0 : readInputLine( getShortcutLabel() + "_residual: CUSTOM ARG=" + getShortcutLabel() + "_dist," + getShortcutLabel() + "_eigsum2 FUNC=sqrt(x-y) PERIODIC=NO");
353 5 : } else if( pdb.getPositions().size()>0 ) {
354 12 : readInputLine( getShortcutLabel() + "_residual: CUSTOM ARG=" + getShortcutLabel() + "_at.dist," + getShortcutLabel() + "_eigsum2 FUNC=sqrt(x-y) PERIODIC=NO");
355 2 : } else if( theargs.size()>0 ) {
356 8 : readInputLine( getShortcutLabel() + "_residual: CUSTOM ARG=" + getShortcutLabel() + "_argdist," + getShortcutLabel() + "_eigsum2 FUNC=sqrt(x-y) PERIODIC=NO");
357 : }
358 147 : }
359 :
360 : }
361 : }
362 :
363 :
|
