Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-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 "Colvar.h"
23 : #include "core/Atoms.h"
24 : #include "core/PlumedMain.h"
25 : #include "ActionRegister.h"
26 : #include "tools/PDB.h"
27 : #include "tools/RMSD.h"
28 : #include "tools/Tools.h"
29 : #include <memory>
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace colvar {
35 :
36 4 : class PCARMSD : public Colvar {
37 :
38 : std::unique_ptr<PLMD::RMSD> rmsd;
39 : bool squared;
40 : bool nopbc;
41 : std::vector< std::vector<Vector> > eigenvectors;
42 : std::vector<PDB> pdbv;
43 : std::vector<string> pca_names;
44 : public:
45 : explicit PCARMSD(const ActionOptions&);
46 : virtual void calculate();
47 : static void registerKeywords(Keywords& keys);
48 : };
49 :
50 :
51 : using namespace std;
52 :
53 : //+PLUMEDOC DCOLVAR PCARMSD
54 : /*
55 : Calculate the PCA components ( see \cite Sutto:2010 and \cite spiwok ) for a number of provided eigenvectors and an average structure. Performs optimal alignment at every step and reports the rmsd so you know if you are far or close from the average structure.
56 : It takes the average structure and eigenvectors in form of a pdb.
57 : Note that beta and occupancy values in the pdb are neglected and all the weights are placed to 1 (differently from the RMSD colvar for example)
58 :
59 : \par Examples
60 :
61 : \plumedfile
62 : PCARMSD AVERAGE=file.pdb EIGENVECTORS=eigenvectors.pdb
63 : \endplumedfile
64 :
65 : The input is taken so to be compatible with the output you get from g_covar utility of gromacs (suitably adapted to have a pdb input format).
66 : The reference configuration (file.pdb) will thus be in a file that looks something like this:
67 :
68 : \auxfile{file.pdb}
69 : TITLE Average structure
70 : MODEL 1
71 : ATOM 1 CL ALA 1 1.042 -3.070 0.946 1.00 0.00
72 : ATOM 5 CLP ALA 1 0.416 -2.033 0.132 1.00 0.00
73 : ATOM 6 OL ALA 1 0.415 -2.082 -0.976 1.00 0.00
74 : ATOM 7 NL ALA 1 -0.134 -1.045 0.677 1.00 0.00
75 : ATOM 9 CA ALA 1 -0.774 0.053 0.003 1.00 0.00
76 : TER
77 : ENDMDL
78 : \endauxfile
79 :
80 : while the eigenvectors will be in a pdb file (eigenvectors.pdb) that looks something like this:
81 :
82 : \auxfile{eigenvectors.pdb}
83 : TITLE frame t= -1.000
84 : MODEL 1
85 : ATOM 1 CL ALA 1 1.194 -2.988 0.724 1.00 0.00
86 : ATOM 5 CLP ALA 1 -0.996 0.042 0.144 1.00 0.00
87 : ATOM 6 OL ALA 1 -1.246 -0.178 -0.886 1.00 0.00
88 : ATOM 7 NL ALA 1 -2.296 0.272 0.934 1.00 0.00
89 : ATOM 9 CA ALA 1 -0.436 2.292 0.814 1.00 0.00
90 : TER
91 : ENDMDL
92 : TITLE frame t= 0.000
93 : MODEL 1
94 : ATOM 1 CL ALA 1 1.042 -3.070 0.946 1.00 0.00
95 : ATOM 5 CLP ALA 1 -0.774 0.053 0.003 1.00 0.00
96 : ATOM 6 OL ALA 1 -0.849 -0.166 -1.034 1.00 0.00
97 : ATOM 7 NL ALA 1 -2.176 0.260 0.563 1.00 0.00
98 : ATOM 9 CA ALA 1 0.314 1.825 0.962 1.00 0.00
99 : TER
100 : ENDMDL
101 :
102 : \endauxfile
103 :
104 : */
105 : //+ENDPLUMEDOC
106 :
107 7360 : PLUMED_REGISTER_ACTION(PCARMSD,"PCARMSD")
108 :
109 3 : void PCARMSD::registerKeywords(Keywords& keys) {
110 3 : Colvar::registerKeywords(keys);
111 12 : keys.add("compulsory","AVERAGE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
112 12 : keys.add("compulsory","EIGENVECTORS","a file in pdb format containing the reference structure and the atoms involved in the CV.");
113 : //useCustomisableComponents(keys);
114 12 : keys.addOutputComponent("eig","default","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
115 12 : keys.addOutputComponent("residual","default","the distance of the present configuration from the configuration supplied as AVERAGE in terms of mean squared displacement after optimal alignment ");
116 9 : keys.addFlag("SQUARED-ROOT",false," This should be set if you want RMSD instead of mean squared displacement ");
117 9 : keys.addFlag("SQUARED_ROOT",false," Same as SQUARED-ROOT");
118 3 : }
119 :
120 2 : PCARMSD::PCARMSD(const ActionOptions&ao):
121 : PLUMED_COLVAR_INIT(ao),
122 : squared(true),
123 6 : nopbc(false)
124 : {
125 : string f_average;
126 4 : parse("AVERAGE",f_average);
127 : string type;
128 2 : type.assign("OPTIMAL");
129 : string f_eigenvectors;
130 4 : parse("EIGENVECTORS",f_eigenvectors);
131 4 : bool sq; parseFlag("SQUARED-ROOT",sq);
132 4 : if(!sq) parseFlag("SQUARED_ROOT",sq);
133 2 : if (sq) { squared=false; }
134 4 : parseFlag("NOPBC",nopbc);
135 2 : checkRead();
136 :
137 4 : PDB pdb;
138 :
139 : // read everything in ang and transform to nm if we are not in natural units
140 4 : if( !pdb.read(f_average,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
141 0 : error("missing input file " + f_average );
142 :
143 2 : rmsd.reset( new RMSD() );
144 : bool remove_com=true;
145 : bool normalize_weights=true;
146 : // here align and displace are a simple vector of ones
147 48 : std::vector<double> align; align=pdb.getOccupancy(); for(unsigned i=0; i<align.size(); i++) {align[i]=1.;} ;
148 48 : std::vector<double> displace; displace=pdb.getBeta(); for(unsigned i=0; i<displace.size(); i++) {displace[i]=1.;} ;
149 : // reset again to reimpose unifrom weights (safe to disable this)
150 2 : rmsd->set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
151 2 : requestAtoms( pdb.getAtomNumbers() );
152 :
153 6 : addComponentWithDerivatives("residual"); componentIsNotPeriodic("residual");
154 :
155 4 : log.printf(" average from file %s\n",f_average.c_str());
156 4 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
157 2 : log.printf(" with indices : ");
158 70 : for(unsigned i=0; i<pdb.getAtomNumbers().size(); ++i) {
159 22 : if(i%25==0) log<<"\n";
160 44 : log.printf("%d ",pdb.getAtomNumbers()[i].serial());
161 : }
162 2 : log.printf("\n");
163 4 : log.printf(" method for alignment : %s \n",type.c_str() );
164 2 : if(nopbc) log.printf(" without periodic boundary conditions\n");
165 1 : else log.printf(" using periodic boundary conditions\n");
166 :
167 6 : log<<" Bibliography "<<plumed.cite("Spiwok, Lipovova and Kralova, JPCB, 111, 3073 (2007) ");
168 6 : log<<" "<<plumed.cite( "Sutto, D'Abramo, Gervasio, JCTC, 6, 3640 (2010)");
169 :
170 : // now get the eigenvectors
171 : // open the file
172 2 : FILE* fp=fopen(f_eigenvectors.c_str(),"r");
173 : std::vector<AtomNumber> aaa;
174 : unsigned neigenvects;
175 2 : neigenvects=0;
176 2 : if (fp!=NULL)
177 : {
178 4 : log<<" Opening the eigenvectors file "<<f_eigenvectors.c_str()<<"\n";
179 : bool do_read=true;
180 72 : while (do_read) {
181 142 : PDB mypdb;
182 : // check the units for reading this file: how can they make sense?
183 144 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
184 72 : if(do_read) {
185 70 : neigenvects++;
186 140 : if(mypdb.getAtomNumbers().size()==0) error("number of atoms in a frame should be more than zero");
187 70 : unsigned nat=mypdb.getAtomNumbers().size();
188 140 : if(nat!=mypdb.getAtomNumbers().size()) error("frames should have the same number of atoms");
189 70 : if(aaa.empty()) aaa=mypdb.getAtomNumbers();
190 140 : if(aaa!=mypdb.getAtomNumbers()) error("frames should contain same atoms in same order");
191 140 : log<<" Found eigenvector: "<<neigenvects<<" containing "<<mypdb.getAtomNumbers().size()<<" atoms\n";
192 70 : pdbv.push_back(mypdb);
193 70 : eigenvectors.push_back(mypdb.getPositions());
194 : } else {break ;}
195 : }
196 2 : fclose (fp);
197 4 : log<<" Found total "<<neigenvects<< " eigenvectors in the file "<<f_eigenvectors.c_str()<<" \n";
198 2 : if(neigenvects==0) error("at least one eigenvector is expected");
199 : }
200 : // the components
201 142 : for(unsigned i=0; i<neigenvects; i++) {
202 70 : std::string num; Tools::convert( i, num );
203 210 : string name; name=string("eig-")+num;
204 70 : pca_names.push_back(name);
205 70 : addComponentWithDerivatives(name); componentIsNotPeriodic(name);
206 : }
207 2 : turnOnDerivatives();
208 :
209 2 : }
210 :
211 : // calculator
212 1092 : void PCARMSD::calculate() {
213 1092 : if(!nopbc) makeWhole();
214 1092 : Tensor rotation,invrotation;
215 : Matrix<std::vector<Vector> > drotdpos(3,3);
216 : std::vector<Vector> alignedpos;
217 : std::vector<Vector> centeredpos;
218 : std::vector<Vector> centeredref;
219 : std::vector<Vector> ddistdpos;
220 2184 : double r=rmsd->calc_PCAelements( getPositions(), ddistdpos, rotation, drotdpos, alignedpos,centeredpos, centeredref,squared);
221 1092 : invrotation=rotation.transpose();
222 :
223 2184 : Value* verr=getPntrToComponent("residual");
224 : verr->set(r);
225 25116 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
226 24024 : setAtomsDerivatives (verr,iat,ddistdpos[iat]);
227 : }
228 :
229 : std::vector< Vector > der;
230 1092 : der.resize(getNumberOfAtoms());
231 :
232 :
233 116844 : for(unsigned i=0; i<eigenvectors.size(); i++) {
234 76440 : Value* value=getPntrToComponent(pca_names[i].c_str());
235 : double val; val=0.;
236 879060 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
237 1261260 : val+=dotProduct(alignedpos[iat]-centeredref[iat],eigenvectors[i][iat]); der[iat].zero();
238 : }
239 : value->set(val);
240 : // here the loop is reversed to better suit the structure of the derivative of the rotation matrix
241 : double tmp1;
242 267540 : for(unsigned a=0; a<3; a++) {
243 802620 : for(unsigned b=0; b<3; b++) {
244 : tmp1=0.;
245 7911540 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
246 11351340 : tmp1+=centeredpos[n][b]*eigenvectors[i][n][a];
247 : }
248 7911540 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
249 11351340 : der[iat]+=drotdpos[a][b][iat]*tmp1;
250 : }
251 : }
252 : }
253 38220 : Vector v1;
254 879060 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
255 1261260 : v1+=(1./getNumberOfAtoms())*matmul(invrotation,eigenvectors[i][n]);
256 : }
257 879060 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
258 1261260 : der[iat]+=matmul(invrotation,eigenvectors[i][iat])-v1;
259 840840 : setAtomsDerivatives (value,iat,der[iat]);
260 : }
261 : }
262 :
263 40404 : for(int i=0; i<getNumberOfComponents(); ++i) setBoxDerivativesNoPbc( getPntrToComponent(i) );
264 :
265 1092 : }
266 :
267 : }
268 5517 : }
269 :
270 :
271 :
|
