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