Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017 of Pipolo Silvio and Fabio Pietrucci.
3 :
4 : The piv module is free software: you can redistribute it and/or modify
5 : it under the terms of the GNU Lesser General Public License as published by
6 : the Free Software Foundation, either version 3 of the License, or
7 : (at your option) any later version.
8 :
9 : The piv module is distributed in the hope that it will be useful,
10 : but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : GNU Lesser General Public License for more details.
13 :
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : #include "colvar/Colvar.h"
18 : #include "colvar/ActionRegister.h"
19 : #include "core/PlumedMain.h"
20 : #include "core/ActionWithVirtualAtom.h"
21 : #include "tools/NeighborList.h"
22 : #include "tools/SwitchingFunction.h"
23 : #include "tools/PDB.h"
24 : #include "tools/Pbc.h"
25 : #include "tools/Stopwatch.h"
26 : #include "core/ActionSet.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 : #include <iostream>
31 :
32 : using namespace std;
33 :
34 : namespace PLMD {
35 : namespace piv {
36 :
37 : //+PLUMEDOC PIVMOD_COLVAR PIV
38 : /*
39 : Calculates the PIV-distance.
40 :
41 : PIV distance is the squared Cartesian distance between the PIV \cite gallet2013structural \cite pipolo2017navigating
42 : associated to the configuration of the system during the dynamics and a reference configuration provided
43 : as input (PDB file format).
44 : PIV can be used together with \ref FUNCPATHMSD to define a path in the PIV space.
45 :
46 : \par Examples
47 :
48 : The following example calculates PIV-distances from three reference configurations in Ref1.pdb, Ref2.pdb and Ref3.pdb
49 : and prints the results in a file named colvar.
50 : Three atoms (PIVATOMS=3) with names (pdb file) A B and C are used to construct the PIV and all PIV blocks (AA, BB, CC, AB, AC, BC) are considered.
51 : SFACTOR is a scaling factor that multiplies the contribution to the PIV-distance given by the single PIV block.
52 : NLIST sets the use of neighbor lists for calculating atom-atom distances.
53 : The SWITCH keyword specifies the parameters of the switching function that transforms atom-atom distances.
54 : SORT=1 means that the PIV block elements are sorted (SORT=0 no sorting.)
55 : Values for SORT, SFACTOR and the neighbor list parameters have to be specified for each block.
56 : The order is the following: AA,BB,CC,AB,AC,BC. If ONLYDIRECT (ONLYCROSS) is used the order is AA,BB,CC (AB,AC,BC).
57 : The sorting operation within each PIV block is performed using the counting sort algorithm, PRECISION specifies the size of the counting array.
58 :
59 : \plumedfile
60 : PIV ...
61 : LABEL=Pivd1
62 : PRECISION=1000
63 : NLIST
64 : REF_FILE=Ref1.pdb
65 : PIVATOMS=3
66 : ATOMTYPES=A,B,C
67 : SFACTOR=0.3,0.5,1.0,0.2,0.2,0.2
68 : SORT=1,1,1,1,1,1
69 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
70 : SWITCH2={RATIONAL R_0=0.4 MM=10 NN=5}
71 : SWITCH3={RATIONAL R_0=0.4 MM=10 NN=5}
72 : SWITCH4={RATIONAL R_0=0.5 MM=12 NN=6}
73 : SWITCH5={RATIONAL R_0=0.5 MM=12 NN=6}
74 : SWITCH6={RATIONAL R_0=0.5 MM=12 NN=6}
75 : NL_CUTOFF=0.8,0.6,0.6,0.7,0.7,0.7
76 : NL_STRIDE=10,10,10,10,10,10
77 : NL_SKIN=0.1,0.1,0.1,0.1,0.1,0.1
78 : ... PIV
79 : PIV ...
80 : LABEL=Pivd2
81 : PRECISION=1000
82 : NLIST
83 : REF_FILE=Ref2.pdb
84 : PIVATOMS=3
85 : ATOMTYPES=A,B,C
86 : SFACTOR=0.3,0.5,1.0,0.2,0.2,0.2
87 : SORT=1,1,1,1,1,1
88 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
89 : SWITCH2={RATIONAL R_0=0.4 MM=10 NN=5}
90 : SWITCH3={RATIONAL R_0=0.4 MM=10 NN=5}
91 : SWITCH4={RATIONAL R_0=0.5 MM=12 NN=6}
92 : SWITCH5={RATIONAL R_0=0.5 MM=12 NN=6}
93 : SWITCH6={RATIONAL R_0=0.5 MM=12 NN=6}
94 : NL_CUTOFF=0.8,0.6,0.6,0.7,0.7,0.7
95 : NL_STRIDE=10,10,10,10,10,10
96 : NL_SKIN=0.1,0.1,0.1,0.1,0.1,0.1
97 : ... PIV
98 : PIV ...
99 : LABEL=Pivd3
100 : PRECISION=1000
101 : NLIST
102 : REF_FILE=Ref3.pdb
103 : PIVATOMS=3
104 : ATOMTYPES=A,B,C
105 : SFACTOR=0.3,0.5,1.0,0.2,0.2,0.2
106 : SORT=1,1,1,1,1,1
107 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
108 : SWITCH2={RATIONAL R_0=0.4 MM=10 NN=5}
109 : SWITCH3={RATIONAL R_0=0.4 MM=10 NN=5}
110 : SWITCH4={RATIONAL R_0=0.5 MM=12 NN=6}
111 : SWITCH5={RATIONAL R_0=0.5 MM=12 NN=6}
112 : SWITCH6={RATIONAL R_0=0.5 MM=12 NN=6}
113 : NL_CUTOFF=0.8,0.6,0.6,0.7,0.7,0.7
114 : NL_STRIDE=10,10,10,10,10,10
115 : NL_SKIN=0.1,0.1,0.1,0.1,0.1,0.1
116 : ... PIV
117 :
118 : PRINT ARG=Pivd1,Pivd2,Pivd3 FILE=colvar
119 : \endplumedfile
120 :
121 : WARNING:
122 : Both the "CRYST" and "ATOM" lines of the PDB files must conform precisely to the official pdb format, including the width of each alphanumerical field:
123 :
124 : \verbatim
125 : CRYST1 31.028 36.957 23.143 89.93 92.31 89.99 P 1 1
126 : ATOM 1 OW1 wate 1 15.630 19.750 1.520 1.00 0.00
127 : \endverbatim
128 :
129 : In each pdb frame, atoms must be numbered in the same order and with the same element symbol as in the input of the MD program.
130 :
131 : The following example calculates the PIV-distances from two reference configurations Ref1.pdb and Ref2.pdb
132 : and uses PIV-distances to define a Path Collective Variable (\ref FUNCPATHMSD) with only two references (Ref1.pdb and Ref2.pdb).
133 : With the VOLUME keyword one scales the atom-atom distances by the cubic root of the ratio between the specified value and the box volume of the initial step of the trajectory file.
134 :
135 : \plumedfile
136 : PIV ...
137 : LABEL=c1
138 : PRECISION=1000
139 : VOLUME=12.15
140 : NLIST
141 : REF_FILE=Ref1.pdb
142 : PIVATOMS=2
143 : ATOMTYPES=A,B
144 : ONLYDIRECT
145 : SFACTOR=1.0,0.2
146 : SORT=1,1
147 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
148 : SWITCH2={RATIONAL R_0=0.5 MM=10 NN=5}
149 : NL_CUTOFF=1.2,1.2
150 : NL_STRIDE=10,10
151 : NL_SKIN=0.1,0.1
152 : ... PIV
153 : PIV ...
154 : LABEL=c2
155 : PRECISION=1000
156 : VOLUME=12.15
157 : NLIST
158 : REF_FILE=Ref2.pdb
159 : PIVATOMS=2
160 : ATOMTYPES=A,B
161 : ONLYDIRECT
162 : SFACTOR=1.0,0.2
163 : SORT=1,1
164 : SWITCH1={RATIONAL R_0=0.6 MM=12 NN=4}
165 : SWITCH2={RATIONAL R_0=0.4 MM=10 NN=5}
166 : NL_CUTOFF=1.2,1.2
167 : NL_STRIDE=10,10
168 : NL_SKIN=0.1,0.1
169 : ... PIV
170 :
171 : p1: FUNCPATHMSD ARG=c1,c2 LAMBDA=0.180338
172 : METAD ARG=p1.s,p1.z SIGMA=0.01,0.2 HEIGHT=0.8 PACE=500 LABEL=res
173 : PRINT ARG=c1,c2,p1.s,p1.z,res.bias STRIDE=500 FILE=colvar FMT=%15.6f
174 : \endplumedfile
175 :
176 : When using PIV please cite \cite pipolo2017navigating .
177 :
178 : (See also \ref PRINT)
179 :
180 : */
181 : //+ENDPLUMEDOC
182 :
183 : class PIV : public Colvar {
184 : private:
185 : bool pbc, serial, timer;
186 : ForwardDecl<Stopwatch> stopwatch_fwd;
187 : Stopwatch& stopwatch=*stopwatch_fwd;
188 : int updatePIV;
189 : size_t Nprec;
190 : unsigned Natm,Nlist,NLsize;
191 : double Fvol,Vol0,m_PIVdistance;
192 : std::string ref_file;
193 : std::unique_ptr<NeighborList> nlall;
194 : std::vector<SwitchingFunction> sfs;
195 : std::vector<std:: vector<double> > rPIV;
196 : std::vector<double> scaling,r00;
197 : std::vector<double> nl_skin;
198 : std::vector<double> fmass;
199 : std::vector<bool> dosort;
200 : std::vector<Vector> compos;
201 : std::vector<string> sw;
202 : std::vector<std::unique_ptr<NeighborList>> nl;
203 : std::vector<std::unique_ptr<NeighborList>> nlcom;
204 : std::vector<Vector> m_deriv;
205 : Tensor m_virial;
206 : bool Svol,cross,direct,doneigh,test,CompDer,com;
207 :
208 : /// Local structure, used to store data that should be
209 : /// shared across multiple PIV instances
210 6 : struct SharedData {
211 : int prev_stp=-1;
212 : int init_stp=1;
213 : std:: vector<std:: vector<Vector> > prev_pos;
214 : std:: vector<std:: vector<double> > cPIV;
215 : std:: vector<std:: vector<int> > Atom0;
216 : std:: vector<std:: vector<int> > Atom1;
217 : };
218 : /// Owning pointer. Will only be allocated by the first PIV instance
219 : std::unique_ptr<SharedData> sharedData_unique;
220 : /// Raw pointer. Will have the same value for all PIV instances
221 : SharedData* sharedData=nullptr;
222 :
223 : public:
224 : static void registerKeywords( Keywords& keys );
225 : explicit PIV(const ActionOptions&);
226 : // active methods:
227 : virtual void calculate();
228 0 : void checkFieldsAllowed() {}
229 : };
230 :
231 13809 : PLUMED_REGISTER_ACTION(PIV,"PIV")
232 :
233 16 : void PIV::registerKeywords( Keywords& keys ) {
234 16 : Colvar::registerKeywords( keys );
235 32 : keys.add("numbered","SWITCH","The switching functions parameter."
236 : "You should specify a Switching function for all PIV blocks."
237 : "Details of the various switching "
238 : "functions you can use are provided on \\ref switchingfunction.");
239 32 : keys.add("compulsory","PRECISION","the precision for approximating reals with integers in sorting.");
240 32 : keys.add("compulsory","REF_FILE","PDB file name that contains the \\f$i\\f$th reference structure.");
241 32 : keys.add("compulsory","PIVATOMS","Number of atoms to use for PIV.");
242 32 : keys.add("compulsory","SORT","Whether to sort or not the PIV block.");
243 32 : keys.add("compulsory","ATOMTYPES","The atom types to use for PIV.");
244 32 : keys.add("optional","SFACTOR","Scale the PIV-distance by such block-specific factor");
245 32 : keys.add("optional","VOLUME","Scale atom-atom distances by the cubic root of the cell volume. The input volume is used to scale the R_0 value of the switching function. ");
246 32 : keys.add("optional","UPDATEPIV","Frequency (in steps) at which the PIV is updated.");
247 32 : keys.addFlag("TEST",false,"Print the actual and reference PIV and exit");
248 32 : keys.addFlag("COM",false,"Use centers of mass of groups of atoms instead of atoms as specified in the Pdb file");
249 32 : keys.addFlag("ONLYCROSS",false,"Use only cross-terms (A-B, A-C, B-C, ...) in PIV");
250 32 : keys.addFlag("ONLYDIRECT",false,"Use only direct-terms (A-A, B-B, C-C, ...) in PIV");
251 32 : keys.addFlag("DERIVATIVES",false,"Activate the calculation of the PIV for every class (needed for numerical derivatives).");
252 32 : keys.addFlag("NLIST",false,"Use a neighbor list for distance calculations.");
253 32 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
254 32 : keys.addFlag("TIMER",false,"Perform timing analysis on heavy loops.");
255 32 : keys.add("optional","NL_CUTOFF","Neighbor lists cutoff.");
256 32 : keys.add("optional","NL_STRIDE","Update neighbor lists every NL_STRIDE steps.");
257 32 : keys.add("optional","NL_SKIN","The maximum atom displacement tolerated for the neighbor lists update.");
258 32 : keys.reset_style("SWITCH","compulsory");
259 16 : }
260 :
261 12 : PIV::PIV(const ActionOptions&ao):
262 : PLUMED_COLVAR_INIT(ao),
263 12 : pbc(true),
264 12 : serial(false),
265 12 : timer(false),
266 12 : updatePIV(1),
267 12 : Nprec(1000),
268 12 : Natm(1),
269 12 : Nlist(1),
270 12 : NLsize(1),
271 12 : Fvol(1.),
272 12 : Vol0(0.),
273 12 : m_PIVdistance(0.),
274 12 : rPIV(std:: vector<std:: vector<double> >(Nlist)),
275 12 : scaling(std:: vector<double>(Nlist)),
276 12 : r00(std:: vector<double>(Nlist)),
277 12 : nl_skin(std:: vector<double>(Nlist)),
278 12 : fmass(std:: vector<double>(Nlist)),
279 12 : dosort(std:: vector<bool>(Nlist)),
280 12 : compos(std:: vector<Vector>(NLsize)),
281 12 : sw(std:: vector<string>(Nlist)),
282 12 : nl(Nlist),
283 12 : nlcom(NLsize),
284 12 : m_deriv(std:: vector<Vector>(1)),
285 12 : Svol(false),
286 12 : cross(true),
287 12 : direct(true),
288 12 : doneigh(false),
289 12 : test(false),
290 12 : CompDer(false),
291 24 : com(false) {
292 12 : log << "Starting PIV Constructor\n";
293 :
294 : {
295 : // look for another PIV instance previously allocated
296 12 : auto* previous=plumed.getActionSet().selectLatest<PIV*>(this);
297 :
298 : // Uncommenting the following line, it is possible to force
299 : // a separate object per instance of the PIB object.
300 : // Results are unaffected, but performance is worse.
301 : // I think this is the expected behavior. GB
302 : // previous=nullptr;
303 :
304 12 : if(!previous) {
305 : // if not found, allocate the shared data struct
306 12 : sharedData_unique=Tools::make_unique<SharedData>();
307 : // then set the raw pointer
308 6 : sharedData=sharedData_unique.get();
309 : } else {
310 : // if found, use the previous raw pointer
311 6 : sharedData=previous->sharedData;
312 6 : log << "(a previous PIV action was found)\n";
313 : }
314 : }
315 :
316 : // Precision on the real-to-integer transformation for the sorting
317 12 : parse("PRECISION",Nprec);
318 12 : if(Nprec<2) {
319 0 : error("Precision must be => 2");
320 : }
321 :
322 : // PBC
323 12 : bool nopbc=!pbc;
324 12 : parseFlag("NOPBC",nopbc);
325 12 : pbc=!nopbc;
326 12 : if(pbc) {
327 12 : log << "Using Periodic Boundary Conditions\n";
328 : } else {
329 0 : log << "Isolated System (NO PBC)\n";
330 : }
331 :
332 : // SERIAL/PARALLEL
333 12 : parseFlag("SERIAL",serial);
334 12 : if(serial) {
335 0 : log << "Serial PIV construction\n";
336 : } else {
337 12 : log << "Parallel PIV construction\n";
338 : }
339 :
340 : // Derivatives
341 12 : parseFlag("DERIVATIVES",CompDer);
342 12 : if(CompDer) {
343 6 : log << "Computing Derivatives\n";
344 : }
345 :
346 : // Timing
347 12 : parseFlag("TIMER",timer);
348 12 : if(timer) {
349 1 : log << "Timing analysis\n";
350 1 : stopwatch.start();
351 1 : stopwatch.pause();
352 : }
353 :
354 : // Test
355 12 : parseFlag("TEST",test);
356 :
357 : // UPDATEPIV
358 24 : if(keywords.exists("UPDATEPIV")) {
359 24 : parse("UPDATEPIV",updatePIV);
360 : }
361 :
362 : // Test
363 12 : parseFlag("COM",com);
364 12 : if(com) {
365 0 : log << "Building PIV using COMs\n";
366 : }
367 :
368 : // Volume Scaling
369 12 : parse("VOLUME",Vol0);
370 12 : if (Vol0>0) {
371 12 : Svol=true;
372 : }
373 :
374 : // PIV direct and cross blocks
375 12 : bool oc=false,od=false;
376 12 : parseFlag("ONLYCROSS",oc);
377 12 : parseFlag("ONLYDIRECT",od);
378 12 : if (oc&&od) {
379 0 : error("ONLYCROSS and ONLYDIRECT are incompatible options!");
380 : }
381 12 : if(oc) {
382 4 : direct=false;
383 4 : log << "Using only CROSS-PIV blocks\n";
384 : }
385 12 : if(od) {
386 4 : cross=false;
387 4 : log << "Using only DIRECT-PIV blocks\n";
388 : }
389 :
390 : // Atoms for PIV
391 12 : parse("PIVATOMS",Natm);
392 12 : std:: vector<string> atype(Natm);
393 12 : parseVector("ATOMTYPES",atype);
394 : //if(atype.size()!=getNumberOfArguments() && atype.size()!=0) error("not enough values for ATOMTYPES");
395 :
396 : // Reference PDB file
397 12 : parse("REF_FILE",ref_file);
398 12 : PDB mypdb;
399 12 : FILE* fp=fopen(ref_file.c_str(),"r");
400 12 : if (fp!=NULL) {
401 12 : log<<"Opening PDB file with reference frame: "<<ref_file.c_str()<<"\n";
402 24 : mypdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
403 12 : fclose (fp);
404 : } else {
405 0 : error("Error in reference PDB file");
406 : }
407 :
408 : // Build COM/Atom lists of AtomNumbers (this might be done in PBC.cpp)
409 : // Atomlist or Plist used to build pair lists
410 12 : std:: vector<std:: vector<AtomNumber> > Plist(Natm);
411 : // Atomlist used to build list of atoms for each COM
412 12 : std:: vector<std:: vector<AtomNumber> > comatm(1);
413 : // NLsize is the number of atoms in the pdb cell
414 12 : NLsize=mypdb.getAtomNumbers().size();
415 : // In the following P stands for Point (either an Atom or a COM)
416 : unsigned resnum=0;
417 : // Presind (array size: number of residues) contains the contains the residue number
418 : // this is because the residue numbers may not always be ordered from 1 to resnum
419 : std:: vector<unsigned> Presind;
420 : // Build Presind
421 19404 : for (unsigned i=0; i<mypdb.getAtomNumbers().size(); i++) {
422 19392 : unsigned rind=mypdb.getResidueNumber(mypdb.getAtomNumbers()[i]);
423 : bool oldres=false;
424 7705008 : for (unsigned j=0; j<Presind.size(); j++) {
425 7685616 : if(rind==Presind[j]) {
426 : oldres=true;
427 : }
428 : }
429 19392 : if(!oldres) {
430 4848 : Presind.push_back(rind);
431 : }
432 : }
433 12 : resnum=Presind.size();
434 :
435 : // Pind0 is the atom/COM used in Nlists (for COM Pind0 is the first atom in the pdb belonging to that COM)
436 : unsigned Pind0size;
437 12 : if(com) {
438 : Pind0size=resnum;
439 : } else {
440 12 : Pind0size=NLsize;
441 : }
442 12 : std:: vector<unsigned> Pind0(Pind0size);
443 : // If COM resize important arrays
444 12 : comatm.resize(NLsize);
445 12 : if(com) {
446 0 : nlcom.resize(NLsize);
447 0 : compos.resize(NLsize);
448 0 : fmass.resize(NLsize,0.);
449 : }
450 12 : log << "Total COM/Atoms: " << Natm*resnum << " \n";
451 : // Build lists of Atoms/COMs for NLists
452 : // comatm filled also for non_COM calculation for analysis purposes
453 36 : for (unsigned j=0; j<Natm; j++) {
454 : unsigned oind;
455 38808 : for (unsigned i=0; i<Pind0.size(); i++) {
456 38784 : Pind0[i]=0;
457 : }
458 38808 : for (unsigned i=0; i<mypdb.getAtomNumbers().size(); i++) {
459 : // Residue/Atom AtomNumber: used to build NL for COMS/Atoms pairs.
460 38784 : AtomNumber anum=mypdb.getAtomNumbers()[i];
461 : // ResidueName/Atomname associated to atom
462 38784 : string rname=mypdb.getResidueName(anum);
463 38784 : string aname=mypdb.getAtomName(anum);
464 : // Index associated to residue/atom: used to separate COM-lists
465 38784 : unsigned rind=mypdb.getResidueNumber(anum);
466 : unsigned aind=anum.index();
467 : // This builds lists for NL
468 : string Pname;
469 : unsigned Pind;
470 38784 : if(com) {
471 : Pname=rname;
472 0 : for(unsigned l=0; l<resnum; l++) {
473 0 : if(rind==Presind[l]) {
474 : Pind=l;
475 : }
476 : }
477 : } else {
478 : Pname=aname;
479 : Pind=aind;
480 : }
481 38784 : if(Pname==atype[j]) {
482 14544 : if(Pind0[Pind]==0) {
483 : // adding the atomnumber to the atom/COM list for pairs
484 14544 : Plist[j].push_back(anum);
485 14544 : Pind0[Pind]=aind+1;
486 : oind=Pind;
487 : }
488 : // adding the atomnumber to list of atoms for every COM/Atoms
489 14544 : comatm[Pind0[Pind]-1].push_back(anum);
490 : }
491 : }
492 : // Output Lists
493 24 : log << " Groups of type " << j << ": " << Plist[j].size() << " \n";
494 : string gname;
495 : unsigned gsize;
496 24 : if(com) {
497 0 : gname=mypdb.getResidueName(comatm[Pind0[oind]-1][0]);
498 0 : gsize=comatm[Pind0[oind]-1].size();
499 : } else {
500 48 : gname=mypdb.getAtomName(comatm[Pind0[oind]-1][0]);
501 : gsize=1;
502 : }
503 24 : log.printf(" %6s %3s %13s %10i %6s\n", "type ", gname.c_str()," containing ",gsize," atoms");
504 : }
505 :
506 : // This is to build the list with all the atoms
507 : std:: vector<AtomNumber> listall;
508 19404 : for (unsigned i=0; i<mypdb.getAtomNumbers().size(); i++) {
509 19392 : listall.push_back(mypdb.getAtomNumbers()[i]);
510 : }
511 :
512 : // PIV blocks and Neighbour Lists
513 12 : Nlist=0;
514 : // Direct adds the A-A ad B-B blocks (N)
515 12 : if(direct) {
516 8 : Nlist=Nlist+unsigned(Natm);
517 : }
518 : // Cross adds the A-B blocks (N*(N-1)/2)
519 12 : if(cross) {
520 8 : Nlist=Nlist+unsigned(double(Natm*(Natm-1))/2.);
521 : }
522 : // Resize vectors according to Nlist
523 12 : rPIV.resize(Nlist);
524 :
525 : // PIV scaled option
526 12 : scaling.resize(Nlist);
527 36 : for(unsigned j=0; j<Nlist; j++) {
528 24 : scaling[j]=1.;
529 : }
530 24 : if(keywords.exists("SFACTOR")) {
531 24 : parseVector("SFACTOR",scaling);
532 : //if(scaling.size()!=getNumberOfArguments() && scaling.size()!=0) error("not enough values for SFACTOR");
533 : }
534 : // Neighbour Lists option
535 12 : parseFlag("NLIST",doneigh);
536 12 : nl.resize(Nlist);
537 12 : nl_skin.resize(Nlist);
538 12 : if(doneigh) {
539 12 : std:: vector<double> nl_cut(Nlist,0.);
540 12 : std:: vector<int> nl_st(Nlist,0);
541 12 : parseVector("NL_CUTOFF",nl_cut);
542 : //if(nl_cut.size()!=getNumberOfArguments() && nl_cut.size()!=0) error("not enough values for NL_CUTOFF");
543 12 : parseVector("NL_STRIDE",nl_st);
544 : //if(nl_st.size()!=getNumberOfArguments() && nl_st.size()!=0) error("not enough values for NL_STRIDE");
545 12 : parseVector("NL_SKIN",nl_skin);
546 : //if(nl_skin.size()!=getNumberOfArguments() && nl_skin.size()!=0) error("not enough values for NL_SKIN");
547 36 : for (unsigned j=0; j<Nlist; j++) {
548 24 : if(nl_cut[j]<=0.0) {
549 0 : error("NL_CUTOFF should be explicitly specified and positive");
550 : }
551 24 : if(nl_st[j]<=0) {
552 0 : error("NL_STRIDE should be explicitly specified and positive");
553 : }
554 24 : if(nl_skin[j]<=0.) {
555 0 : error("NL_SKIN should be explicitly specified and positive");
556 : }
557 24 : nl_cut[j]=nl_cut[j]+nl_skin[j];
558 : }
559 12 : log << "Creating Neighbor Lists \n";
560 : // WARNING: is nl_cut meaningful here?
561 24 : nlall= Tools::make_unique<NeighborList>(listall,true,pbc,getPbc(),comm,nl_cut[0],nl_st[0]);
562 12 : if(com) {
563 : //Build lists of Atoms for every COM
564 0 : for (unsigned i=0; i<compos.size(); i++) {
565 : // WARNING: is nl_cut meaningful here?
566 0 : nlcom[i] = Tools::make_unique<NeighborList>(comatm[i],true,pbc,getPbc(),comm,nl_cut[0],nl_st[0]);
567 : }
568 : }
569 : unsigned ncnt=0;
570 : // Direct blocks AA, BB, CC, ...
571 12 : if(direct) {
572 24 : for (unsigned j=0; j<Natm; j++) {
573 16 : nl[ncnt]= Tools::make_unique<NeighborList>(Plist[j],true,pbc,getPbc(),comm,nl_cut[j],nl_st[j]);
574 16 : ncnt+=1;
575 : }
576 : }
577 : // Cross blocks AB, AC, BC, ...
578 12 : if(cross) {
579 24 : for (unsigned j=0; j<Natm; j++) {
580 24 : for (unsigned i=j+1; i<Natm; i++) {
581 16 : nl[ncnt]= Tools::make_unique<NeighborList>(Plist[i],Plist[j],true,false,pbc,getPbc(),comm,nl_cut[ncnt],nl_st[ncnt]);
582 8 : ncnt+=1;
583 : }
584 : }
585 : }
586 : } else {
587 0 : log << "WARNING: Neighbor List not activated this has not been tested!! \n";
588 0 : nlall= Tools::make_unique<NeighborList>(listall,true,pbc,getPbc(),comm);
589 0 : for (unsigned j=0; j<Nlist; j++) {
590 0 : nl[j]= Tools::make_unique<NeighborList>(Plist[j],Plist[j],true,true,pbc,getPbc(),comm);
591 : }
592 : }
593 : // Output Nlist
594 12 : log << "Total Nlists: " << Nlist << " \n";
595 36 : for (unsigned j=0; j<Nlist; j++) {
596 24 : log << " list " << j+1 << " size " << nl[j]->size() << " \n";
597 : }
598 : // Calculate COM masses once and for all from lists
599 12 : if(com) {
600 0 : for(unsigned j=0; j<compos.size(); j++) {
601 : double commass=0.;
602 0 : for(unsigned i=0; i<nlcom[j]->getFullAtomList().size(); i++) {
603 0 : unsigned andx=nlcom[j]->getFullAtomList()[i].index();
604 0 : commass+=mypdb.getOccupancy()[andx];
605 : }
606 0 : for(unsigned i=0; i<nlcom[j]->getFullAtomList().size(); i++) {
607 0 : unsigned andx=nlcom[j]->getFullAtomList()[i].index();
608 0 : if(commass>0.) {
609 0 : fmass[andx]=mypdb.getOccupancy()[andx]/commass;
610 : } else {
611 0 : fmass[andx]=1.;
612 : }
613 : }
614 : }
615 : }
616 :
617 : // Sorting
618 12 : dosort.resize(Nlist);
619 12 : std:: vector<int> ynsort(Nlist);
620 12 : parseVector("SORT",ynsort);
621 36 : for (unsigned i=0; i<Nlist; i++) {
622 24 : if(ynsort[i]==0||CompDer) {
623 : dosort[i]=false;
624 : } else {
625 : dosort[i]=true;
626 : }
627 : }
628 :
629 : //build box vectors and correct for pbc
630 12 : log << "Building the box from PDB data ... \n";
631 12 : Tensor Box=mypdb.getBoxVec();
632 12 : log << " Done! A,B,C vectors in Cartesian space: \n";
633 12 : log.printf(" A: %12.6f%12.6f%12.6f\n", Box[0][0],Box[0][1],Box[0][2]);
634 12 : log.printf(" B: %12.6f%12.6f%12.6f\n", Box[1][0],Box[1][1],Box[1][2]);
635 12 : log.printf(" C: %12.6f%12.6f%12.6f\n", Box[2][0],Box[2][1],Box[2][2]);
636 12 : log << "Changing the PBC according to the new box \n";
637 12 : Pbc mypbc;
638 12 : mypbc.setBox(Box);
639 12 : log << "The box volume is " << mypbc.getBox().determinant() << " \n";
640 :
641 : //Compute scaling factor
642 12 : if(Svol) {
643 12 : Fvol=cbrt(Vol0/mypbc.getBox().determinant());
644 12 : log << "Scaling atom distances by " << Fvol << " \n";
645 : } else {
646 0 : log << "Using unscaled atom distances \n";
647 : }
648 :
649 12 : r00.resize(Nlist);
650 12 : sw.resize(Nlist);
651 36 : for (unsigned j=0; j<Nlist; j++) {
652 48 : if( !parseNumbered( "SWITCH", j+1, sw[j] ) ) {
653 : break;
654 : }
655 : }
656 12 : if(CompDer) {
657 : // Set switching function parameters here only if computing derivatives
658 : // now set at the beginning of the dynamics to solve the r0 issue
659 6 : log << "Switching Function Parameters \n";
660 6 : sfs.resize(Nlist);
661 : std::string errors;
662 18 : for (unsigned j=0; j<Nlist; j++) {
663 12 : if(Svol) {
664 : double r0;
665 : std::string old_r0;
666 12 : vector<string> data=Tools::getWords(sw[j]);
667 : data.erase(data.begin());
668 12 : Tools::parse(data,"R_0",old_r0);
669 12 : Tools::convert(old_r0,r0);
670 12 : r0*=Fvol;
671 : std::string new_r0;
672 12 : Tools::convert(r0,new_r0);
673 12 : std::size_t pos = sw[j].find("R_0");
674 12 : sw[j].replace(pos+4,old_r0.size(),new_r0);
675 12 : }
676 12 : sfs[j].set(sw[j],errors);
677 : std::string num;
678 12 : Tools::convert(j+1, num);
679 12 : if( errors.length()!=0 ) {
680 0 : error("problem reading SWITCH" + num + " keyword : " + errors );
681 : }
682 12 : r00[j]=sfs[j].get_r0();
683 12 : log << " Swf: " << j << " r0=" << (sfs[j].description()).c_str() << " \n";
684 : }
685 : }
686 :
687 : // build COMs from positions if requested
688 12 : if(com) {
689 0 : for(unsigned j=0; j<compos.size(); j++) {
690 0 : compos[j][0]=0.;
691 0 : compos[j][1]=0.;
692 0 : compos[j][2]=0.;
693 0 : for(unsigned i=0; i<nlcom[j]->getFullAtomList().size(); i++) {
694 0 : unsigned andx=nlcom[j]->getFullAtomList()[i].index();
695 0 : compos[j]+=fmass[andx]*mypdb.getPositions()[andx];
696 : }
697 : }
698 : }
699 : // build the rPIV distances (transformation and sorting is done afterwards)
700 12 : if(CompDer) {
701 6 : log << " PIV | block | Size | Zeros | Ones |" << " \n";
702 : }
703 36 : for(unsigned j=0; j<Nlist; j++) {
704 11516328 : for(unsigned i=0; i<nl[j]->size(); i++) {
705 11516304 : unsigned i0=(nl[j]->getClosePairAtomNumber(i).first).index();
706 11516304 : unsigned i1=(nl[j]->getClosePairAtomNumber(i).second).index();
707 : //calculate/get COM position of centers i0 and i1
708 11516304 : Vector Pos0,Pos1;
709 11516304 : if(com) {
710 : //if(pbc) makeWhole();
711 0 : Pos0=compos[i0];
712 0 : Pos1=compos[i1];
713 : } else {
714 11516304 : Pos0=mypdb.getPositions()[i0];
715 11516304 : Pos1=mypdb.getPositions()[i1];
716 : }
717 11516304 : Vector ddist;
718 11516304 : if(pbc) {
719 11516304 : ddist=mypbc.distance(Pos0,Pos1);
720 : } else {
721 0 : ddist=delta(Pos0,Pos1);
722 : }
723 11516304 : double df=0.;
724 : // Transformation and sorting done at the first timestep to solve the r0 definition issue
725 11516304 : if(CompDer) {
726 1104 : rPIV[j].push_back(sfs[j].calculate(ddist.modulo()*Fvol, df));
727 : } else {
728 11515200 : rPIV[j].push_back(ddist.modulo()*Fvol);
729 : }
730 : }
731 24 : if(CompDer) {
732 12 : if(dosort[j]) {
733 0 : std::sort(rPIV[j].begin(),rPIV[j].end());
734 : }
735 : int lmt0=0;
736 : int lmt1=0;
737 1116 : for(unsigned i=0; i<rPIV[j].size(); i++) {
738 1104 : if(int(rPIV[j][i]*double(Nprec-1))==0) {
739 0 : lmt0+=1;
740 : }
741 1104 : if(int(rPIV[j][i]*double(Nprec-1))==1) {
742 0 : lmt1+=1;
743 : }
744 : }
745 12 : log.printf(" |%10i|%15zu|%15i|%15i|\n", j, rPIV[j].size(), lmt0, lmt1);
746 : }
747 : }
748 :
749 12 : checkRead();
750 : // From the plumed manual on how to build-up a new Colvar
751 12 : addValueWithDerivatives();
752 12 : requestAtoms(nlall->getFullAtomList());
753 12 : setNotPeriodic();
754 : // getValue()->setPeridodicity(false);
755 : // set size of derivative vector
756 12 : m_deriv.resize(getNumberOfAtoms());
757 12 : }
758 :
759 327 : void PIV::calculate() {
760 :
761 : // Local variables
762 :
763 327 : if(sharedData_unique) {
764 : // This is executed by the first PIV instance.
765 : // We initialize variables with the correct Nlist.
766 321 : sharedData_unique->prev_pos.resize(Nlist);
767 321 : sharedData_unique->cPIV.resize(Nlist);
768 321 : sharedData_unique->Atom0.resize(Nlist);
769 321 : sharedData_unique->Atom1.resize(Nlist);
770 : }
771 :
772 : // create references to minimize the impact of the code below
773 327 : auto & prev_stp(sharedData->prev_stp);
774 : auto & init_stp(sharedData->init_stp);
775 : auto & prev_pos(sharedData->prev_pos);
776 : auto & cPIV(sharedData->cPIV);
777 : auto & Atom0(sharedData->Atom0);
778 : auto & Atom1(sharedData->Atom1);
779 :
780 327 : std:: vector<std:: vector<int> > A0(Nprec);
781 327 : std:: vector<std:: vector<int> > A1(Nprec);
782 : size_t stride=1;
783 : unsigned rank=0;
784 :
785 327 : if(!serial) {
786 327 : stride=comm.Get_size();
787 327 : rank=comm.Get_rank();
788 : } else {
789 : stride=1;
790 : rank=0;
791 : }
792 :
793 : // Transform (and sort) the rPIV before starting the dynamics
794 327 : if (((prev_stp==-1) || (init_stp==1)) &&!CompDer) {
795 6 : if(prev_stp!=-1) {
796 3 : init_stp=0;
797 : }
798 : // Calculate the volume scaling factor
799 6 : if(Svol) {
800 6 : Fvol=cbrt(Vol0/getBox().determinant());
801 : }
802 : //Set switching function parameters
803 6 : log << "\n";
804 6 : log << "REFERENCE PDB # " << prev_stp+2 << " \n";
805 : // Set switching function parameters here only if computing derivatives
806 : // now set at the beginning of the dynamics to solve the r0 issue
807 6 : log << "Switching Function Parameters \n";
808 6 : sfs.resize(Nlist);
809 : std::string errors;
810 18 : for (unsigned j=0; j<Nlist; j++) {
811 12 : if(Svol) {
812 : double r0;
813 : std::string old_r0;
814 12 : vector<string> data=Tools::getWords(sw[j]);
815 : data.erase(data.begin());
816 12 : Tools::parse(data,"R_0",old_r0);
817 12 : Tools::convert(old_r0,r0);
818 12 : r0*=Fvol;
819 : std::string new_r0;
820 12 : Tools::convert(r0,new_r0);
821 12 : std::size_t pos = sw[j].find("R_0");
822 12 : sw[j].replace(pos+4,old_r0.size(),new_r0);
823 12 : }
824 12 : sfs[j].set(sw[j],errors);
825 : std::string num;
826 12 : Tools::convert(j+1, num);
827 12 : if( errors.length()!=0 ) {
828 0 : error("problem reading SWITCH" + num + " keyword : " + errors );
829 : }
830 12 : r00[j]=sfs[j].get_r0();
831 12 : log << " Swf: " << j << " r0=" << (sfs[j].description()).c_str() << " \n";
832 : }
833 : //Transform and sort
834 6 : log << "Building Reference PIV Vector \n";
835 6 : log << " PIV | block | Size | Zeros | Ones |" << " \n";
836 6 : double df=0.;
837 18 : for (unsigned j=0; j<Nlist; j++) {
838 11515212 : for (unsigned i=0; i<rPIV[j].size(); i++) {
839 11515200 : rPIV[j][i]=sfs[j].calculate(rPIV[j][i], df);
840 : }
841 12 : if(dosort[j]) {
842 12 : std::sort(rPIV[j].begin(),rPIV[j].end());
843 : }
844 : int lmt0=0;
845 : int lmt1=0;
846 11515212 : for(unsigned i=0; i<rPIV[j].size(); i++) {
847 11515200 : if(int(rPIV[j][i]*double(Nprec-1))==0) {
848 26 : lmt0+=1;
849 : }
850 11515200 : if(int(rPIV[j][i]*double(Nprec-1))==1) {
851 63358 : lmt1+=1;
852 : }
853 : }
854 12 : log.printf(" |%10i|%15zu|%15i|%15i|\n", j, rPIV[j].size(), lmt0, lmt1);
855 : }
856 6 : log << "\n";
857 : }
858 : // Do the sorting only once per timestep to avoid building the PIV N times for N rPIV PDB structures!
859 327 : if ((getStep()>prev_stp&&getStep()%updatePIV==0)||CompDer) {
860 324 : if (CompDer) {
861 321 : log << " Step " << getStep() << " Computing Derivatives NON-SORTED PIV \n";
862 : }
863 : //
864 : // build COMs from positions if requested
865 324 : if(com) {
866 0 : if(pbc) {
867 0 : makeWhole();
868 : }
869 0 : for(unsigned j=0; j<compos.size(); j++) {
870 0 : compos[j][0]=0.;
871 0 : compos[j][1]=0.;
872 0 : compos[j][2]=0.;
873 0 : for(unsigned i=0; i<nlcom[j]->getFullAtomList().size(); i++) {
874 0 : unsigned andx=nlcom[j]->getFullAtomList()[i].index();
875 0 : compos[j]+=fmass[andx]*getPosition(andx);
876 : }
877 : }
878 : }
879 : // update neighbor lists when an atom moves out of the Neighbor list skin
880 324 : if (doneigh) {
881 : bool doupdate=false;
882 : // For the first step build previous positions = actual positions
883 324 : if (prev_stp==-1) {
884 6 : bool docom=com;
885 18 : for (unsigned j=0; j<Nlist; j++) {
886 9708 : for (unsigned i=0; i<nl[j]->getFullAtomList().size(); i++) {
887 9696 : Vector Pos;
888 9696 : if(docom) {
889 0 : Pos=compos[i];
890 : } else {
891 9696 : Pos=getPosition(nl[j]->getFullAtomList()[i].index());
892 : }
893 9696 : prev_pos[j].push_back(Pos);
894 : }
895 : }
896 : doupdate=true;
897 : }
898 : // Decide whether to update lists based on atom displacement, every stride
899 324 : std:: vector<std:: vector<Vector> > tmp_pos(Nlist);
900 324 : if (getStep() % nlall->getStride() ==0) {
901 324 : bool docom=com;
902 972 : for (unsigned j=0; j<Nlist; j++) {
903 20520 : for (unsigned i=0; i<nl[j]->getFullAtomList().size(); i++) {
904 19872 : Vector Pos;
905 19872 : if(docom) {
906 0 : Pos=compos[i];
907 : } else {
908 19872 : Pos=getPosition(nl[j]->getFullAtomList()[i].index());
909 : }
910 19872 : tmp_pos[j].push_back(Pos);
911 19872 : if (pbcDistance(tmp_pos[j][i],prev_pos[j][i]).modulo()>=nl_skin[j]) {
912 : doupdate=true;
913 : }
914 : }
915 : }
916 : }
917 : // Update Nlists if needed
918 324 : if (doupdate==true) {
919 18 : for (unsigned j=0; j<Nlist; j++) {
920 9708 : for (unsigned i=0; i<nl[j]->getFullAtomList().size(); i++) {
921 9696 : prev_pos[j][i]=tmp_pos[j][i];
922 : }
923 12 : nl[j]->update(prev_pos[j]);
924 12 : log << " Step " << getStep() << " Neighbour lists updated " << nl[j]->size() << " \n";
925 : }
926 : }
927 324 : }
928 : // Calculate the volume scaling factor
929 324 : if(Svol) {
930 324 : Fvol=cbrt(Vol0/getBox().determinant());
931 : }
932 324 : Vector ddist;
933 : // Global to local variables
934 324 : bool doserial=serial;
935 : // Build "Nlist" PIV blocks
936 972 : for(unsigned j=0; j<Nlist; j++) {
937 648 : if(dosort[j]) {
938 : // from global to local variables to speedup the for loop with if statements
939 6 : bool docom=com;
940 6 : bool dopbc=pbc;
941 : // Vectors collecting occupancies: OrdVec one rank, OrdVecAll all ranks
942 6 : std:: vector<int> OrdVec(Nprec,0);
943 6 : cPIV[j].resize(0);
944 6 : Atom0[j].resize(0);
945 6 : Atom1[j].resize(0);
946 : // Building distances for the PIV vector at time t
947 6 : if(timer) {
948 6 : stopwatch.start("1 Build cPIV");
949 : }
950 5757606 : for(unsigned i=rank; i<nl[j]->size(); i+=stride) {
951 5757600 : unsigned i0=(nl[j]->getClosePairAtomNumber(i).first).index();
952 5757600 : unsigned i1=(nl[j]->getClosePairAtomNumber(i).second).index();
953 5757600 : Vector Pos0,Pos1;
954 5757600 : if(docom) {
955 0 : Pos0=compos[i0];
956 0 : Pos1=compos[i1];
957 : } else {
958 5757600 : Pos0=getPosition(i0);
959 5757600 : Pos1=getPosition(i1);
960 : }
961 5757600 : if(dopbc) {
962 5757600 : ddist=pbcDistance(Pos0,Pos1);
963 : } else {
964 0 : ddist=delta(Pos0,Pos1);
965 : }
966 5757600 : double df=0.;
967 : //Integer sorting ... faster!
968 : //Transforming distances with the Switching function + real to integer transformation
969 5757600 : int Vint=int(sfs[j].calculate(ddist.modulo()*Fvol, df)*double(Nprec-1)+0.5);
970 : //Integer transformed distance values as index of the Ordering Vector OrdVec
971 5757600 : OrdVec[Vint]+=1;
972 : //Keeps track of atom indices for force and virial calculations
973 5757600 : A0[Vint].push_back(i0);
974 5757600 : A1[Vint].push_back(i1);
975 : }
976 6 : if(timer) {
977 6 : stopwatch.stop("1 Build cPIV");
978 : }
979 6 : if(timer) {
980 6 : stopwatch.start("2 Sort cPIV");
981 : }
982 6 : if(!doserial && comm.initialized()) {
983 : // Vectors keeping track of the dimension and the starting-position of the rank-specific pair vector in the big pair vector.
984 0 : std:: vector<int> Vdim(stride,0);
985 0 : std:: vector<int> Vpos(stride,0);
986 : // Vectors collecting occupancies: OrdVec one rank, OrdVecAll all ranks
987 0 : std:: vector<int> OrdVecAll(stride*Nprec);
988 : // Big vectors containing all Atom indexes for every occupancy (Atom0O(Nprec,n) and Atom1O(Nprec,n) matrices in one vector)
989 : std:: vector<int> Atom0F;
990 : std:: vector<int> Atom1F;
991 : // Vector used to reconstruct arrays
992 0 : std:: vector<unsigned> k(stride,0);
993 : // Zeros might be many, this slows down a lot due to MPI communication
994 : // Avoid passing the zeros (i=1) for atom indices
995 0 : for(unsigned i=1; i<Nprec; i++) {
996 : // Building long vectors with all atom indexes for occupancies ordered from i=1 to i=Nprec-1
997 : // Can this be avoided ???
998 0 : Atom0F.insert(Atom0F.end(),A0[i].begin(),A0[i].end());
999 0 : Atom1F.insert(Atom1F.end(),A1[i].begin(),A1[i].end());
1000 0 : A0[i].resize(0);
1001 0 : A1[i].resize(0);
1002 : }
1003 : // Resize partial arrays to fill up for the next PIV block
1004 0 : A0[0].resize(0);
1005 0 : A1[0].resize(0);
1006 0 : A0[Nprec-1].resize(0);
1007 0 : A1[Nprec-1].resize(0);
1008 : // Avoid passing the zeros (i=1) for atom indices
1009 0 : OrdVec[0]=0;
1010 0 : OrdVec[Nprec-1]=0;
1011 :
1012 : // Wait for all ranks before communication of Vectors
1013 0 : comm.Barrier();
1014 :
1015 : // pass the array sizes before passing the arrays
1016 0 : int dim=Atom0F.size();
1017 : // Vdim and Vpos keep track of the dimension and the starting-position of the rank-specific pair vector in the big pair vector.
1018 0 : comm.Allgather(&dim,1,&Vdim[0],1);
1019 :
1020 : // TO BE IMPROVED: the following may be done by the rank 0 (now every rank does it)
1021 : int Fdim=0;
1022 0 : for(unsigned i=1; i<stride; i++) {
1023 0 : Vpos[i]=Vpos[i-1]+Vdim[i-1];
1024 0 : Fdim+=Vdim[i];
1025 : }
1026 0 : Fdim+=Vdim[0];
1027 : // build big vectors for atom pairs on all ranks for all ranks
1028 0 : std:: vector<int> Atom0FAll(Fdim);
1029 0 : std:: vector<int> Atom1FAll(Fdim);
1030 : // TO BE IMPROVED: Allgathers may be substituted by gathers by proc 0
1031 : // Moreover vectors are gathered head-to-tail and assembled later-on in a serial step.
1032 : // Gather the full Ordering Vector (occupancies). This is what we need to build the PIV
1033 0 : comm.Allgather(&OrdVec[0],Nprec,&OrdVecAll[0],Nprec);
1034 : // Gather the vectors of atom pairs to keep track of the idexes for the forces
1035 0 : comm.Allgatherv(Atom0F.data(),Atom0F.size(),&Atom0FAll[0],&Vdim[0],&Vpos[0]);
1036 0 : comm.Allgatherv(Atom1F.data(),Atom1F.size(),&Atom1FAll[0],&Vdim[0],&Vpos[0]);
1037 :
1038 : // Reconstruct the full vectors from collections of Allgathered parts (this is a serial step)
1039 : // This is the tricky serial step, to assemble together PIV and atom-pair info from head-tail big vectors
1040 : // Loop before on l and then on i would be better but the allgather should be modified
1041 : // Loop on blocks
1042 : //for(unsigned m=0;m<Nlist;m++) {
1043 : // Loop on Ordering Vector size excluding zeros (i=1)
1044 0 : if(timer) {
1045 0 : stopwatch.stop("2 Sort cPIV");
1046 : }
1047 0 : if(timer) {
1048 0 : stopwatch.start("3 Reconstruct cPIV");
1049 : }
1050 0 : for(unsigned i=1; i<Nprec; i++) {
1051 : // Loop on the ranks
1052 0 : for(unsigned l=0; l<stride; l++) {
1053 : // Loop on the number of head-to-tail pieces
1054 0 : for(unsigned m=0; m<OrdVecAll[i+l*Nprec]; m++) {
1055 : // cPIV is the current PIV at time t
1056 0 : cPIV[j].push_back(double(i)/double(Nprec-1));
1057 0 : Atom0[j].push_back(Atom0FAll[k[l]+Vpos[l]]);
1058 0 : Atom1[j].push_back(Atom1FAll[k[l]+Vpos[l]]);
1059 0 : k[l]+=1;
1060 : }
1061 : }
1062 : }
1063 0 : if(timer) {
1064 0 : stopwatch.stop("3 Reconstruct cPIV");
1065 : }
1066 : } else {
1067 6000000 : for(unsigned i=1; i<Nprec; i++) {
1068 11757594 : for(unsigned m=0; m<OrdVec[i]; m++) {
1069 5757600 : cPIV[j].push_back(double(i)/double(Nprec-1));
1070 5757600 : Atom0[j].push_back(A0[i][m]);
1071 5757600 : Atom1[j].push_back(A1[i][m]);
1072 : }
1073 : }
1074 : }
1075 : }
1076 : }
1077 : }
1078 327 : Vector distance;
1079 327 : double dfunc=0.;
1080 : // Calculate volume scaling factor
1081 327 : if(Svol) {
1082 327 : Fvol=cbrt(Vol0/getBox().determinant());
1083 : }
1084 :
1085 : // This test may be run by specifying the TEST keyword as input, it pritnts rPIV and cPIV and quits
1086 327 : if(test) {
1087 : unsigned limit=0;
1088 0 : for(unsigned j=0; j<Nlist; j++) {
1089 0 : if(dosort[j]) {
1090 0 : limit = cPIV[j].size();
1091 : } else {
1092 0 : limit = rPIV[j].size();
1093 : }
1094 0 : log.printf("PIV Block: %6i %12s %6i \n", j, " Size:", limit);
1095 0 : log.printf("%6s%6s%12s%12s%36s\n"," i"," j", " c-PIV "," r-PIV "," i-j distance vector ");
1096 0 : for(unsigned i=0; i<limit; i++) {
1097 : unsigned i0=0;
1098 : unsigned i1=0;
1099 0 : if(dosort[j]) {
1100 0 : i0=Atom0[j][i];
1101 0 : i1=Atom1[j][i];
1102 : } else {
1103 0 : i0=(nl[j]->getClosePairAtomNumber(i).first).index();
1104 0 : i1=(nl[j]->getClosePairAtomNumber(i).second).index();
1105 : }
1106 0 : Vector Pos0,Pos1;
1107 0 : if(com) {
1108 0 : Pos0=compos[i0];
1109 0 : Pos1=compos[i1];
1110 : } else {
1111 0 : Pos0=getPosition(i0);
1112 0 : Pos1=getPosition(i1);
1113 : }
1114 0 : if(pbc) {
1115 0 : distance=pbcDistance(Pos0,Pos1);
1116 : } else {
1117 0 : distance=delta(Pos0,Pos1);
1118 : }
1119 0 : dfunc=0.;
1120 : double cP,rP;
1121 0 : if(dosort[j]) {
1122 0 : cP = cPIV[j][i];
1123 0 : rP = rPIV[j][rPIV[j].size()-cPIV[j].size()+i];
1124 : } else {
1125 0 : double dm=distance.modulo();
1126 0 : cP = sfs[j].calculate(dm*Fvol, dfunc);
1127 0 : rP = rPIV[j][i];
1128 : }
1129 0 : log.printf("%6i%6i%12.6f%12.6f%12.6f%12.6f%12.6f\n",i0,i1,cP,rP,distance[0],distance[1],distance[2]);
1130 : }
1131 : }
1132 0 : log.printf("This was a test, now exit \n");
1133 0 : exit();
1134 : }
1135 :
1136 327 : if(timer) {
1137 2 : stopwatch.start("4 Build For Derivatives");
1138 : }
1139 : // non-global variables Nder and Scalevol defined to speedup if structures in cycles
1140 327 : bool Nder=CompDer;
1141 327 : bool Scalevol=Svol;
1142 327 : if(getStep()%updatePIV==0) {
1143 : // set to zero PIVdistance, derivatives and virial when they are calculated
1144 29799 : for(unsigned j=0; j<m_deriv.size(); j++) {
1145 117888 : for(unsigned k=0; k<3; k++) {
1146 88416 : m_deriv[j][k]=0.;
1147 : }
1148 : }
1149 1308 : for(unsigned j=0; j<3; j++) {
1150 3924 : for(unsigned k=0; k<3; k++) {
1151 2943 : m_virial[j][k]=0.;
1152 : }
1153 : }
1154 327 : m_PIVdistance=0.;
1155 : // Re-compute atomic distances for derivatives and compute PIV-PIV distance
1156 981 : for(unsigned j=0; j<Nlist; j++) {
1157 : unsigned limit=0;
1158 : // dosorting definition is to speedup if structure in cycles with non-global variables
1159 654 : bool dosorting=dosort[j];
1160 654 : bool docom=com;
1161 654 : bool dopbc=pbc;
1162 654 : if(dosorting) {
1163 12 : limit = cPIV[j].size();
1164 : } else {
1165 642 : limit = rPIV[j].size();
1166 : }
1167 11574918 : for(unsigned i=rank; i<limit; i+=stride) {
1168 : unsigned i0=0;
1169 : unsigned i1=0;
1170 11574264 : if(dosorting) {
1171 11515200 : i0=Atom0[j][i];
1172 11515200 : i1=Atom1[j][i];
1173 : } else {
1174 59064 : i0=(nl[j]->getClosePairAtomNumber(i).first).index();
1175 59064 : i1=(nl[j]->getClosePairAtomNumber(i).second).index();
1176 : }
1177 11574264 : Vector Pos0,Pos1;
1178 11574264 : if(docom) {
1179 0 : Pos0=compos[i0];
1180 0 : Pos1=compos[i1];
1181 : } else {
1182 11574264 : Pos0=getPosition(i0);
1183 11574264 : Pos1=getPosition(i1);
1184 : }
1185 11574264 : if(dopbc) {
1186 11574264 : distance=pbcDistance(Pos0,Pos1);
1187 : } else {
1188 0 : distance=delta(Pos0,Pos1);
1189 : }
1190 11574264 : dfunc=0.;
1191 : // this is needed for dfunc and dervatives
1192 11574264 : double dm=distance.modulo();
1193 11574264 : double tPIV = sfs[j].calculate(dm*Fvol, dfunc);
1194 : // PIV distance
1195 : double coord=0.;
1196 11574264 : if(!dosorting||Nder) {
1197 59064 : coord = tPIV - rPIV[j][i];
1198 : } else {
1199 11515200 : coord = cPIV[j][i] - rPIV[j][rPIV[j].size()-cPIV[j].size()+i];
1200 : }
1201 : // Calculate derivatives, virial, and variable=sum_j (scaling[j] *(cPIV-rPIV)_j^2)
1202 : // WARNING: dfunc=dswf/(Fvol*dm) (this may change in future Plumed versions)
1203 11574264 : double tmp = 2.*scaling[j]*coord*Fvol*Fvol*dfunc;
1204 11574264 : Vector tmpder = tmp*distance;
1205 : // 0.5*(x_i-x_k)*f_ik (force on atom k due to atom i)
1206 11574264 : if(docom) {
1207 0 : Vector dist;
1208 0 : for(unsigned k=0; k<nlcom[i0]->getFullAtomList().size(); k++) {
1209 0 : unsigned x0=nlcom[i0]->getFullAtomList()[k].index();
1210 0 : m_deriv[x0] -= tmpder*fmass[x0];
1211 0 : for(unsigned l=0; l<3; l++) {
1212 0 : dist[l]=0.;
1213 : }
1214 0 : Vector P0=getPosition(x0);
1215 0 : for(unsigned l=0; l<nlcom[i0]->getFullAtomList().size(); l++) {
1216 0 : unsigned x1=nlcom[i0]->getFullAtomList()[l].index();
1217 0 : Vector P1=getPosition(x1);
1218 0 : if(dopbc) {
1219 0 : dist+=pbcDistance(P0,P1);
1220 : } else {
1221 0 : dist+=delta(P0,P1);
1222 : }
1223 : }
1224 0 : for(unsigned l=0; l<nlcom[i1]->getFullAtomList().size(); l++) {
1225 0 : unsigned x1=nlcom[i1]->getFullAtomList()[l].index();
1226 0 : Vector P1=getPosition(x1);
1227 0 : if(dopbc) {
1228 0 : dist+=pbcDistance(P0,P1);
1229 : } else {
1230 0 : dist+=delta(P0,P1);
1231 : }
1232 : }
1233 0 : m_virial -= 0.25*fmass[x0]*Tensor(dist,tmpder);
1234 : }
1235 0 : for(unsigned k=0; k<nlcom[i1]->getFullAtomList().size(); k++) {
1236 0 : unsigned x1=nlcom[i1]->getFullAtomList()[k].index();
1237 0 : m_deriv[x1] += tmpder*fmass[x1];
1238 0 : for(unsigned l=0; l<3; l++) {
1239 0 : dist[l]=0.;
1240 : }
1241 0 : Vector P1=getPosition(x1);
1242 0 : for(unsigned l=0; l<nlcom[i1]->getFullAtomList().size(); l++) {
1243 0 : unsigned x0=nlcom[i1]->getFullAtomList()[l].index();
1244 0 : Vector P0=getPosition(x0);
1245 0 : if(dopbc) {
1246 0 : dist+=pbcDistance(P1,P0);
1247 : } else {
1248 0 : dist+=delta(P1,P0);
1249 : }
1250 : }
1251 0 : for(unsigned l=0; l<nlcom[i0]->getFullAtomList().size(); l++) {
1252 0 : unsigned x0=nlcom[i0]->getFullAtomList()[l].index();
1253 0 : Vector P0=getPosition(x0);
1254 0 : if(dopbc) {
1255 0 : dist+=pbcDistance(P1,P0);
1256 : } else {
1257 0 : dist+=delta(P1,P0);
1258 : }
1259 : }
1260 0 : m_virial += 0.25*fmass[x1]*Tensor(dist,tmpder);
1261 : }
1262 : } else {
1263 11574264 : m_deriv[i0] -= tmpder;
1264 11574264 : m_deriv[i1] += tmpder;
1265 11574264 : m_virial -= tmp*Tensor(distance,distance);
1266 : }
1267 11574264 : if(Scalevol) {
1268 11574264 : m_virial+=1./3.*tmp*dm*dm*Tensor::identity();
1269 : }
1270 11574264 : m_PIVdistance += scaling[j]*coord*coord;
1271 : }
1272 : }
1273 :
1274 327 : if (!serial && comm.initialized()) {
1275 0 : comm.Barrier();
1276 0 : comm.Sum(&m_PIVdistance,1);
1277 0 : if(!m_deriv.empty()) {
1278 0 : comm.Sum(&m_deriv[0][0],3*m_deriv.size());
1279 : }
1280 0 : comm.Sum(&m_virial[0][0],9);
1281 : }
1282 : }
1283 327 : prev_stp=getStep();
1284 :
1285 : //Timing
1286 327 : if(timer) {
1287 2 : stopwatch.stop("4 Build For Derivatives");
1288 : }
1289 327 : if(timer) {
1290 1 : log.printf("Timings for action %s with label %s \n", getName().c_str(), getLabel().c_str() );
1291 1 : log<<stopwatch;
1292 : }
1293 :
1294 : // Update derivatives, virial, and variable (PIV-distance^2)
1295 29799 : for(unsigned i=0; i<m_deriv.size(); ++i) {
1296 29472 : setAtomsDerivatives(i,m_deriv[i]);
1297 : }
1298 327 : setValue (m_PIVdistance);
1299 327 : setBoxDerivatives (m_virial);
1300 327 : }
1301 : //Close Namespaces at the very beginning
1302 : }
1303 : }
1304 :
|
