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