Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2019 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 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace colvar {
34 :
35 : class PCARMSD : public Colvar {
36 :
37 : PLMD::RMSD* rmsd;
38 : bool squared;
39 : std::vector< std::vector<Vector> > eigenvectors;
40 : std::vector<PDB> pdbv;
41 : std::vector<string> pca_names;
42 : public:
43 : explicit PCARMSD(const ActionOptions&);
44 : ~PCARMSD();
45 : virtual void calculate();
46 : static void registerKeywords(Keywords& keys);
47 : };
48 :
49 :
50 : using namespace std;
51 :
52 : //+PLUMEDOC DCOLVAR PCARMSD
53 : /*
54 : Calculate the PCA components ( see \cite Sutto:2010 and \cite spiwok ) for a number of provided eigenvectors and an average structure. Performs optimal alignment at every step and reports the rmsd so you know if you are far or close from the average structure.
55 : It takes the average structure and eigenvectors in form of a pdb.
56 : 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)
57 :
58 : \par Examples
59 :
60 : \plumedfile
61 : PCARMSD AVERAGE=file.pdb EIGENVECTORS=eigenvectors.pdb
62 : \endplumedfile
63 :
64 : 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).
65 :
66 : */
67 : //+ENDPLUMEDOC
68 :
69 4822 : PLUMED_REGISTER_ACTION(PCARMSD,"PCARMSD")
70 :
71 2 : void PCARMSD::registerKeywords(Keywords& keys) {
72 2 : Colvar::registerKeywords(keys);
73 2 : keys.add("compulsory","AVERAGE","a file in pdb format containing the reference structure and the atoms involved in the CV.");
74 2 : keys.add("compulsory","EIGENVECTORS","a file in pdb format containing the reference structure and the atoms involved in the CV.");
75 : //useCustomisableComponents(keys);
76 2 : keys.addOutputComponent("eig","default","the projections on each eigenvalue are stored on values labeled eig-1, eig-2, ...");
77 2 : keys.addOutputComponent("residual","default","the distance of the present configuration from the configuration supplied as AVERAGE in terms of MSD after optimal alignment ");
78 2 : keys.addFlag("SQUARED-ROOT",false," This should be setted if you want RMSD instead of MSD ");
79 2 : keys.addFlag("SQUARED_ROOT",false," Same as SQUARED-ROOT");
80 2 : }
81 :
82 1 : PCARMSD::PCARMSD(const ActionOptions&ao):
83 1 : PLUMED_COLVAR_INIT(ao),squared(true)
84 : {
85 1 : string f_average;
86 1 : parse("AVERAGE",f_average);
87 2 : string type;
88 1 : type.assign("OPTIMAL");
89 2 : string f_eigenvectors;
90 1 : parse("EIGENVECTORS",f_eigenvectors);
91 1 : bool sq; parseFlag("SQUARED-ROOT",sq);
92 1 : if(!sq) parseFlag("SQUARED_ROOT",sq);
93 1 : if (sq) { squared=false; }
94 1 : checkRead();
95 :
96 2 : PDB pdb;
97 :
98 : // read everything in ang and transform to nm if we are not in natural units
99 1 : if( !pdb.read(f_average,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength()) )
100 0 : error("missing input file " + f_average );
101 :
102 1 : rmsd = new RMSD();
103 1 : bool remove_com=true;
104 1 : bool normalize_weights=true;
105 : // here align and displace are a simple vector of ones
106 2 : std::vector<double> align; align=pdb.getOccupancy(); for(unsigned i=0; i<align.size(); i++) {align[i]=1.;} ;
107 2 : std::vector<double> displace; displace=pdb.getBeta(); for(unsigned i=0; i<displace.size(); i++) {displace[i]=1.;} ;
108 : // reset again to reimpose unifrom weights (safe to disable this)
109 1 : rmsd->set(align,displace,pdb.getPositions(),type,remove_com,normalize_weights);
110 1 : requestAtoms( pdb.getAtomNumbers() );
111 :
112 1 : addComponentWithDerivatives("residual"); componentIsNotPeriodic("residual");
113 :
114 1 : log.printf(" average from file %s\n",f_average.c_str());
115 1 : log.printf(" which contains %d atoms\n",getNumberOfAtoms());
116 1 : log.printf(" method for alignment : %s \n",type.c_str() );
117 :
118 1 : log<<" Bibliography "<<plumed.cite("Spiwok, Lipovova and Kralova, JPCB, 111, 3073 (2007) ");
119 1 : log<<" "<<plumed.cite( "Sutto, D'Abramo, Gervasio, JCTC, 6, 3640 (2010)");
120 :
121 : // now get the eigenvectors
122 : // open the file
123 1 : FILE* fp=fopen(f_eigenvectors.c_str(),"r");
124 2 : std::vector<AtomNumber> aaa;
125 : unsigned neigenvects;
126 1 : neigenvects=0;
127 1 : if (fp!=NULL)
128 : {
129 1 : log<<" Opening the eigenvectors file "<<f_eigenvectors.c_str()<<"\n";
130 1 : bool do_read=true;
131 37 : while (do_read) {
132 36 : PDB mypdb;
133 : // check the units for reading this file: how can they make sense?
134 36 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
135 36 : if(do_read) {
136 35 : neigenvects++;
137 35 : if(mypdb.getAtomNumbers().size()==0) error("number of atoms in a frame should be more than zero");
138 35 : unsigned nat=mypdb.getAtomNumbers().size();
139 35 : if(nat!=mypdb.getAtomNumbers().size()) error("frames should have the same number of atoms");
140 35 : if(aaa.empty()) aaa=mypdb.getAtomNumbers();
141 35 : if(aaa!=mypdb.getAtomNumbers()) error("frames should contain same atoms in same order");
142 35 : log<<" Found eigenvector: "<<neigenvects<<" containing "<<mypdb.getAtomNumbers().size()<<" atoms\n";
143 35 : pdbv.push_back(mypdb);
144 35 : eigenvectors.push_back(mypdb.getPositions());
145 1 : } else {break ;}
146 35 : }
147 1 : fclose (fp);
148 1 : log<<" Found total "<<neigenvects<< " eigenvectors in the file "<<f_eigenvectors.c_str()<<" \n";
149 1 : if(neigenvects==0) error("at least one eigenvector is expected");
150 : }
151 : // the components
152 36 : for(unsigned i=0; i<neigenvects; i++) {
153 35 : std::string num; Tools::convert( i, num );
154 70 : string name; name=string("eig-")+num;
155 35 : pca_names.push_back(name);
156 35 : addComponentWithDerivatives(name); componentIsNotPeriodic(name);
157 35 : }
158 2 : turnOnDerivatives();
159 :
160 1 : }
161 :
162 4 : PCARMSD::~PCARMSD() {
163 1 : delete rmsd;
164 3 : }
165 :
166 :
167 : // calculator
168 546 : void PCARMSD::calculate() {
169 546 : Tensor rotation,invrotation;
170 546 : Matrix<std::vector<Vector> > drotdpos(3,3);
171 1092 : std::vector<Vector> alignedpos;
172 1092 : std::vector<Vector> centeredpos;
173 1092 : std::vector<Vector> centeredref;
174 1092 : std::vector<Vector> ddistdpos;
175 546 : double r=rmsd->calc_PCAelements( getPositions(), ddistdpos, rotation, drotdpos, alignedpos,centeredpos, centeredref,squared);
176 546 : invrotation=rotation.transpose();
177 :
178 546 : Value* verr=getPntrToComponent("residual");
179 546 : verr->set(r);
180 6552 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
181 6006 : setAtomsDerivatives (verr,iat,ddistdpos[iat]);
182 : }
183 :
184 1092 : std::vector< Vector > der;
185 546 : der.resize(getNumberOfAtoms());
186 :
187 :
188 19656 : for(unsigned i=0; i<eigenvectors.size(); i++) {
189 19110 : Value* value=getPntrToComponent(pca_names[i].c_str());
190 19110 : double val; val=0.;
191 229320 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
192 210210 : val+=dotProduct(alignedpos[iat]-centeredref[iat],eigenvectors[i][iat]); der[iat].zero();
193 : }
194 19110 : value->set(val);
195 : // here the loop is reversed to better suit the structure of the derivative of the rotation matrix
196 : double tmp1;
197 76440 : for(unsigned a=0; a<3; a++) {
198 229320 : for(unsigned b=0; b<3; b++) {
199 171990 : tmp1=0.;
200 2063880 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
201 1891890 : tmp1+=centeredpos[n][b]*eigenvectors[i][n][a];
202 : }
203 2063880 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
204 1891890 : der[iat]+=drotdpos[a][b][iat]*tmp1;
205 : }
206 : }
207 : }
208 19110 : Vector v1;
209 229320 : for(unsigned n=0; n<getNumberOfAtoms(); n++) {
210 210210 : v1+=(1./getNumberOfAtoms())*matmul(invrotation,eigenvectors[i][n]);
211 : }
212 229320 : for(unsigned iat=0; iat<getNumberOfAtoms(); iat++) {
213 210210 : der[iat]+=matmul(invrotation,eigenvectors[i][iat])-v1;
214 210210 : setAtomsDerivatives (value,iat,der[iat]);
215 : }
216 : }
217 :
218 1092 : for(int i=0; i<getNumberOfComponents(); ++i) setBoxDerivativesNoPbc( getPntrToComponent(i) );
219 :
220 546 : }
221 :
222 : }
223 4821 : }
224 :
225 :
226 :
|
