Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 "PathMSDBase.h"
23 : #include "Colvar.h"
24 : #include "ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/Atoms.h"
27 : #include "tools/PDB.h"
28 : #include "tools/RMSD.h"
29 : #include "tools/Tools.h"
30 : #include <cmath>
31 :
32 : using namespace std;
33 :
34 : namespace PLMD {
35 : namespace colvar {
36 :
37 24 : void PathMSDBase::registerKeywords(Keywords& keys) {
38 24 : Colvar::registerKeywords(keys);
39 24 : keys.remove("NOPBC");
40 24 : keys.add("compulsory","LAMBDA","the lambda parameter is needed for smoothing, is in the units of plumed");
41 24 : keys.add("compulsory","REFERENCE","the pdb is needed to provide the various milestones");
42 24 : keys.add("optional","NEIGH_SIZE","size of the neighbor list");
43 24 : keys.add("optional","NEIGH_STRIDE","how often the neighbor list needs to be calculated in time units");
44 24 : }
45 :
46 22 : PathMSDBase::PathMSDBase(const ActionOptions&ao):
47 : PLUMED_COLVAR_INIT(ao),
48 : neigh_size(-1),
49 : neigh_stride(-1),
50 22 : nframes(0)
51 : {
52 22 : parse("LAMBDA",lambda);
53 22 : parse("NEIGH_SIZE",neigh_size);
54 22 : parse("NEIGH_STRIDE",neigh_stride);
55 22 : parse("REFERENCE",reference);
56 :
57 : // open the file
58 22 : FILE* fp=fopen(reference.c_str(),"r");
59 22 : std::vector<AtomNumber> aaa;
60 22 : if (fp!=NULL)
61 : {
62 22 : log<<"Opening reference file "<<reference.c_str()<<"\n";
63 22 : bool do_read=true;
64 992 : while (do_read) {
65 970 : PDB mypdb;
66 1918 : RMSD mymsd;
67 970 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
68 970 : if(do_read) {
69 948 : nframes++;
70 948 : if(mypdb.getAtomNumbers().size()==0) error("number of atoms in a frame should be more than zero");
71 948 : unsigned nat=mypdb.getAtomNumbers().size();
72 948 : if(nat!=mypdb.getAtomNumbers().size()) error("frames should have the same number of atoms");
73 948 : if(aaa.empty()) aaa=mypdb.getAtomNumbers();
74 948 : if(aaa!=mypdb.getAtomNumbers()) error("frames should contain same atoms in same order");
75 948 : log<<"Found PDB: "<<nframes<<" containing "<<mypdb.getAtomNumbers().size()<<" atoms\n";
76 948 : pdbv.push_back(mypdb);
77 : // requestAtoms(mypdb.getAtomNumbers()); // is done in non base classes
78 948 : derivs_s.resize(mypdb.getAtomNumbers().size());
79 948 : derivs_z.resize(mypdb.getAtomNumbers().size());
80 948 : mymsd.set(mypdb,"OPTIMAL");
81 948 : msdv.push_back(mymsd); // the vector that stores the frames
82 : //log<<mypdb;
83 22 : } else {break ;}
84 948 : }
85 22 : fclose (fp);
86 22 : log<<"Found TOTAL "<<nframes<< " PDB in the file "<<reference.c_str()<<" \n";
87 22 : if(nframes==0) error("at least one frame expected");
88 : }
89 22 : if(neigh_stride>0 || neigh_size>0) {
90 13 : if(neigh_size>int(nframes)) {
91 0 : log.printf(" List size required ( %d ) is too large: resizing to the maximum number of frames required: %u \n",neigh_size,nframes);
92 0 : neigh_size=nframes;
93 : }
94 13 : log.printf(" Neighbor list enabled: \n");
95 13 : log.printf(" size : %d elements\n",neigh_size);
96 13 : log.printf(" stride : %d timesteps \n",neigh_stride);
97 : } else {
98 9 : log.printf(" Neighbor list NOT enabled \n");
99 22 : }
100 :
101 22 : }
102 :
103 10036 : void PathMSDBase::calculate() {
104 :
105 10036 : if(neigh_size>0 && getExchangeStep()) error("Neighbor lists for this collective variable are not compatible with replica exchange, sorry for that!");
106 :
107 : //log.printf("NOW CALCULATE! \n");
108 :
109 :
110 : // resize the list to full
111 10036 : if(imgVec.empty()) { // this is the signal that means: recalculate all
112 6616 : imgVec.resize(nframes);
113 284512 : for(unsigned i=0; i<nframes; i++) {
114 277896 : imgVec[i].property=indexvec[i];
115 277896 : imgVec[i].index=i;
116 : }
117 : }
118 :
119 : // THIS IS THE HEAVY PART (RMSD STUFF)
120 10036 : unsigned stride=comm.Get_size();
121 10036 : unsigned rank=comm.Get_rank();
122 10036 : unsigned nat=pdbv[0].size();
123 10036 : plumed_assert(nat>0);
124 10036 : plumed_assert(nframes>0);
125 10036 : plumed_assert(imgVec.size()>0);
126 :
127 10036 : std::vector<double> tmp_distances(imgVec.size(),0.0);
128 20072 : std::vector<Vector> tmp_derivs;
129 : // this array is a merge of all tmp_derivs, so as to allow a single comm.Sum below
130 20072 : std::vector<Vector> tmp_derivs2(imgVec.size()*nat);
131 :
132 : // if imgVec.size() is less than nframes, it means that only some msd will be calculated
133 295180 : for(unsigned i=rank; i<imgVec.size(); i+=stride) {
134 : // store temporary local results
135 285144 : tmp_distances[i]=msdv[imgVec[i].index].calculate(getPositions(),tmp_derivs,true);
136 285144 : plumed_assert(tmp_derivs.size()==nat);
137 285144 : for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
138 : }
139 : // reduce over all processors
140 10036 : comm.Sum(tmp_distances);
141 10036 : comm.Sum(tmp_derivs2);
142 : // assign imgVec[i].distance and imgVec[i].distder
143 432772 : for(unsigned i=0; i<imgVec.size(); i++) {
144 422736 : imgVec[i].distance=tmp_distances[i];
145 422736 : imgVec[i].distder.assign(&tmp_derivs2[i*nat],nat+&tmp_derivs2[i*nat]);
146 : }
147 :
148 : // END OF THE HEAVY PART
149 :
150 20072 : vector<Value*> val_s_path;
151 10036 : if(labels.size()>0) {
152 4914 : for(unsigned i=0; i<labels.size(); i++) { val_s_path.push_back(getPntrToComponent(labels[i].c_str()));}
153 : } else {
154 5122 : val_s_path.push_back(getPntrToComponent("sss"));
155 : }
156 10036 : Value* val_z_path=getPntrToComponent("zzz");
157 :
158 20072 : vector<double> s_path(val_s_path.size()); for(unsigned i=0; i<s_path.size(); i++)s_path[i]=0.;
159 10036 : double partition=0.;
160 : double tmp;
161 :
162 : // clean vector
163 10036 : for(unsigned i=0; i< derivs_z.size(); i++) {derivs_z[i].zero();}
164 :
165 : typedef vector< class ImagePath >::iterator imgiter;
166 432772 : for(imgiter it=imgVec.begin(); it!=imgVec.end(); ++it) {
167 422736 : (*it).similarity=exp(-lambda*((*it).distance));
168 : //log<<"DISTANCE "<<(*it).distance<<"\n";
169 1051860 : for(unsigned i=0; i<s_path.size(); i++) {
170 629124 : s_path[i]+=((*it).property[i])*(*it).similarity;
171 : }
172 422736 : partition+=(*it).similarity;
173 : }
174 10036 : for(unsigned i=0; i<s_path.size(); i++) { s_path[i]/=partition; val_s_path[i]->set(s_path[i]) ;}
175 10036 : val_z_path->set(-(1./lambda)*std::log(partition));
176 24986 : for(unsigned j=0; j<s_path.size(); j++) {
177 : // clean up
178 14950 : for(unsigned i=0; i< derivs_s.size(); i++) {derivs_s[i].zero();}
179 : // do the derivative
180 644074 : for(imgiter it=imgVec.begin(); it!=imgVec.end(); ++it) {
181 629124 : double expval=(*it).similarity;
182 629124 : tmp=lambda*expval*(s_path[j]-(*it).property[j])/partition;
183 629124 : for(unsigned i=0; i< derivs_s.size(); i++) { derivs_s[i]+=tmp*(*it).distder[i] ;}
184 629124 : if(j==0) {for(unsigned i=0; i< derivs_z.size(); i++) { derivs_z[i]+=(*it).distder[i]*expval/partition;}}
185 : }
186 209147 : for(unsigned i=0; i< derivs_s.size(); i++) {
187 194197 : setAtomsDerivatives (val_s_path[j],i,derivs_s[i]);
188 194197 : if(j==0) {setAtomsDerivatives (val_z_path,i,derivs_z[i]);}
189 : }
190 : }
191 10036 : for(unsigned i=0; i<val_s_path.size(); ++i) setBoxDerivativesNoPbc(val_s_path[i]);
192 10036 : setBoxDerivativesNoPbc(val_z_path);
193 : //
194 : // here set next round neighbors
195 : //
196 10036 : if (neigh_size>0) {
197 : //if( int(getStep())%int(neigh_stride/getTimeStep())==0 ){
198 : // enforce consistency: the stride is in time steps
199 6607 : if( int(getStep())%int(neigh_stride)==0 ) {
200 :
201 : // next round do it all:empty the vector
202 6607 : imgVec.clear();
203 : }
204 : // time to analyze the results:
205 6607 : if(imgVec.size()==nframes) {
206 : //sort by msd
207 0 : sort(imgVec.begin(), imgVec.end(), imgOrderByDist());
208 : //resize
209 0 : imgVec.resize(neigh_size);
210 : }
211 10036 : }
212 : //log.printf("CALCULATION DONE! \n");
213 10036 : }
214 :
215 : }
216 :
217 2523 : }
|
