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 "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/ActionShortcut.h"
29 : #include "tools/Communicator.h"
30 : #include "tools/Random.h"
31 : #include "tools/Pbc.h"
32 : #include <cstdio>
33 : #include <cstring>
34 : #include <vector>
35 : #include <map>
36 : #include <memory>
37 : #include "tools/Units.h"
38 : #include "tools/PDB.h"
39 : #include "tools/FileBase.h"
40 : #include "tools/IFile.h"
41 : #include "xdrfile/xdrfile_trr.h"
42 : #include "xdrfile/xdrfile_xtc.h"
43 :
44 :
45 : // when using molfile plugin
46 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
47 : #ifndef __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS
48 : /* Use the internal ones. Alternatively:
49 : * ifeq (,$(findstring __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS,$(CPPFLAGS)))
50 : * CPPFLAGS+=-I../molfile
51 : */
52 : #include "molfile/libmolfile_plugin.h"
53 : #include "molfile/molfile_plugin.h"
54 : using namespace PLMD::molfile;
55 : #else
56 : #include <libmolfile_plugin.h>
57 : #include <molfile_plugin.h>
58 : #endif
59 : #endif
60 :
61 : namespace PLMD {
62 : namespace cltools {
63 :
64 : //+PLUMEDOC TOOLS driver-float
65 : /*
66 : Equivalent to driver, but using single precision reals.
67 :
68 : The purpose of this tool is just to test what PLUMED does when linked from
69 : a single precision code. The documentation is identical to that for \ref driver
70 :
71 : \par Examples
72 :
73 : \verbatim
74 : plumed driver-float --plumed plumed.dat --ixyz trajectory.xyz
75 : \endverbatim
76 :
77 : See also examples in \ref driver
78 :
79 : */
80 : //+ENDPLUMEDOC
81 : //
82 :
83 :
84 : //+PLUMEDOC TOOLS driver
85 : /*
86 : driver is a tool that allows one to to use plumed to post-process an existing trajectory.
87 :
88 : The input to driver is specified using the command line arguments described below.
89 :
90 : In addition, you can use the special \subpage READ command inside your plumed input
91 : to read in colvar files that were generated during your MD simulation. The values
92 : read in can then be treated like calculated colvars.
93 :
94 : \warning
95 : Notice that by default the driver has no knowledge about the masses and charges
96 : of your atoms! Thus, if you want to compute quantities depending charges (e.g. \ref DHENERGY)
97 : or masses (e.g. \ref COM) you should pass the proper information to the driver.
98 : You can do it either with the --pdb option or with the --mc option. The latter
99 : will read a file produced by \ref DUMPMASSCHARGE .
100 :
101 :
102 : \par Examples
103 :
104 : The following command tells plumed to post process the trajectory contained in `trajectory.xyz`
105 : by performing the actions described in the input file `plumed.dat`. If an action that takes the
106 : stride keyword is given a stride equal to \f$n\f$ then it will be performed only on every \f$n\f$th
107 : frames in the trajectory file.
108 : \verbatim
109 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz
110 : \endverbatim
111 :
112 : Notice that `xyz` files are expected to be in internal PLUMED units, that is by default nm.
113 : You can change this behavior by using the `--length-units` option:
114 : \verbatim
115 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --length-units A
116 : \endverbatim
117 : The strings accepted by the `--length-units` options are the same ones accepted by the \ref UNITS action.
118 : Other file formats typically have their default coordinates (e.g., `gro` files are always in nm)
119 : and it thus should not be necessary to use the `--length-units` option. Additionally,
120 : consider that the units used by the `driver` might be different by the units used in the PLUMED input
121 : file `plumed.dat`. For instance consider the command:
122 : \verbatim
123 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --length-units A
124 : \endverbatim
125 : where `plumed.dat` is
126 : \plumedfile
127 : # no explicit UNITS action here
128 : d: DISTANCE ATOMS=1,2
129 : PRINT ARG=d FILE=colvar
130 : \endplumedfile
131 : In this case, the driver reads the `xyz` file assuming it to contain coordinates in Angstrom units.
132 : However, the resulting `colvar` file contains a distance expressed in nm.
133 :
134 : The following command tells plumed to post process the trajectory contained in trajectory.xyz.
135 : by performing the actions described in the input file plumed.dat.
136 : \verbatim
137 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --trajectory-stride 100 --timestep 0.001
138 : \endverbatim
139 : Here though
140 : `--trajectory-stride` is set equal to the frequency with which frames were output during the trajectory
141 : and the `--timestep` is equal to the simulation timestep. As such the `STRIDE` parameters in the `plumed.dat`
142 : files are referred to the original timestep and any files output resemble those that would have been generated
143 : had we run the calculation we are running with driver when the MD simulation was running.
144 :
145 : PLUMED can read xyz files (in PLUMED units) and gro files (in nm). In addition,
146 : PLUMED includes by default support for a
147 : subset of the trajectory file formats supported by VMD, e.g. xtc and dcd:
148 :
149 : \verbatim
150 : plumed driver --plumed plumed.dat --pdb diala.pdb --mf_xtc traj.xtc --trajectory-stride 100 --timestep 0.001
151 : \endverbatim
152 :
153 : where `--mf_` prefixes the extension of one of the accepted molfile plugin format.
154 : If PLUMED has been \ref Installation "installed" with full molfile support, other formats will be available.
155 : Just type `plumed driver --help` to see which plugins are available.
156 :
157 : Molfile plugin require periodic cell to be triangular (i.e. first vector oriented along x and
158 : second vector in xy plane). This is true for many MD codes. However, it could be false
159 : if you rotate the coordinates in your trajectory before reading them in the driver.
160 : Also notice that some formats (e.g. amber crd) do not specify atom number. In this case you can use
161 : the `--natoms` option:
162 : \verbatim
163 : plumed driver --plumed plumed.dat --imf_crd trajectory.crd --natoms 128
164 : \endverbatim
165 :
166 : Check the available molfile plugins and limitations at [this link](http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/).
167 :
168 : Additionally, you can use the xdrfile implementation of xtc and trr. To this aim, just
169 : download and install properly the xdrfile library (see [this link](http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library)).
170 : If the xdrfile library is installed properly the PLUMED configure script should be able to
171 : detect it and enable it.
172 : Notice that the xdrfile implementation of xtc and trr
173 : is more robust than the molfile one, since it provides support for generic cell shapes.
174 : In addition, it allows \ref DUMPATOMS to write compressed xtc files.
175 :
176 : \par Multiple replicas
177 :
178 : When PLUMED is compiled with MPI support, you can emulate a multi-simulation setup with the `driver` by providing the `--multi`
179 : option with the appropriate number of ranks. This allows to use the \ref special-replica-syntax and to analyze multiple
180 : trajectories (see \ref trieste-5-replicas). PLUMED will automatically append a numbered suffix to output files
181 : (e.g. `COLVAR.0`, `COLVAR.1`, …). Similarly, each replica will search for the corresponding suffixed input file (e.g. `traj.0.xtc`, …)
182 : or default to the unsuffixed one.
183 :
184 : */
185 : //+ENDPLUMEDOC
186 : //
187 :
188 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
189 : static std::vector<molfile_plugin_t *> plugins;
190 : static std::map <std::string, unsigned> pluginmap;
191 82710 : static int register_cb(void *v, vmdplugin_t *p) {
192 : //const char *key = p->name;
193 165420 : const auto ret = pluginmap.insert ( std::pair<std::string,unsigned>(std::string(p->name),plugins.size()) );
194 82710 : if (ret.second==false) {
195 : //cerr<<"MOLFILE: found duplicate plugin for "<<key<<" : not inserted "<<endl;
196 : } else {
197 : //cerr<<"MOLFILE: loading plugin "<<key<<" number "<<plugins.size()-1<<endl;
198 82710 : plugins.push_back(reinterpret_cast<molfile_plugin_t *>(p));
199 : }
200 82710 : return VMDPLUGIN_SUCCESS;
201 : }
202 : #endif
203 :
204 : template<typename real>
205 : class Driver : public CLTool {
206 : public:
207 : static void registerKeywords( Keywords& keys );
208 : explicit Driver(const CLToolOptions& co );
209 : int main(FILE* in,FILE*out,Communicator& pc) override;
210 : void evaluateNumericalDerivatives( const long long int& step, PlumedMain& p, const std::vector<real>& coordinates,
211 : const std::vector<real>& masses, const std::vector<real>& charges,
212 : std::vector<real>& cell, const double& base, std::vector<real>& numder );
213 : std::string description()const override;
214 : };
215 :
216 : template<typename real>
217 9190 : void Driver<real>::registerKeywords( Keywords& keys ) {
218 9190 : CLTool::registerKeywords( keys );
219 : keys.isDriver();
220 18380 : keys.addFlag("--help-debug",false,"print special options that can be used to create regtests");
221 18380 : keys.add("compulsory","--plumed","plumed.dat","specify the name of the plumed input file");
222 18380 : keys.add("compulsory","--timestep","1.0","the timestep that was used in the calculation that produced this trajectory in picoseconds");
223 18380 : keys.add("compulsory","--trajectory-stride","1","the frequency with which frames were output to this trajectory during the simulation"
224 : " (0 means that the number of the step is read from the trajectory file,"
225 : " currently working only for xtc/trr files read with --ixtc/--trr)"
226 : );
227 18380 : keys.add("compulsory","--multi","0","set number of replicas for multi environment (needs MPI)");
228 18380 : keys.addFlag("--noatoms",false,"don't read in a trajectory. Just use colvar files as specified in plumed.dat");
229 18380 : keys.addFlag("--parse-only",false,"read the plumed input file and stop");
230 18380 : keys.addFlag("--restart",false,"makes driver behave as if restarting");
231 18380 : keys.add("atoms","--ixyz","the trajectory in xyz format");
232 18380 : keys.add("atoms","--igro","the trajectory in gro format");
233 18380 : keys.add("atoms","--idlp4","the trajectory in DL_POLY_4 format");
234 18380 : keys.add("atoms","--ixtc","the trajectory in xtc format (xdrfile implementation)");
235 18380 : keys.add("atoms","--itrr","the trajectory in trr format (xdrfile implementation)");
236 18380 : 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");
237 18380 : keys.add("optional","--length-units","units for length, either as a string or a number");
238 18380 : keys.add("optional","--mass-units","units for mass in pdb and mc file, either as a string or a number");
239 18380 : keys.add("optional","--charge-units","units for charge in pdb and mc file, either as a string or a number");
240 18380 : keys.add("optional","--kt","set \\f$k_B T\\f$, it will not be necessary to specify temperature in input file");
241 18380 : keys.add("optional","--dump-forces","dump the forces on a file");
242 18380 : keys.add("optional","--dump-forces-fmt","( default=%%f ) the format to use to dump the forces");
243 18380 : keys.addFlag("--dump-full-virial",false,"with --dump-forces, it dumps the 9 components of the virial");
244 18380 : keys.add("optional","--pdb","provides a pdb with masses and charges");
245 18380 : keys.add("optional","--mc","provides a file with masses and charges as produced with DUMPMASSCHARGE");
246 18380 : keys.add("optional","--box","comma-separated box dimensions (3 for orthorhombic, 9 for generic)");
247 18380 : keys.add("optional","--natoms","provides number of atoms - only used if file format does not contain number of atoms");
248 18380 : keys.add("optional","--initial-step","provides a number for the initial step, default is 0");
249 18380 : keys.add("optional","--debug-forces","output a file containing the forces due to the bias evaluated using numerical derivatives "
250 : "and using the analytical derivatives implemented in plumed");
251 18380 : keys.add("hidden","--debug-float","[yes/no] turns on the single precision version (to check float interface)");
252 18380 : keys.add("hidden","--debug-dd","[yes/no] use a fake domain decomposition");
253 18380 : keys.add("hidden","--debug-pd","[yes/no] use a fake particle decomposition");
254 18380 : keys.add("hidden","--debug-grex","use a fake gromacs-like replica exchange, specify exchange stride");
255 18380 : keys.add("hidden","--debug-grex-log","log file for debug=grex");
256 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
257 9190 : MOLFILE_INIT_ALL
258 9190 : MOLFILE_REGISTER_ALL(NULL, register_cb)
259 91900 : for(unsigned i=0; i<plugins.size(); i++) {
260 165420 : std::string kk="--mf_"+std::string(plugins[i]->name);
261 165420 : std::string mm=" molfile: the trajectory in "+std::string(plugins[i]->name)+" format " ;
262 165420 : keys.add("atoms",kk,mm);
263 : }
264 : #endif
265 9190 : }
266 : template<typename real>
267 839 : Driver<real>::Driver(const CLToolOptions& co ):
268 839 : CLTool(co) {
269 839 : inputdata=commandline;
270 839 : }
271 : template<typename real>
272 4 : std::string Driver<real>::description()const {
273 4 : return "analyze trajectories with plumed";
274 : }
275 :
276 : template<typename real>
277 831 : int Driver<real>::main(FILE* in,FILE*out,Communicator& pc) {
278 :
279 831 : Units units;
280 831 : PDB pdb;
281 :
282 : // Parse everything
283 : bool printhelpdebug;
284 831 : parseFlag("--help-debug",printhelpdebug);
285 831 : if( printhelpdebug ) {
286 : std::fprintf(out,"%s",
287 : "Additional options for debug (only to be used in regtest):\n"
288 : " [--debug-float yes] : turns on the single precision version (to check float interface)\n"
289 : " [--debug-dd yes] : use a fake domain decomposition\n"
290 : " [--debug-pd yes] : use a fake particle decomposition\n"
291 : );
292 : return 0;
293 : }
294 : // Are we reading trajectory data
295 : bool noatoms;
296 831 : parseFlag("--noatoms",noatoms);
297 : bool parseOnly;
298 1662 : parseFlag("--parse-only",parseOnly);
299 : std::string full_outputfile;
300 831 : parse("--shortcut-ofile",full_outputfile);
301 : bool restart;
302 1662 : parseFlag("--restart",restart);
303 :
304 : std::string fakein;
305 : bool debug_float=false;
306 : fakein="";
307 1662 : if(parse("--debug-float",fakein)) {
308 0 : if(fakein=="yes") {
309 : debug_float=true;
310 0 : } else if(fakein=="no") {
311 : debug_float=false;
312 : } else {
313 0 : error("--debug-float should have argument yes or no");
314 : }
315 : }
316 :
317 : if(debug_float && sizeof(real)!=sizeof(float)) {
318 0 : auto cl=cltoolRegister().create(CLToolOptions("driver-float"));
319 : cl->setInputData(this->getInputData());
320 0 : int ret=cl->main(in,out,pc);
321 : return ret;
322 0 : }
323 :
324 : bool debug_pd=false;
325 : fakein="";
326 1662 : if(parse("--debug-pd",fakein)) {
327 12 : if(fakein=="yes") {
328 : debug_pd=true;
329 0 : } else if(fakein=="no") {
330 : debug_pd=false;
331 : } else {
332 0 : error("--debug-pd should have argument yes or no");
333 : }
334 : }
335 : if(debug_pd) {
336 : std::fprintf(out,"DEBUGGING PARTICLE DECOMPOSITION\n");
337 : }
338 :
339 : bool debug_dd=false;
340 : fakein="";
341 1662 : if(parse("--debug-dd",fakein)) {
342 60 : if(fakein=="yes") {
343 : debug_dd=true;
344 0 : } else if(fakein=="no") {
345 : debug_dd=false;
346 : } else {
347 0 : error("--debug-dd should have argument yes or no");
348 : }
349 : }
350 : if(debug_dd) {
351 : std::fprintf(out,"DEBUGGING DOMAIN DECOMPOSITION\n");
352 : }
353 :
354 831 : if( debug_pd || debug_dd ) {
355 72 : if(noatoms) {
356 0 : error("cannot debug without atoms");
357 : }
358 : }
359 :
360 : // set up for multi replica driver:
361 831 : int multi=0;
362 831 : parse("--multi",multi);
363 831 : Communicator intracomm;
364 831 : Communicator intercomm;
365 831 : if(multi) {
366 174 : int ntot=pc.Get_size();
367 174 : int nintra=ntot/multi;
368 174 : if(multi*nintra!=ntot) {
369 0 : error("invalid number of processes for multi environment");
370 : }
371 174 : pc.Split(pc.Get_rank()/nintra,pc.Get_rank(),intracomm);
372 174 : pc.Split(pc.Get_rank()%nintra,pc.Get_rank(),intercomm);
373 : } else {
374 657 : intracomm.Set_comm(pc.Get_comm());
375 : }
376 :
377 : // set up for debug replica exchange:
378 831 : bool debug_grex=parse("--debug-grex",fakein);
379 831 : int grex_stride=0;
380 : FILE*grex_log=NULL;
381 : // call fclose when fp goes out of scope
382 993 : auto deleter=[](FILE* f) {
383 : if(f) {
384 993 : std::fclose(f);
385 : }
386 : };
387 : std::unique_ptr<FILE,decltype(deleter)> grex_log_deleter(grex_log,deleter);
388 :
389 831 : if(debug_grex) {
390 30 : if(noatoms) {
391 0 : error("must have atoms to debug_grex");
392 : }
393 30 : if(multi<2) {
394 0 : error("--debug_grex needs --multi with at least two replicas");
395 : }
396 30 : Tools::convert(fakein,grex_stride);
397 : std::string n;
398 30 : Tools::convert(intercomm.Get_rank(),n);
399 : std::string file;
400 60 : parse("--debug-grex-log",file);
401 30 : if(file.length()>0) {
402 60 : file+="."+n;
403 30 : grex_log=std::fopen(file.c_str(),"w");
404 : grex_log_deleter.reset(grex_log);
405 : }
406 : }
407 :
408 : // Read the plumed input file name
409 : std::string plumedFile;
410 831 : parse("--plumed",plumedFile);
411 : // the timestep
412 : double t;
413 831 : parse("--timestep",t);
414 831 : real timestep=real(t);
415 : // the stride
416 : unsigned stride;
417 831 : parse("--trajectory-stride",stride);
418 : // are we writing forces
419 831 : std::string dumpforces(""), debugforces(""), dumpforcesFmt("%f");;
420 831 : bool dumpfullvirial=false;
421 831 : if(!noatoms) {
422 777 : parse("--dump-forces",dumpforces);
423 1554 : parse("--debug-forces",debugforces);
424 : }
425 1204 : if(dumpforces!="" || debugforces!="" ) {
426 936 : parse("--dump-forces-fmt",dumpforcesFmt);
427 : }
428 831 : if(dumpforces!="") {
429 916 : parseFlag("--dump-full-virial",dumpfullvirial);
430 : }
431 831 : if( debugforces!="" && (debug_dd || debug_pd) ) {
432 0 : error("cannot debug forces and domain/particle decomposition at same time");
433 : }
434 0 : if( debugforces!="" && sizeof(real)!=sizeof(double) ) {
435 0 : error("cannot debug forces in single precision mode");
436 : }
437 :
438 831 : real kt=-1.0;
439 1662 : parse("--kt",kt);
440 : std::string trajectory_fmt;
441 :
442 : bool use_molfile=false;
443 : molfile_plugin_t *api=NULL;
444 :
445 : // Read in an xyz file
446 831 : std::string trajectoryFile(""), pdbfile(""), mcfile("");
447 : bool pbc_cli_given=false;
448 831 : std::vector<double> pbc_cli_box(9,0.0);
449 831 : int command_line_natoms=-1;
450 :
451 831 : if(!noatoms) {
452 : std::string traj_xyz;
453 1554 : parse("--ixyz",traj_xyz);
454 : std::string traj_gro;
455 1554 : parse("--igro",traj_gro);
456 : std::string traj_dlp4;
457 1554 : parse("--idlp4",traj_dlp4);
458 : std::string traj_xtc;
459 : std::string traj_trr;
460 777 : parse("--ixtc",traj_xtc);
461 777 : parse("--itrr",traj_trr);
462 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
463 7770 : for(unsigned i=0; i<plugins.size(); i++) {
464 13986 : std::string molfile_key="--mf_"+std::string(plugins[i]->name);
465 : std::string traj_molfile;
466 6993 : parse(molfile_key,traj_molfile);
467 6993 : if(traj_molfile.length()>0) {
468 267 : std::fprintf(out,"\nDRIVER: Found molfile format trajectory %s with name %s\n",plugins[i]->name,traj_molfile.c_str());
469 : trajectoryFile=traj_molfile;
470 534 : trajectory_fmt=std::string(plugins[i]->name);
471 : use_molfile=true;
472 267 : api = plugins[i];
473 : }
474 : }
475 : #endif
476 : {
477 : // check that only one fmt is specified
478 : int nn=0;
479 777 : if(traj_xyz.length()>0) {
480 : nn++;
481 : }
482 777 : if(traj_gro.length()>0) {
483 112 : nn++;
484 : }
485 777 : if(traj_dlp4.length()>0) {
486 2 : nn++;
487 : }
488 777 : if(traj_xtc.length()>0) {
489 2 : nn++;
490 : }
491 777 : if(traj_trr.length()>0) {
492 3 : nn++;
493 : }
494 777 : if(nn>1) {
495 0 : std::fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
496 : return 1;
497 : }
498 : }
499 777 : if(traj_xyz.length()>0 && trajectoryFile.length()==0) {
500 : trajectoryFile=traj_xyz;
501 : trajectory_fmt="xyz";
502 : }
503 777 : if(traj_gro.length()>0 && trajectoryFile.length()==0) {
504 : trajectoryFile=traj_gro;
505 : trajectory_fmt="gro";
506 : }
507 777 : if(traj_dlp4.length()>0 && trajectoryFile.length()==0) {
508 : trajectoryFile=traj_dlp4;
509 : trajectory_fmt="dlp4";
510 : }
511 777 : if(traj_xtc.length()>0 && trajectoryFile.length()==0) {
512 : trajectoryFile=traj_xtc;
513 : trajectory_fmt="xdr-xtc";
514 : }
515 777 : if(traj_trr.length()>0 && trajectoryFile.length()==0) {
516 : trajectoryFile=traj_trr;
517 : trajectory_fmt="xdr-trr";
518 : }
519 777 : if(trajectoryFile.length()==0&&!parseOnly) {
520 0 : std::fprintf(stderr,"ERROR: missing trajectory data\n");
521 : return 1;
522 : }
523 777 : std::string lengthUnits("");
524 1554 : parse("--length-units",lengthUnits);
525 777 : if(lengthUnits.length()>0) {
526 16 : units.setLength(lengthUnits);
527 : }
528 777 : std::string chargeUnits("");
529 1554 : parse("--charge-units",chargeUnits);
530 777 : if(chargeUnits.length()>0) {
531 1 : units.setCharge(chargeUnits);
532 : }
533 777 : std::string massUnits("");
534 1554 : parse("--mass-units",massUnits);
535 777 : if(massUnits.length()>0) {
536 1 : units.setMass(massUnits);
537 : }
538 :
539 1554 : parse("--pdb",pdbfile);
540 777 : if(pdbfile.length()>0) {
541 75 : bool check=pdb.read(pdbfile,false,1.0);
542 75 : if(!check) {
543 0 : error("error reading pdb file");
544 : }
545 : }
546 :
547 1554 : parse("--mc",mcfile);
548 :
549 : std::string pbc_cli_list;
550 1554 : parse("--box",pbc_cli_list);
551 777 : if(pbc_cli_list.length()>0) {
552 : pbc_cli_given=true;
553 35 : std::vector<std::string> words=Tools::getWords(pbc_cli_list,",");
554 35 : if(words.size()==3) {
555 140 : for(int i=0; i<3; i++) {
556 105 : std::sscanf(words[i].c_str(),"%100lf",&(pbc_cli_box[4*i]));
557 : }
558 0 : } else if(words.size()==9) {
559 0 : for(int i=0; i<9; i++) {
560 0 : std::sscanf(words[i].c_str(),"%100lf",&(pbc_cli_box[i]));
561 : }
562 : } else {
563 0 : std::string msg="ERROR: cannot parse command-line box "+pbc_cli_list;
564 0 : std::fprintf(stderr,"%s\n",msg.c_str());
565 : return 1;
566 : }
567 :
568 35 : }
569 :
570 1554 : parse("--natoms",command_line_natoms);
571 :
572 : }
573 :
574 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
575 267 : auto mf_deleter=[api](void* h_in) {
576 267 : if(h_in) {
577 267 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
578 267 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) {
579 420 : lck=Tools::molfile_lock();
580 : }
581 267 : api->close_file_read(h_in);
582 267 : }
583 : };
584 : void *h_in=NULL;
585 : std::unique_ptr<void,decltype(mf_deleter)> h_in_deleter(h_in,mf_deleter);
586 :
587 : molfile_timestep_t ts_in; // this is the structure that has the timestep
588 : // a std::vector<float> with the same scope as ts_in
589 : // it is necessary in order to store the pointer to ts_in.coords
590 : std::vector<float> ts_in_coords;
591 831 : ts_in.coords=ts_in_coords.data();
592 831 : ts_in.velocities=NULL;
593 831 : ts_in.A=-1; // we use this to check whether cell is provided or not
594 : #endif
595 :
596 :
597 :
598 831 : if(debug_dd && debug_pd) {
599 0 : error("cannot use debug-dd and debug-pd at the same time");
600 : }
601 831 : if(debug_pd || debug_dd) {
602 72 : if( !Communicator::initialized() ) {
603 0 : error("needs mpi for debug-pd");
604 : }
605 : }
606 :
607 831 : PlumedMain p;
608 831 : if( parseOnly ) {
609 0 : p.activateParseOnlyMode();
610 : }
611 1662 : p.cmd("setRealPrecision",(int)sizeof(real));
612 : int checknatoms=-1;
613 831 : long long int step=0;
614 831 : parse("--initial-step",step);
615 :
616 831 : if(restart) {
617 4 : p.cmd("setRestart",1);
618 : }
619 :
620 831 : if(Communicator::initialized()) {
621 314 : if(multi) {
622 174 : if(intracomm.Get_rank()==0) {
623 378 : p.cmd("GREX setMPIIntercomm",&intercomm.Get_comm());
624 : }
625 522 : p.cmd("GREX setMPIIntracomm",&intracomm.Get_comm());
626 348 : p.cmd("GREX init");
627 : }
628 942 : p.cmd("setMPIComm",&intracomm.Get_comm());
629 : }
630 1662 : p.cmd("setLog",out);
631 1662 : p.cmd("setMDLengthUnits",units.getLength());
632 1662 : p.cmd("setMDChargeUnits",units.getCharge());
633 1662 : p.cmd("setMDMassUnits",units.getMass());
634 1662 : p.cmd("setMDEngine","driver");
635 1662 : p.cmd("setTimestep",timestep);
636 831 : if( !parseOnly || full_outputfile.length()==0 ) {
637 1662 : p.cmd("setPlumedDat",plumedFile.c_str());
638 : }
639 :
640 : int natoms;
641 831 : int lvl=0;
642 831 : int pb=1;
643 :
644 831 : if(parseOnly) {
645 0 : if(command_line_natoms<0) {
646 0 : error("--parseOnly requires setting the number of atoms with --natoms");
647 : }
648 0 : natoms=command_line_natoms;
649 : }
650 :
651 :
652 : FILE* fp=NULL;
653 : FILE* fp_forces=NULL;
654 831 : OFile fp_dforces;
655 :
656 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
657 : std::unique_ptr<FILE,decltype(deleter)> fp_forces_deleter(fp_forces,deleter);
658 :
659 5 : auto xdr_deleter=[](xdrfile::XDRFILE* xd) {
660 : if(xd) {
661 5 : xdrfile::xdrfile_close(xd);
662 : }
663 : };
664 :
665 : xdrfile::XDRFILE* xd=NULL;
666 :
667 : std::unique_ptr<xdrfile::XDRFILE,decltype(xdr_deleter)> xd_deleter(xd,xdr_deleter);
668 :
669 831 : if(!noatoms&&!parseOnly) {
670 777 : if (trajectoryFile=="-") {
671 : fp=in;
672 : } else {
673 777 : if(multi) {
674 : std::string n;
675 174 : Tools::convert(intercomm.Get_rank(),n);
676 348 : std::string testfile=FileBase::appendSuffix(trajectoryFile,"."+n);
677 174 : FILE* tmp_fp=std::fopen(testfile.c_str(),"r");
678 : // no exceptions here
679 174 : if(tmp_fp) {
680 135 : std::fclose(tmp_fp);
681 : trajectoryFile=testfile;
682 : }
683 : }
684 777 : if(use_molfile==true) {
685 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
686 267 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
687 267 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) {
688 420 : lck=Tools::molfile_lock();
689 : }
690 267 : h_in = api->open_file_read(trajectoryFile.c_str(), trajectory_fmt.c_str(), &natoms);
691 : h_in_deleter.reset(h_in);
692 267 : if(natoms==MOLFILE_NUMATOMS_UNKNOWN) {
693 2 : if(command_line_natoms>=0) {
694 2 : natoms=command_line_natoms;
695 : } else {
696 0 : error("this file format does not provide number of atoms; use --natoms on the command line");
697 : }
698 : }
699 267 : ts_in_coords.resize(3*natoms);
700 267 : ts_in.coords = ts_in_coords.data();
701 : #endif
702 1285 : } else if(trajectory_fmt=="xdr-xtc" || trajectory_fmt=="xdr-trr") {
703 5 : xd=xdrfile::xdrfile_open(trajectoryFile.c_str(),"r");
704 : xd_deleter.reset(xd);
705 5 : if(!xd) {
706 0 : std::string msg="ERROR: Error opening trajectory file "+trajectoryFile;
707 0 : std::fprintf(stderr,"%s\n",msg.c_str());
708 : return 1;
709 : }
710 5 : if(trajectory_fmt=="xdr-xtc") {
711 2 : xdrfile::read_xtc_natoms(&trajectoryFile[0],&natoms);
712 : }
713 5 : if(trajectory_fmt=="xdr-trr") {
714 3 : xdrfile::read_trr_natoms(&trajectoryFile[0],&natoms);
715 : }
716 : } else {
717 505 : fp=std::fopen(trajectoryFile.c_str(),"r");
718 : fp_deleter.reset(fp);
719 505 : if(!fp) {
720 0 : std::string msg="ERROR: Error opening trajectory file "+trajectoryFile;
721 0 : std::fprintf(stderr,"%s\n",msg.c_str());
722 : return 1;
723 : }
724 : }
725 : }
726 777 : if(dumpforces.length()>0) {
727 458 : if(Communicator::initialized() && pc.Get_size()>1) {
728 : std::string n;
729 232 : Tools::convert(pc.Get_rank(),n);
730 464 : dumpforces+="."+n;
731 : }
732 458 : fp_forces=std::fopen(dumpforces.c_str(),"w");
733 : fp_forces_deleter.reset(fp_forces);
734 : }
735 777 : if(debugforces.length()>0) {
736 11 : if(Communicator::initialized() && pc.Get_size()>1) {
737 : std::string n;
738 6 : Tools::convert(pc.Get_rank(),n);
739 12 : debugforces+="."+n;
740 : }
741 11 : fp_dforces.open(debugforces);
742 : }
743 : }
744 :
745 : std::string line;
746 : std::vector<real> coordinates;
747 : std::vector<real> forces;
748 : std::vector<real> masses;
749 : std::vector<real> charges;
750 : std::vector<real> cell;
751 : std::vector<real> virial;
752 : std::vector<real> numder;
753 :
754 : // variables to test particle decomposition
755 : int pd_nlocal=0;
756 : int pd_start=0;
757 : // variables to test random decomposition (=domain decomposition)
758 : std::vector<int> dd_gatindex;
759 : std::vector<int> dd_g2l;
760 : std::vector<real> dd_masses;
761 : std::vector<real> dd_charges;
762 : std::vector<real> dd_forces;
763 : std::vector<real> dd_coordinates;
764 : int dd_nlocal=0;
765 : // random stream to choose decompositions
766 831 : Random rnd;
767 :
768 831 : if(trajectory_fmt=="dlp4") {
769 2 : if(!Tools::getline(fp,line)) {
770 0 : error("error reading title");
771 : }
772 2 : if(!Tools::getline(fp,line)) {
773 0 : error("error reading atoms");
774 : }
775 2 : std::sscanf(line.c_str(),"%d %d %d",&lvl,&pb,&natoms);
776 :
777 : }
778 : bool lstep=true;
779 246402 : while(true) {
780 247233 : if(!noatoms&&!parseOnly) {
781 46893 : if(use_molfile==true) {
782 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
783 20887 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
784 20887 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) {
785 41454 : lck=Tools::molfile_lock();
786 : }
787 : int rc;
788 20887 : rc = api->read_next_timestep(h_in, natoms, &ts_in);
789 20887 : if(rc==MOLFILE_EOF) {
790 : break;
791 : }
792 : #endif
793 49524 : } else if(trajectory_fmt=="xyz" || trajectory_fmt=="gro" || trajectory_fmt=="dlp4") {
794 25937 : if(!Tools::getline(fp,line)) {
795 : break;
796 : }
797 : }
798 : }
799 : bool first_step=false;
800 246465 : if(!noatoms&&!parseOnly) {
801 74059 : if(use_molfile==false && (trajectory_fmt=="xyz" || trajectory_fmt=="gro")) {
802 25416 : if(trajectory_fmt=="gro")
803 2340 : if(!Tools::getline(fp,line)) {
804 0 : error("premature end of trajectory file");
805 : }
806 25416 : std::sscanf(line.c_str(),"%100d",&natoms);
807 : }
808 71630 : if(use_molfile==false && trajectory_fmt=="dlp4") {
809 : char xa[9];
810 : int xb,xc,xd;
811 : double t;
812 20 : std::sscanf(line.c_str(),"%8s %lld %d %d %d %lf",xa,&step,&xb,&xc,&xd,&t);
813 20 : if (lstep) {
814 4 : p.cmd("setTimestep",real(t));
815 : lstep = false;
816 : }
817 : }
818 : }
819 246465 : if(checknatoms<0 && !noatoms) {
820 777 : pd_nlocal=natoms;
821 : pd_start=0;
822 : first_step=true;
823 777 : masses.assign(natoms,std::numeric_limits<real>::quiet_NaN());
824 777 : charges.assign(natoms,std::numeric_limits<real>::quiet_NaN());
825 : //case pdb: structure
826 777 : if(pdbfile.length()>0) {
827 28092 : for(unsigned i=0; i<pdb.size(); ++i) {
828 28017 : AtomNumber an=pdb.getAtomNumbers()[i];
829 : unsigned index=an.index();
830 28017 : if( index>=unsigned(natoms) ) {
831 0 : error("atom index in pdb exceeds the number of atoms in trajectory");
832 : }
833 28017 : masses[index]=pdb.getOccupancy()[i];
834 28017 : charges[index]=pdb.getBeta()[i];
835 : }
836 : }
837 777 : if(mcfile.length()>0) {
838 5 : IFile ifile;
839 5 : ifile.open(mcfile);
840 : int index;
841 : double mass;
842 : double charge;
843 1138 : while(ifile.scanField("index",index).scanField("mass",mass).scanField("charge",charge).scanField()) {
844 564 : masses[index]=mass;
845 564 : charges[index]=charge;
846 : }
847 5 : }
848 245688 : } else if( checknatoms<0 && noatoms ) {
849 54 : natoms=0;
850 : }
851 246465 : if( checknatoms<0 ) {
852 831 : if(kt>=0) {
853 6 : p.cmd("setKbT",kt);
854 : }
855 831 : checknatoms=natoms;
856 1662 : p.cmd("setNatoms",natoms);
857 1662 : p.cmd("init");
858 : // Check if we have been asked to output the long version of the input and if there are shortcuts
859 831 : if( parseOnly && full_outputfile.length()>0 ) {
860 :
861 : // Read in the plumed input file and store what is in there
862 : std::map<std::string,std::vector<std::string> > data;
863 0 : IFile ifile;
864 0 : ifile.open(plumedFile);
865 : std::vector<std::string> words;
866 0 : while( Tools::getParsedLine(ifile,words) && !p.getEndPlumed() ) {
867 0 : p.readInputWords(words);
868 : Action* aa=p.getActionSet()[p.getActionSet().size()-1].get();
869 0 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(aa);
870 0 : if( av && aa->getDefaultString().length()>0 ) {
871 : std::vector<std::string> def;
872 0 : def.push_back( "defaults " + aa->getDefaultString() );
873 0 : data[ aa->getLabel() ] = def;
874 0 : }
875 0 : ActionShortcut* as=dynamic_cast<ActionShortcut*>( aa );
876 0 : if( as ) {
877 0 : if( aa->getDefaultString().length()>0 ) {
878 : std::vector<std::string> def;
879 0 : def.push_back( "defaults " + aa->getDefaultString() );
880 0 : data[ as->getShortcutLabel() ] = def;
881 0 : }
882 0 : if( data.find( as->getShortcutLabel() )!=data.end() ) {
883 0 : std::vector<std::string> shortcut_commands = as->getSavedInputLines();
884 0 : for(unsigned i=0; i<shortcut_commands.size(); ++i) {
885 0 : data[ as->getShortcutLabel() ].push_back( shortcut_commands[i] );
886 : }
887 0 : } else {
888 0 : data[ as->getShortcutLabel() ] = as->getSavedInputLines();
889 : }
890 : }
891 : }
892 0 : ifile.close();
893 : // Only output the full version of the input file if there are shortcuts
894 0 : if( data.size()>0 ) {
895 0 : OFile long_file;
896 0 : long_file.open( full_outputfile );
897 0 : long_file.printf("{\n");
898 : bool firstpass=true;
899 0 : for(auto& x : data ) {
900 0 : if( !firstpass ) {
901 0 : long_file.printf(" },\n");
902 : }
903 0 : long_file.printf(" \"%s\" : {\n", x.first.c_str() );
904 0 : plumed_assert( x.second.size()>0 );
905 : unsigned sstart=0;
906 0 : if( x.second[0].find("defaults")!=std::string::npos ) {
907 : sstart=1;
908 0 : long_file.printf(" \"defaults\" : \"%s\"", x.second[0].substr( 9 ).c_str() );
909 0 : if( x.second.size()>1 ) {
910 0 : long_file.printf(",\n");
911 : } else {
912 0 : long_file.printf("\n");
913 : }
914 : }
915 0 : if( x.second.size()>sstart ) {
916 0 : long_file.printf(" \"expansion\" : \"%s", x.second[sstart].c_str() );
917 0 : for(unsigned j=sstart+1; j<x.second.size(); ++j) {
918 0 : long_file.printf("\\n%s", x.second[j].c_str() );
919 : }
920 0 : long_file.printf("\"\n");
921 : }
922 : firstpass=false;
923 : }
924 0 : long_file.printf(" }\n}\n");
925 0 : long_file.close();
926 0 : }
927 0 : }
928 831 : if(parseOnly) {
929 : break;
930 : }
931 : }
932 246465 : if(checknatoms!=natoms) {
933 : std::string stepstr;
934 0 : Tools::convert(step,stepstr);
935 0 : error("number of atoms in frame " + stepstr + " does not match number of atoms in first frame");
936 : }
937 :
938 246465 : coordinates.assign(3*natoms,real(0.0));
939 246465 : forces.assign(3*natoms,real(0.0));
940 246465 : cell.assign(9,real(0.0));
941 246465 : virial.assign(9,real(0.0));
942 :
943 246465 : if( first_step || rnd.U01()>0.5) {
944 124895 : if(debug_pd) {
945 152 : int npe=intracomm.Get_size();
946 152 : std::vector<int> loc(npe,0);
947 152 : std::vector<int> start(npe,0);
948 608 : for(int i=0; i<npe-1; i++) {
949 456 : int cc=(natoms*2*rnd.U01())/npe;
950 456 : if(start[i]+cc>natoms) {
951 16 : cc=natoms-start[i];
952 : }
953 456 : loc[i]=cc;
954 456 : start[i+1]=start[i]+loc[i];
955 : }
956 152 : loc[npe-1]=natoms-start[npe-1];
957 152 : intracomm.Bcast(loc,0);
958 152 : intracomm.Bcast(start,0);
959 152 : pd_nlocal=loc[intracomm.Get_rank()];
960 152 : pd_start=start[intracomm.Get_rank()];
961 152 : if(intracomm.Get_rank()==0) {
962 : std::fprintf(out,"\nDRIVER: Reassigning particle decomposition\n");
963 : std::fprintf(out,"DRIVER: ");
964 190 : for(int i=0; i<npe; i++) {
965 152 : std::fprintf(out,"%d ",loc[i]);
966 : }
967 : printf("\n");
968 : std::fprintf(out,"DRIVER: ");
969 190 : for(int i=0; i<npe; i++) {
970 152 : std::fprintf(out,"%d ",start[i]);
971 : }
972 : printf("\n");
973 : }
974 304 : p.cmd("setAtomsNlocal",pd_nlocal);
975 304 : p.cmd("setAtomsContiguous",pd_start);
976 124743 : } else if(debug_dd) {
977 956 : int npe=intracomm.Get_size();
978 956 : int rank=intracomm.Get_rank();
979 956 : dd_charges.assign(natoms,0.0);
980 956 : dd_masses.assign(natoms,0.0);
981 956 : dd_gatindex.assign(natoms,-1);
982 956 : dd_g2l.assign(natoms,-1);
983 956 : dd_coordinates.assign(3*natoms,0.0);
984 956 : dd_forces.assign(3*natoms,0.0);
985 : dd_nlocal=0;
986 53786 : for(int i=0; i<natoms; ++i) {
987 52830 : double r=rnd.U01()*npe;
988 : int n;
989 112376 : for(n=0; n<npe; n++)
990 112376 : if(n+1>r) {
991 : break;
992 : }
993 52830 : plumed_assert(n<npe);
994 52830 : if(n==rank) {
995 19827 : dd_gatindex[dd_nlocal]=i;
996 19827 : dd_g2l[i]=dd_nlocal;
997 19827 : dd_charges[dd_nlocal]=charges[i];
998 19827 : dd_masses[dd_nlocal]=masses[i];
999 19827 : dd_nlocal++;
1000 : }
1001 : }
1002 956 : if(intracomm.Get_rank()==0) {
1003 : std::fprintf(out,"\nDRIVER: Reassigning domain decomposition\n");
1004 : }
1005 1912 : p.cmd("setAtomsNlocal",dd_nlocal);
1006 1912 : p.cmd("setAtomsGatindex",&dd_gatindex[0],dd_nlocal);
1007 : }
1008 : }
1009 :
1010 246465 : int plumedStopCondition=0;
1011 246465 : if(!noatoms) {
1012 46125 : if(use_molfile) {
1013 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
1014 20620 : if(pbc_cli_given==false) {
1015 20580 : if(ts_in.A>0.0) { // this is negative if molfile does not provide box
1016 : // info on the cell: convert using pbcset.tcl from pbctools in vmd distribution
1017 20575 : real cosBC=cos(real(ts_in.alpha)*pi/180.);
1018 : //double sinBC=std::sin(ts_in.alpha*pi/180.);
1019 20575 : real cosAC=std::cos(real(ts_in.beta)*pi/180.);
1020 20575 : real cosAB=std::cos(real(ts_in.gamma)*pi/180.);
1021 20575 : real sinAB=std::sin(real(ts_in.gamma)*pi/180.);
1022 20575 : real Ax=real(ts_in.A);
1023 20575 : real Bx=real(ts_in.B)*cosAB;
1024 20575 : real By=real(ts_in.B)*sinAB;
1025 20575 : real Cx=real(ts_in.C)*cosAC;
1026 20575 : real Cy=(real(ts_in.C)*real(ts_in.B)*cosBC-Cx*Bx)/By;
1027 20575 : real Cz=std::sqrt(real(ts_in.C)*real(ts_in.C)-Cx*Cx-Cy*Cy);
1028 20575 : cell[0]=Ax/10.;
1029 20575 : cell[1]=0.;
1030 20575 : cell[2]=0.;
1031 20575 : cell[3]=Bx/10.;
1032 20575 : cell[4]=By/10.;
1033 20575 : cell[5]=0.;
1034 20575 : cell[6]=Cx/10.;
1035 20575 : cell[7]=Cy/10.;
1036 20575 : cell[8]=Cz/10.;
1037 : } else {
1038 5 : cell[0]=0.0;
1039 5 : cell[1]=0.0;
1040 5 : cell[2]=0.0;
1041 5 : cell[3]=0.0;
1042 5 : cell[4]=0.0;
1043 5 : cell[5]=0.0;
1044 5 : cell[6]=0.0;
1045 5 : cell[7]=0.0;
1046 5 : cell[8]=0.0;
1047 : }
1048 : } else {
1049 400 : for(unsigned i=0; i<9; i++) {
1050 360 : cell[i]=pbc_cli_box[i];
1051 : }
1052 : }
1053 : // info on coords
1054 : // the order is xyzxyz...
1055 63325219 : for(int i=0; i<3*natoms; i++) {
1056 63304599 : coordinates[i]=real(ts_in.coords[i])/real(10.); //convert to nm
1057 : //cerr<<"COOR "<<coordinates[i]<<endl;
1058 : }
1059 : #endif
1060 50980 : } else if(trajectory_fmt=="xdr-xtc" || trajectory_fmt=="xdr-trr") {
1061 : int localstep;
1062 : float time;
1063 : xdrfile::matrix box;
1064 : // here we cannot use a std::vector<rvec> since it does not compile.
1065 : // we thus use a std::unique_ptr<rvec[]>
1066 69 : auto pos=Tools::make_unique<xdrfile::rvec[]>(natoms);
1067 : float prec,lambda;
1068 : int ret=xdrfile::exdrOK;
1069 69 : if(trajectory_fmt=="xdr-xtc") {
1070 30 : ret=xdrfile::read_xtc(xd,natoms,&localstep,&time,box,pos.get(),&prec);
1071 : }
1072 69 : if(trajectory_fmt=="xdr-trr") {
1073 39 : ret=xdrfile::read_trr(xd,natoms,&localstep,&time,&lambda,box,pos.get(),NULL,NULL);
1074 : }
1075 69 : if(stride==0) {
1076 30 : step=localstep;
1077 : }
1078 69 : if(ret==xdrfile::exdrENDOFFILE) {
1079 : break;
1080 : }
1081 67 : if(ret!=xdrfile::exdrOK) {
1082 : break;
1083 : }
1084 256 : for(unsigned i=0; i<3; i++)
1085 768 : for(unsigned j=0; j<3; j++) {
1086 576 : cell[3*i+j]=box[i][j];
1087 : }
1088 2792 : for(int i=0; i<natoms; i++)
1089 10912 : for(unsigned j=0; j<3; j++) {
1090 8184 : coordinates[3*i+j]=real(pos[i][j]);
1091 : }
1092 : } else {
1093 25436 : if(trajectory_fmt=="xyz") {
1094 23076 : if(!Tools::getline(fp,line)) {
1095 0 : error("premature end of trajectory file");
1096 : }
1097 :
1098 23076 : std::vector<double> celld(9,0.0);
1099 23076 : if(pbc_cli_given==false) {
1100 : std::vector<std::string> words;
1101 45740 : words=Tools::getWords(line);
1102 22870 : if(words.size()==3) {
1103 22297 : std::sscanf(line.c_str(),"%100lf %100lf %100lf",&celld[0],&celld[4],&celld[8]);
1104 573 : } else if(words.size()==9) {
1105 573 : std::sscanf(line.c_str(),"%100lf %100lf %100lf %100lf %100lf %100lf %100lf %100lf %100lf",
1106 : &celld[0], &celld[1], &celld[2],
1107 : &celld[3], &celld[4], &celld[5],
1108 : &celld[6], &celld[7], &celld[8]);
1109 : } else {
1110 0 : error("needed box in second line of xyz file");
1111 : }
1112 22870 : } else { // from command line
1113 206 : celld=pbc_cli_box;
1114 : }
1115 230760 : for(unsigned i=0; i<9; i++) {
1116 207684 : cell[i]=real(celld[i]);
1117 : }
1118 : }
1119 25436 : if(trajectory_fmt=="dlp4") {
1120 20 : std::vector<double> celld(9,0.0);
1121 20 : if(pbc_cli_given==false) {
1122 20 : if(!Tools::getline(fp,line)) {
1123 0 : error("error reading vector a of cell");
1124 : }
1125 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[0],&celld[1],&celld[2]);
1126 20 : if(!Tools::getline(fp,line)) {
1127 0 : error("error reading vector b of cell");
1128 : }
1129 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[3],&celld[4],&celld[5]);
1130 20 : if(!Tools::getline(fp,line)) {
1131 0 : error("error reading vector c of cell");
1132 : }
1133 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[6],&celld[7],&celld[8]);
1134 : } else {
1135 0 : celld=pbc_cli_box;
1136 : }
1137 200 : for(auto i=0; i<9; i++) {
1138 180 : cell[i]=real(celld[i])*0.1;
1139 : }
1140 : }
1141 : int ddist=0;
1142 : // Read coordinates
1143 2483105 : for(int i=0; i<natoms; i++) {
1144 2457669 : bool ok=Tools::getline(fp,line);
1145 2457669 : if(!ok) {
1146 0 : error("premature end of trajectory file");
1147 : }
1148 : double cc[3];
1149 2457669 : if(trajectory_fmt=="xyz") {
1150 : char dummy[1000];
1151 1744555 : int ret=std::sscanf(line.c_str(),"%999s %100lf %100lf %100lf",dummy,&cc[0],&cc[1],&cc[2]);
1152 1744555 : if(ret!=4) {
1153 0 : error("cannot read line"+line);
1154 : }
1155 713114 : } else if(trajectory_fmt=="gro") {
1156 : // do the gromacs way
1157 712474 : if(!i) {
1158 : //
1159 : // calculate the distance between dots (as in gromacs gmxlib/confio.c, routine get_w_conf )
1160 : //
1161 : const char *p1, *p2, *p3;
1162 : p1 = std::strchr(line.c_str(), '.');
1163 2340 : if (p1 == NULL) {
1164 0 : error("seems there are no coordinates in the gro file");
1165 : }
1166 2340 : p2 = std::strchr(&p1[1], '.');
1167 2340 : if (p2 == NULL) {
1168 0 : error("seems there is only one coordinates in the gro file");
1169 : }
1170 2340 : ddist = p2 - p1;
1171 2340 : p3 = std::strchr(&p2[1], '.');
1172 2340 : if (p3 == NULL) {
1173 0 : error("seems there are only two coordinates in the gro file");
1174 : }
1175 2340 : if (p3 - p2 != ddist) {
1176 0 : error("not uniform spacing in fields in the gro file");
1177 : }
1178 : }
1179 712474 : Tools::convert(line.substr(20,ddist),cc[0]);
1180 712474 : Tools::convert(line.substr(20+ddist,ddist),cc[1]);
1181 1424948 : Tools::convert(line.substr(20+ddist+ddist,ddist),cc[2]);
1182 640 : } else if(trajectory_fmt=="dlp4") {
1183 : char dummy[9];
1184 : int idummy;
1185 : double m,c;
1186 640 : std::sscanf(line.c_str(),"%8s %d %lf %lf",dummy,&idummy,&m,&c);
1187 640 : masses[i]=real(m);
1188 640 : charges[i]=real(c);
1189 640 : if(!Tools::getline(fp,line)) {
1190 0 : error("error reading coordinates");
1191 : }
1192 640 : std::sscanf(line.c_str(),"%lf %lf %lf",&cc[0],&cc[1],&cc[2]);
1193 640 : cc[0]*=0.1;
1194 640 : cc[1]*=0.1;
1195 640 : cc[2]*=0.1;
1196 640 : if(lvl>0) {
1197 640 : if(!Tools::getline(fp,line)) {
1198 0 : error("error skipping velocities");
1199 : }
1200 : }
1201 640 : if(lvl>1) {
1202 640 : if(!Tools::getline(fp,line)) {
1203 0 : error("error skipping forces");
1204 : }
1205 : }
1206 : } else {
1207 0 : plumed_error();
1208 : }
1209 2457669 : if(!debug_pd || ( i>=pd_start && i<pd_start+pd_nlocal) ) {
1210 2438229 : coordinates[3*i]=real(cc[0]);
1211 2438229 : coordinates[3*i+1]=real(cc[1]);
1212 2438229 : coordinates[3*i+2]=real(cc[2]);
1213 : }
1214 : }
1215 25436 : if(trajectory_fmt=="gro") {
1216 2340 : if(!Tools::getline(fp,line)) {
1217 0 : error("premature end of trajectory file");
1218 : }
1219 2340 : std::vector<std::string> words=Tools::getWords(line);
1220 2340 : if(words.size()<3) {
1221 0 : error("cannot understand box format");
1222 : }
1223 2340 : Tools::convert(words[0],cell[0]);
1224 2340 : Tools::convert(words[1],cell[4]);
1225 2340 : Tools::convert(words[2],cell[8]);
1226 2340 : if(words.size()>3) {
1227 510 : Tools::convert(words[3],cell[1]);
1228 : }
1229 2340 : if(words.size()>4) {
1230 510 : Tools::convert(words[4],cell[2]);
1231 : }
1232 2340 : if(words.size()>5) {
1233 510 : Tools::convert(words[5],cell[3]);
1234 : }
1235 2340 : if(words.size()>6) {
1236 510 : Tools::convert(words[6],cell[5]);
1237 : }
1238 2340 : if(words.size()>7) {
1239 510 : Tools::convert(words[7],cell[6]);
1240 : }
1241 2340 : if(words.size()>8) {
1242 510 : Tools::convert(words[8],cell[7]);
1243 : }
1244 2340 : }
1245 :
1246 : }
1247 :
1248 92240 : p.cmd("setStepLongLong",step);
1249 92240 : p.cmd("setStopFlag",&plumedStopCondition);
1250 :
1251 46120 : if(debug_dd) {
1252 38944 : for(int i=0; i<dd_nlocal; ++i) {
1253 37156 : int kk=dd_gatindex[i];
1254 37156 : dd_coordinates[3*i+0]=coordinates[3*kk+0];
1255 37156 : dd_coordinates[3*i+1]=coordinates[3*kk+1];
1256 37156 : dd_coordinates[3*i+2]=coordinates[3*kk+2];
1257 : }
1258 3576 : p.cmd("setForces",&dd_forces[0],3*dd_nlocal);
1259 3576 : p.cmd("setPositions",&dd_coordinates[0],3*dd_nlocal);
1260 3576 : p.cmd("setMasses",&dd_masses[0],dd_nlocal);
1261 3576 : p.cmd("setCharges",&dd_charges[0],dd_nlocal);
1262 : } else {
1263 : // this is required to avoid troubles when the last domain
1264 : // contains zero atoms
1265 : // Basically, for empty domains we pass null pointers
1266 : #define fix_pd(xx) (pd_nlocal!=0?&xx:NULL)
1267 132996 : p.cmd("setForces",fix_pd(forces[3*pd_start]),3*pd_nlocal);
1268 132996 : p.cmd("setPositions",fix_pd(coordinates[3*pd_start]),3*pd_nlocal);
1269 132996 : p.cmd("setMasses",fix_pd(masses[pd_start]),pd_nlocal);
1270 132996 : p.cmd("setCharges",fix_pd(charges[pd_start]),pd_nlocal);
1271 : }
1272 92240 : p.cmd("setBox",cell.data(),9);
1273 92240 : p.cmd("setVirial",virial.data(),9);
1274 : } else {
1275 400680 : p.cmd("setStepLongLong",step);
1276 400680 : p.cmd("setStopFlag",&plumedStopCondition);
1277 : }
1278 492920 : p.cmd("calc");
1279 246460 : if(debugforces.length()>0) {
1280 25 : virial.assign(9,real(0.0));
1281 25 : forces.assign(3*natoms,real(0.0));
1282 50 : p.cmd("prepareCalc");
1283 50 : p.cmd("performCalcNoUpdate");
1284 : }
1285 :
1286 : // this is necessary as only processor zero is adding to the virial:
1287 246460 : intracomm.Bcast(virial,0);
1288 246460 : if(debug_pd) {
1289 240 : intracomm.Sum(forces);
1290 : }
1291 246460 : if(debug_dd) {
1292 38944 : for(int i=0; i<dd_nlocal; i++) {
1293 37156 : forces[3*dd_gatindex[i]+0]=dd_forces[3*i+0];
1294 37156 : forces[3*dd_gatindex[i]+1]=dd_forces[3*i+1];
1295 37156 : forces[3*dd_gatindex[i]+2]=dd_forces[3*i+2];
1296 : }
1297 1788 : dd_forces.assign(3*natoms,0.0);
1298 1788 : intracomm.Sum(forces);
1299 : }
1300 246460 : if(debug_grex &&step%grex_stride==0) {
1301 228 : p.cmd("GREX savePositions");
1302 114 : if(intracomm.Get_rank()>0) {
1303 114 : p.cmd("GREX prepare");
1304 : } else {
1305 57 : int r=intercomm.Get_rank();
1306 57 : int n=intercomm.Get_size();
1307 57 : int partner=r+(2*((r+step/grex_stride)%2))-1;
1308 : if(partner<0) {
1309 : partner=0;
1310 : }
1311 57 : if(partner>=n) {
1312 8 : partner=n-1;
1313 : }
1314 114 : p.cmd("GREX setPartner",partner);
1315 114 : p.cmd("GREX calculate");
1316 114 : p.cmd("GREX shareAllDeltaBias");
1317 228 : for(int i=0; i<n; i++) {
1318 : std::string s;
1319 171 : Tools::convert(i,s);
1320 171 : real a=std::numeric_limits<real>::quiet_NaN();
1321 342 : s="GREX getDeltaBias "+s;
1322 171 : p.cmd(s,&a);
1323 171 : if(grex_log) {
1324 171 : std::fprintf(grex_log," %f",a);
1325 : }
1326 : }
1327 57 : if(grex_log) {
1328 : std::fprintf(grex_log,"\n");
1329 : }
1330 : }
1331 : }
1332 :
1333 :
1334 246460 : if(fp_forces) {
1335 24547 : std::fprintf(fp_forces,"%d\n",natoms);
1336 49094 : std::string fmtv=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
1337 49094 : std::string fmt=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
1338 24547 : if(dumpfullvirial) {
1339 350 : 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]);
1340 : } else {
1341 24197 : std::fprintf(fp_forces,fmt.c_str(),virial[0],virial[4],virial[8]);
1342 : }
1343 24547 : fmt="X "+fmt;
1344 2480254 : for(int i=0; i<natoms; i++) {
1345 2455707 : std::fprintf(fp_forces,fmt.c_str(),forces[3*i],forces[3*i+1],forces[3*i+2]);
1346 : }
1347 : }
1348 246460 : if(debugforces.length()>0) {
1349 : // Now call the routine to work out the derivatives numerically
1350 25 : numder.assign(3*natoms+9,real(0.0));
1351 25 : real base=0;
1352 50 : p.cmd("getBias",&base);
1353 25 : if( std::abs(base)<epsilon ) {
1354 : printf("WARNING: bias for configuration appears to be zero so debugging forces is trivial");
1355 : }
1356 25 : evaluateNumericalDerivatives( step, p, coordinates, masses, charges, cell, base, numder );
1357 :
1358 : // And output everything to a file
1359 25 : fp_dforces.fmtField(" " + dumpforcesFmt);
1360 1795 : for(int i=0; i<3*natoms; ++i) {
1361 1770 : fp_dforces.printField("parameter",(int)i);
1362 3540 : fp_dforces.printField("analytical",forces[i]);
1363 1770 : fp_dforces.printField("numerical",-numder[i]);
1364 1770 : fp_dforces.printField();
1365 : }
1366 : // And print the virial
1367 250 : for(int i=0; i<9; ++i) {
1368 225 : fp_dforces.printField("parameter",3*natoms+i );
1369 225 : fp_dforces.printField("analytical",virial[i] );
1370 225 : fp_dforces.printField("numerical",-numder[3*natoms+i]);
1371 225 : fp_dforces.printField();
1372 : }
1373 : }
1374 :
1375 246460 : if(plumedStopCondition) {
1376 : break;
1377 : }
1378 :
1379 246402 : step+=stride;
1380 : }
1381 831 : if(!parseOnly) {
1382 1662 : p.cmd("runFinalJobs");
1383 : }
1384 :
1385 : return 0;
1386 3324 : }
1387 :
1388 : template<typename real>
1389 25 : void Driver<real>::evaluateNumericalDerivatives( const long long int& step, PlumedMain& p, const std::vector<real>& coordinates,
1390 : const std::vector<real>& masses, const std::vector<real>& charges,
1391 : std::vector<real>& cell, const double& base, std::vector<real>& numder ) {
1392 :
1393 25 : int natoms = coordinates.size() / 3;
1394 : real delta = std::sqrt(epsilon);
1395 25 : std::vector<Vector> pos(natoms);
1396 25 : real bias=0;
1397 25 : std::vector<real> fake_forces( 3*natoms ), fake_virial(9);
1398 615 : for(int i=0; i<natoms; ++i) {
1399 2360 : for(unsigned j=0; j<3; ++j) {
1400 1770 : pos[i][j]=coordinates[3*i+j];
1401 : }
1402 : }
1403 :
1404 615 : for(int i=0; i<natoms; ++i) {
1405 2360 : for(unsigned j=0; j<3; ++j) {
1406 1770 : pos[i][j]=pos[i][j]+delta;
1407 3540 : p.cmd("setStepLongLong",step);
1408 3540 : p.cmd("setPositions",&pos[0][0],3*natoms);
1409 3540 : p.cmd("setForces",&fake_forces[0],3*natoms);
1410 3540 : p.cmd("setMasses",&masses[0],natoms);
1411 3540 : p.cmd("setCharges",&charges[0],natoms);
1412 3540 : p.cmd("setBox",&cell[0],9);
1413 3540 : p.cmd("setVirial",&fake_virial[0],9);
1414 3540 : p.cmd("prepareCalc");
1415 3540 : p.cmd("performCalcNoUpdate");
1416 3540 : p.cmd("getBias",&bias);
1417 1770 : pos[i][j]=coordinates[3*i+j];
1418 1770 : numder[3*i+j] = (bias - base) / delta;
1419 : }
1420 : }
1421 :
1422 : // Create the cell
1423 25 : Tensor box( cell[0], cell[1], cell[2], cell[3], cell[4], cell[5], cell[6], cell[7], cell[8] );
1424 : // And the virial
1425 25 : Pbc pbc;
1426 25 : pbc.setBox( box );
1427 25 : Tensor nvirial;
1428 100 : for(unsigned i=0; i<3; i++)
1429 300 : for(unsigned k=0; k<3; k++) {
1430 225 : double arg0=box(i,k);
1431 5535 : for(int j=0; j<natoms; ++j) {
1432 5310 : pos[j]=pbc.realToScaled( pos[j] );
1433 : }
1434 225 : cell[3*i+k]=box(i,k)=box(i,k)+delta;
1435 225 : pbc.setBox(box);
1436 5535 : for(int j=0; j<natoms; j++) {
1437 5310 : pos[j]=pbc.scaledToReal( pos[j] );
1438 : }
1439 450 : p.cmd("setStepLongLong",step);
1440 450 : p.cmd("setPositions",&pos[0][0],3*natoms);
1441 450 : p.cmd("setForces",&fake_forces[0],3*natoms);
1442 450 : p.cmd("setMasses",&masses[0],natoms);
1443 450 : p.cmd("setCharges",&charges[0],natoms);
1444 450 : p.cmd("setBox",&cell[0],9);
1445 450 : p.cmd("setVirial",&fake_virial[0],9);
1446 450 : p.cmd("prepareCalc");
1447 450 : p.cmd("performCalcNoUpdate");
1448 450 : p.cmd("getBias",&bias);
1449 225 : cell[3*i+k]=box(i,k)=arg0;
1450 225 : pbc.setBox(box);
1451 5535 : for(int j=0; j<natoms; j++)
1452 21240 : for(unsigned n=0; n<3; ++n) {
1453 15930 : pos[j][n]=coordinates[3*j+n];
1454 : }
1455 225 : nvirial(i,k) = ( bias - base ) / delta;
1456 : }
1457 25 : nvirial=-matmul(box.transpose(),nvirial);
1458 100 : for(unsigned i=0; i<3; i++)
1459 300 : for(unsigned k=0; k<3; k++) {
1460 225 : numder[3*natoms+3*i+k] = nvirial(i,k);
1461 : }
1462 :
1463 25 : }
1464 :
1465 : }
1466 : }
|
