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 : Information about this method can be found in the reference papers in the bibliography below. An example input is provided below:
51 :
52 : ```plumed
53 : #SETTINGS INPUTFILES=regtest/trajectories/pca/average.pdb,regtest/trajectories/pca/eigenvec.pdb
54 : PCARMSD ...
55 : AVERAGE=regtest/trajectories/pca/average.pdb
56 : EIGENVECTORS=regtest/trajectories/pca/eigenvec.pdb
57 : ...
58 : ```
59 :
60 : This input performs optimal alignment at every step and reports the rmsd so you know if you are far or close from the average structure.
61 : It takes the average structure and eigenvectors in form of a pdb.
62 : 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)
63 :
64 : */
65 : //+ENDPLUMEDOC
66 :
67 : PLUMED_REGISTER_ACTION(PCARMSD,"PCARMSD")
68 :
69 7 : void PCARMSD::registerKeywords(Keywords& keys) {
70 7 : Colvar::registerKeywords(keys);
71 7 : keys.add("compulsory","AVERAGE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
72 7 : keys.add("compulsory","EIGENVECTORS","a file in pdb format containing the reference structure and the atoms involved in the CV.");
73 14 : keys.addOutputComponent("eig","default","scalar","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
74 14 : 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 ");
75 7 : keys.addFlag("SQUARED_ROOT",false," This should be set if you want RMSD instead of mean squared displacement ");
76 7 : keys.addDOI("10.1021/ct100413b");
77 7 : keys.addDOI("10.1021/jp068587c");
78 7 : }
79 :
80 3 : PCARMSD::PCARMSD(const ActionOptions&ao):
81 : PLUMED_COLVAR_INIT(ao),
82 3 : squared(true),
83 3 : nopbc(false) {
84 : std::string f_average;
85 7 : parse("AVERAGE",f_average);
86 : std::string type;
87 3 : type.assign("OPTIMAL");
88 : std::string f_eigenvectors;
89 3 : parse("EIGENVECTORS",f_eigenvectors);
90 : bool sq;
91 3 : parseFlag("SQUARED_ROOT",sq);
92 3 : if (sq) {
93 0 : squared=false;
94 : }
95 4 : parseFlag("NOPBC",nopbc);
96 3 : checkRead();
97 :
98 3 : PDB pdb;
99 :
100 : // read everything in ang and transform to nm if we are not in natural units
101 3 : if( !pdb.read(f_average,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
102 0 : error("missing input file " + f_average );
103 : }
104 :
105 3 : rmsd=Tools::make_unique<RMSD>();
106 : bool remove_com=true;
107 : bool normalize_weights=true;
108 : // here align and displace are a simple vector of ones
109 3 : std::vector<double> align;
110 3 : align=pdb.getOccupancy();
111 37 : for(unsigned i=0; i<align.size(); i++) {
112 34 : align[i]=1.;
113 : } ;
114 3 : std::vector<double> displace;
115 3 : displace=pdb.getBeta();
116 37 : for(unsigned i=0; i<displace.size(); i++) {
117 34 : displace[i]=1.;
118 : } ;
119 : // reset again to reimpose unifrom weights (safe to disable this)
120 3 : rmsd->set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
121 3 : requestAtoms( pdb.getAtomNumbers() );
122 :
123 6 : addComponentWithDerivatives("residual");
124 3 : componentIsNotPeriodic("residual");
125 :
126 3 : log.printf(" average from file %s\n",f_average.c_str());
127 3 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
128 3 : log.printf(" with indices : ");
129 37 : for(unsigned i=0; i<pdb.getAtomNumbers().size(); ++i) {
130 34 : if(i%25==0) {
131 3 : log<<"\n";
132 : }
133 34 : log.printf("%d ",pdb.getAtomNumbers()[i].serial());
134 : }
135 3 : log.printf("\n");
136 3 : log.printf(" method for alignment : %s \n",type.c_str() );
137 3 : if(nopbc) {
138 1 : log.printf(" without periodic boundary conditions\n");
139 : } else {
140 2 : log.printf(" using periodic boundary conditions\n");
141 : }
142 :
143 3 : log<<" Bibliography "<<plumed.cite("Spiwok, Lipovova and Kralova, JPCB, 111, 3073 (2007) ");
144 3 : log<<" "<<plumed.cite( "Sutto, D'Abramo, Gervasio, JCTC, 6, 3640 (2010)");
145 :
146 : unsigned neigenvects;
147 3 : neigenvects=0;
148 : // now get the eigenvectors
149 : // open the file
150 3 : if (FILE* fp=this->fopen(f_eigenvectors.c_str(),"r")) {
151 : // call fclose when exiting this block
152 3 : auto deleter=[this](FILE* f) {
153 3 : this->fclose(f);
154 3 : };
155 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
156 :
157 3 : std::vector<AtomNumber> aaa;
158 3 : log<<" Opening the eigenvectors file "<<f_eigenvectors.c_str()<<"\n";
159 : bool do_read=true;
160 : unsigned nat=0;
161 73 : while (do_read) {
162 73 : PDB mypdb;
163 : // check the units for reading this file: how can they make sense?
164 73 : do_read=mypdb.readFromFilepointer(fp,usingNaturalUnits(),0.1/getUnits().getLength());
165 73 : if(do_read) {
166 71 : neigenvects++;
167 71 : if(mypdb.getAtomNumbers().size()==0) {
168 0 : error("number of atoms in a frame should be more than zero");
169 : }
170 71 : if(nat==0) {
171 3 : nat=mypdb.getAtomNumbers().size();
172 : }
173 71 : if(nat!=mypdb.getAtomNumbers().size()) {
174 0 : error("frames should have the same number of atoms");
175 : }
176 71 : if(aaa.empty()) {
177 3 : aaa=mypdb.getAtomNumbers();
178 : }
179 71 : if(aaa!=mypdb.getAtomNumbers()) {
180 0 : error("frames should contain same atoms in same order");
181 : }
182 71 : if(mypdb.getAtomNumbers().size()!=getNumberOfAtoms()) {
183 5 : error("eigenvector " + std::to_string(neigenvects) + " contains " + std::to_string(mypdb.getAtomNumbers().size()) +
184 2 : " atoms but the AVERAGE structure contains " + std::to_string(getNumberOfAtoms()));
185 : }
186 70 : log<<" Found eigenvector: "<<neigenvects<<" containing "<<mypdb.getAtomNumbers().size()<<" atoms\n";
187 70 : pdbv.push_back(mypdb);
188 70 : eigenvectors.push_back(mypdb.getPositions());
189 : } else {
190 : break ;
191 : }
192 73 : }
193 2 : log<<" Found total "<<neigenvects<< " eigenvectors in the file "<<f_eigenvectors.c_str()<<" \n";
194 2 : if(neigenvects==0) {
195 0 : error("at least one eigenvector is expected");
196 : }
197 4 : }
198 : // the components
199 72 : for(unsigned i=0; i<neigenvects; i++) {
200 : std::string num;
201 70 : Tools::convert( i, num );
202 70 : std::string compName=std::string("eig-")+num;
203 70 : pca_names.push_back(compName);
204 70 : addComponentWithDerivatives(compName);
205 70 : componentIsNotPeriodic(compName);
206 : }
207 2 : turnOnDerivatives();
208 :
209 13 : }
210 :
211 : // calculator
212 1092 : void PCARMSD::calculate() {
213 1092 : if(!nopbc) {
214 546 : makeWhole();
215 : }
216 1092 : Tensor rotation,invrotation;
217 1092 : Matrix<std::vector<Vector> > drotdpos(3,3);
218 1092 : std::vector<Vector> alignedpos;
219 1092 : std::vector<Vector> centeredpos;
220 1092 : std::vector<Vector> centeredref;
221 1092 : std::vector<Vector> ddistdpos;
222 1092 : double r=rmsd->calc_PCAelements( getPositions(), ddistdpos, rotation, drotdpos, alignedpos,centeredpos, centeredref,squared);
223 1092 : invrotation=rotation.transpose();
224 :
225 2184 : Value* verr=getPntrToComponent("residual");
226 : verr->set(r);
227 13104 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
228 12012 : setAtomsDerivatives (verr,iat,ddistdpos[iat]);
229 : }
230 :
231 1092 : std::vector< Vector > der;
232 1092 : der.resize(getNumberOfAtoms());
233 :
234 :
235 39312 : for(unsigned i=0; i<eigenvectors.size(); i++) {
236 38220 : Value* value=getPntrToComponent(pca_names[i].c_str());
237 : double val;
238 : val=0.;
239 458640 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
240 420420 : val+=dotProduct(alignedpos[iat]-centeredref[iat],eigenvectors[i][iat]);
241 : der[iat].zero();
242 : }
243 : value->set(val);
244 : // here the loop is reversed to better suit the structure of the derivative of the rotation matrix
245 : double tmp1;
246 152880 : for(unsigned a=0; a<3; a++) {
247 458640 : for(unsigned b=0; b<3; b++) {
248 : tmp1=0.;
249 4127760 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
250 3783780 : tmp1+=centeredpos[n][b]*eigenvectors[i][n][a];
251 : }
252 4127760 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
253 3783780 : der[iat]+=drotdpos[a][b][iat]*tmp1;
254 : }
255 : }
256 : }
257 : Vector v1;
258 458640 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
259 420420 : v1+=(1./getNumberOfAtoms())*matmul(invrotation,eigenvectors[i][n]);
260 : }
261 458640 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
262 420420 : der[iat]+=matmul(invrotation,eigenvectors[i][iat])-v1;
263 420420 : setAtomsDerivatives (value,iat,der[iat]);
264 : }
265 : }
266 :
267 40404 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
268 39312 : setBoxDerivativesNoPbc( getPntrToComponent(i) );
269 : }
270 :
271 2184 : }
272 :
273 : }
274 : }
275 :
276 :
277 :
|
