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