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