Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "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 27 : void PathMSDBase::registerKeywords(Keywords& keys) {
38 27 : Colvar::registerKeywords(keys);
39 108 : keys.add("compulsory","LAMBDA","the lambda parameter is needed for smoothing, is in the units of plumed");
40 108 : keys.add("compulsory","REFERENCE","the pdb is needed to provide the various milestones");
41 108 : keys.add("optional","NEIGH_SIZE","size of the neighbor list");
42 108 : keys.add("optional","NEIGH_STRIDE","how often the neighbor list needs to be calculated in time units");
43 108 : keys.add("optional", "EPSILON", "(default=-1) the maximum distance between the close and the current structure, the positive value turn on the close structure method");
44 108 : keys.add("optional", "LOG_CLOSE", "(default=0) value 1 enables logging regarding the close structure");
45 108 : keys.add("optional", "DEBUG_CLOSE", "(default=0) value 1 enables extensive debugging info regarding the close structure, the simulation will run much slower");
46 27 : }
47 :
48 25 : PathMSDBase::PathMSDBase(const ActionOptions&ao):
49 : PLUMED_COLVAR_INIT(ao),
50 : nopbc(false),
51 : neigh_size(-1),
52 : neigh_stride(-1),
53 : epsilonClose(-1),
54 : debugClose(0),
55 : logClose(0),
56 : computeRefClose(false),
57 100 : nframes(0)
58 : {
59 50 : parse("LAMBDA",lambda);
60 50 : parse("NEIGH_SIZE",neigh_size);
61 50 : parse("NEIGH_STRIDE",neigh_stride);
62 50 : parse("REFERENCE",reference);
63 50 : parse("EPSILON", epsilonClose);
64 50 : parse("LOG_CLOSE", logClose);
65 50 : parse("DEBUG_CLOSE", debugClose);
66 50 : parseFlag("NOPBC",nopbc);
67 :
68 : // open the file
69 25 : FILE* fp=fopen(reference.c_str(),"r");
70 : std::vector<AtomNumber> aaa;
71 25 : if (fp!=NULL)
72 : {
73 25 : log<<"Opening reference file "<<reference.c_str()<<"\n";
74 : bool do_read=true;
75 1132 : while (do_read) {
76 1107 : PDB mypdb;
77 2189 : RMSD mymsd;
78 2214 : do_read=mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
79 1107 : if(do_read) {
80 1082 : nframes++;
81 2164 : if(mypdb.getAtomNumbers().size()==0) error("number of atoms in a frame should be more than zero");
82 1082 : unsigned nat=mypdb.getAtomNumbers().size();
83 2164 : if(nat!=mypdb.getAtomNumbers().size()) error("frames should have the same number of atoms");
84 1082 : if(aaa.empty()) {
85 25 : aaa=mypdb.getAtomNumbers();
86 25 : log.printf(" found %z atoms in input \n",aaa.size());
87 25 : log.printf(" with indices : ");
88 667 : for(unsigned i=0; i<aaa.size(); ++i) {
89 321 : if(i%25==0) log<<"\n";
90 321 : log.printf("%d ",aaa[i].serial());
91 : }
92 25 : log.printf("\n");
93 : }
94 2164 : if(aaa!=mypdb.getAtomNumbers()) error("frames should contain same atoms in same order");
95 2164 : log<<"Found PDB: "<<nframes<<" containing "<<mypdb.getAtomNumbers().size()<<" atoms\n";
96 1082 : pdbv.push_back(mypdb);
97 2164 : derivs_s.resize(mypdb.getAtomNumbers().size());
98 2164 : derivs_z.resize(mypdb.getAtomNumbers().size());
99 2164 : mymsd.set(mypdb,"OPTIMAL");
100 1082 : msdv.push_back(mymsd); // the vector that stores the frames
101 : } else {break ;}
102 1082 : }
103 25 : fclose (fp);
104 25 : log<<"Found TOTAL "<<nframes<< " PDB in the file "<<reference.c_str()<<" \n";
105 25 : if(nframes==0) error("at least one frame expected");
106 : //set up rmsdRefClose, initialize it to the first structure loaded from reference file
107 75 : rmsdPosClose.set(pdbv[0], "OPTIMAL");
108 25 : firstPosClose = true;
109 : }
110 25 : if(neigh_stride>0 || neigh_size>0) {
111 14 : if(neigh_size>int(nframes)) {
112 0 : log.printf(" List size required ( %d ) is too large: resizing to the maximum number of frames required: %u \n",neigh_size,nframes);
113 0 : neigh_size=nframes;
114 : }
115 14 : log.printf(" Neighbor list enabled: \n");
116 14 : log.printf(" size : %d elements\n",neigh_size);
117 14 : log.printf(" stride : %d timesteps \n",neigh_stride);
118 : } else {
119 11 : log.printf(" Neighbor list NOT enabled \n");
120 : }
121 25 : if (epsilonClose > 0) {
122 2 : log.printf(" Computing with the close structure, epsilon = %lf\n", epsilonClose);
123 6 : log << " Bibliography " << plumed.cite("Pazurikova J, Krenek A, Spiwok V, Simkova M J. Chem. Phys. 146, 115101 (2017)") << "\n";
124 : }
125 : else {
126 23 : debugClose = 0;
127 23 : logClose = 0;
128 : }
129 25 : if (debugClose)
130 2 : log.printf(" Extensive debug info regarding close structure turned on\n");
131 :
132 25 : rotationRefClose.resize(nframes);
133 50 : savedIndices = vector<unsigned>(nframes);
134 :
135 25 : if(nopbc) log.printf(" without periodic boundary conditions\n");
136 24 : else log.printf(" using periodic boundary conditions\n");
137 :
138 25 : }
139 :
140 75 : PathMSDBase::~PathMSDBase() {
141 25 : }
142 :
143 11179 : void PathMSDBase::calculate() {
144 :
145 11179 : if(neigh_size>0 && getExchangeStep()) error("Neighbor lists for this collective variable are not compatible with replica exchange, sorry for that!");
146 :
147 : //log.printf("NOW CALCULATE! \n");
148 :
149 11179 : if(!nopbc) makeWhole();
150 :
151 :
152 : // resize the list to full
153 11179 : if(imgVec.empty()) { // this is the signal that means: recalculate all
154 7164 : imgVec.resize(nframes);
155 : #pragma omp simd
156 300920 : for(unsigned i=0; i<nframes; i++) {
157 601840 : imgVec[i].property=indexvec[i];
158 300920 : imgVec[i].index=i;
159 : }
160 : }
161 :
162 : // THIS IS THE HEAVY PART (RMSD STUFF)
163 11179 : unsigned stride=comm.Get_size();
164 11179 : unsigned rank=comm.Get_rank();
165 11179 : unsigned nat=pdbv[0].size();
166 11179 : plumed_assert(nat>0);
167 11179 : plumed_assert(nframes>0);
168 11179 : plumed_assert(imgVec.size()>0);
169 :
170 11179 : std::vector<Tensor> tmp_rotationRefClose(nframes);
171 :
172 11179 : if (epsilonClose > 0) {
173 : //compute rmsd between positions and close structure, save rotation matrix, drotation_drr01
174 1092 : double posclose = rmsdPosClose.calc_Rot_DRotDRr01(getPositions(), rotationPosClose, drotationPosCloseDrr01, true);
175 : //if we compute for the first time or the existing close structure is too far from current structure
176 1092 : if (firstPosClose || (posclose > epsilonClose)) {
177 : //set the current structure as close one for a few next steps
178 16 : if (logClose)
179 16 : log << "PLUMED_CLOSE: new close structure, rmsd pos close " << posclose << "\n";
180 16 : rmsdPosClose.clear();
181 16 : rmsdPosClose.setReference(getPositions());
182 : //as this is a new close structure, we need to save the rotation matrices fitted to the reference structures
183 : // and we need to accurately recalculate for all reference structures
184 16 : computeRefClose = true;
185 16 : imgVec.resize(nframes);
186 672 : for(unsigned i=0; i<nframes; i++) {
187 1344 : imgVec[i].property=indexvec[i];
188 672 : imgVec[i].index=i;
189 : }
190 16 : firstPosClose = false;
191 : }
192 : else {
193 : //the current structure is pretty close to the close structure, so we use saved rotation matrices to decrease the complexity of rmsd comuptation
194 1076 : if (debugClose)
195 1076 : log << "PLUMED-CLOSE: old close structure, rmsd pos close " << posclose << "\n";
196 1076 : computeRefClose = false;
197 : }
198 : }
199 :
200 22358 : std::vector<double> tmp_distances(imgVec.size(),0.0);
201 : std::vector<Vector> tmp_derivs;
202 : // this array is a merge of all tmp_derivs, so as to allow a single comm.Sum below
203 11179 : std::vector<Vector> tmp_derivs2(imgVec.size()*nat);
204 :
205 : // if imgVec.size() is less than nframes, it means that only some msd will be calculated
206 11179 : if (epsilonClose > 0) {
207 1092 : if (computeRefClose) {
208 : //recompute rotation matrices accurately
209 1360 : for(unsigned i=rank; i<imgVec.size(); i+=stride) {
210 2016 : tmp_distances[i] = msdv[imgVec[i].index].calc_Rot(getPositions(), tmp_derivs, tmp_rotationRefClose[imgVec[i].index], true);
211 672 : plumed_assert(tmp_derivs.size()==nat);
212 : #pragma omp simd
213 26208 : for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
214 : }
215 : }
216 : else {
217 : //approximate distance with saved rotation matrices
218 92536 : for(unsigned i=rank; i<imgVec.size(); i+=stride) {
219 180768 : tmp_distances[i] = msdv[imgVec[i].index].calculateWithCloseStructure(getPositions(), tmp_derivs, rotationPosClose, rotationRefClose[imgVec[i].index], drotationPosCloseDrr01, true);
220 45192 : plumed_assert(tmp_derivs.size()==nat);
221 : #pragma omp simd
222 1762488 : for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
223 45192 : if (debugClose) {
224 90384 : double withclose = tmp_distances[i];
225 45192 : RMSD opt;
226 90384 : opt.setType("OPTIMAL");
227 180768 : opt.setReference(msdv[imgVec[i].index].getReference());
228 : vector<Vector> ders;
229 45192 : double withoutclose = opt.calculate(getPositions(), ders, true);
230 45192 : float difference = fabs(withoutclose-withclose);
231 45192 : log<<"PLUMED-CLOSE: difference original "<<withoutclose;
232 135576 : log<<" - with close "<<withclose<<" = "<<difference<<", step "<<getStep()<<", i "<<i<<" imgVec[i].index "<<imgVec[i].index<<"\n";
233 : }
234 : }
235 : }
236 : }
237 : else {
238 : // store temporary local results
239 585475 : for(unsigned i=rank; i<imgVec.size(); i+=stride) {
240 863082 : tmp_distances[i]=msdv[imgVec[i].index].calculate(getPositions(),tmp_derivs,true);
241 287694 : plumed_assert(tmp_derivs.size()==nat);
242 : #pragma omp simd
243 11189466 : for(unsigned j=0; j<nat; j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
244 : }
245 : }
246 :
247 : // reduce over all processors
248 11179 : comm.Sum(tmp_distances);
249 11179 : comm.Sum(tmp_derivs2);
250 11179 : if (epsilonClose > 0 && computeRefClose) {
251 16 : comm.Sum(tmp_rotationRefClose);
252 672 : for (unsigned i=0; i<nframes; i++) {
253 1344 : rotationRefClose[i] = tmp_rotationRefClose[i];
254 : }
255 : }
256 : // assign imgVec[i].distance and imgVec[i].distder
257 953479 : for(unsigned i=0; i<imgVec.size(); i++) {
258 471150 : imgVec[i].distance=tmp_distances[i];
259 1413450 : imgVec[i].distder.assign(&tmp_derivs2[i*nat],nat+&tmp_derivs2[i*nat]);
260 : }
261 :
262 : // END OF THE HEAVY PART
263 :
264 : vector<Value*> val_s_path;
265 11179 : if(labels.size()>0) {
266 54054 : for(unsigned i=0; i<labels.size(); i++) { val_s_path.push_back(getPntrToComponent(labels[i].c_str()));}
267 : } else {
268 15519 : val_s_path.push_back(getPntrToComponent("sss"));
269 : }
270 22358 : Value* val_z_path=getPntrToComponent("zzz");
271 :
272 56728 : vector<double> s_path(val_s_path.size()); for(unsigned i=0; i<s_path.size(); i++)s_path[i]=0.;
273 : double partition=0.;
274 : double tmp;
275 :
276 : // clean vector
277 301425 : for(unsigned i=0; i< derivs_z.size(); i++) {derivs_z[i].zero();}
278 :
279 493508 : for(auto & it : imgVec) {
280 471150 : it.similarity=exp(-lambda*(it.distance));
281 : //log<<"DISTANCE "<<(*it).distance<<"\n";
282 2389104 : for(unsigned i=0; i<s_path.size(); i++) {
283 1446804 : s_path[i]+=(it.property[i])*it.similarity;
284 : }
285 471150 : partition+=it.similarity;
286 : }
287 97104 : for(unsigned i=0; i<s_path.size(); i++) { s_path[i]/=partition; val_s_path[i]->set(s_path[i]) ;}
288 11179 : val_z_path->set(-(1./lambda)*std::log(partition));
289 45549 : for(unsigned j=0; j<s_path.size(); j++) {
290 : // clean up
291 : #pragma omp simd
292 463587 : for(unsigned i=0; i< derivs_s.size(); i++) {derivs_s[i].zero();}
293 : // do the derivative
294 757772 : for(const auto & it : imgVec) {
295 723402 : double expval=it.similarity;
296 2170206 : tmp=lambda*expval*(s_path[j]-it.property[j])/partition;
297 : #pragma omp simd
298 29628882 : for(unsigned i=0; i< derivs_s.size(); i++) { derivs_s[i]+=tmp*it.distder[i] ;}
299 723402 : if(j==0) {
300 : #pragma omp simd
301 18815400 : for(unsigned i=0; i< derivs_z.size(); i++) { derivs_z[i]+=it.distder[i]*expval/partition;}
302 : }
303 : }
304 463587 : for(unsigned i=0; i< derivs_s.size(); i++) {
305 446402 : setAtomsDerivatives (val_s_path[j],i,derivs_s[i]);
306 368324 : if(j==0) {setAtomsDerivatives (val_z_path,i,derivs_z[i]);}
307 : }
308 : }
309 45549 : for(unsigned i=0; i<val_s_path.size(); ++i) setBoxDerivativesNoPbc(val_s_path[i]);
310 11179 : setBoxDerivativesNoPbc(val_z_path);
311 : //
312 : // here set next round neighbors
313 : //
314 11179 : if (neigh_size>0) {
315 : //if( int(getStep())%int(neigh_stride/getTimeStep())==0 ){
316 : // enforce consistency: the stride is in time steps
317 7153 : if( int(getStep())%int(neigh_stride)==0 ) {
318 :
319 : // next round do it all:empty the vector
320 : imgVec.clear();
321 : }
322 : // time to analyze the results:
323 7153 : if(imgVec.size()==nframes) {
324 : //sort by msd
325 0 : sort(imgVec.begin(), imgVec.end(), imgOrderByDist());
326 : //resize
327 0 : imgVec.resize(neigh_size);
328 : }
329 : }
330 : //log.printf("CALCULATION DONE! \n");
331 11179 : }
332 :
333 : }
334 :
335 5874 : }
|
