Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "CLTool.h"
23 : #include "core/CLToolRegister.h"
24 : #include "tools/Tools.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "core/ActionWithValue.h"
28 : #include "core/ActionWithVirtualAtom.h"
29 : #include "core/ActionShortcut.h"
30 : #include "tools/Communicator.h"
31 : #include "tools/Random.h"
32 : #include "tools/Pbc.h"
33 : #include <cstdio>
34 : #include <cstring>
35 : #include <vector>
36 : #include <map>
37 : #include <memory>
38 : #include "tools/Units.h"
39 : #include "tools/PDB.h"
40 : #include "tools/FileBase.h"
41 : #include "tools/IFile.h"
42 : #include "tools/TrajectoryParser.h"
43 :
44 : namespace PLMD {
45 : namespace cltools {
46 :
47 : //+PLUMEDOC TOOLS driver-float
48 : /*
49 : A tool that does the same as driver, but using single precision reals.
50 :
51 : We recommend always using [driver](driver.md) to postprocess your trajectories. However,
52 : if you want to test what PLUMED does when linked from a single precision code you can use
53 : this version of driver. The documentation for this command line tool is identical to that for
54 : [driver](driver.md).
55 :
56 : ## Examples
57 :
58 : ```plumed
59 : plumed driver-float --plumed plumed.dat --ixyz trajectory.xyz
60 : ```
61 :
62 : For more examples see the documentation for [driver](driver.md).
63 :
64 : */
65 : //+ENDPLUMEDOC
66 : //
67 :
68 :
69 : //+PLUMEDOC TOOLS driver
70 : /*
71 : driver is a tool that allows one to use plumed to post-process an existing trajectory.
72 :
73 : The input to driver is specified using the command line arguments described below.
74 :
75 : In addition, you can use the special [READ](READ.md) command inside your plumed input
76 : to read in colvar files that were generated during your MD simulation. The values
77 : read in from the file can then be treated like calculated colvars.
78 :
79 : !!! warning "masses and chages"
80 :
81 : Notice that by default the driver has no knowledge about the masses and charges
82 : of your atoms! Thus, if you want to compute quantities that depend on charges (e.g. [DHENERGY](DHENERGY.md))
83 : or masses (e.g. [COM](COM.md)) you need to pass masses and charges to driver.
84 : You can do this by either using --pdb option or the --mc option. The latter
85 : will read a file produced by [DUMPMASSCHARGE](DUMPMASSCHARGE.md).
86 :
87 :
88 : ## Examples
89 :
90 : The following command tells plumed to post process the trajectory contained in `trajectory.xyz`
91 : by performing the actions described in the input file `plumed.dat`. If an action that takes the
92 : stride keyword is given a stride equal to $n$ then it will be performed only on every $n$th
93 : frame in the trajectory file.
94 :
95 : ```plumed
96 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz
97 : ```
98 :
99 : Notice that `xyz` files are expected to be in the internal PLUMED unit of nm.
100 : You can change this behavior by using the `--length-units` option as shown below
101 :
102 : ```plumed
103 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --length-units A
104 : ```
105 :
106 : The strings accepted by the `--length-units` options are the same ones accepted by the [UNITS](UNITS.md) action.
107 : Other file formats typically have their default coordinates (e.g., `gro` files are always in nm)
108 : and it thus should not be necessary to use the `--length-units` option. Additionally,
109 : consider that the units used by the `driver` might be different by the units used in the PLUMED input
110 : file `plumed.dat`. For instance consider the command:
111 :
112 : ```plumed
113 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --length-units A
114 : ```
115 :
116 : which uses the following `plumed.dat` file
117 :
118 : ```plumed
119 : # no explicit UNITS action here
120 : d: DISTANCE ATOMS=1,2
121 : PRINT ARG=d FILE=colvar
122 : ```
123 :
124 : In this case, the driver reads the `xyz` file assuming it to contain coordinates in Angstroms.
125 : However, the resulting `colvar` file contains a distance expressed in nm.
126 :
127 : The following command tells plumed to post process the trajectory contained in trajectory.xyz.
128 : by performing the actions described in the input file plumed.dat.
129 :
130 : ```plumed
131 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --trajectory-stride 100 --timestep 0.001
132 : ```
133 :
134 : Here though
135 : `--trajectory-stride` is set equal to the frequency with which frames were output during the trajectory
136 : and the `--timestep` is equal to the simulation timestep. As such the `STRIDE` parameters in the `plumed.dat`
137 : files are referred to the original timestep and any files output resemble those that would have been generated
138 : had we run the calculation we are running with driver when the MD simulation was running.
139 :
140 : PLUMED can read xyz files (in PLUMED units) and gro files (in nm). In addition,
141 : PLUMED includes support for a
142 : subset of the trajectory file formats supported by VMD, e.g. xtc and dcd:
143 :
144 : ```plumed
145 : plumed driver --plumed plumed.dat --pdb diala.pdb --mf_xtc traj.xtc --trajectory-stride 100 --timestep 0.001
146 : ```
147 :
148 : where `--mf_` prefixes the extension of one of the accepted molfile plugin format.
149 : If PLUMED has been installed with full molfile support, other formats will be available.
150 : Just type `plumed driver --help` to see which file formats you can use.
151 :
152 : Molfile plugin requires the periodic cell to be triangular (i.e. the first vector oriented along x and
153 : second vector in xy plane). This is true for many MD codes. However, it could be false
154 : if you rotate the coordinates in your trajectory before reading them in the driver.
155 : Also notice that some formats (e.g. amber crd) do not specify the total number of atoms. In this case you can use
156 : the `--natoms` option as shown below
157 :
158 : ```plumed
159 : plumed driver --plumed plumed.dat --mf_crd trajectory.crd --natoms 128
160 : ```
161 :
162 : Check the available molfile plugins and limitations at [this link](http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/).
163 :
164 : You can also use the xdrfile implementation of xtc and trr with plumed driver. If you want to do so you just
165 : download and install properly the xdrfile library (see [this link](http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library)).
166 : If the xdrfile library is installed properly the PLUMED configure script should be able to
167 : detect it and enable it.
168 : Notice that the xdrfile implementation of xtc and trr
169 : is more robust than the molfile one, since it provides support for generic cell shapes.
170 : In addition, if you install xdrfile you can then use the [DUMPATOMS](DUMPATOMS.md) command to write compressed xtc files.
171 :
172 : ## Multiple replicas
173 :
174 : When PLUMED is compiled with MPI support, you can emulate a multi-simulation setup with `driver` by providing the `--multi`
175 : option with the appropriate number of ranks. This allows you to use the ref special-replica-syntax that is discussed [here](parsing.md) to analyze multiple
176 : trajectories (see [this tutorial](https://www.plumed-tutorials.org/lessons/21/005/data/NAVIGATION.html)). PLUMED will also automatically append a numbered suffix to output files
177 : (e.g. `COLVAR.0`, `COLVAR.1`, …) as discussed [here](Files.md). Similarly, each replica will search for the corresponding suffixed input file (e.g. `traj.0.xtc`, …)
178 : or default to the unsuffixed one.
179 :
180 : */
181 : //+ENDPLUMEDOC
182 : //
183 :
184 : template<typename real>
185 : class Driver : public CLTool {
186 : public:
187 : static void registerKeywords( Keywords& keys );
188 : explicit Driver(const CLToolOptions& co );
189 : int main(FILE* in,FILE*out,Communicator& pc) override;
190 : void evaluateNumericalDerivatives( const long long int& step, PlumedMain& p, const std::vector<real>& coordinates,
191 : const std::vector<real>& masses, const std::vector<real>& charges,
192 : std::vector<real>& cell, const double& base, std::vector<real>& numder );
193 : std::string description()const override;
194 : };
195 :
196 : template<typename real>
197 14642 : void Driver<real>::registerKeywords( Keywords& keys ) {
198 14642 : CLTool::registerKeywords( keys );
199 : keys.isDriver();
200 14642 : keys.addFlag("--help-debug",false,"print special options that can be used to create regtests");
201 14642 : keys.add("compulsory","--plumed","plumed.dat","specify the name of the plumed input file");
202 14642 : keys.add("compulsory","--timestep","1.0","the timestep that was used in the calculation that produced this trajectory in picoseconds");
203 14642 : keys.add("compulsory","--trajectory-stride","1","the frequency with which frames were output to this trajectory during the simulation"
204 : " (0 means that the number of the step is read from the trajectory file,"
205 : " currently working only for xtc/trr files read with --ixtc/--trr)"
206 : );
207 14642 : keys.add("compulsory","--multi","0","set number of replicas for multi environment (needs MPI)");
208 14642 : keys.addFlag("--noatoms",false,"don't read in a trajectory. Just use colvar files as specified in plumed.dat");
209 14642 : keys.addFlag("--parse-only",false,"read the plumed input file and stop");
210 14642 : keys.addFlag("--restart",false,"makes driver behave as if restarting");
211 14642 : TrajectoryParser::registerKeywords(keys);
212 14642 : keys.add("optional","--shortcut-ofile","the name of the file to output info on the way shortcuts have been expanded. If there are no shortcuts in your input file nothing is output");
213 14642 : keys.add("optional","--valuedict-ofile","output a dictionary giving information about each value in the input file");
214 14642 : keys.add("optional","--length-units","units for length, either as a string or a number");
215 14642 : keys.add("optional","--mass-units","units for mass in pdb and mc file, either as a string or a number");
216 14642 : keys.add("optional","--charge-units","units for charge in pdb and mc file, either as a string or a number");
217 14642 : keys.add("optional","--kt","set kT it will not be necessary to specify temperature in input file");
218 14642 : keys.add("optional","--dump-forces","dump the forces on a file");
219 14642 : keys.add("optional","--dump-forces-fmt","( default=%%f ) the format to use to dump the forces");
220 14642 : keys.addFlag("--dump-full-virial",false,"with --dump-forces, it dumps the 9 components of the virial");
221 14642 : keys.add("optional","--pdb","provides a pdb with masses and charges");
222 14642 : keys.add("optional","--mc","provides a file with masses and charges as produced with DUMPMASSCHARGE");
223 14642 : keys.add("optional","--box","comma-separated box dimensions (3 for orthorhombic, 9 for generic)");
224 14642 : keys.add("optional","--natoms","provides number of atoms - only used if file format does not contain number of atoms");
225 14642 : keys.add("optional","--initial-step","provides a number for the initial step, default is 0");
226 14642 : keys.add("optional","--debug-forces","output a file containing the forces due to the bias evaluated using numerical derivatives "
227 : "and using the analytical derivatives implemented in plumed");
228 14642 : keys.add("hidden","--debug-float","[yes/no] turns on the single precision version (to check float interface)");
229 14642 : keys.add("hidden","--debug-dd","[yes/no] use a fake domain decomposition");
230 14642 : keys.add("hidden","--debug-pd","[yes/no] use a fake particle decomposition");
231 14642 : keys.add("hidden","--debug-grex","use a fake gromacs-like replica exchange, specify exchange stride");
232 14642 : keys.add("hidden","--debug-grex-log","log file for debug=grex");
233 14642 : }
234 : template<typename real>
235 1345 : Driver<real>::Driver(const CLToolOptions& co ):
236 1345 : CLTool(co) {
237 1345 : inputdata=inputType::commandline;
238 1345 : }
239 : template<typename real>
240 7 : std::string Driver<real>::description()const {
241 7 : return "analyze trajectories with plumed";
242 : }
243 :
244 : template<typename real>
245 1331 : int Driver<real>::main(FILE* in,FILE*out,Communicator& pc) {
246 :
247 1331 : Units units;
248 1331 : PDB pdb;
249 :
250 : // Parse everything
251 : bool printhelpdebug;
252 1331 : parseFlag("--help-debug",printhelpdebug);
253 1331 : if( printhelpdebug ) {
254 : std::fprintf(out,"%s",
255 : "Additional options for debug (only to be used in regtest):\n"
256 : " [--debug-float yes] : turns on the single precision version (to check float interface)\n"
257 : " [--debug-dd yes] : use a fake domain decomposition\n"
258 : " [--debug-pd yes] : use a fake particle decomposition\n"
259 : );
260 : return 0;
261 : }
262 : // Are we reading trajectory data
263 : bool noatoms;
264 1331 : parseFlag("--noatoms",noatoms);
265 : bool parseOnly;
266 2662 : parseFlag("--parse-only",parseOnly);
267 : std::string full_outputfile;
268 2732 : parse("--shortcut-ofile",full_outputfile);
269 : std::string valuedict_file;
270 1331 : parse("--valuedict-ofile",valuedict_file);
271 : bool restart;
272 2732 : parseFlag("--restart",restart);
273 :
274 : std::string fakein;
275 : bool debug_float=false;
276 : fakein="";
277 2662 : if(parse("--debug-float",fakein)) {
278 2 : if(fakein=="yes") {
279 : debug_float=true;
280 0 : } else if(fakein=="no") {
281 : debug_float=false;
282 : } else {
283 0 : error("--debug-float should have argument yes or no");
284 : }
285 : }
286 :
287 : if(debug_float && sizeof(real)!=sizeof(float)) {
288 2 : auto cl=cltoolRegister().create(CLToolOptions("driver-float"));
289 : cl->setInputData(this->getInputData());
290 1 : int ret=cl->main(in,out,pc);
291 : return ret;
292 1 : }
293 :
294 : bool debug_pd=false;
295 : fakein="";
296 2660 : if(parse("--debug-pd",fakein)) {
297 12 : if(fakein=="yes") {
298 : debug_pd=true;
299 0 : } else if(fakein=="no") {
300 : debug_pd=false;
301 : } else {
302 0 : error("--debug-pd should have argument yes or no");
303 : }
304 : }
305 : if(debug_pd) {
306 : std::fprintf(out,"DEBUGGING PARTICLE DECOMPOSITION\n");
307 : }
308 :
309 : bool debug_dd=false;
310 : fakein="";
311 2660 : if(parse("--debug-dd",fakein)) {
312 92 : if(fakein=="yes") {
313 : debug_dd=true;
314 0 : } else if(fakein=="no") {
315 : debug_dd=false;
316 : } else {
317 0 : error("--debug-dd should have argument yes or no");
318 : }
319 : }
320 : if(debug_dd) {
321 : std::fprintf(out,"DEBUGGING DOMAIN DECOMPOSITION\n");
322 : }
323 :
324 1330 : if( debug_pd || debug_dd ) {
325 104 : if(noatoms) {
326 0 : error("cannot debug without atoms");
327 : }
328 : }
329 :
330 : // set up for multi replica driver:
331 1330 : int multi=0;
332 1400 : parse("--multi",multi);
333 1330 : Communicator intracomm;
334 1330 : Communicator intercomm;
335 1330 : if(multi) {
336 300 : int ntot=pc.Get_size();
337 300 : int nintra=ntot/multi;
338 300 : if(multi*nintra!=ntot) {
339 0 : error("invalid number of processes for multi environment");
340 : }
341 300 : pc.Split(pc.Get_rank()/nintra,pc.Get_rank(),intracomm);
342 300 : pc.Split(pc.Get_rank()%nintra,pc.Get_rank(),intercomm);
343 : } else {
344 1030 : intracomm.Set_comm(pc.Get_comm());
345 : }
346 :
347 : // set up for debug replica exchange:
348 1330 : bool debug_grex=parse("--debug-grex",fakein);
349 1330 : int grex_stride=0;
350 : FILE*grex_log=NULL;
351 : // call fclose when fp goes out of scope
352 724 : auto deleter=[](auto f) {
353 : if(f) {
354 724 : std::fclose(f);
355 : }
356 : };
357 : std::unique_ptr<FILE,decltype(deleter)> grex_log_deleter(grex_log,deleter);
358 :
359 1330 : if(debug_grex) {
360 30 : if(noatoms) {
361 0 : error("must have atoms to debug_grex");
362 : }
363 30 : if(multi<2) {
364 0 : error("--debug_grex needs --multi with at least two replicas");
365 : }
366 30 : if(!Tools::convertNoexcept(fakein,grex_stride) || grex_stride<=0) {
367 70 : error("--debug-grex should be a positive integer");
368 : }
369 : std::string n;
370 30 : Tools::convert(intercomm.Get_rank(),n);
371 : std::string file;
372 60 : parse("--debug-grex-log",file);
373 30 : if(file.length()>0) {
374 30 : file+="."+n;
375 30 : grex_log=std::fopen(file.c_str(),"w");
376 : grex_log_deleter.reset(grex_log);
377 : }
378 : }
379 :
380 : // Read the plumed input file name
381 : std::string plumedFile;
382 1330 : parse("--plumed",plumedFile);
383 : // the timestep
384 : double t;
385 1330 : parse("--timestep",t);
386 1330 : real timestep=real(t);
387 : // the stride
388 : unsigned stride;
389 1330 : parse("--trajectory-stride",stride);
390 : // are we writing forces
391 1540 : std::string dumpforces(""), debugforces(""), dumpforcesFmt("%f");;
392 1330 : bool dumpfullvirial=false;
393 1330 : if(!noatoms) {
394 1273 : parse("--dump-forces",dumpforces);
395 2546 : parse("--debug-forces",debugforces);
396 : }
397 1330 : if(dumpforces!="" || debugforces!="" ) {
398 1414 : parse("--dump-forces-fmt",dumpforcesFmt);
399 : }
400 1330 : if(dumpforces!="") {
401 1390 : parseFlag("--dump-full-virial",dumpfullvirial);
402 : }
403 1330 : if( debugforces!="" && (debug_dd || debug_pd) ) {
404 0 : error("cannot debug forces and domain/particle decomposition at same time");
405 : }
406 1 : if( debugforces!="" && sizeof(real)!=sizeof(double) ) {
407 0 : error("cannot debug forces in single precision mode");
408 : }
409 :
410 1330 : real kt=-1.0;
411 2730 : parse("--kt",kt);
412 : std::string trajectory_fmt;
413 :
414 : bool use_molfile=false;
415 :
416 : // Read in an xyz file
417 1540 : std::string trajectoryFile(""), pdbfile(""), mcfile("");
418 : bool pbc_cli_given=false;
419 1400 : std::vector<double> pbc_cli_box(9,0.0);
420 1330 : int command_line_natoms=-1;
421 :
422 1330 : TrajectoryParser parser;
423 1330 : if(!noatoms) {
424 : int nn=0;
425 : trajectory_fmt="";
426 8911 : for (const auto & trj_type : TrajectoryParser::trajectoryOptions()) {
427 : std::string tmp;
428 12730 : parse("--i"+trj_type, tmp);
429 6365 : if (tmp.length()>0) {
430 : trajectory_fmt=trj_type;
431 898 : ++nn;
432 : trajectoryFile=tmp;
433 : }
434 : }
435 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
436 : {
437 1273 : auto plugins_names=TrajectoryParser::getMolfilePluginsnames() ;
438 12730 : for(unsigned i=0; i<plugins_names.size(); i++) {
439 11457 : std::string molfile_key="--mf_"+plugins_names[i];
440 : std::string traj_molfile;
441 11457 : parse(molfile_key,traj_molfile);
442 11457 : if(traj_molfile.length()>0) {
443 358 : ++nn;
444 : std::fprintf(out,"\nDRIVER: Found molfile format trajectory %s with name %s\n",plugins_names[i].c_str(),traj_molfile.c_str());
445 : trajectoryFile=traj_molfile;
446 : trajectory_fmt=plugins_names[i];
447 : use_molfile=true;
448 : }
449 : }
450 1273 : }
451 : #endif
452 : {
453 :
454 1273 : if(nn>1) {
455 0 : std::fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
456 1 : return 1;
457 : }
458 : }
459 :
460 1273 : if(trajectoryFile.length()==0&&!parseOnly) {
461 0 : std::fprintf(stderr,"ERROR: missing trajectory data\n");
462 : return 1;
463 : }
464 1273 : std::string lengthUnits("");
465 2546 : parse("--length-units",lengthUnits);
466 1273 : if(lengthUnits.length()>0) {
467 31 : units.setLength(lengthUnits);
468 : }
469 1274 : std::string chargeUnits("");
470 2546 : parse("--charge-units",chargeUnits);
471 1273 : if(chargeUnits.length()>0) {
472 1 : units.setCharge(chargeUnits);
473 : }
474 1274 : std::string massUnits("");
475 2546 : parse("--mass-units",massUnits);
476 1273 : if(massUnits.length()>0) {
477 1 : units.setMass(massUnits);
478 : }
479 :
480 2546 : parse("--pdb",pdbfile);
481 1273 : if(pdbfile.length()>0) {
482 80 : bool check=pdb.read(pdbfile,false,1.0);
483 80 : if(!check) {
484 0 : error("error reading pdb file");
485 : }
486 : }
487 :
488 2547 : parse("--mc",mcfile);
489 :
490 : std::string pbc_cli_list;
491 2546 : parse("--box",pbc_cli_list);
492 1273 : if(pbc_cli_list.length()>0) {
493 : pbc_cli_given=true;
494 52 : std::vector<std::string> words=Tools::getWords(pbc_cli_list,",");
495 52 : if(words.size()==3 || words.size()==9) {
496 221 : for(unsigned i=0; i<words.size(); i++) {
497 : // with 3 numbers, they are the diagonal elements
498 171 : if(!Tools::convertNoexcept(words[i],pbc_cli_box[words.size()==3?4*i:i])) {
499 1 : error("cannot parse command-line box "+pbc_cli_list);
500 : }
501 : }
502 : } else {
503 1 : std::string msg="ERROR: cannot parse command-line box "+pbc_cli_list;
504 1 : std::fprintf(stderr,"%s\n",msg.c_str());
505 : return 1;
506 : }
507 :
508 52 : }
509 :
510 2543 : parse("--natoms",command_line_natoms);
511 :
512 : } // if(!noatoms)
513 :
514 1328 : if(debug_dd && debug_pd) {
515 0 : error("cannot use debug-dd and debug-pd at the same time");
516 : }
517 1328 : if(debug_pd || debug_dd) {
518 104 : if( !Communicator::initialized() ) {
519 0 : error("needs mpi for debug-pd");
520 : }
521 : }
522 :
523 1328 : PlumedMain p;
524 1328 : if( parseOnly ) {
525 17 : p.activateParseOnlyMode();
526 : }
527 1328 : p.cmd("setRealPrecision",(int)sizeof(real));
528 : int checknatoms=-1;
529 1328 : long long int step=0;
530 1328 : parse("--initial-step",step);
531 :
532 1328 : if(restart) {
533 2 : p.cmd("setRestart",1);
534 : }
535 :
536 1328 : if(Communicator::initialized()) {
537 483 : if(multi) {
538 300 : if(intracomm.Get_rank()==0) {
539 412 : p.cmd("GREX setMPIIntercomm",&intercomm.Get_comm());
540 : }
541 600 : p.cmd("GREX setMPIIntracomm",&intracomm.Get_comm());
542 300 : p.cmd("GREX init");
543 : }
544 966 : p.cmd("setMPIComm",&intracomm.Get_comm());
545 : }
546 1328 : p.cmd("setLog",out);
547 : // the units should be passed with the same precision as the other real numbers
548 1328 : p.cmd("setMDLengthUnits",real(units.getLength()));
549 1328 : p.cmd("setMDChargeUnits",real(units.getCharge()));
550 1328 : p.cmd("setMDMassUnits",real(units.getMass()));
551 1328 : p.cmd("setMDEngine","driver");
552 1328 : p.cmd("setTimestep",timestep);
553 1328 : if( !parseOnly || full_outputfile.length()==0 ) {
554 1323 : p.cmd("setPlumedDat",plumedFile.c_str());
555 : }
556 :
557 : int natoms=0;
558 :
559 1328 : if(parseOnly) {
560 17 : if(command_line_natoms<0) {
561 0 : error("--parseOnly requires setting the number of atoms with --natoms");
562 : }
563 : natoms=command_line_natoms;
564 : }
565 :
566 :
567 : FILE* fp_forces=NULL;
568 1328 : OFile fp_dforces;
569 :
570 : std::unique_ptr<FILE,decltype(deleter)> fp_forces_deleter(fp_forces,deleter);
571 :
572 1328 : if(!noatoms&&!parseOnly) {
573 1254 : if (trajectoryFile=="-") {
574 1 : auto errormessage=parser.init(trajectory_fmt, in);
575 0 : if(errormessage) {
576 0 : error(*errormessage);
577 : }
578 : } else {
579 1253 : if(multi) {
580 : std::string n;
581 300 : Tools::convert(intercomm.Get_rank(),n);
582 600 : std::string testfile=FileBase::appendSuffix(trajectoryFile,"."+n);
583 300 : FILE* tmp_fp=std::fopen(testfile.c_str(),"r");
584 : // no exceptions here
585 300 : if(tmp_fp) {
586 235 : std::fclose(tmp_fp);
587 : trajectoryFile=testfile;
588 : }
589 : }
590 1253 : auto errormessage=parser.init(trajectory_fmt,trajectoryFile,use_molfile,command_line_natoms);
591 1253 : if(errormessage) {
592 1 : error(*errormessage);
593 : }
594 1252 : natoms = parser.nOfAtoms();
595 : }
596 1252 : if(dumpforces.length()>0) {
597 695 : if(Communicator::initialized() && pc.Get_size()>1) {
598 : std::string n;
599 282 : Tools::convert(pc.Get_rank(),n);
600 564 : dumpforces+="."+n;
601 : }
602 695 : fp_forces=std::fopen(dumpforces.c_str(),"w");
603 695 : if(!fp_forces) {
604 71 : error("cannot open file "+dumpforces+" for writing");
605 : }
606 : fp_forces_deleter.reset(fp_forces);
607 : }
608 1251 : if(debugforces.length()>0) {
609 19 : if(Communicator::initialized() && pc.Get_size()>1) {
610 : std::string n;
611 6 : Tools::convert(pc.Get_rank(),n);
612 12 : debugforces+="."+n;
613 : }
614 19 : fp_dforces.open(debugforces);
615 : }
616 : }
617 :
618 : std::string line;
619 1325 : std::vector<real> coordinates;
620 1325 : std::vector<real> forces;
621 1325 : std::vector<real> masses;
622 1325 : std::vector<real> charges;
623 1325 : std::vector<real> cell;
624 1325 : std::vector<real> virial;
625 1325 : std::vector<real> numder;
626 :
627 : // variables to test particle decomposition
628 : int pd_nlocal=0;
629 : int pd_start=0;
630 : // variables to test random decomposition (=domain decomposition)
631 1325 : std::vector<int> dd_gatindex;
632 1325 : std::vector<int> dd_g2l;
633 1325 : std::vector<real> dd_masses;
634 1325 : std::vector<real> dd_charges;
635 1325 : std::vector<real> dd_forces;
636 1325 : std::vector<real> dd_coordinates;
637 : int dd_nlocal=0;
638 : // random stream to choose decompositions
639 1325 : Random rnd;
640 :
641 :
642 : bool lstep=true;
643 277988 : while(true) {
644 : bool first_step=false;
645 279313 : if(!noatoms&&!parseOnly) {
646 :
647 78711 : real timeStep=-1.0;
648 :
649 78711 : auto errormessage=parser.readHeader(step,timeStep);
650 :
651 78711 : if (errormessage) {
652 1182 : if (*errormessage =="EOF") {
653 : break;
654 : } else {
655 0 : error(*errormessage);
656 : }
657 : }
658 77529 : if (lstep && timeStep>0.0) {
659 2 : p.cmd("setTimestep",real(timeStep));
660 : lstep = false;
661 : }
662 77529 : natoms=parser.nOfAtoms();
663 : }
664 278131 : if(checknatoms<0 && !noatoms) {
665 : pd_nlocal=natoms;
666 : pd_start=0;
667 : first_step=true;
668 1268 : masses.assign(natoms,std::numeric_limits<real>::quiet_NaN());
669 1268 : charges.assign(natoms,std::numeric_limits<real>::quiet_NaN());
670 : //case pdb: structure
671 1268 : if(pdbfile.length()>0) {
672 150997 : for(unsigned i=0; i<pdb.size(); ++i) {
673 150917 : AtomNumber an=pdb.getAtomNumbers()[i];
674 : unsigned index=an.index();
675 150917 : if( index>=unsigned(natoms) ) {
676 66 : error("atom index in pdb exceeds the number of atoms in trajectory");
677 : }
678 150917 : masses[index]=pdb.getOccupancy()[i];
679 150917 : charges[index]=pdb.getBeta()[i];
680 : }
681 : }
682 1268 : if(mcfile.length()>0) {
683 17 : IFile ifile;
684 17 : ifile.open(mcfile);
685 : int index;
686 : double mass;
687 : double charge;
688 495756 : while(ifile.scanField("index",index).scanField("mass",mass).scanField("charge",charge).scanField()) {
689 247862 : if( index<0 || index>=natoms ) {
690 2 : error("atom index in file " + mcfile + " is not between 0 and the number of atoms in trajectory minus one");
691 : }
692 247861 : masses[index]=mass;
693 247861 : charges[index]=charge;
694 : }
695 17 : }
696 : } else if( checknatoms<0 && noatoms ) {
697 : natoms=0;
698 : }
699 :
700 278130 : if( checknatoms<0 ) {
701 1324 : if(kt>=0) {
702 5 : p.cmd("setKbT",kt);
703 : }
704 : checknatoms=natoms;
705 1324 : p.cmd("setNatoms",natoms);
706 1324 : p.cmd("init");
707 : // Check if we have been asked to output the long version of the input and if there are shortcuts
708 1264 : if( parseOnly && full_outputfile.length()>0 ) {
709 :
710 : // Read in the plumed input file and store what is in there
711 : std::map<std::string,std::vector<std::string> > data;
712 5 : IFile ifile;
713 5 : ifile.open(plumedFile);
714 5 : std::vector<std::string> words;
715 34 : while( Tools::getParsedLine(ifile,words) && !p.getEndPlumed() ) {
716 24 : p.readInputWords(words,false);
717 : Action* aa=p.getActionSet()[p.getActionSet().size()-1].get();
718 24 : plumed_assert(aa); // needed for following calls, see #1046
719 24 : ActionWithValue* av=aa->castToActionWithValue();
720 32 : if( av && aa->getDefaultString().length()>0 ) {
721 3 : std::vector<std::string> def;
722 6 : def.push_back( "defaults " + aa->getDefaultString() );
723 3 : data[ aa->getLabel() ] = def;
724 3 : }
725 24 : ActionShortcut* as=aa->castToActionShortcut();
726 24 : if( as ) {
727 20 : if( aa->getDefaultString().length()>0 ) {
728 5 : std::vector<std::string> def;
729 5 : def.push_back( "defaults " + aa->getDefaultString() );
730 5 : data[ as->getShortcutLabel() ] = def;
731 5 : }
732 20 : const auto & shortcut_commands = as->getSavedInputLines();
733 20 : if( shortcut_commands.size()==0 ) {
734 14 : continue;
735 : }
736 12 : if( data.find( as->getShortcutLabel() )!=data.end() ) {
737 16 : for(unsigned i=0; i<shortcut_commands.size(); ++i) {
738 14 : data[ as->getShortcutLabel() ].push_back( shortcut_commands[i] );
739 : }
740 : } else {
741 4 : data[ as->getShortcutLabel() ] = shortcut_commands;
742 : }
743 : }
744 : }
745 5 : ifile.close();
746 : // Only output the full version of the input file if there are shortcuts
747 5 : if( data.size()>0 && intracomm.Get_rank()==0 && intercomm.Get_rank()==0 ) {
748 5 : OFile long_file;
749 5 : long_file.open( full_outputfile );
750 5 : long_file.printf("{\n");
751 : bool firstpass=true;
752 17 : for(auto& x : data ) {
753 12 : if( !firstpass ) {
754 7 : long_file.printf(" },\n");
755 : }
756 12 : long_file.printf(" \"%s\" : {\n", x.first.c_str() );
757 12 : plumed_assert( x.second.size()>0 );
758 : unsigned sstart=0;
759 12 : if( x.second[0].find("defaults")!=std::string::npos ) {
760 : sstart=1;
761 16 : long_file.printf(" \"defaults\" : \"%s\"", x.second[0].substr( 9 ).c_str() );
762 8 : if( x.second.size()>1 ) {
763 2 : long_file.printf(",\n");
764 : } else {
765 6 : long_file.printf("\n");
766 : }
767 : }
768 12 : if( x.second.size()>sstart ) {
769 : {
770 6 : std::string cmnd=Tools::convertRegexForJson(x.second[sstart]);
771 6 : long_file.printf(" \"expansion\" : \"%s", cmnd.c_str() );
772 : }
773 36 : for(unsigned j=sstart+1; j<x.second.size(); ++j) {
774 : //reinitializing it to RVO the string
775 30 : std::string cmnd=Tools::convertRegexForJson(x.second[j]);
776 30 : long_file.printf("\\n%s", cmnd.c_str() );
777 : }
778 6 : long_file.printf("\"\n");
779 : }
780 : firstpass=false;
781 : }
782 5 : long_file.printf(" }\n}\n");
783 5 : long_file.close();
784 5 : }
785 5 : }
786 1264 : if( valuedict_file.length()>0 && intracomm.Get_rank()==0 && intercomm.Get_rank()==0 ) {
787 5 : OFile valuefile;
788 5 : valuefile.open( valuedict_file );
789 5 : valuefile.printf("{\n");
790 : bool firsta=true;
791 149 : for(const auto & pp : p.getActionSet()) {
792 144 : if( pp.get()->getName()=="CENTER" ) {
793 3 : continue ;
794 : }
795 141 : ActionWithVirtualAtom* avv=dynamic_cast<ActionWithVirtualAtom*>(pp.get());
796 141 : if( avv || pp.get()->getName()=="GROUP" || pp.get()->getName()=="DENSITY" ) {
797 : Action* actionBase=pp.get();
798 6 : if( firsta ) {
799 0 : valuefile.printf(" \"%s\" : {\n \"action\" : \"%s\"",
800 : actionBase->getLabel().c_str(),
801 : actionBase->getName().c_str() );
802 : firsta=false;
803 : } else {
804 6 : valuefile.printf(",\n \"%s\" : {\n \"action\" : \"%s\"",
805 : actionBase->getLabel().c_str(),
806 : actionBase->getName().c_str() );
807 : }
808 6 : if( avv ) {
809 3 : valuefile.printf(",\n \"%s\" : { \"type\": \"atoms\","
810 : " \"description\": \"virtual atom calculated by %s action\" }",
811 : avv->getLabel().c_str(),
812 3 : avv->getName().c_str() );
813 : } else {
814 3 : valuefile.printf(",\n \"%s\" : { \"type\": \"atoms\","
815 : " \"description\": \"indices of atoms specified in GROUP\" }",
816 : actionBase->getLabel().c_str() );
817 : }
818 6 : valuefile.printf("\n }");
819 6 : continue;
820 6 : }
821 135 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(pp.get());
822 135 : if( av && av->getNumberOfComponents()>0 ) {
823 89 : Keywords keys;
824 89 : p.getKeywordsForAction( av->getName(), keys );
825 89 : if( firsta ) {
826 10 : valuefile.printf(" \"%s\" : {\n \"action\" : \"%s\"",
827 5 : av->getLabel().c_str(),
828 : keys.getDisplayName().c_str() );
829 : firsta=false;
830 : } else {
831 168 : valuefile.printf(",\n \"%s\" : {\n \"action\" : \"%s\"",
832 84 : av->getLabel().c_str(),
833 : keys.getDisplayName().c_str() );
834 : }
835 207 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
836 118 : Value* myval = av->copyOutput(i);
837 118 : std::string compname = myval->getName(), description;
838 118 : if( av->getLabel()==compname ) {
839 150 : description = keys.getOutputComponentDescription(".#!value");
840 : } else {
841 43 : std::size_t dot=compname.find(av->getLabel() + ".");
842 43 : std::string cname = compname.substr(dot + av->getLabel().length() + 1);
843 86 : description = av->getOutputComponentDescription( cname, keys );
844 : }
845 118 : if( description.find("\\")!=std::string::npos ) {
846 0 : error("found invalid backslash character in documentation for component " + compname + " in action " + av->getName() + " with label " + av->getLabel() );
847 : }
848 236 : valuefile.printf(",\n \"%s\" : { \"type\": \"%s\", \"description\": \"%s\" }", myval->getName().c_str(), myval->getValueType().c_str(), description.c_str() );
849 : }
850 89 : valuefile.printf("\n }");
851 89 : }
852 135 : ActionShortcut* as=pp->castToActionShortcut();
853 135 : if( as ) {
854 41 : const auto & cnames = as->getSavedOutputs();
855 41 : if( cnames.size()==0 ) {
856 34 : continue ;
857 : }
858 :
859 7 : if( firsta ) {
860 0 : valuefile.printf(" \"shortcut_%s\" : {\n \"action\" : \"%s\"", as->getShortcutLabel().c_str(), as->getName().c_str() );
861 : firsta=false;
862 : } else {
863 7 : valuefile.printf(",\n \"shortcut_%s\" : {\n \"action\" : \"%s\"", as->getShortcutLabel().c_str(), as->getName().c_str() );
864 : }
865 7 : Keywords keys;
866 7 : p.getKeywordsForAction( as->getName(), keys );
867 18 : for(unsigned i=0; i<cnames.size(); ++i) {
868 11 : ActionWithValue* av2=p.getActionSet().selectWithLabel<ActionWithValue*>( cnames[i] );
869 11 : if( !av2 ) {
870 0 : plumed_merror("could not find value created by shortcut with name " + cnames[i] );
871 : }
872 11 : if( av2->getNumberOfComponents()==1 ) {
873 11 : Value* myval = av2->copyOutput(0);
874 11 : std::string compname = myval->getName(), description;
875 11 : if( compname==as->getShortcutLabel() ) {
876 14 : description = keys.getOutputComponentDescription(".#!value");
877 : } else {
878 4 : std::size_t shortcutPos=compname.find(as->getShortcutLabel());
879 8 : description = keys.getOutputComponentDescription(
880 4 : compname.substr(shortcutPos+as->getShortcutLabel().length()+1) );
881 : }
882 11 : if( description.find("\\")!=std::string::npos ) {
883 0 : error("found invalid backslash character in documentation for component " + compname + " in action " + as->getName() + " with label " + as->getLabel() );
884 : }
885 22 : valuefile.printf(",\n \"%s\" : { \"type\": \"%s\", \"description\": \"%s\" }",
886 : myval->getName().c_str(),
887 : myval->getValueType().c_str(),
888 : description.c_str() );
889 : } else {
890 0 : for(unsigned j=0; j<av2->getNumberOfComponents(); ++j) {
891 0 : Value* myval = av2->copyOutput(j);
892 0 : std::string compname = myval->getName(), description;
893 0 : if( av2->getLabel()==compname ) {
894 0 : plumed_merror("should not be outputting description of value from action when using shortcuts");
895 : } else {
896 0 : std::size_t dot=compname.find(av2->getLabel() + ".");
897 0 : std::string cname = compname.substr(dot+av2->getLabel().length() + 1);
898 0 : description = av2->getOutputComponentDescription( cname, keys );
899 : }
900 0 : if( description.find("\\")!=std::string::npos ) {
901 0 : error("found invalid backslash character in documentation for component " + compname + " in action " + av2->getName() + " with label " + av2->getLabel() );
902 : }
903 0 : valuefile.printf(",\n \"%s\" : { \"type\": \"%s\", \"description\": \"%s\" }", myval->getName().c_str(), myval->getValueType().c_str(), description.c_str() );
904 : }
905 : }
906 : }
907 7 : valuefile.printf("\n }");
908 7 : }
909 : }
910 5 : valuefile.printf("\n}\n");
911 5 : valuefile.close();
912 5 : }
913 1264 : if(parseOnly) {
914 : break;
915 : }
916 : }
917 278064 : if(checknatoms!=natoms) {
918 : std::string stepstr;
919 0 : Tools::convert(step,stepstr);
920 0 : error("number of atoms in frame " + stepstr + " does not match number of atoms in first frame");
921 : }
922 :
923 278064 : coordinates.assign(3*natoms,real(0.0));
924 278064 : forces.assign(3*natoms,real(0.0));
925 278064 : cell.assign(9,real(0.0));
926 278064 : virial.assign(9,real(0.0));
927 :
928 278064 : if( first_step || rnd.U01()>0.5) {
929 141084 : if(debug_pd) {
930 152 : int npe=intracomm.Get_size();
931 152 : std::vector<int> loc(npe,0);
932 152 : std::vector<int> start(npe,0);
933 608 : for(int i=0; i<npe-1; i++) {
934 456 : int cc=(natoms*2*rnd.U01())/npe;
935 456 : if(start[i]+cc>natoms) {
936 16 : cc=natoms-start[i];
937 : }
938 456 : loc[i]=cc;
939 456 : start[i+1]=start[i]+loc[i];
940 : }
941 152 : loc[npe-1]=natoms-start[npe-1];
942 152 : intracomm.Bcast(loc,0);
943 152 : intracomm.Bcast(start,0);
944 152 : pd_nlocal=loc[intracomm.Get_rank()];
945 152 : pd_start=start[intracomm.Get_rank()];
946 152 : if(intracomm.Get_rank()==0) {
947 : std::fprintf(out,"\nDRIVER: Reassigning particle decomposition\n");
948 : std::fprintf(out,"DRIVER: ");
949 190 : for(int i=0; i<npe; i++) {
950 152 : std::fprintf(out,"%d ",loc[i]);
951 : }
952 : std::fprintf(out,"\n");
953 : std::fprintf(out,"DRIVER: ");
954 190 : for(int i=0; i<npe; i++) {
955 152 : std::fprintf(out,"%d ",start[i]);
956 : }
957 : std::fprintf(out,"\n");
958 : }
959 152 : p.cmd("setAtomsNlocal",pd_nlocal);
960 152 : p.cmd("setAtomsContiguous",pd_start);
961 141084 : } else if(debug_dd) {
962 1100 : int npe=intracomm.Get_size();
963 1100 : int rank=intracomm.Get_rank();
964 1100 : dd_charges.assign(natoms,0.0);
965 1100 : dd_masses.assign(natoms,0.0);
966 1100 : dd_gatindex.assign(natoms,-1);
967 1100 : dd_g2l.assign(natoms,-1);
968 1100 : dd_coordinates.assign(3*natoms,0.0);
969 1100 : dd_forces.assign(3*natoms,0.0);
970 : dd_nlocal=0;
971 56098 : for(int i=0; i<natoms; ++i) {
972 54998 : double r=rnd.U01()*npe;
973 : int n;
974 115728 : for(n=0; n<npe; n++)
975 115728 : if(n+1>r) {
976 : break;
977 : }
978 54998 : plumed_assert(n<npe);
979 54998 : if(n==rank) {
980 20911 : dd_gatindex[dd_nlocal]=i;
981 20911 : dd_g2l[i]=dd_nlocal;
982 20911 : dd_charges[dd_nlocal]=charges[i];
983 20911 : dd_masses[dd_nlocal]=masses[i];
984 20911 : dd_nlocal++;
985 : }
986 : }
987 1100 : if(intracomm.Get_rank()==0) {
988 : std::fprintf(out,"\nDRIVER: Reassigning domain decomposition\n");
989 : }
990 1100 : p.cmd("setAtomsNlocal",dd_nlocal);
991 1100 : p.cmd("setAtomsGatindex",&dd_gatindex[0], {dd_nlocal});
992 : }
993 : }
994 :
995 278064 : int plumedStopCondition=0;
996 278064 : if(!noatoms) {
997 77479 : auto errormessage=parser.readAtoms(
998 : stride,
999 : pbc_cli_given,
1000 : debug_pd,
1001 : pd_start,
1002 : pd_nlocal,
1003 : step,
1004 : masses.data(),
1005 : charges.data(),
1006 : coordinates.data(),
1007 : cell.data()
1008 : );
1009 :
1010 77479 : if (errormessage) {
1011 11 : if (*errormessage =="EOF") {
1012 : break;
1013 : } else {
1014 0 : error(*errormessage);
1015 : }
1016 : }
1017 :
1018 77468 : if(pbc_cli_given) {
1019 4270 : for(unsigned i=0; i<9; i++) {
1020 3843 : cell[i]=real(pbc_cli_box[i]);
1021 : }
1022 : }
1023 77468 : p.cmd("setStepLongLong",step);
1024 77468 : p.cmd("setStopFlag",&plumedStopCondition);
1025 :
1026 77468 : if(debug_dd) {
1027 40880 : for(int i=0; i<dd_nlocal; ++i) {
1028 38856 : int kk=dd_gatindex[i];
1029 38856 : dd_coordinates[3*i+0]=coordinates[3*kk+0];
1030 38856 : dd_coordinates[3*i+1]=coordinates[3*kk+1];
1031 38856 : dd_coordinates[3*i+2]=coordinates[3*kk+2];
1032 : // masses and charges might have been read in readAtoms (e.g., dlp4)
1033 38856 : dd_masses[i]=masses[kk];
1034 38856 : dd_charges[i]=charges[kk];
1035 : }
1036 2024 : p.cmd("setForces",&dd_forces[0], {dd_nlocal,3});
1037 2024 : p.cmd("setPositions",&dd_coordinates[0], {dd_nlocal,3});
1038 2024 : p.cmd("setMasses",&dd_masses[0], {dd_nlocal});
1039 2024 : p.cmd("setCharges",&dd_charges[0], {dd_nlocal});
1040 : } else {
1041 : // this is required to avoid troubles when the last domain
1042 : // contains zero atoms
1043 : // Basically, for empty domains we pass null pointers
1044 : #define fix_pd(xx) (pd_nlocal!=0?&xx:NULL)
1045 150888 : p.cmd("setForces",fix_pd(forces[3*pd_start]), {pd_nlocal,3});
1046 150888 : p.cmd("setPositions",fix_pd(coordinates[3*pd_start]), {pd_nlocal,3});
1047 150888 : p.cmd("setMasses",fix_pd(masses[pd_start]), {pd_nlocal});
1048 150888 : p.cmd("setCharges",fix_pd(charges[pd_start]), {pd_nlocal});
1049 : }
1050 77468 : p.cmd("setBox",cell.data(), {3,3});
1051 77468 : p.cmd("setVirial",virial.data(), {3,3});
1052 : } else {
1053 200585 : p.cmd("setStepLongLong",step);
1054 200585 : p.cmd("setStopFlag",&plumedStopCondition);
1055 : }
1056 278053 : p.cmd("calc");
1057 278049 : if(debugforces.length()>0) {
1058 163 : virial.assign(9,real(0.0));
1059 163 : forces.assign(3*natoms,real(0.0));
1060 163 : p.cmd("prepareCalc");
1061 163 : p.cmd("performCalcNoUpdate");
1062 : }
1063 :
1064 : // this is necessary as only processor zero is adding to the virial:
1065 278049 : intracomm.Bcast(virial,0);
1066 278049 : if(debug_pd) {
1067 240 : intracomm.Sum(forces);
1068 : }
1069 278049 : if(debug_dd) {
1070 40880 : for(int i=0; i<dd_nlocal; i++) {
1071 38856 : forces[3*dd_gatindex[i]+0]=dd_forces[3*i+0];
1072 38856 : forces[3*dd_gatindex[i]+1]=dd_forces[3*i+1];
1073 38856 : forces[3*dd_gatindex[i]+2]=dd_forces[3*i+2];
1074 : }
1075 2024 : dd_forces.assign(3*natoms,0.0);
1076 2024 : intracomm.Sum(forces);
1077 : }
1078 278049 : if(debug_grex &&step%grex_stride==0) {
1079 114 : p.cmd("GREX savePositions");
1080 114 : if(intracomm.Get_rank()>0) {
1081 57 : p.cmd("GREX prepare");
1082 : } else {
1083 57 : int r=intercomm.Get_rank();
1084 57 : int n=intercomm.Get_size();
1085 57 : int partner=r+(2*((r+step/grex_stride)%2))-1;
1086 57 : if(partner<0) {
1087 : partner=0;
1088 : }
1089 57 : if(partner>=n) {
1090 8 : partner=n-1;
1091 : }
1092 57 : p.cmd("GREX setPartner",partner);
1093 57 : p.cmd("GREX calculate");
1094 57 : p.cmd("GREX shareAllDeltaBias");
1095 228 : for(int i=0; i<n; i++) {
1096 : std::string s;
1097 171 : Tools::convert(i,s);
1098 171 : real a=std::numeric_limits<real>::quiet_NaN();
1099 342 : s="GREX getDeltaBias "+s;
1100 171 : p.cmd(s,&a);
1101 171 : if(grex_log) {
1102 171 : std::fprintf(grex_log," %f",a);
1103 : }
1104 : }
1105 57 : if(grex_log) {
1106 : std::fprintf(grex_log,"\n");
1107 : }
1108 : }
1109 : }
1110 :
1111 :
1112 278049 : if(fp_forces) {
1113 : std::fprintf(fp_forces,"%d\n",natoms);
1114 321880 : std::string fmtv=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
1115 80470 : std::string fmt=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
1116 40235 : if(dumpfullvirial) {
1117 441 : std::fprintf(fp_forces,fmtv.c_str(),virial[0],virial[1],virial[2],virial[3],virial[4],virial[5],virial[6],virial[7],virial[8]);
1118 : } else {
1119 39794 : std::fprintf(fp_forces,fmt.c_str(),virial[0],virial[4],virial[8]);
1120 : }
1121 40235 : fmt="X "+fmt;
1122 5088484 : for(int i=0; i<natoms; i++) {
1123 5048249 : std::fprintf(fp_forces,fmt.c_str(),forces[3*i],forces[3*i+1],forces[3*i+2]);
1124 : }
1125 : }
1126 278049 : if(debugforces.length()>0) {
1127 : // Now call the routine to work out the derivatives numerically
1128 163 : numder.assign(3*natoms+9,real(0.0));
1129 163 : real base=0;
1130 163 : p.cmd("getBias",&base);
1131 163 : if( std::abs(base)<epsilon ) {
1132 : std::fprintf(out,"WARNING: bias for configuration appears to be zero so debugging forces is trivial\n");
1133 : }
1134 163 : evaluateNumericalDerivatives( step, p, coordinates, masses, charges, cell, base, numder );
1135 :
1136 : // And output everything to a file
1137 163 : fp_dforces.fmtField(" " + dumpforcesFmt);
1138 9772 : for(int i=0; i<3*natoms; ++i) {
1139 9609 : fp_dforces.printField("parameter",(int)i);
1140 19218 : fp_dforces.printField("analytical",forces[i]);
1141 9609 : fp_dforces.printField("numerical",-numder[i]);
1142 9609 : fp_dforces.printField();
1143 : }
1144 : // And print the virial
1145 1630 : for(int i=0; i<9; ++i) {
1146 1467 : fp_dforces.printField("parameter",3*natoms+i );
1147 1467 : fp_dforces.printField("analytical",virial[i] );
1148 1467 : fp_dforces.printField("numerical",-numder[3*natoms+i]);
1149 1467 : fp_dforces.printField();
1150 : }
1151 : }
1152 :
1153 278049 : if(plumedStopCondition) {
1154 : break;
1155 : }
1156 :
1157 277988 : step+=stride;
1158 : }
1159 1260 : if(!parseOnly) {
1160 1319 : p.cmd("runFinalJobs");
1161 : }
1162 :
1163 : return 0;
1164 5690 : }
1165 :
1166 : template<typename real>
1167 163 : void Driver<real>::evaluateNumericalDerivatives( const long long int& step, PlumedMain& p, const std::vector<real>& coordinates,
1168 : const std::vector<real>& masses, const std::vector<real>& charges,
1169 : std::vector<real>& cell, const double& base, std::vector<real>& numder ) {
1170 :
1171 163 : int natoms = coordinates.size() / 3;
1172 : real delta = std::sqrt(epsilon);
1173 163 : std::vector<Vector> pos(natoms);
1174 163 : real bias=0;
1175 163 : std::vector<real> fake_forces( 3*natoms ), fake_virial(9);
1176 3366 : for(int i=0; i<natoms; ++i) {
1177 12812 : for(unsigned j=0; j<3; ++j) {
1178 9609 : pos[i][j]=coordinates[3*i+j];
1179 : }
1180 : }
1181 :
1182 3366 : for(int i=0; i<natoms; ++i) {
1183 12812 : for(unsigned j=0; j<3; ++j) {
1184 9609 : pos[i][j]=pos[i][j]+delta;
1185 9609 : p.cmd("setStepLongLong",step);
1186 9609 : p.cmd("setPositions",&pos[0][0], {natoms,3});
1187 9609 : p.cmd("setForces",&fake_forces[0], {natoms,3});
1188 9609 : p.cmd("setMasses",&masses[0], {natoms});
1189 9609 : p.cmd("setCharges",&charges[0], {natoms});
1190 9609 : p.cmd("setBox",&cell[0], {3,3});
1191 9609 : p.cmd("setVirial",&fake_virial[0], {3,3});
1192 9609 : p.cmd("prepareCalc");
1193 9609 : p.cmd("performCalcNoUpdate");
1194 9609 : p.cmd("getBias",&bias);
1195 9609 : pos[i][j]=coordinates[3*i+j];
1196 9609 : numder[3*i+j] = (bias - base) / delta;
1197 : }
1198 : }
1199 :
1200 : // Create the cell
1201 163 : Tensor box( cell[0], cell[1], cell[2], cell[3], cell[4], cell[5], cell[6], cell[7], cell[8] );
1202 : // And the virial
1203 163 : Pbc pbc;
1204 163 : pbc.setBox( box );
1205 163 : Tensor nvirial;
1206 652 : for(unsigned i=0; i<3; i++)
1207 1956 : for(unsigned k=0; k<3; k++) {
1208 1467 : double arg0=box(i,k);
1209 30294 : for(int j=0; j<natoms; ++j) {
1210 28827 : pos[j]=pbc.realToScaled( pos[j] );
1211 : }
1212 1467 : cell[3*i+k]=box(i,k)=box(i,k)+delta;
1213 1467 : pbc.setBox(box);
1214 30294 : for(int j=0; j<natoms; j++) {
1215 28827 : pos[j]=pbc.scaledToReal( pos[j] );
1216 : }
1217 1467 : p.cmd("setStepLongLong",step);
1218 1467 : p.cmd("setPositions",&pos[0][0], {natoms,3});
1219 1467 : p.cmd("setForces",&fake_forces[0], {natoms,3});
1220 1467 : p.cmd("setMasses",&masses[0], {natoms});
1221 1467 : p.cmd("setCharges",&charges[0], {natoms});
1222 1467 : p.cmd("setBox",&cell[0], {3,3});
1223 1467 : p.cmd("setVirial",&fake_virial[0], {3,3});
1224 1467 : p.cmd("prepareCalc");
1225 1467 : p.cmd("performCalcNoUpdate");
1226 1467 : p.cmd("getBias",&bias);
1227 1467 : cell[3*i+k]=box(i,k)=arg0;
1228 1467 : pbc.setBox(box);
1229 30294 : for(int j=0; j<natoms; j++)
1230 115308 : for(unsigned n=0; n<3; ++n) {
1231 86481 : pos[j][n]=coordinates[3*j+n];
1232 : }
1233 1467 : nvirial(i,k) = ( bias - base ) / delta;
1234 : }
1235 163 : nvirial=-matmul(box.transpose(),nvirial);
1236 652 : for(unsigned i=0; i<3; i++)
1237 1956 : for(unsigned k=0; k<3; k++) {
1238 1467 : numder[3*natoms+3*i+k] = nvirial(i,k);
1239 : }
1240 :
1241 163 : }
1242 :
1243 : }
1244 : }
|
