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