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