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