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 "tools/Communicator.h"
27 : #include "tools/Random.h"
28 : #include "tools/Pbc.h"
29 : #include <cstdio>
30 : #include <cstring>
31 : #include <vector>
32 : #include <map>
33 : #include <memory>
34 : #include "tools/Units.h"
35 : #include "tools/PDB.h"
36 : #include "tools/FileBase.h"
37 : #include "tools/IFile.h"
38 : #include "xdrfile/xdrfile_trr.h"
39 : #include "xdrfile/xdrfile_xtc.h"
40 :
41 :
42 : // when using molfile plugin
43 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
44 : #ifndef __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS
45 : /* Use the internal ones. Alternatively:
46 : * ifeq (,$(findstring __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS,$(CPPFLAGS)))
47 : * CPPFLAGS+=-I../molfile
48 : */
49 : #include "molfile/libmolfile_plugin.h"
50 : #include "molfile/molfile_plugin.h"
51 : using namespace PLMD::molfile;
52 : #else
53 : #include <libmolfile_plugin.h>
54 : #include <molfile_plugin.h>
55 : #endif
56 : #endif
57 :
58 : namespace PLMD {
59 : namespace cltools {
60 :
61 : //+PLUMEDOC TOOLS driver-float
62 : /*
63 : Equivalent to \ref driver, but using single precision reals.
64 :
65 : The purpose of this tool is just to test what PLUMED does when linked from
66 : a single precision code.
67 :
68 : \par Examples
69 :
70 : \verbatim
71 : plumed driver-float --plumed plumed.dat --ixyz trajectory.xyz
72 : \endverbatim
73 :
74 : See also examples in \ref driver
75 :
76 : */
77 : //+ENDPLUMEDOC
78 : //
79 :
80 :
81 : //+PLUMEDOC TOOLS driver
82 : /*
83 : driver is a tool that allows one to to use plumed to post-process an existing trajectory.
84 :
85 : The input to driver is specified using the command line arguments described below.
86 :
87 : In addition, you can use the special \subpage READ command inside your plumed input
88 : to read in colvar files that were generated during your MD simulation. The values
89 : read in can then be treated like calculated colvars.
90 :
91 : \warning
92 : Notice that by default the driver has no knowledge about the masses and charges
93 : of your atoms! Thus, if you want to compute quantities depending charges (e.g. \ref DHENERGY)
94 : or masses (e.g. \ref COM) you should pass the proper information to the driver.
95 : You can do it either with the --pdb option or with the --mc option. The latter
96 : will read a file produced by \ref DUMPMASSCHARGE .
97 :
98 :
99 : \par Examples
100 :
101 : The following command tells plumed to post process the trajectory contained in `trajectory.xyz`
102 : by performing the actions described in the input file `plumed.dat`. If an action that takes the
103 : stride keyword is given a stride equal to \f$n\f$ then it will be performed only on every \f$n\f$th
104 : frames in the trajectory file.
105 : \verbatim
106 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz
107 : \endverbatim
108 :
109 : Notice that `xyz` files are expected to be in internal PLUMED units, that is by default nm.
110 : You can change this behavior by using the `--length-units` option:
111 : \verbatim
112 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --length-units A
113 : \endverbatim
114 : The strings accepted by the `--length-units` options are the same ones accepted by the \ref UNITS action.
115 : Other file formats typically have their default coordinates (e.g., `gro` files are always in nm)
116 : and it thus should not be necessary to use the `--length-units` option. Additionally,
117 : consider that the units used by the `driver` might be different by the units used in the PLUMED input
118 : file `plumed.dat`. For instance consider the command:
119 : \verbatim
120 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --length-units A
121 : \endverbatim
122 : where `plumed.dat` is
123 : \plumedfile
124 : # no explicit UNITS action here
125 : d: DISTANCE ATOMS=1,2
126 : PRINT ARG=d FILE=colvar
127 : \endplumedfile
128 : In this case, the driver reads the `xyz` file assuming it to contain coordinates in Angstrom units.
129 : However, the resulting `colvar` file contains a distance expressed in nm.
130 :
131 : The following command tells plumed to post process the trajectory contained in trajectory.xyz.
132 : by performing the actions described in the input file plumed.dat.
133 : \verbatim
134 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --trajectory-stride 100 --timestep 0.001
135 : \endverbatim
136 : Here though
137 : `--trajectory-stride` is set equal to the frequency with which frames were output during the trajectory
138 : and the `--timestep` is equal to the simulation timestep. As such the `STRIDE` parameters in the `plumed.dat`
139 : files are referred to the original timestep and any files output resemble those that would have been generated
140 : had we run the calculation we are running with driver when the MD simulation was running.
141 :
142 : PLUMED can read xyz files (in PLUMED units) and gro files (in nm). In addition,
143 : PLUMED includes by default support for a
144 : subset of the trajectory file formats supported by VMD, e.g. xtc and dcd:
145 :
146 : \verbatim
147 : plumed driver --plumed plumed.dat --pdb diala.pdb --mf_xtc traj.xtc --trajectory-stride 100 --timestep 0.001
148 : \endverbatim
149 :
150 : where `--mf_` prefixes the extension of one of the accepted molfile plugin format.
151 : If PLUMED has been \ref Installation "installed" with full molfile support, other formats will be available.
152 : Just type `plumed driver --help` to see which plugins are available.
153 :
154 : Molfile plugin require periodic cell to be triangular (i.e. first vector oriented along x and
155 : second vector in xy plane). This is true for many MD codes. However, it could be false
156 : if you rotate the coordinates in your trajectory before reading them in the driver.
157 : Also notice that some formats (e.g. amber crd) do not specify atom number. In this case you can use
158 : the `--natoms` option:
159 : \verbatim
160 : plumed driver --plumed plumed.dat --imf_crd trajectory.crd --natoms 128
161 : \endverbatim
162 :
163 : Check the available molfile plugins and limitations at [this link](http://www.ks.uiuc.edu/Research/vmd/plugins/molfile/).
164 :
165 : Additionally, you can use the xdrfile implementation of xtc and trr. To this aim, just
166 : download and install properly the xdrfile library (see [this link](http://www.gromacs.org/Developer_Zone/Programming_Guide/XTC_Library)).
167 : If the xdrfile library is installed properly the PLUMED configure script should be able to
168 : detect it and enable it.
169 : Notice that the xdrfile implementation of xtc and trr
170 : is more robust than the molfile one, since it provides support for generic cell shapes.
171 : In addition, it allows \ref DUMPATOMS to write compressed xtc files.
172 :
173 :
174 : */
175 : //+ENDPLUMEDOC
176 : //
177 :
178 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
179 : static std::vector<molfile_plugin_t *> plugins;
180 : static std::map <std::string, unsigned> pluginmap;
181 62514 : static int register_cb(void *v, vmdplugin_t *p) {
182 : //const char *key = p->name;
183 125028 : const auto ret = pluginmap.insert ( std::pair<std::string,unsigned>(std::string(p->name),plugins.size()) );
184 62514 : if (ret.second==false) {
185 : //cerr<<"MOLFILE: found duplicate plugin for "<<key<<" : not inserted "<<endl;
186 : } else {
187 : //cerr<<"MOLFILE: loading plugin "<<key<<" number "<<plugins.size()-1<<endl;
188 62514 : plugins.push_back(reinterpret_cast<molfile_plugin_t *>(p));
189 : }
190 62514 : return VMDPLUGIN_SUCCESS;
191 : }
192 : #endif
193 :
194 : template<typename real>
195 : class Driver : public CLTool {
196 : public:
197 : static void registerKeywords( Keywords& keys );
198 : explicit Driver(const CLToolOptions& co );
199 : int main(FILE* in,FILE*out,Communicator& pc) override;
200 : void evaluateNumericalDerivatives( const long int& step, PlumedMain& p, const std::vector<real>& coordinates,
201 : const std::vector<real>& masses, const std::vector<real>& charges,
202 : std::vector<real>& cell, const double& base, std::vector<real>& numder );
203 : std::string description()const override;
204 : };
205 :
206 : template<typename real>
207 6946 : void Driver<real>::registerKeywords( Keywords& keys ) {
208 6946 : CLTool::registerKeywords( keys ); keys.isDriver();
209 13892 : keys.addFlag("--help-debug",false,"print special options that can be used to create regtests");
210 13892 : keys.add("compulsory","--plumed","plumed.dat","specify the name of the plumed input file");
211 13892 : keys.add("compulsory","--timestep","1.0","the timestep that was used in the calculation that produced this trajectory in picoseconds");
212 20838 : keys.add("compulsory","--trajectory-stride","1","the frequency with which frames were output to this trajectory during the simulation"
213 : " (0 means that the number of the step is read from the trajectory file,"
214 : " currently working only for xtc/trr files read with --ixtc/--trr)"
215 : );
216 13892 : keys.add("compulsory","--multi","0","set number of replicas for multi environment (needs MPI)");
217 13892 : keys.addFlag("--noatoms",false,"don't read in a trajectory. Just use colvar files as specified in plumed.dat");
218 13892 : keys.addFlag("--parse-only",false,"read the plumed input file and stop");
219 13892 : keys.addFlag("--restart",false,"makes driver behave as if restarting");
220 13892 : keys.add("atoms","--ixyz","the trajectory in xyz format");
221 13892 : keys.add("atoms","--igro","the trajectory in gro format");
222 13892 : keys.add("atoms","--idlp4","the trajectory in DL_POLY_4 format");
223 13892 : keys.add("atoms","--ixtc","the trajectory in xtc format (xdrfile implementation)");
224 13892 : keys.add("atoms","--itrr","the trajectory in trr format (xdrfile implementation)");
225 13892 : keys.add("optional","--length-units","units for length, either as a string or a number");
226 13892 : keys.add("optional","--mass-units","units for mass in pdb and mc file, either as a string or a number");
227 13892 : keys.add("optional","--charge-units","units for charge in pdb and mc file, either as a string or a number");
228 13892 : keys.add("optional","--kt","set \\f$k_B T\\f$, it will not be necessary to specify temperature in input file");
229 13892 : keys.add("optional","--dump-forces","dump the forces on a file");
230 20838 : keys.add("optional","--dump-forces-fmt","( default=%%f ) the format to use to dump the forces");
231 13892 : keys.addFlag("--dump-full-virial",false,"with --dump-forces, it dumps the 9 components of the virial");
232 13892 : keys.add("optional","--pdb","provides a pdb with masses and charges");
233 13892 : keys.add("optional","--mc","provides a file with masses and charges as produced with DUMPMASSCHARGE");
234 13892 : keys.add("optional","--box","comma-separated box dimensions (3 for orthorhombic, 9 for generic)");
235 13892 : keys.add("optional","--natoms","provides number of atoms - only used if file format does not contain number of atoms");
236 13892 : keys.add("optional","--initial-step","provides a number for the initial step, default is 0");
237 13892 : keys.add("optional","--debug-forces","output a file containing the forces due to the bias evaluated using numerical derivatives "
238 : "and using the analytical derivatives implemented in plumed");
239 13892 : keys.add("hidden","--debug-float","[yes/no] turns on the single precision version (to check float interface)");
240 13892 : keys.add("hidden","--debug-dd","[yes/no] use a fake domain decomposition");
241 13892 : keys.add("hidden","--debug-pd","[yes/no] use a fake particle decomposition");
242 13892 : keys.add("hidden","--debug-grex","use a fake gromacs-like replica exchange, specify exchange stride");
243 20838 : keys.add("hidden","--debug-grex-log","log file for debug=grex");
244 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
245 6946 : MOLFILE_INIT_ALL
246 6946 : MOLFILE_REGISTER_ALL(NULL, register_cb)
247 69460 : for(unsigned i=0; i<plugins.size(); i++) {
248 125028 : std::string kk="--mf_"+std::string(plugins[i]->name);
249 125028 : std::string mm=" molfile: the trajectory in "+std::string(plugins[i]->name)+" format " ;
250 125028 : keys.add("atoms",kk,mm);
251 : }
252 : #endif
253 6946 : }
254 : template<typename real>
255 720 : Driver<real>::Driver(const CLToolOptions& co ):
256 720 : CLTool(co)
257 : {
258 720 : inputdata=commandline;
259 720 : }
260 : template<typename real>
261 4 : std::string Driver<real>::description()const { return "analyze trajectories with plumed"; }
262 :
263 : template<typename real>
264 712 : int Driver<real>::main(FILE* in,FILE*out,Communicator& pc) {
265 :
266 712 : Units units;
267 712 : PDB pdb;
268 :
269 : // Parse everything
270 712 : bool printhelpdebug; parseFlag("--help-debug",printhelpdebug);
271 712 : if( printhelpdebug ) {
272 : std::fprintf(out,"%s",
273 : "Additional options for debug (only to be used in regtest):\n"
274 : " [--debug-float yes] : turns on the single precision version (to check float interface)\n"
275 : " [--debug-dd yes] : use a fake domain decomposition\n"
276 : " [--debug-pd yes] : use a fake particle decomposition\n"
277 : );
278 : return 0;
279 : }
280 : // Are we reading trajectory data
281 712 : bool noatoms; parseFlag("--noatoms",noatoms);
282 712 : bool parseOnly; parseFlag("--parse-only",parseOnly);
283 1424 : bool restart; parseFlag("--restart",restart);
284 :
285 : std::string fakein;
286 : bool debug_float=false;
287 : fakein="";
288 1424 : if(parse("--debug-float",fakein)) {
289 0 : if(fakein=="yes") debug_float=true;
290 0 : else if(fakein=="no") debug_float=false;
291 0 : else error("--debug-float should have argument yes or no");
292 : }
293 :
294 : if(debug_float && sizeof(real)!=sizeof(float)) {
295 0 : auto cl=cltoolRegister().create(CLToolOptions("driver-float"));
296 : cl->setInputData(this->getInputData());
297 0 : int ret=cl->main(in,out,pc);
298 : return ret;
299 0 : }
300 :
301 : bool debug_pd=false;
302 : fakein="";
303 1424 : if(parse("--debug-pd",fakein)) {
304 12 : if(fakein=="yes") debug_pd=true;
305 0 : else if(fakein=="no") debug_pd=false;
306 0 : else error("--debug-pd should have argument yes or no");
307 : }
308 : if(debug_pd) std::fprintf(out,"DEBUGGING PARTICLE DECOMPOSITION\n");
309 :
310 : bool debug_dd=false;
311 : fakein="";
312 1424 : if(parse("--debug-dd",fakein)) {
313 60 : if(fakein=="yes") debug_dd=true;
314 0 : else if(fakein=="no") debug_dd=false;
315 0 : else error("--debug-dd should have argument yes or no");
316 : }
317 : if(debug_dd) std::fprintf(out,"DEBUGGING DOMAIN DECOMPOSITION\n");
318 :
319 712 : if( debug_pd || debug_dd ) {
320 72 : if(noatoms) error("cannot debug without atoms");
321 : }
322 :
323 : // set up for multi replica driver:
324 712 : int multi=0;
325 712 : parse("--multi",multi);
326 712 : Communicator intracomm;
327 712 : Communicator intercomm;
328 712 : if(multi) {
329 170 : int ntot=pc.Get_size();
330 170 : int nintra=ntot/multi;
331 170 : if(multi*nintra!=ntot) error("invalid number of processes for multi environment");
332 170 : pc.Split(pc.Get_rank()/nintra,pc.Get_rank(),intracomm);
333 170 : pc.Split(pc.Get_rank()%nintra,pc.Get_rank(),intercomm);
334 : } else {
335 542 : intracomm.Set_comm(pc.Get_comm());
336 : }
337 :
338 : // set up for debug replica exchange:
339 712 : bool debug_grex=parse("--debug-grex",fakein);
340 712 : int grex_stride=0;
341 : FILE*grex_log=NULL;
342 712 : if(debug_grex) {
343 30 : if(noatoms) error("must have atoms to debug_grex");
344 30 : if(multi<2) error("--debug_grex needs --multi with at least two replicas");
345 30 : Tools::convert(fakein,grex_stride);
346 30 : std::string n; Tools::convert(intercomm.Get_rank(),n);
347 : std::string file;
348 60 : parse("--debug-grex-log",file);
349 30 : if(file.length()>0) {
350 60 : file+="."+n;
351 30 : grex_log=fopen(file.c_str(),"w");
352 : }
353 : }
354 :
355 : // Read the plumed input file name
356 712 : std::string plumedFile; parse("--plumed",plumedFile);
357 : // the timestep
358 712 : double t; parse("--timestep",t);
359 712 : real timestep=real(t);
360 : // the stride
361 712 : unsigned stride; parse("--trajectory-stride",stride);
362 : // are we writing forces
363 712 : std::string dumpforces(""), debugforces(""), dumpforcesFmt("%f");;
364 712 : bool dumpfullvirial=false;
365 712 : if(!noatoms) {
366 658 : parse("--dump-forces",dumpforces);
367 1316 : parse("--debug-forces",debugforces);
368 : }
369 1434 : if(dumpforces!="" || debugforces!="" ) parse("--dump-forces-fmt",dumpforcesFmt);
370 1148 : if(dumpforces!="") parseFlag("--dump-full-virial",dumpfullvirial);
371 712 : if( debugforces!="" && (debug_dd || debug_pd) ) error("cannot debug forces and domain/particle decomposition at same time");
372 0 : if( debugforces!="" && sizeof(real)!=sizeof(double) ) error("cannot debug forces in single precision mode");
373 :
374 712 : real kt=-1.0;
375 1424 : parse("--kt",kt);
376 : std::string trajectory_fmt;
377 :
378 : bool use_molfile=false;
379 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
380 : molfile_plugin_t *api=NULL;
381 : void *h_in=NULL;
382 : molfile_timestep_t ts_in; // this is the structure that has the timestep
383 : // a std::unique_ptr<float> with the same scope as ts_in
384 : // it is necessary in order to store the pointer to ts_in.coords
385 : std::unique_ptr<float[]> ts_in_coords;
386 712 : ts_in.coords=ts_in_coords.get();
387 712 : ts_in.velocities=NULL;
388 712 : ts_in.A=-1; // we use this to check whether cell is provided or not
389 : #endif
390 :
391 : // Read in an xyz file
392 712 : std::string trajectoryFile(""), pdbfile(""), mcfile("");
393 712 : bool pbc_cli_given=false; std::vector<double> pbc_cli_box(9,0.0);
394 712 : int command_line_natoms=-1;
395 :
396 712 : if(!noatoms) {
397 1316 : std::string traj_xyz; parse("--ixyz",traj_xyz);
398 1316 : std::string traj_gro; parse("--igro",traj_gro);
399 1316 : std::string traj_dlp4; parse("--idlp4",traj_dlp4);
400 : std::string traj_xtc;
401 : std::string traj_trr;
402 658 : parse("--ixtc",traj_xtc);
403 658 : parse("--itrr",traj_trr);
404 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
405 6580 : for(unsigned i=0; i<plugins.size(); i++) {
406 11844 : std::string molfile_key="--mf_"+std::string(plugins[i]->name);
407 : std::string traj_molfile;
408 5922 : parse(molfile_key,traj_molfile);
409 5922 : if(traj_molfile.length()>0) {
410 159 : std::fprintf(out,"\nDRIVER: Found molfile format trajectory %s with name %s\n",plugins[i]->name,traj_molfile.c_str());
411 : trajectoryFile=traj_molfile;
412 318 : trajectory_fmt=std::string(plugins[i]->name);
413 : use_molfile=true;
414 159 : api = plugins[i];
415 : }
416 : }
417 : #endif
418 : { // check that only one fmt is specified
419 : int nn=0;
420 658 : if(traj_xyz.length()>0) nn++;
421 658 : if(traj_gro.length()>0) nn++;
422 658 : if(traj_dlp4.length()>0) nn++;
423 658 : if(traj_xtc.length()>0) nn++;
424 658 : if(traj_trr.length()>0) nn++;
425 658 : if(nn>1) {
426 0 : std::fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
427 0 : if(grex_log)fclose(grex_log);
428 0 : return 1;
429 : }
430 : }
431 658 : if(traj_xyz.length()>0 && trajectoryFile.length()==0) {
432 : trajectoryFile=traj_xyz;
433 : trajectory_fmt="xyz";
434 : }
435 658 : if(traj_gro.length()>0 && trajectoryFile.length()==0) {
436 : trajectoryFile=traj_gro;
437 : trajectory_fmt="gro";
438 : }
439 658 : if(traj_dlp4.length()>0 && trajectoryFile.length()==0) {
440 : trajectoryFile=traj_dlp4;
441 : trajectory_fmt="dlp4";
442 : }
443 658 : if(traj_xtc.length()>0 && trajectoryFile.length()==0) {
444 : trajectoryFile=traj_xtc;
445 : trajectory_fmt="xdr-xtc";
446 : }
447 658 : if(traj_trr.length()>0 && trajectoryFile.length()==0) {
448 : trajectoryFile=traj_trr;
449 : trajectory_fmt="xdr-trr";
450 : }
451 658 : if(trajectoryFile.length()==0&&!parseOnly) {
452 0 : std::fprintf(stderr,"ERROR: missing trajectory data\n");
453 0 : if(grex_log)fclose(grex_log);
454 0 : return 1;
455 : }
456 1316 : std::string lengthUnits(""); parse("--length-units",lengthUnits);
457 658 : if(lengthUnits.length()>0) units.setLength(lengthUnits);
458 1316 : std::string chargeUnits(""); parse("--charge-units",chargeUnits);
459 658 : if(chargeUnits.length()>0) units.setCharge(chargeUnits);
460 1316 : std::string massUnits(""); parse("--mass-units",massUnits);
461 658 : if(massUnits.length()>0) units.setMass(massUnits);
462 :
463 1316 : parse("--pdb",pdbfile);
464 658 : if(pdbfile.length()>0) {
465 75 : bool check=pdb.read(pdbfile,false,1.0);
466 75 : if(!check) error("error reading pdb file");
467 : }
468 :
469 1316 : parse("--mc",mcfile);
470 :
471 1316 : std::string pbc_cli_list; parse("--box",pbc_cli_list);
472 658 : if(pbc_cli_list.length()>0) {
473 : pbc_cli_given=true;
474 30 : std::vector<std::string> words=Tools::getWords(pbc_cli_list,",");
475 30 : if(words.size()==3) {
476 120 : for(int i=0; i<3; i++) std::sscanf(words[i].c_str(),"%100lf",&(pbc_cli_box[4*i]));
477 0 : } else if(words.size()==9) {
478 0 : for(int i=0; i<9; i++) std::sscanf(words[i].c_str(),"%100lf",&(pbc_cli_box[i]));
479 : } else {
480 0 : std::string msg="ERROR: cannot parse command-line box "+pbc_cli_list;
481 0 : std::fprintf(stderr,"%s\n",msg.c_str());
482 : return 1;
483 : }
484 :
485 30 : }
486 :
487 1316 : parse("--natoms",command_line_natoms);
488 :
489 : }
490 :
491 :
492 712 : if(debug_dd && debug_pd) error("cannot use debug-dd and debug-pd at the same time");
493 712 : if(debug_pd || debug_dd) {
494 72 : if( !Communicator::initialized() ) error("needs mpi for debug-pd");
495 : }
496 :
497 712 : PlumedMain p;
498 1424 : p.cmd("setRealPrecision",(int)sizeof(real));
499 : int checknatoms=-1;
500 712 : long int step=0;
501 712 : parse("--initial-step",step);
502 :
503 714 : if(restart) p.cmd("setRestart",1);
504 :
505 712 : if(Communicator::initialized()) {
506 303 : if(multi) {
507 414 : if(intracomm.Get_rank()==0) p.cmd("GREX setMPIIntercomm",&intercomm.Get_comm());
508 510 : p.cmd("GREX setMPIIntracomm",&intracomm.Get_comm());
509 340 : p.cmd("GREX init");
510 : }
511 909 : p.cmd("setMPIComm",&intracomm.Get_comm());
512 : }
513 1424 : p.cmd("setMDLengthUnits",units.getLength());
514 1424 : p.cmd("setMDChargeUnits",units.getCharge());
515 1424 : p.cmd("setMDMassUnits",units.getMass());
516 1424 : p.cmd("setMDEngine","driver");
517 1424 : p.cmd("setTimestep",timestep);
518 1424 : p.cmd("setPlumedDat",plumedFile.c_str());
519 1424 : p.cmd("setLog",out);
520 :
521 : int natoms;
522 712 : int lvl=0;
523 712 : int pb=1;
524 :
525 712 : if(parseOnly) {
526 0 : if(command_line_natoms<0) error("--parseOnly requires setting the number of atoms with --natoms");
527 0 : natoms=command_line_natoms;
528 : }
529 :
530 :
531 712 : FILE* fp=NULL; FILE* fp_forces=NULL; OFile fp_dforces;
532 : xdrfile::XDRFILE* xd=NULL;
533 712 : if(!noatoms&&!parseOnly) {
534 658 : if (trajectoryFile=="-")
535 : fp=in;
536 : else {
537 658 : if(multi) {
538 : std::string n;
539 170 : Tools::convert(intercomm.Get_rank(),n);
540 340 : std::string testfile=FileBase::appendSuffix(trajectoryFile,"."+n);
541 170 : FILE* tmp_fp=fopen(testfile.c_str(),"r");
542 170 : if(tmp_fp) { fclose(tmp_fp); trajectoryFile=testfile.c_str();}
543 : }
544 658 : if(use_molfile==true) {
545 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
546 159 : h_in = api->open_file_read(trajectoryFile.c_str(), trajectory_fmt.c_str(), &natoms);
547 159 : if(natoms==MOLFILE_NUMATOMS_UNKNOWN) {
548 2 : if(command_line_natoms>=0) natoms=command_line_natoms;
549 0 : else error("this file format does not provide number of atoms; use --natoms on the command line");
550 : }
551 159 : ts_in_coords=Tools::make_unique<float[]>(3*natoms);
552 159 : ts_in.coords = ts_in_coords.get();
553 : #endif
554 996 : } else if(trajectory_fmt=="xdr-xtc" || trajectory_fmt=="xdr-trr") {
555 5 : xd=xdrfile::xdrfile_open(trajectoryFile.c_str(),"r");
556 5 : if(!xd) {
557 0 : std::string msg="ERROR: Error opening trajectory file "+trajectoryFile;
558 0 : std::fprintf(stderr,"%s\n",msg.c_str());
559 : return 1;
560 : }
561 5 : if(trajectory_fmt=="xdr-xtc") xdrfile::read_xtc_natoms(&trajectoryFile[0],&natoms);
562 5 : if(trajectory_fmt=="xdr-trr") xdrfile::read_trr_natoms(&trajectoryFile[0],&natoms);
563 : } else {
564 494 : fp=fopen(trajectoryFile.c_str(),"r");
565 494 : if(!fp) {
566 0 : std::string msg="ERROR: Error opening trajectory file "+trajectoryFile;
567 0 : std::fprintf(stderr,"%s\n",msg.c_str());
568 : return 1;
569 : }
570 : }
571 : }
572 658 : if(dumpforces.length()>0) {
573 436 : if(Communicator::initialized() && pc.Get_size()>1) {
574 : std::string n;
575 226 : Tools::convert(pc.Get_rank(),n);
576 452 : dumpforces+="."+n;
577 : }
578 436 : fp_forces=fopen(dumpforces.c_str(),"w");
579 : }
580 658 : if(debugforces.length()>0) {
581 11 : if(Communicator::initialized() && pc.Get_size()>1) {
582 : std::string n;
583 6 : Tools::convert(pc.Get_rank(),n);
584 12 : debugforces+="."+n;
585 : }
586 11 : fp_dforces.open(debugforces);
587 : }
588 : }
589 :
590 : std::string line;
591 : std::vector<real> coordinates;
592 : std::vector<real> forces;
593 : std::vector<real> masses;
594 : std::vector<real> charges;
595 : std::vector<real> cell;
596 : std::vector<real> virial;
597 : std::vector<real> numder;
598 :
599 : // variables to test particle decomposition
600 : int pd_nlocal=0;
601 : int pd_start=0;
602 : // variables to test random decomposition (=domain decomposition)
603 : std::vector<int> dd_gatindex;
604 : std::vector<int> dd_g2l;
605 : std::vector<real> dd_masses;
606 : std::vector<real> dd_charges;
607 : std::vector<real> dd_forces;
608 : std::vector<real> dd_coordinates;
609 : int dd_nlocal=0;
610 : // random stream to choose decompositions
611 712 : Random rnd;
612 :
613 712 : if(trajectory_fmt=="dlp4") {
614 2 : if(!Tools::getline(fp,line)) error("error reading title");
615 2 : if(!Tools::getline(fp,line)) error("error reading atoms");
616 2 : std::sscanf(line.c_str(),"%d %d %d",&lvl,&pb,&natoms);
617 :
618 : }
619 : bool lstep=true;
620 235752 : while(true) {
621 236464 : if(!noatoms&&!parseOnly) {
622 36124 : if(use_molfile==true) {
623 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
624 : int rc;
625 10204 : rc = api->read_next_timestep(h_in, natoms, &ts_in);
626 10204 : if(rc==MOLFILE_EOF) {
627 : break;
628 : }
629 : #endif
630 28551 : } else if(trajectory_fmt=="xyz" || trajectory_fmt=="gro" || trajectory_fmt=="dlp4") {
631 25851 : if(!Tools::getline(fp,line)) break;
632 : }
633 : }
634 : bool first_step=false;
635 235815 : if(!noatoms&&!parseOnly) {
636 63334 : if(use_molfile==false && (trajectory_fmt=="xyz" || trajectory_fmt=="gro")) {
637 25341 : if(trajectory_fmt=="gro") if(!Tools::getline(fp,line)) error("premature end of trajectory file");
638 25341 : std::sscanf(line.c_str(),"%100d",&natoms);
639 : }
640 60905 : if(use_molfile==false && trajectory_fmt=="dlp4") {
641 : char xa[9];
642 : int xb,xc,xd;
643 : double t;
644 20 : std::sscanf(line.c_str(),"%8s %ld %d %d %d %lf",xa,&step,&xb,&xc,&xd,&t);
645 20 : if (lstep) {
646 4 : p.cmd("setTimestep",real(t));
647 : lstep = false;
648 : }
649 : }
650 : }
651 235815 : if(checknatoms<0 && !noatoms) {
652 658 : pd_nlocal=natoms;
653 : pd_start=0;
654 : first_step=true;
655 658 : masses.assign(natoms,std::numeric_limits<real>::quiet_NaN());
656 658 : charges.assign(natoms,std::numeric_limits<real>::quiet_NaN());
657 : //case pdb: structure
658 658 : if(pdbfile.length()>0) {
659 28092 : for(unsigned i=0; i<pdb.size(); ++i) {
660 28017 : AtomNumber an=pdb.getAtomNumbers()[i];
661 : unsigned index=an.index();
662 28017 : if( index>=unsigned(natoms) ) error("atom index in pdb exceeds the number of atoms in trajectory");
663 28017 : masses[index]=pdb.getOccupancy()[i];
664 28017 : charges[index]=pdb.getBeta()[i];
665 : }
666 : }
667 658 : if(mcfile.length()>0) {
668 5 : IFile ifile;
669 5 : ifile.open(mcfile);
670 : int index; double mass; double charge;
671 1138 : while(ifile.scanField("index",index).scanField("mass",mass).scanField("charge",charge).scanField()) {
672 564 : masses[index]=mass;
673 564 : charges[index]=charge;
674 : }
675 5 : }
676 235157 : } else if( checknatoms<0 && noatoms ) {
677 54 : natoms=0;
678 : }
679 235815 : if( checknatoms<0 ) {
680 712 : if(kt>=0) {
681 6 : p.cmd("setKbT",kt);
682 : }
683 712 : checknatoms=natoms;
684 1424 : p.cmd("setNatoms",natoms);
685 1424 : p.cmd("init");
686 712 : if(parseOnly) break;
687 : }
688 235815 : if(checknatoms!=natoms) {
689 0 : std::string stepstr; Tools::convert(step,stepstr);
690 0 : error("number of atoms in frame " + stepstr + " does not match number of atoms in first frame");
691 : }
692 :
693 235815 : coordinates.assign(3*natoms,real(0.0));
694 235815 : forces.assign(3*natoms,real(0.0));
695 235815 : cell.assign(9,real(0.0));
696 235815 : virial.assign(9,real(0.0));
697 :
698 235815 : if( first_step || rnd.U01()>0.5) {
699 119166 : if(debug_pd) {
700 152 : int npe=intracomm.Get_size();
701 152 : std::vector<int> loc(npe,0);
702 152 : std::vector<int> start(npe,0);
703 608 : for(int i=0; i<npe-1; i++) {
704 456 : int cc=(natoms*2*rnd.U01())/npe;
705 456 : if(start[i]+cc>natoms) cc=natoms-start[i];
706 456 : loc[i]=cc;
707 456 : start[i+1]=start[i]+loc[i];
708 : }
709 152 : loc[npe-1]=natoms-start[npe-1];
710 152 : intracomm.Bcast(loc,0);
711 152 : intracomm.Bcast(start,0);
712 152 : pd_nlocal=loc[intracomm.Get_rank()];
713 152 : pd_start=start[intracomm.Get_rank()];
714 152 : if(intracomm.Get_rank()==0) {
715 : std::fprintf(out,"\nDRIVER: Reassigning particle decomposition\n");
716 190 : std::fprintf(out,"DRIVER: "); for(int i=0; i<npe; i++) std::fprintf(out,"%d ",loc[i]); printf("\n");
717 190 : std::fprintf(out,"DRIVER: "); for(int i=0; i<npe; i++) std::fprintf(out,"%d ",start[i]); printf("\n");
718 : }
719 304 : p.cmd("setAtomsNlocal",pd_nlocal);
720 304 : p.cmd("setAtomsContiguous",pd_start);
721 119014 : } else if(debug_dd) {
722 956 : int npe=intracomm.Get_size();
723 956 : int rank=intracomm.Get_rank();
724 956 : dd_charges.assign(natoms,0.0);
725 956 : dd_masses.assign(natoms,0.0);
726 956 : dd_gatindex.assign(natoms,-1);
727 956 : dd_g2l.assign(natoms,-1);
728 956 : dd_coordinates.assign(3*natoms,0.0);
729 956 : dd_forces.assign(3*natoms,0.0);
730 : dd_nlocal=0;
731 53786 : for(int i=0; i<natoms; ++i) {
732 52830 : double r=rnd.U01()*npe;
733 112376 : int n; for(n=0; n<npe; n++) if(n+1>r)break;
734 52830 : plumed_assert(n<npe);
735 52830 : if(n==rank) {
736 19827 : dd_gatindex[dd_nlocal]=i;
737 19827 : dd_g2l[i]=dd_nlocal;
738 19827 : dd_charges[dd_nlocal]=charges[i];
739 19827 : dd_masses[dd_nlocal]=masses[i];
740 19827 : dd_nlocal++;
741 : }
742 : }
743 956 : if(intracomm.Get_rank()==0) {
744 : std::fprintf(out,"\nDRIVER: Reassigning domain decomposition\n");
745 : }
746 1912 : p.cmd("setAtomsNlocal",dd_nlocal);
747 1912 : p.cmd("setAtomsGatindex",&dd_gatindex[0],dd_nlocal);
748 : }
749 : }
750 :
751 235815 : int plumedStopCondition=0;
752 235815 : if(!noatoms) {
753 35475 : if(use_molfile) {
754 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
755 10045 : if(pbc_cli_given==false) {
756 10024 : if(ts_in.A>0.0) { // this is negative if molfile does not provide box
757 : // info on the cell: convert using pbcset.tcl from pbctools in vmd distribution
758 10021 : real cosBC=cos(real(ts_in.alpha)*pi/180.);
759 : //double sinBC=std::sin(ts_in.alpha*pi/180.);
760 10021 : real cosAC=std::cos(real(ts_in.beta)*pi/180.);
761 10021 : real cosAB=std::cos(real(ts_in.gamma)*pi/180.);
762 10021 : real sinAB=std::sin(real(ts_in.gamma)*pi/180.);
763 10021 : real Ax=real(ts_in.A);
764 10021 : real Bx=real(ts_in.B)*cosAB;
765 10021 : real By=real(ts_in.B)*sinAB;
766 10021 : real Cx=real(ts_in.C)*cosAC;
767 10021 : real Cy=(real(ts_in.C)*real(ts_in.B)*cosBC-Cx*Bx)/By;
768 10021 : real Cz=std::sqrt(real(ts_in.C)*real(ts_in.C)-Cx*Cx-Cy*Cy);
769 10021 : cell[0]=Ax/10.; cell[1]=0.; cell[2]=0.;
770 10021 : cell[3]=Bx/10.; cell[4]=By/10.; cell[5]=0.;
771 10021 : cell[6]=Cx/10.; cell[7]=Cy/10.; cell[8]=Cz/10.;
772 : } else {
773 3 : cell[0]=0.0; cell[1]=0.0; cell[2]=0.0;
774 3 : cell[3]=0.0; cell[4]=0.0; cell[5]=0.0;
775 3 : cell[6]=0.0; cell[7]=0.0; cell[8]=0.0;
776 : }
777 : } else {
778 210 : for(unsigned i=0; i<9; i++)cell[i]=pbc_cli_box[i];
779 : }
780 : // info on coords
781 : // the order is xyzxyz...
782 5795473 : for(int i=0; i<3*natoms; i++) {
783 5785428 : coordinates[i]=real(ts_in.coords[i])/real(10.); //convert to nm
784 : //cerr<<"COOR "<<coordinates[i]<<endl;
785 : }
786 : #endif
787 50830 : } else if(trajectory_fmt=="xdr-xtc" || trajectory_fmt=="xdr-trr") {
788 : int localstep;
789 : float time;
790 : xdrfile::matrix box;
791 69 : auto pos=Tools::make_unique<xdrfile::rvec[]>(natoms);
792 : float prec,lambda;
793 : int ret=xdrfile::exdrOK;
794 69 : if(trajectory_fmt=="xdr-xtc") ret=xdrfile::read_xtc(xd,natoms,&localstep,&time,box,pos.get(),&prec);
795 69 : if(trajectory_fmt=="xdr-trr") ret=xdrfile::read_trr(xd,natoms,&localstep,&time,&lambda,box,pos.get(),NULL,NULL);
796 69 : if(stride==0) step=localstep;
797 69 : if(ret==xdrfile::exdrENDOFFILE) break;
798 67 : if(ret!=xdrfile::exdrOK) break;
799 832 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) cell[3*i+j]=box[i][j];
800 10976 : for(int i=0; i<natoms; i++) for(unsigned j=0; j<3; j++)
801 8184 : coordinates[3*i+j]=real(pos[i][j]);
802 : } else {
803 25361 : if(trajectory_fmt=="xyz") {
804 23001 : if(!Tools::getline(fp,line)) error("premature end of trajectory file");
805 :
806 23001 : std::vector<double> celld(9,0.0);
807 23001 : if(pbc_cli_given==false) {
808 : std::vector<std::string> words;
809 45590 : words=Tools::getWords(line);
810 22795 : if(words.size()==3) {
811 22222 : std::sscanf(line.c_str(),"%100lf %100lf %100lf",&celld[0],&celld[4],&celld[8]);
812 573 : } else if(words.size()==9) {
813 573 : std::sscanf(line.c_str(),"%100lf %100lf %100lf %100lf %100lf %100lf %100lf %100lf %100lf",
814 : &celld[0], &celld[1], &celld[2],
815 : &celld[3], &celld[4], &celld[5],
816 : &celld[6], &celld[7], &celld[8]);
817 0 : } else error("needed box in second line of xyz file");
818 22795 : } else { // from command line
819 206 : celld=pbc_cli_box;
820 : }
821 230010 : for(unsigned i=0; i<9; i++)cell[i]=real(celld[i]);
822 : }
823 25361 : if(trajectory_fmt=="dlp4") {
824 20 : std::vector<double> celld(9,0.0);
825 20 : if(pbc_cli_given==false) {
826 20 : if(!Tools::getline(fp,line)) error("error reading vector a of cell");
827 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[0],&celld[1],&celld[2]);
828 20 : if(!Tools::getline(fp,line)) error("error reading vector b of cell");
829 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[3],&celld[4],&celld[5]);
830 20 : if(!Tools::getline(fp,line)) error("error reading vector c of cell");
831 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[6],&celld[7],&celld[8]);
832 : } else {
833 0 : celld=pbc_cli_box;
834 : }
835 200 : for(auto i=0; i<9; i++)cell[i]=real(celld[i])*0.1;
836 : }
837 : int ddist=0;
838 : // Read coordinates
839 2432085 : for(int i=0; i<natoms; i++) {
840 2406724 : bool ok=Tools::getline(fp,line);
841 2406724 : if(!ok) error("premature end of trajectory file");
842 : double cc[3];
843 2406724 : if(trajectory_fmt=="xyz") {
844 : char dummy[1000];
845 1693610 : int ret=std::sscanf(line.c_str(),"%999s %100lf %100lf %100lf",dummy,&cc[0],&cc[1],&cc[2]);
846 1693610 : if(ret!=4) error("cannot read line"+line);
847 713114 : } else if(trajectory_fmt=="gro") {
848 : // do the gromacs way
849 712474 : if(!i) {
850 : //
851 : // calculate the distance between dots (as in gromacs gmxlib/confio.c, routine get_w_conf )
852 : //
853 : const char *p1, *p2, *p3;
854 : p1 = std::strchr(line.c_str(), '.');
855 2340 : if (p1 == NULL) error("seems there are no coordinates in the gro file");
856 2340 : p2 = std::strchr(&p1[1], '.');
857 2340 : if (p2 == NULL) error("seems there is only one coordinates in the gro file");
858 2340 : ddist = p2 - p1;
859 2340 : p3 = std::strchr(&p2[1], '.');
860 2340 : if (p3 == NULL)error("seems there are only two coordinates in the gro file");
861 2340 : if (p3 - p2 != ddist)error("not uniform spacing in fields in the gro file");
862 : }
863 712474 : Tools::convert(line.substr(20,ddist),cc[0]);
864 712474 : Tools::convert(line.substr(20+ddist,ddist),cc[1]);
865 1424948 : Tools::convert(line.substr(20+ddist+ddist,ddist),cc[2]);
866 640 : } else if(trajectory_fmt=="dlp4") {
867 : char dummy[9];
868 : int idummy;
869 : double m,c;
870 640 : std::sscanf(line.c_str(),"%8s %d %lf %lf",dummy,&idummy,&m,&c);
871 640 : masses[i]=real(m);
872 640 : charges[i]=real(c);
873 640 : if(!Tools::getline(fp,line)) error("error reading coordinates");
874 640 : std::sscanf(line.c_str(),"%lf %lf %lf",&cc[0],&cc[1],&cc[2]);
875 640 : cc[0]*=0.1;
876 640 : cc[1]*=0.1;
877 640 : cc[2]*=0.1;
878 640 : if(lvl>0) {
879 640 : if(!Tools::getline(fp,line)) error("error skipping velocities");
880 : }
881 640 : if(lvl>1) {
882 640 : if(!Tools::getline(fp,line)) error("error skipping forces");
883 : }
884 0 : } else plumed_error();
885 2406724 : if(!debug_pd || ( i>=pd_start && i<pd_start+pd_nlocal) ) {
886 2387284 : coordinates[3*i]=real(cc[0]);
887 2387284 : coordinates[3*i+1]=real(cc[1]);
888 2387284 : coordinates[3*i+2]=real(cc[2]);
889 : }
890 : }
891 25361 : if(trajectory_fmt=="gro") {
892 2340 : if(!Tools::getline(fp,line)) error("premature end of trajectory file");
893 2340 : std::vector<std::string> words=Tools::getWords(line);
894 2340 : if(words.size()<3) error("cannot understand box format");
895 2340 : Tools::convert(words[0],cell[0]);
896 2340 : Tools::convert(words[1],cell[4]);
897 2340 : Tools::convert(words[2],cell[8]);
898 2340 : if(words.size()>3) Tools::convert(words[3],cell[1]);
899 2340 : if(words.size()>4) Tools::convert(words[4],cell[2]);
900 2340 : if(words.size()>5) Tools::convert(words[5],cell[3]);
901 2340 : if(words.size()>6) Tools::convert(words[6],cell[5]);
902 2340 : if(words.size()>7) Tools::convert(words[7],cell[6]);
903 2340 : if(words.size()>8) Tools::convert(words[8],cell[7]);
904 2340 : }
905 :
906 : }
907 :
908 70940 : p.cmd("setStepLong",step);
909 70940 : p.cmd("setStopFlag",&plumedStopCondition);
910 :
911 35470 : if(debug_dd) {
912 38944 : for(int i=0; i<dd_nlocal; ++i) {
913 37156 : int kk=dd_gatindex[i];
914 37156 : dd_coordinates[3*i+0]=coordinates[3*kk+0];
915 37156 : dd_coordinates[3*i+1]=coordinates[3*kk+1];
916 37156 : dd_coordinates[3*i+2]=coordinates[3*kk+2];
917 : }
918 3576 : p.cmd("setForces",&dd_forces[0],3*dd_nlocal);
919 3576 : p.cmd("setPositions",&dd_coordinates[0],3*dd_nlocal);
920 3576 : p.cmd("setMasses",&dd_masses[0],dd_nlocal);
921 3576 : p.cmd("setCharges",&dd_charges[0],dd_nlocal);
922 : } else {
923 : // this is required to avoid troubles when the last domain
924 : // contains zero atoms
925 : // Basically, for empty domains we pass null pointers
926 : #define fix_pd(xx) (pd_nlocal!=0?&xx:NULL)
927 101046 : p.cmd("setForces",fix_pd(forces[3*pd_start]),3*pd_nlocal);
928 101046 : p.cmd("setPositions",fix_pd(coordinates[3*pd_start]),3*pd_nlocal);
929 101046 : p.cmd("setMasses",fix_pd(masses[pd_start]),pd_nlocal);
930 101046 : p.cmd("setCharges",fix_pd(charges[pd_start]),pd_nlocal);
931 : }
932 70940 : p.cmd("setBox",cell.data(),9);
933 70940 : p.cmd("setVirial",virial.data(),9);
934 : } else {
935 400680 : p.cmd("setStepLong",step);
936 400680 : p.cmd("setStopFlag",&plumedStopCondition);
937 : }
938 471620 : p.cmd("calc");
939 235810 : if(debugforces.length()>0) {
940 25 : virial.assign(9,real(0.0));
941 25 : forces.assign(3*natoms,real(0.0));
942 50 : p.cmd("prepareCalc");
943 50 : p.cmd("performCalcNoUpdate");
944 : }
945 :
946 : // this is necessary as only processor zero is adding to the virial:
947 235810 : intracomm.Bcast(virial,0);
948 235810 : if(debug_pd) intracomm.Sum(forces);
949 235810 : if(debug_dd) {
950 38944 : for(int i=0; i<dd_nlocal; i++) {
951 37156 : forces[3*dd_gatindex[i]+0]=dd_forces[3*i+0];
952 37156 : forces[3*dd_gatindex[i]+1]=dd_forces[3*i+1];
953 37156 : forces[3*dd_gatindex[i]+2]=dd_forces[3*i+2];
954 : }
955 1788 : dd_forces.assign(3*natoms,0.0);
956 1788 : intracomm.Sum(forces);
957 : }
958 235810 : if(debug_grex &&step%grex_stride==0) {
959 228 : p.cmd("GREX savePositions");
960 114 : if(intracomm.Get_rank()>0) {
961 114 : p.cmd("GREX prepare");
962 : } else {
963 57 : int r=intercomm.Get_rank();
964 57 : int n=intercomm.Get_size();
965 57 : int partner=r+(2*((r+step/grex_stride)%2))-1;
966 : if(partner<0)partner=0;
967 57 : if(partner>=n) partner=n-1;
968 114 : p.cmd("GREX setPartner",partner);
969 114 : p.cmd("GREX calculate");
970 114 : p.cmd("GREX shareAllDeltaBias");
971 228 : for(int i=0; i<n; i++) {
972 171 : std::string s; Tools::convert(i,s);
973 342 : real a=std::numeric_limits<real>::quiet_NaN(); s="GREX getDeltaBias "+s; p.cmd(s,&a);
974 171 : if(grex_log) std::fprintf(grex_log," %f",a);
975 : }
976 57 : if(grex_log) std::fprintf(grex_log,"\n");
977 : }
978 : }
979 :
980 :
981 235810 : if(fp_forces) {
982 22249 : std::fprintf(fp_forces,"%d\n",natoms);
983 44498 : std::string fmtv=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
984 44498 : std::string fmt=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
985 22249 : if(dumpfullvirial) {
986 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]);
987 : } else {
988 21899 : std::fprintf(fp_forces,fmt.c_str(),virial[0],virial[4],virial[8]);
989 : }
990 22249 : fmt="X "+fmt;
991 2153450 : for(int i=0; i<natoms; i++)
992 2131201 : std::fprintf(fp_forces,fmt.c_str(),forces[3*i],forces[3*i+1],forces[3*i+2]);
993 : }
994 235810 : if(debugforces.length()>0) {
995 : // Now call the routine to work out the derivatives numerically
996 25 : numder.assign(3*natoms+9,real(0.0)); real base=0;
997 50 : p.cmd("getBias",&base);
998 25 : if( std::abs(base)<epsilon ) printf("WARNING: bias for configuration appears to be zero so debugging forces is trivial");
999 25 : evaluateNumericalDerivatives( step, p, coordinates, masses, charges, cell, base, numder );
1000 :
1001 : // And output everything to a file
1002 25 : fp_dforces.fmtField(" " + dumpforcesFmt);
1003 1795 : for(int i=0; i<3*natoms; ++i) {
1004 1770 : fp_dforces.printField("parameter",(int)i);
1005 3540 : fp_dforces.printField("analytical",forces[i]);
1006 1770 : fp_dforces.printField("numerical",-numder[i]);
1007 1770 : fp_dforces.printField();
1008 : }
1009 : // And print the virial
1010 250 : for(int i=0; i<9; ++i) {
1011 225 : fp_dforces.printField("parameter",3*natoms+i );
1012 225 : fp_dforces.printField("analytical",virial[i] );
1013 225 : fp_dforces.printField("numerical",-numder[3*natoms+i]);
1014 225 : fp_dforces.printField();
1015 : }
1016 : }
1017 :
1018 235810 : if(plumedStopCondition) break;
1019 :
1020 235752 : step+=stride;
1021 : }
1022 1424 : if(!parseOnly) p.cmd("runFinalJobs");
1023 :
1024 712 : if(fp_forces) fclose(fp_forces);
1025 712 : if(debugforces.length()>0) fp_dforces.close();
1026 712 : if(fp && fp!=in)fclose(fp);
1027 712 : if(xd) xdrfile::xdrfile_close(xd);
1028 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
1029 712 : if(h_in) api->close_file_read(h_in);
1030 : #endif
1031 712 : if(grex_log) fclose(grex_log);
1032 :
1033 : return 0;
1034 2136 : }
1035 :
1036 : template<typename real>
1037 25 : void Driver<real>::evaluateNumericalDerivatives( const long int& step, PlumedMain& p, const std::vector<real>& coordinates,
1038 : const std::vector<real>& masses, const std::vector<real>& charges,
1039 : std::vector<real>& cell, const double& base, std::vector<real>& numder ) {
1040 :
1041 25 : int natoms = coordinates.size() / 3; real delta = std::sqrt(epsilon);
1042 25 : std::vector<Vector> pos(natoms); real bias=0;
1043 25 : std::vector<real> fake_forces( 3*natoms ), fake_virial(9);
1044 615 : for(int i=0; i<natoms; ++i) {
1045 2360 : for(unsigned j=0; j<3; ++j) pos[i][j]=coordinates[3*i+j];
1046 : }
1047 :
1048 615 : for(int i=0; i<natoms; ++i) {
1049 2360 : for(unsigned j=0; j<3; ++j) {
1050 1770 : pos[i][j]=pos[i][j]+delta;
1051 3540 : p.cmd("setStepLong",step);
1052 3540 : p.cmd("setPositions",&pos[0][0],3*natoms);
1053 3540 : p.cmd("setForces",&fake_forces[0],3*natoms);
1054 3540 : p.cmd("setMasses",&masses[0],natoms);
1055 3540 : p.cmd("setCharges",&charges[0],natoms);
1056 3540 : p.cmd("setBox",&cell[0],9);
1057 3540 : p.cmd("setVirial",&fake_virial[0],9);
1058 3540 : p.cmd("prepareCalc");
1059 3540 : p.cmd("performCalcNoUpdate");
1060 3540 : p.cmd("getBias",&bias);
1061 1770 : pos[i][j]=coordinates[3*i+j];
1062 1770 : numder[3*i+j] = (bias - base) / delta;
1063 : }
1064 : }
1065 :
1066 : // Create the cell
1067 25 : Tensor box( cell[0], cell[1], cell[2], cell[3], cell[4], cell[5], cell[6], cell[7], cell[8] );
1068 : // And the virial
1069 25 : Pbc pbc; pbc.setBox( box ); Tensor nvirial;
1070 325 : for(unsigned i=0; i<3; i++) for(unsigned k=0; k<3; k++) {
1071 225 : double arg0=box(i,k);
1072 5535 : for(int j=0; j<natoms; ++j) pos[j]=pbc.realToScaled( pos[j] );
1073 225 : cell[3*i+k]=box(i,k)=box(i,k)+delta; pbc.setBox(box);
1074 5535 : for(int j=0; j<natoms; j++) pos[j]=pbc.scaledToReal( pos[j] );
1075 450 : p.cmd("setStepLong",step);
1076 450 : p.cmd("setPositions",&pos[0][0],3*natoms);
1077 450 : p.cmd("setForces",&fake_forces[0],3*natoms);
1078 450 : p.cmd("setMasses",&masses[0],natoms);
1079 450 : p.cmd("setCharges",&charges[0],natoms);
1080 450 : p.cmd("setBox",&cell[0],9);
1081 450 : p.cmd("setVirial",&fake_virial[0],9);
1082 450 : p.cmd("prepareCalc");
1083 450 : p.cmd("performCalcNoUpdate");
1084 450 : p.cmd("getBias",&bias);
1085 225 : cell[3*i+k]=box(i,k)=arg0; pbc.setBox(box);
1086 21465 : for(int j=0; j<natoms; j++) for(unsigned n=0; n<3; ++n) pos[j][n]=coordinates[3*j+n];
1087 225 : nvirial(i,k) = ( bias - base ) / delta;
1088 : }
1089 25 : nvirial=-matmul(box.transpose(),nvirial);
1090 325 : for(unsigned i=0; i<3; i++) for(unsigned k=0; k<3; k++) numder[3*natoms+3*i+k] = nvirial(i,k);
1091 :
1092 50 : }
1093 :
1094 : }
1095 : }
|
