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