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","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
109 8 : 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 ");
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 : 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 : bool sq;
124 2 : parseFlag("SQUARED_ROOT",sq);
125 2 : if (sq) {
126 0 : squared=false;
127 : }
128 2 : parseFlag("NOPBC",nopbc);
129 2 : checkRead();
130 :
131 2 : PDB pdb;
132 :
133 : // read everything in ang and transform to nm if we are not in natural units
134 2 : if( !pdb.read(f_average,usingNaturalUnits(),0.1/getUnits().getLength()) ) {
135 0 : error("missing input file " + f_average );
136 : }
137 :
138 2 : rmsd=Tools::make_unique<RMSD>();
139 : bool remove_com=true;
140 : bool normalize_weights=true;
141 : // here align and displace are a simple vector of ones
142 : std::vector<double> align;
143 2 : align=pdb.getOccupancy();
144 24 : for(unsigned i=0; i<align.size(); i++) {
145 22 : align[i]=1.;
146 : } ;
147 : std::vector<double> displace;
148 2 : displace=pdb.getBeta();
149 24 : for(unsigned i=0; i<displace.size(); i++) {
150 22 : displace[i]=1.;
151 : } ;
152 : // reset again to reimpose unifrom weights (safe to disable this)
153 2 : rmsd->set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
154 2 : requestAtoms( pdb.getAtomNumbers() );
155 :
156 4 : addComponentWithDerivatives("residual");
157 2 : componentIsNotPeriodic("residual");
158 :
159 2 : log.printf(" average from file %s\n",f_average.c_str());
160 2 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
161 2 : log.printf(" with indices : ");
162 24 : for(unsigned i=0; i<pdb.getAtomNumbers().size(); ++i) {
163 22 : if(i%25==0) {
164 2 : log<<"\n";
165 : }
166 22 : log.printf("%d ",pdb.getAtomNumbers()[i].serial());
167 : }
168 2 : log.printf("\n");
169 2 : log.printf(" method for alignment : %s \n",type.c_str() );
170 2 : if(nopbc) {
171 1 : log.printf(" without periodic boundary conditions\n");
172 : } else {
173 1 : log.printf(" using periodic boundary conditions\n");
174 : }
175 :
176 4 : log<<" Bibliography "<<plumed.cite("Spiwok, Lipovova and Kralova, JPCB, 111, 3073 (2007) ");
177 4 : log<<" "<<plumed.cite( "Sutto, D'Abramo, Gervasio, JCTC, 6, 3640 (2010)");
178 :
179 : unsigned neigenvects;
180 2 : neigenvects=0;
181 : // now get the eigenvectors
182 : // open the file
183 2 : if (FILE* fp=this->fopen(f_eigenvectors.c_str(),"r")) {
184 : // call fclose when exiting this block
185 2 : auto deleter=[this](FILE* f) {
186 2 : this->fclose(f);
187 2 : };
188 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
189 :
190 : std::vector<AtomNumber> aaa;
191 2 : log<<" Opening the eigenvectors file "<<f_eigenvectors.c_str()<<"\n";
192 : bool do_read=true;
193 : unsigned nat=0;
194 72 : while (do_read) {
195 72 : PDB mypdb;
196 : // check the units for reading this file: how can they make sense?
197 72 : do_read=mypdb.readFromFilepointer(fp,usingNaturalUnits(),0.1/getUnits().getLength());
198 72 : if(do_read) {
199 70 : neigenvects++;
200 70 : if(mypdb.getAtomNumbers().size()==0) {
201 0 : error("number of atoms in a frame should be more than zero");
202 : }
203 70 : if(nat==0) {
204 2 : nat=mypdb.getAtomNumbers().size();
205 : }
206 70 : if(nat!=mypdb.getAtomNumbers().size()) {
207 0 : error("frames should have the same number of atoms");
208 : }
209 70 : if(aaa.empty()) {
210 2 : aaa=mypdb.getAtomNumbers();
211 : }
212 140 : if(aaa!=mypdb.getAtomNumbers()) {
213 0 : error("frames should contain same atoms in same order");
214 : }
215 70 : log<<" Found eigenvector: "<<neigenvects<<" containing "<<mypdb.getAtomNumbers().size()<<" atoms\n";
216 70 : pdbv.push_back(mypdb);
217 70 : eigenvectors.push_back(mypdb.getPositions());
218 : } else {
219 : break ;
220 : }
221 72 : }
222 2 : log<<" Found total "<<neigenvects<< " eigenvectors in the file "<<f_eigenvectors.c_str()<<" \n";
223 2 : if(neigenvects==0) {
224 0 : error("at least one eigenvector is expected");
225 : }
226 2 : }
227 : // the components
228 72 : for(unsigned i=0; i<neigenvects; i++) {
229 : std::string num;
230 70 : Tools::convert( i, num );
231 : std::string name;
232 140 : name=std::string("eig-")+num;
233 70 : pca_names.push_back(name);
234 70 : addComponentWithDerivatives(name);
235 70 : componentIsNotPeriodic(name);
236 : }
237 2 : turnOnDerivatives();
238 :
239 4 : }
240 :
241 : // calculator
242 1092 : void PCARMSD::calculate() {
243 1092 : if(!nopbc) {
244 546 : makeWhole();
245 : }
246 1092 : Tensor rotation,invrotation;
247 : Matrix<std::vector<Vector> > drotdpos(3,3);
248 : std::vector<Vector> alignedpos;
249 : std::vector<Vector> centeredpos;
250 : std::vector<Vector> centeredref;
251 : std::vector<Vector> ddistdpos;
252 1092 : double r=rmsd->calc_PCAelements( getPositions(), ddistdpos, rotation, drotdpos, alignedpos,centeredpos, centeredref,squared);
253 1092 : invrotation=rotation.transpose();
254 :
255 2184 : Value* verr=getPntrToComponent("residual");
256 : verr->set(r);
257 13104 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
258 12012 : setAtomsDerivatives (verr,iat,ddistdpos[iat]);
259 : }
260 :
261 : std::vector< Vector > der;
262 1092 : der.resize(getNumberOfAtoms());
263 :
264 :
265 39312 : for(unsigned i=0; i<eigenvectors.size(); i++) {
266 38220 : Value* value=getPntrToComponent(pca_names[i].c_str());
267 : double val;
268 : val=0.;
269 458640 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
270 420420 : val+=dotProduct(alignedpos[iat]-centeredref[iat],eigenvectors[i][iat]);
271 420420 : der[iat].zero();
272 : }
273 : value->set(val);
274 : // here the loop is reversed to better suit the structure of the derivative of the rotation matrix
275 : double tmp1;
276 152880 : for(unsigned a=0; a<3; a++) {
277 458640 : for(unsigned b=0; b<3; b++) {
278 : tmp1=0.;
279 4127760 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
280 3783780 : tmp1+=centeredpos[n][b]*eigenvectors[i][n][a];
281 : }
282 4127760 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
283 3783780 : der[iat]+=drotdpos[a][b][iat]*tmp1;
284 : }
285 : }
286 : }
287 38220 : Vector v1;
288 458640 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
289 420420 : v1+=(1./getNumberOfAtoms())*matmul(invrotation,eigenvectors[i][n]);
290 : }
291 458640 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
292 420420 : der[iat]+=matmul(invrotation,eigenvectors[i][iat])-v1;
293 420420 : setAtomsDerivatives (value,iat,der[iat]);
294 : }
295 : }
296 :
297 40404 : for(int i=0; i<getNumberOfComponents(); ++i) {
298 39312 : setBoxDerivativesNoPbc( getPntrToComponent(i) );
299 : }
300 :
301 1092 : }
302 :
303 : }
304 : }
305 :
306 :
307 :
|
