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