Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "CLTool.h"
23 : #include "core/CLToolRegister.h"
24 : #include "tools/Tools.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/ActionSet.h"
27 : #include "core/ActionWithValue.h"
28 : #include "core/ActionWithVirtualAtom.h"
29 : #include "core/ActionShortcut.h"
30 : #include "tools/Communicator.h"
31 : #include "tools/Random.h"
32 : #include "tools/Pbc.h"
33 : #include <cstdio>
34 : #include <cstring>
35 : #include <vector>
36 : #include <map>
37 : #include <memory>
38 : #include "tools/Units.h"
39 : #include "tools/PDB.h"
40 : #include "tools/FileBase.h"
41 : #include "tools/IFile.h"
42 : #include "xdrfile/xdrfile_trr.h"
43 : #include "xdrfile/xdrfile_xtc.h"
44 :
45 :
46 : // when using molfile plugin
47 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
48 : #ifndef __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS
49 : /* Use the internal ones. Alternatively:
50 : * ifeq (,$(findstring __PLUMED_HAS_EXTERNAL_MOLFILE_PLUGINS,$(CPPFLAGS)))
51 : * CPPFLAGS+=-I../molfile
52 : */
53 : #include "molfile/libmolfile_plugin.h"
54 : #include "molfile/molfile_plugin.h"
55 : using namespace PLMD::molfile;
56 : #else
57 : #include <libmolfile_plugin.h>
58 : #include <molfile_plugin.h>
59 : #endif
60 : #endif
61 :
62 : namespace PLMD {
63 : namespace cltools {
64 :
65 : //+PLUMEDOC TOOLS driver-float
66 : /*
67 : Equivalent to 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. The documentation is identical to that for \ref driver
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 : static std::vector<molfile_plugin_t *> plugins;
184 : static std::map <std::string, unsigned> pluginmap;
185 95688 : static int register_cb(void *v, vmdplugin_t *p) {
186 : //const char *key = p->name;
187 191376 : const auto ret = pluginmap.insert ( std::pair<std::string,unsigned>(std::string(p->name),plugins.size()) );
188 95688 : 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 95688 : plugins.push_back(reinterpret_cast<molfile_plugin_t *>(p));
193 : }
194 95688 : return VMDPLUGIN_SUCCESS;
195 : }
196 : #endif
197 :
198 : template<typename real>
199 : 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 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 : std::string description()const override;
208 : };
209 :
210 : template<typename real>
211 10632 : void Driver<real>::registerKeywords( Keywords& keys ) {
212 10632 : CLTool::registerKeywords( keys ); keys.isDriver();
213 21264 : keys.addFlag("--help-debug",false,"print special options that can be used to create regtests");
214 21264 : keys.add("compulsory","--plumed","plumed.dat","specify the name of the plumed input file");
215 21264 : keys.add("compulsory","--timestep","1.0","the timestep that was used in the calculation that produced this trajectory in picoseconds");
216 21264 : keys.add("compulsory","--trajectory-stride","1","the frequency with which frames were output to this trajectory during the simulation"
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 : );
220 21264 : keys.add("compulsory","--multi","0","set number of replicas for multi environment (needs MPI)");
221 21264 : keys.addFlag("--noatoms",false,"don't read in a trajectory. Just use colvar files as specified in plumed.dat");
222 21264 : keys.addFlag("--parse-only",false,"read the plumed input file and stop");
223 21264 : keys.addFlag("--restart",false,"makes driver behave as if restarting");
224 21264 : keys.add("atoms","--ixyz","the trajectory in xyz format");
225 21264 : keys.add("atoms","--igro","the trajectory in gro format");
226 21264 : keys.add("atoms","--idlp4","the trajectory in DL_POLY_4 format");
227 21264 : keys.add("atoms","--ixtc","the trajectory in xtc format (xdrfile implementation)");
228 21264 : keys.add("atoms","--itrr","the trajectory in trr format (xdrfile implementation)");
229 21264 : keys.add("optional","--shortcut-ofile","the name of the file to output info on the way shortcuts have been expanded. If there are no shortcuts in your input file nothing is output");
230 21264 : keys.add("optional","--valuedict-ofile","output a dictionary giving information about each value in the input file");
231 21264 : keys.add("optional","--length-units","units for length, either as a string or a number");
232 21264 : keys.add("optional","--mass-units","units for mass in pdb and mc file, either as a string or a number");
233 21264 : keys.add("optional","--charge-units","units for charge in pdb and mc file, either as a string or a number");
234 21264 : keys.add("optional","--kt","set \\f$k_B T\\f$, it will not be necessary to specify temperature in input file");
235 21264 : keys.add("optional","--dump-forces","dump the forces on a file");
236 21264 : keys.add("optional","--dump-forces-fmt","( default=%%f ) the format to use to dump the forces");
237 21264 : keys.addFlag("--dump-full-virial",false,"with --dump-forces, it dumps the 9 components of the virial");
238 21264 : keys.add("optional","--pdb","provides a pdb with masses and charges");
239 21264 : keys.add("optional","--mc","provides a file with masses and charges as produced with DUMPMASSCHARGE");
240 21264 : keys.add("optional","--box","comma-separated box dimensions (3 for orthorhombic, 9 for generic)");
241 21264 : keys.add("optional","--natoms","provides number of atoms - only used if file format does not contain number of atoms");
242 21264 : keys.add("optional","--initial-step","provides a number for the initial step, default is 0");
243 21264 : keys.add("optional","--debug-forces","output a file containing the forces due to the bias evaluated using numerical derivatives "
244 : "and using the analytical derivatives implemented in plumed");
245 21264 : keys.add("hidden","--debug-float","[yes/no] turns on the single precision version (to check float interface)");
246 21264 : keys.add("hidden","--debug-dd","[yes/no] use a fake domain decomposition");
247 21264 : keys.add("hidden","--debug-pd","[yes/no] use a fake particle decomposition");
248 21264 : keys.add("hidden","--debug-grex","use a fake gromacs-like replica exchange, specify exchange stride");
249 21264 : keys.add("hidden","--debug-grex-log","log file for debug=grex");
250 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
251 10632 : MOLFILE_INIT_ALL
252 10632 : MOLFILE_REGISTER_ALL(NULL, register_cb)
253 106320 : for(unsigned i=0; i<plugins.size(); i++) {
254 191376 : std::string kk="--mf_"+std::string(plugins[i]->name);
255 191376 : std::string mm=" molfile: the trajectory in "+std::string(plugins[i]->name)+" format " ;
256 191376 : keys.add("atoms",kk,mm);
257 : }
258 : #endif
259 10632 : }
260 : template<typename real>
261 950 : Driver<real>::Driver(const CLToolOptions& co ):
262 950 : CLTool(co)
263 : {
264 950 : inputdata=commandline;
265 950 : }
266 : template<typename real>
267 4 : std::string Driver<real>::description()const { return "analyze trajectories with plumed"; }
268 :
269 : template<typename real>
270 942 : int Driver<real>::main(FILE* in,FILE*out,Communicator& pc) {
271 :
272 942 : Units units;
273 942 : PDB pdb;
274 :
275 : // Parse everything
276 942 : bool printhelpdebug; parseFlag("--help-debug",printhelpdebug);
277 942 : if( printhelpdebug ) {
278 : std::fprintf(out,"%s",
279 : "Additional options for debug (only to be used in regtest):\n"
280 : " [--debug-float yes] : turns on the single precision version (to check float interface)\n"
281 : " [--debug-dd yes] : use a fake domain decomposition\n"
282 : " [--debug-pd yes] : use a fake particle decomposition\n"
283 : );
284 : return 0;
285 : }
286 : // Are we reading trajectory data
287 942 : bool noatoms; parseFlag("--noatoms",noatoms);
288 1884 : bool parseOnly; parseFlag("--parse-only",parseOnly);
289 1884 : std::string full_outputfile; parse("--shortcut-ofile",full_outputfile);
290 942 : std::string valuedict_file; parse("--valuedict-ofile",valuedict_file);
291 1884 : bool restart; parseFlag("--restart",restart);
292 :
293 : std::string fakein;
294 : bool debug_float=false;
295 : fakein="";
296 1884 : if(parse("--debug-float",fakein)) {
297 0 : if(fakein=="yes") debug_float=true;
298 0 : else if(fakein=="no") debug_float=false;
299 0 : else error("--debug-float should have argument yes or no");
300 : }
301 :
302 : if(debug_float && sizeof(real)!=sizeof(float)) {
303 0 : auto cl=cltoolRegister().create(CLToolOptions("driver-float"));
304 : cl->setInputData(this->getInputData());
305 0 : int ret=cl->main(in,out,pc);
306 : return ret;
307 0 : }
308 :
309 : bool debug_pd=false;
310 : fakein="";
311 1884 : if(parse("--debug-pd",fakein)) {
312 12 : if(fakein=="yes") debug_pd=true;
313 0 : else if(fakein=="no") debug_pd=false;
314 0 : else error("--debug-pd should have argument yes or no");
315 : }
316 : if(debug_pd) std::fprintf(out,"DEBUGGING PARTICLE DECOMPOSITION\n");
317 :
318 : bool debug_dd=false;
319 : fakein="";
320 1884 : if(parse("--debug-dd",fakein)) {
321 60 : if(fakein=="yes") debug_dd=true;
322 0 : else if(fakein=="no") debug_dd=false;
323 0 : else error("--debug-dd should have argument yes or no");
324 : }
325 : if(debug_dd) std::fprintf(out,"DEBUGGING DOMAIN DECOMPOSITION\n");
326 :
327 942 : if( debug_pd || debug_dd ) {
328 72 : if(noatoms) error("cannot debug without atoms");
329 : }
330 :
331 : // set up for multi replica driver:
332 942 : int multi=0;
333 942 : parse("--multi",multi);
334 942 : Communicator intracomm;
335 942 : Communicator intercomm;
336 942 : if(multi) {
337 170 : int ntot=pc.Get_size();
338 170 : int nintra=ntot/multi;
339 170 : if(multi*nintra!=ntot) error("invalid number of processes for multi environment");
340 170 : pc.Split(pc.Get_rank()/nintra,pc.Get_rank(),intracomm);
341 170 : pc.Split(pc.Get_rank()%nintra,pc.Get_rank(),intercomm);
342 : } else {
343 772 : intracomm.Set_comm(pc.Get_comm());
344 : }
345 :
346 : // set up for debug replica exchange:
347 942 : bool debug_grex=parse("--debug-grex",fakein);
348 942 : int grex_stride=0;
349 : FILE*grex_log=NULL;
350 : // call fclose when fp goes out of scope
351 1184 : auto deleter=[](auto f) { if(f) std::fclose(f); };
352 : std::unique_ptr<FILE,decltype(deleter)> grex_log_deleter(grex_log,deleter);
353 :
354 942 : if(debug_grex) {
355 30 : if(noatoms) error("must have atoms to debug_grex");
356 30 : if(multi<2) error("--debug_grex needs --multi with at least two replicas");
357 30 : Tools::convert(fakein,grex_stride);
358 30 : std::string n; Tools::convert(intercomm.Get_rank(),n);
359 : std::string file;
360 60 : parse("--debug-grex-log",file);
361 30 : if(file.length()>0) {
362 60 : file+="."+n;
363 30 : grex_log=std::fopen(file.c_str(),"w");
364 : grex_log_deleter.reset(grex_log);
365 : }
366 : }
367 :
368 : // Read the plumed input file name
369 942 : std::string plumedFile; parse("--plumed",plumedFile);
370 : // the timestep
371 942 : double t; parse("--timestep",t);
372 942 : real timestep=real(t);
373 : // the stride
374 942 : unsigned stride; parse("--trajectory-stride",stride);
375 : // are we writing forces
376 942 : std::string dumpforces(""), debugforces(""), dumpforcesFmt("%f");;
377 942 : bool dumpfullvirial=false;
378 942 : if(!noatoms) {
379 886 : parse("--dump-forces",dumpforces);
380 1772 : parse("--debug-forces",debugforces);
381 : }
382 1894 : if(dumpforces!="" || debugforces!="" ) parse("--dump-forces-fmt",dumpforcesFmt);
383 1495 : if(dumpforces!="") parseFlag("--dump-full-virial",dumpfullvirial);
384 942 : if( debugforces!="" && (debug_dd || debug_pd) ) error("cannot debug forces and domain/particle decomposition at same time");
385 0 : if( debugforces!="" && sizeof(real)!=sizeof(double) ) error("cannot debug forces in single precision mode");
386 :
387 942 : real kt=-1.0;
388 1884 : parse("--kt",kt);
389 : std::string trajectory_fmt;
390 :
391 : bool use_molfile=false;
392 : molfile_plugin_t *api=NULL;
393 :
394 : // Read in an xyz file
395 942 : std::string trajectoryFile(""), pdbfile(""), mcfile("");
396 942 : bool pbc_cli_given=false; std::vector<double> pbc_cli_box(9,0.0);
397 942 : int command_line_natoms=-1;
398 :
399 942 : if(!noatoms) {
400 1772 : std::string traj_xyz; parse("--ixyz",traj_xyz);
401 1772 : std::string traj_gro; parse("--igro",traj_gro);
402 1772 : std::string traj_dlp4; parse("--idlp4",traj_dlp4);
403 : std::string traj_xtc;
404 : std::string traj_trr;
405 886 : parse("--ixtc",traj_xtc);
406 886 : parse("--itrr",traj_trr);
407 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
408 8860 : for(unsigned i=0; i<plugins.size(); i++) {
409 15948 : std::string molfile_key="--mf_"+std::string(plugins[i]->name);
410 : std::string traj_molfile;
411 7974 : parse(molfile_key,traj_molfile);
412 7974 : if(traj_molfile.length()>0) {
413 269 : std::fprintf(out,"\nDRIVER: Found molfile format trajectory %s with name %s\n",plugins[i]->name,traj_molfile.c_str());
414 : trajectoryFile=traj_molfile;
415 538 : trajectory_fmt=std::string(plugins[i]->name);
416 : use_molfile=true;
417 269 : api = plugins[i];
418 : }
419 : }
420 : #endif
421 : { // check that only one fmt is specified
422 : int nn=0;
423 886 : if(traj_xyz.length()>0) nn++;
424 886 : if(traj_gro.length()>0) nn++;
425 886 : if(traj_dlp4.length()>0) nn++;
426 886 : if(traj_xtc.length()>0) nn++;
427 886 : if(traj_trr.length()>0) nn++;
428 886 : if(nn>1) {
429 0 : std::fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
430 : return 1;
431 : }
432 : }
433 886 : if(traj_xyz.length()>0 && trajectoryFile.length()==0) {
434 : trajectoryFile=traj_xyz;
435 : trajectory_fmt="xyz";
436 : }
437 886 : if(traj_gro.length()>0 && trajectoryFile.length()==0) {
438 : trajectoryFile=traj_gro;
439 : trajectory_fmt="gro";
440 : }
441 886 : if(traj_dlp4.length()>0 && trajectoryFile.length()==0) {
442 : trajectoryFile=traj_dlp4;
443 : trajectory_fmt="dlp4";
444 : }
445 886 : if(traj_xtc.length()>0 && trajectoryFile.length()==0) {
446 : trajectoryFile=traj_xtc;
447 : trajectory_fmt="xdr-xtc";
448 : }
449 886 : if(traj_trr.length()>0 && trajectoryFile.length()==0) {
450 : trajectoryFile=traj_trr;
451 : trajectory_fmt="xdr-trr";
452 : }
453 886 : if(trajectoryFile.length()==0&&!parseOnly) {
454 0 : std::fprintf(stderr,"ERROR: missing trajectory data\n");
455 : return 1;
456 : }
457 1772 : std::string lengthUnits(""); parse("--length-units",lengthUnits);
458 886 : if(lengthUnits.length()>0) units.setLength(lengthUnits);
459 1772 : std::string chargeUnits(""); parse("--charge-units",chargeUnits);
460 886 : if(chargeUnits.length()>0) units.setCharge(chargeUnits);
461 1772 : std::string massUnits(""); parse("--mass-units",massUnits);
462 886 : if(massUnits.length()>0) units.setMass(massUnits);
463 :
464 1772 : parse("--pdb",pdbfile);
465 886 : if(pdbfile.length()>0) {
466 79 : bool check=pdb.read(pdbfile,false,1.0);
467 79 : if(!check) error("error reading pdb file");
468 : }
469 :
470 1772 : parse("--mc",mcfile);
471 :
472 1772 : std::string pbc_cli_list; parse("--box",pbc_cli_list);
473 886 : if(pbc_cli_list.length()>0) {
474 : pbc_cli_given=true;
475 38 : std::vector<std::string> words=Tools::getWords(pbc_cli_list,",");
476 38 : if(words.size()==3) {
477 152 : for(int i=0; i<3; i++) std::sscanf(words[i].c_str(),"%100lf",&(pbc_cli_box[4*i]));
478 0 : } else if(words.size()==9) {
479 0 : for(int i=0; i<9; i++) std::sscanf(words[i].c_str(),"%100lf",&(pbc_cli_box[i]));
480 : } else {
481 0 : std::string msg="ERROR: cannot parse command-line box "+pbc_cli_list;
482 0 : std::fprintf(stderr,"%s\n",msg.c_str());
483 : return 1;
484 : }
485 :
486 38 : }
487 :
488 1772 : parse("--natoms",command_line_natoms);
489 :
490 : }
491 :
492 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
493 269 : auto mf_deleter=[api](void* h_in) {
494 269 : if(h_in) {
495 269 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
496 478 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) lck=Tools::molfile_lock();
497 269 : api->close_file_read(h_in);
498 269 : }
499 : };
500 : void *h_in=NULL;
501 : std::unique_ptr<void,decltype(mf_deleter)> h_in_deleter(h_in,mf_deleter);
502 :
503 : molfile_timestep_t ts_in; // this is the structure that has the timestep
504 : // a std::vector<float> with the same scope as ts_in
505 : // it is necessary in order to store the pointer to ts_in.coords
506 : std::vector<float> ts_in_coords;
507 942 : ts_in.coords=ts_in_coords.data();
508 942 : ts_in.velocities=NULL;
509 942 : ts_in.A=-1; // we use this to check whether cell is provided or not
510 : #endif
511 :
512 :
513 :
514 942 : if(debug_dd && debug_pd) error("cannot use debug-dd and debug-pd at the same time");
515 942 : if(debug_pd || debug_dd) {
516 72 : if( !Communicator::initialized() ) error("needs mpi for debug-pd");
517 : }
518 :
519 942 : PlumedMain p; if( parseOnly ) p.activateParseOnlyMode();
520 942 : p.cmd("setRealPrecision",(int)sizeof(real));
521 : int checknatoms=-1;
522 942 : long long int step=0;
523 942 : parse("--initial-step",step);
524 :
525 942 : if(restart) p.cmd("setRestart",1);
526 :
527 942 : if(Communicator::initialized()) {
528 319 : if(multi) {
529 292 : if(intracomm.Get_rank()==0) p.cmd("GREX setMPIIntercomm",&intercomm.Get_comm());
530 340 : p.cmd("GREX setMPIIntracomm",&intracomm.Get_comm());
531 170 : p.cmd("GREX init");
532 : }
533 638 : p.cmd("setMPIComm",&intracomm.Get_comm());
534 : }
535 942 : p.cmd("setMDLengthUnits",units.getLength());
536 942 : p.cmd("setMDChargeUnits",units.getCharge());
537 942 : p.cmd("setMDMassUnits",units.getMass());
538 942 : p.cmd("setMDEngine","driver");
539 942 : p.cmd("setTimestep",timestep);
540 1879 : if( !parseOnly || full_outputfile.length()==0 ) p.cmd("setPlumedDat",plumedFile.c_str());
541 942 : p.cmd("setLog",out);
542 :
543 : int natoms;
544 942 : int lvl=0;
545 942 : int pb=1;
546 :
547 942 : if(parseOnly) {
548 5 : if(command_line_natoms<0) error("--parseOnly requires setting the number of atoms with --natoms");
549 5 : natoms=command_line_natoms;
550 : }
551 :
552 :
553 942 : FILE* fp=NULL; FILE* fp_forces=NULL; OFile fp_dforces;
554 :
555 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
556 : std::unique_ptr<FILE,decltype(deleter)> fp_forces_deleter(fp_forces,deleter);
557 :
558 11 : auto xdr_deleter=[](auto xd) { if(xd) xdrfile::xdrfile_close(xd); };
559 :
560 : xdrfile::XDRFILE* xd=NULL;
561 :
562 : std::unique_ptr<xdrfile::XDRFILE,decltype(xdr_deleter)> xd_deleter(xd,xdr_deleter);
563 :
564 942 : if(!noatoms&&!parseOnly) {
565 881 : if (trajectoryFile=="-")
566 : fp=in;
567 : else {
568 881 : if(multi) {
569 : std::string n;
570 170 : Tools::convert(intercomm.Get_rank(),n);
571 340 : std::string testfile=FileBase::appendSuffix(trajectoryFile,"."+n);
572 170 : FILE* tmp_fp=std::fopen(testfile.c_str(),"r");
573 : // no exceptions here
574 170 : if(tmp_fp) { std::fclose(tmp_fp); trajectoryFile=testfile;}
575 : }
576 881 : if(use_molfile==true) {
577 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
578 269 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
579 478 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) lck=Tools::molfile_lock();
580 269 : h_in = api->open_file_read(trajectoryFile.c_str(), trajectory_fmt.c_str(), &natoms);
581 : h_in_deleter.reset(h_in);
582 269 : if(natoms==MOLFILE_NUMATOMS_UNKNOWN) {
583 2 : if(command_line_natoms>=0) natoms=command_line_natoms;
584 0 : else error("this file format does not provide number of atoms; use --natoms on the command line");
585 : }
586 269 : ts_in_coords.resize(3*natoms);
587 269 : ts_in.coords = ts_in_coords.data();
588 : #endif
589 1491 : } else if(trajectory_fmt=="xdr-xtc" || trajectory_fmt=="xdr-trr") {
590 11 : xd=xdrfile::xdrfile_open(trajectoryFile.c_str(),"r");
591 : xd_deleter.reset(xd);
592 11 : if(!xd) {
593 0 : std::string msg="ERROR: Error opening trajectory file "+trajectoryFile;
594 0 : std::fprintf(stderr,"%s\n",msg.c_str());
595 : return 1;
596 : }
597 11 : if(trajectory_fmt=="xdr-xtc") xdrfile::read_xtc_natoms(&trajectoryFile[0],&natoms);
598 11 : if(trajectory_fmt=="xdr-trr") xdrfile::read_trr_natoms(&trajectoryFile[0],&natoms);
599 : } else {
600 601 : fp=std::fopen(trajectoryFile.c_str(),"r");
601 : fp_deleter.reset(fp);
602 601 : if(!fp) {
603 0 : std::string msg="ERROR: Error opening trajectory file "+trajectoryFile;
604 0 : std::fprintf(stderr,"%s\n",msg.c_str());
605 : return 1;
606 : }
607 : }
608 : }
609 881 : if(dumpforces.length()>0) {
610 553 : if(Communicator::initialized() && pc.Get_size()>1) {
611 : std::string n;
612 236 : Tools::convert(pc.Get_rank(),n);
613 472 : dumpforces+="."+n;
614 : }
615 553 : fp_forces=std::fopen(dumpforces.c_str(),"w");
616 : fp_forces_deleter.reset(fp_forces);
617 : }
618 881 : if(debugforces.length()>0) {
619 15 : if(Communicator::initialized() && pc.Get_size()>1) {
620 : std::string n;
621 6 : Tools::convert(pc.Get_rank(),n);
622 12 : debugforces+="."+n;
623 : }
624 15 : fp_dforces.open(debugforces);
625 : }
626 : }
627 :
628 : std::string line;
629 : std::vector<real> coordinates;
630 : std::vector<real> forces;
631 : std::vector<real> masses;
632 : std::vector<real> charges;
633 : std::vector<real> cell;
634 : std::vector<real> virial;
635 : std::vector<real> numder;
636 :
637 : // variables to test particle decomposition
638 : int pd_nlocal=0;
639 : int pd_start=0;
640 : // variables to test random decomposition (=domain decomposition)
641 : std::vector<int> dd_gatindex;
642 : std::vector<int> dd_g2l;
643 : std::vector<real> dd_masses;
644 : std::vector<real> dd_charges;
645 : std::vector<real> dd_forces;
646 : std::vector<real> dd_coordinates;
647 : int dd_nlocal=0;
648 : // random stream to choose decompositions
649 942 : Random rnd;
650 :
651 942 : if(trajectory_fmt=="dlp4") {
652 2 : if(!Tools::getline(fp,line)) error("error reading title");
653 2 : if(!Tools::getline(fp,line)) error("error reading atoms");
654 2 : std::sscanf(line.c_str(),"%d %d %d",&lvl,&pb,&natoms);
655 :
656 : }
657 : bool lstep=true;
658 258151 : while(true) {
659 259093 : if(!noatoms&&!parseOnly) {
660 58506 : if(use_molfile==true) {
661 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
662 18915 : std::unique_ptr<std::lock_guard<std::mutex>> lck;
663 37673 : if(api->is_reentrant==VMDPLUGIN_THREADUNSAFE) lck=Tools::molfile_lock();
664 : int rc;
665 18915 : rc = api->read_next_timestep(h_in, natoms, &ts_in);
666 18915 : if(rc==MOLFILE_EOF) {
667 : break;
668 : }
669 : #endif
670 71214 : } else if(trajectory_fmt=="xyz" || trajectory_fmt=="gro" || trajectory_fmt=="dlp4") {
671 39486 : if(!Tools::getline(fp,line)) break;
672 : }
673 : }
674 : bool first_step=false;
675 258227 : if(!noatoms&&!parseOnly) {
676 109102 : if(use_molfile==false && (trajectory_fmt=="xyz" || trajectory_fmt=="gro")) {
677 38869 : if(trajectory_fmt=="gro") if(!Tools::getline(fp,line)) error("premature end of trajectory file");
678 38869 : std::sscanf(line.c_str(),"%100d",&natoms);
679 : }
680 96634 : if(use_molfile==false && trajectory_fmt=="dlp4") {
681 : char xa[9];
682 : int xb,xc,xd;
683 : double t;
684 20 : std::sscanf(line.c_str(),"%8s %lld %d %d %d %lf",xa,&step,&xb,&xc,&xd,&t);
685 20 : if (lstep) {
686 2 : p.cmd("setTimestep",real(t));
687 : lstep = false;
688 : }
689 : }
690 : }
691 258227 : if(checknatoms<0 && !noatoms) {
692 886 : pd_nlocal=natoms;
693 : pd_start=0;
694 : first_step=true;
695 886 : masses.assign(natoms,std::numeric_limits<real>::quiet_NaN());
696 886 : charges.assign(natoms,std::numeric_limits<real>::quiet_NaN());
697 : //case pdb: structure
698 886 : if(pdbfile.length()>0) {
699 150980 : for(unsigned i=0; i<pdb.size(); ++i) {
700 150901 : AtomNumber an=pdb.getAtomNumbers()[i];
701 : unsigned index=an.index();
702 150901 : if( index>=unsigned(natoms) ) error("atom index in pdb exceeds the number of atoms in trajectory");
703 150901 : masses[index]=pdb.getOccupancy()[i];
704 150901 : charges[index]=pdb.getBeta()[i];
705 : }
706 : }
707 886 : if(mcfile.length()>0) {
708 12 : IFile ifile;
709 12 : ifile.open(mcfile);
710 : int index; double mass; double charge;
711 495694 : while(ifile.scanField("index",index).scanField("mass",mass).scanField("charge",charge).scanField()) {
712 247835 : masses[index]=mass;
713 247835 : charges[index]=charge;
714 : }
715 12 : }
716 257341 : } else if( checknatoms<0 && noatoms ) {
717 56 : natoms=0;
718 : }
719 258227 : if( checknatoms<0 ) {
720 942 : if(kt>=0) {
721 3 : p.cmd("setKbT",kt);
722 : }
723 942 : checknatoms=natoms;
724 942 : p.cmd("setNatoms",natoms);
725 942 : p.cmd("init");
726 : // Check if we have been asked to output the long version of the input and if there are shortcuts
727 942 : if( parseOnly && full_outputfile.length()>0 ) {
728 :
729 : // Read in the plumed input file and store what is in there
730 : std::map<std::string,std::vector<std::string> > data;
731 5 : IFile ifile; ifile.open(plumedFile); std::vector<std::string> words;
732 29 : while( Tools::getParsedLine(ifile,words) && !p.getEndPlumed() ) {
733 24 : p.readInputWords(words,false);
734 : Action* aa=p.getActionSet()[p.getActionSet().size()-1].get();
735 24 : plumed_assert(aa); // needed for following calls, see #1046
736 24 : ActionWithValue* av=aa->castToActionWithValue();
737 36 : if( av && aa->getDefaultString().length()>0 ) {
738 10 : std::vector<std::string> def; def.push_back( "defaults " + aa->getDefaultString() );
739 5 : data[ aa->getLabel() ] = def;
740 5 : }
741 24 : ActionShortcut* as=aa->castToActionShortcut();
742 24 : if( as ) {
743 18 : if( aa->getDefaultString().length()>0 ) {
744 6 : std::vector<std::string> def; def.push_back( "defaults " + aa->getDefaultString() );
745 3 : data[ as->getShortcutLabel() ] = def;
746 3 : }
747 18 : std::vector<std::string> shortcut_commands = as->getSavedInputLines(); if( shortcut_commands.size()==0 ) continue;
748 12 : if( data.find( as->getShortcutLabel() )!=data.end() ) {
749 16 : for(unsigned i=0; i<shortcut_commands.size(); ++i) data[ as->getShortcutLabel() ].push_back( shortcut_commands[i] );
750 4 : } else data[ as->getShortcutLabel() ] = shortcut_commands;
751 18 : }
752 : }
753 5 : ifile.close();
754 : // Only output the full version of the input file if there are shortcuts
755 5 : if( data.size()>0 ) {
756 5 : OFile long_file; long_file.open( full_outputfile ); long_file.printf("{\n"); bool firstpass=true;
757 17 : for(auto& x : data ) {
758 12 : if( !firstpass ) long_file.printf(" },\n");
759 12 : long_file.printf(" \"%s\" : {\n", x.first.c_str() );
760 12 : plumed_assert( x.second.size()>0 ); unsigned sstart=0;
761 12 : if( x.second[0].find("defaults")!=std::string::npos ) {
762 16 : sstart=1; long_file.printf(" \"defaults\" : \"%s\"", x.second[0].substr( 9 ).c_str() );
763 8 : if( x.second.size()>1 ) long_file.printf(",\n"); else long_file.printf("\n");
764 : }
765 12 : if( x.second.size()>sstart ) {
766 6 : long_file.printf(" \"expansion\" : \"%s", x.second[sstart].c_str() );
767 36 : for(unsigned j=sstart+1; j<x.second.size(); ++j) long_file.printf("\\n%s", x.second[j].c_str() );
768 6 : long_file.printf("\"\n");
769 : }
770 : firstpass=false;
771 : }
772 5 : long_file.printf(" }\n}\n"); long_file.close();
773 5 : }
774 5 : }
775 942 : if( valuedict_file.length()>0 ) {
776 5 : OFile valuefile; valuefile.open( valuedict_file ); valuefile.printf("{\n"); bool firsta=true;
777 149 : for(const auto & pp : p.getActionSet()) {
778 144 : if( pp.get()->getName()=="CENTER" ) continue ;
779 141 : ActionWithVirtualAtom* avv=dynamic_cast<ActionWithVirtualAtom*>(pp.get());
780 414 : if( avv || pp.get()->getName()=="GROUP" || pp.get()->getName()=="DENSITY" ) {
781 : Action* p(pp.get());
782 6 : if( firsta ) { valuefile.printf(" \"%s\" : {\n \"action\" : \"%s\"", p->getLabel().c_str(), p->getName().c_str() ); firsta=false; }
783 6 : else valuefile.printf(",\n \"%s\" : {\n \"action\" : \"%s\"", p->getLabel().c_str(), p->getName().c_str() );
784 6 : if( avv ) valuefile.printf(",\n \"%s\" : { \"type\": \"atoms\", \"description\": \"virtual atom calculated by %s action\" }", avv->getLabel().c_str(), avv->getName().c_str() );
785 3 : else valuefile.printf(",\n \"%s\" : { \"type\": \"atoms\", \"description\": \"indices of atoms specified in GROUP\" }", p->getLabel().c_str() );
786 6 : valuefile.printf("\n }"); continue;
787 6 : }
788 135 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(pp.get());
789 135 : if( av && av->getNumberOfComponents()>0 ) {
790 89 : Keywords keys; p.getKeywordsForAction( av->getName(), keys );
791 94 : if( firsta ) { valuefile.printf(" \"%s\" : {\n \"action\" : \"%s\"", av->getLabel().c_str(), keys.getDisplayName().c_str() ); firsta=false; }
792 168 : else valuefile.printf(",\n \"%s\" : {\n \"action\" : \"%s\"", av->getLabel().c_str(), keys.getDisplayName().c_str() );
793 207 : for(unsigned i=0; i<av->getNumberOfComponents(); ++i) {
794 118 : Value* myval = av->copyOutput(i); std::string compname = myval->getName(), description;
795 118 : if( av->getLabel()==compname ) {
796 150 : description = keys.getOutputComponentDescription(".#!value");
797 : } else {
798 86 : std::size_t dot=compname.find(av->getLabel() + "."); std::string cname = compname.substr(dot + av->getLabel().length() + 1);
799 86 : description = av->getOutputComponentDescription( cname, keys );
800 : }
801 118 : if( description.find("\\")!=std::string::npos ) error("found invalid backslash character in documentation for component " + compname + " in action " + av->getName() + " with label " + av->getLabel() );
802 236 : valuefile.printf(",\n \"%s\" : { \"type\": \"%s\", \"description\": \"%s\" }", myval->getName().c_str(), myval->getValueType().c_str(), description.c_str() );
803 : }
804 89 : valuefile.printf("\n }");
805 89 : }
806 135 : ActionShortcut* as=pp->castToActionShortcut();
807 135 : if( as ) {
808 41 : std::vector<std::string> cnames( as->getSavedOutputs() );
809 41 : if( cnames.size()==0 ) continue ;
810 :
811 7 : if( firsta ) { valuefile.printf(" \"shortcut_%s\" : {\n \"action\" : \"%s\"", as->getShortcutLabel().c_str(), as->getName().c_str() ); firsta=false; }
812 7 : else valuefile.printf(",\n \"shortcut_%s\" : {\n \"action\" : \"%s\"", as->getShortcutLabel().c_str(), as->getName().c_str() );
813 7 : Keywords keys; p.getKeywordsForAction( as->getName(), keys );
814 18 : for(unsigned i=0; i<cnames.size(); ++i) {
815 11 : ActionWithValue* av2=p.getActionSet().selectWithLabel<ActionWithValue*>( cnames[i] );
816 11 : if( !av2 ) plumed_merror("could not find value created by shortcut with name " + cnames[i] );
817 11 : if( av2->getNumberOfComponents()==1 ) {
818 11 : Value* myval = av2->copyOutput(0); std::string compname = myval->getName(), description;
819 18 : if( compname==as->getShortcutLabel() ) description = keys.getOutputComponentDescription(".#!value");
820 12 : else { std::size_t pp=compname.find(as->getShortcutLabel()); description = keys.getOutputComponentDescription( compname.substr(pp+as->getShortcutLabel().length()+1) ); }
821 11 : if( description.find("\\")!=std::string::npos ) error("found invalid backslash character in documentation for component " + compname + " in action " + as->getName() + " with label " + as->getLabel() );
822 22 : valuefile.printf(",\n \"%s\" : { \"type\": \"%s\", \"description\": \"%s\" }", myval->getName().c_str(), myval->getValueType().c_str(), description.c_str() );
823 : } else {
824 0 : for(unsigned j=0; j<av2->getNumberOfComponents(); ++j) {
825 0 : Value* myval = av2->copyOutput(j); std::string compname = myval->getName(), description;
826 0 : if( av2->getLabel()==compname ) {
827 0 : plumed_merror("should not be outputting description of value from action when using shortcuts");
828 : } else {
829 0 : std::size_t dot=compname.find(av2->getLabel() + "."); std::string cname = compname.substr(dot+av2->getLabel().length() + 1);
830 0 : description = av2->getOutputComponentDescription( cname, keys );
831 : }
832 0 : if( description.find("\\")!=std::string::npos ) error("found invalid backslash character in documentation for component " + compname + " in action " + av2->getName() + " with label " + av2->getLabel() );
833 0 : valuefile.printf(",\n \"%s\" : { \"type\": \"%s\", \"description\": \"%s\" }", myval->getName().c_str(), myval->getValueType().c_str(), description.c_str() );
834 : }
835 : }
836 : }
837 7 : valuefile.printf("\n }");
838 41 : }
839 : }
840 5 : valuefile.printf("\n}\n"); valuefile.close();
841 5 : }
842 942 : if(parseOnly) break;
843 : }
844 258222 : if(checknatoms!=natoms) {
845 0 : std::string stepstr; Tools::convert(step,stepstr);
846 0 : error("number of atoms in frame " + stepstr + " does not match number of atoms in first frame");
847 : }
848 :
849 258222 : coordinates.assign(3*natoms,real(0.0));
850 258222 : forces.assign(3*natoms,real(0.0));
851 258222 : cell.assign(9,real(0.0));
852 258222 : virial.assign(9,real(0.0));
853 :
854 258222 : if( first_step || rnd.U01()>0.5) {
855 130891 : if(debug_pd) {
856 152 : int npe=intracomm.Get_size();
857 152 : std::vector<int> loc(npe,0);
858 152 : std::vector<int> start(npe,0);
859 608 : for(int i=0; i<npe-1; i++) {
860 456 : int cc=(natoms*2*rnd.U01())/npe;
861 456 : if(start[i]+cc>natoms) cc=natoms-start[i];
862 456 : loc[i]=cc;
863 456 : start[i+1]=start[i]+loc[i];
864 : }
865 152 : loc[npe-1]=natoms-start[npe-1];
866 152 : intracomm.Bcast(loc,0);
867 152 : intracomm.Bcast(start,0);
868 152 : pd_nlocal=loc[intracomm.Get_rank()];
869 152 : pd_start=start[intracomm.Get_rank()];
870 152 : if(intracomm.Get_rank()==0) {
871 : std::fprintf(out,"\nDRIVER: Reassigning particle decomposition\n");
872 190 : std::fprintf(out,"DRIVER: "); for(int i=0; i<npe; i++) std::fprintf(out,"%d ",loc[i]); printf("\n");
873 190 : std::fprintf(out,"DRIVER: "); for(int i=0; i<npe; i++) std::fprintf(out,"%d ",start[i]); printf("\n");
874 : }
875 152 : p.cmd("setAtomsNlocal",pd_nlocal);
876 152 : p.cmd("setAtomsContiguous",pd_start);
877 130739 : } else if(debug_dd) {
878 956 : int npe=intracomm.Get_size();
879 956 : int rank=intracomm.Get_rank();
880 956 : dd_charges.assign(natoms,0.0);
881 956 : dd_masses.assign(natoms,0.0);
882 956 : dd_gatindex.assign(natoms,-1);
883 956 : dd_g2l.assign(natoms,-1);
884 956 : dd_coordinates.assign(3*natoms,0.0);
885 956 : dd_forces.assign(3*natoms,0.0);
886 : dd_nlocal=0;
887 53786 : for(int i=0; i<natoms; ++i) {
888 52830 : double r=rnd.U01()*npe;
889 112376 : int n; for(n=0; n<npe; n++) if(n+1>r)break;
890 52830 : plumed_assert(n<npe);
891 52830 : if(n==rank) {
892 19827 : dd_gatindex[dd_nlocal]=i;
893 19827 : dd_g2l[i]=dd_nlocal;
894 19827 : dd_charges[dd_nlocal]=charges[i];
895 19827 : dd_masses[dd_nlocal]=masses[i];
896 19827 : dd_nlocal++;
897 : }
898 : }
899 956 : if(intracomm.Get_rank()==0) {
900 : std::fprintf(out,"\nDRIVER: Reassigning domain decomposition\n");
901 : }
902 956 : p.cmd("setAtomsNlocal",dd_nlocal);
903 956 : p.cmd("setAtomsGatindex",&dd_gatindex[0], {dd_nlocal});
904 : }
905 : }
906 :
907 258222 : int plumedStopCondition=0;
908 258222 : if(!noatoms) {
909 57640 : if(use_molfile) {
910 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
911 18646 : if(pbc_cli_given==false) {
912 18622 : if(ts_in.A>0.0) { // this is negative if molfile does not provide box
913 : // info on the cell: convert using pbcset.tcl from pbctools in vmd distribution
914 18609 : real cosBC=cos(real(ts_in.alpha)*pi/180.);
915 : //double sinBC=std::sin(ts_in.alpha*pi/180.);
916 18609 : real cosAC=std::cos(real(ts_in.beta)*pi/180.);
917 18609 : real cosAB=std::cos(real(ts_in.gamma)*pi/180.);
918 18609 : real sinAB=std::sin(real(ts_in.gamma)*pi/180.);
919 18609 : real Ax=real(ts_in.A);
920 18609 : real Bx=real(ts_in.B)*cosAB;
921 18609 : real By=real(ts_in.B)*sinAB;
922 18609 : real Cx=real(ts_in.C)*cosAC;
923 18609 : real Cy=(real(ts_in.C)*real(ts_in.B)*cosBC-Cx*Bx)/By;
924 18609 : real Cz=std::sqrt(real(ts_in.C)*real(ts_in.C)-Cx*Cx-Cy*Cy);
925 18609 : cell[0]=Ax/10.; cell[1]=0.; cell[2]=0.;
926 18609 : cell[3]=Bx/10.; cell[4]=By/10.; cell[5]=0.;
927 18609 : cell[6]=Cx/10.; cell[7]=Cy/10.; cell[8]=Cz/10.;
928 : } else {
929 13 : cell[0]=0.0; cell[1]=0.0; cell[2]=0.0;
930 13 : cell[3]=0.0; cell[4]=0.0; cell[5]=0.0;
931 13 : cell[6]=0.0; cell[7]=0.0; cell[8]=0.0;
932 : }
933 : } else {
934 240 : for(unsigned i=0; i<9; i++)cell[i]=pbc_cli_box[i];
935 : }
936 : // info on coords
937 : // the order is xyzxyz...
938 14247430 : for(int i=0; i<3*natoms; i++) {
939 14228784 : coordinates[i]=real(ts_in.coords[i])/real(10.); //convert to nm
940 : //cerr<<"COOR "<<coordinates[i]<<endl;
941 : }
942 : #endif
943 77958 : } else if(trajectory_fmt=="xdr-xtc" || trajectory_fmt=="xdr-trr") {
944 : int localstep;
945 : float time;
946 : xdrfile::matrix box;
947 : // here we cannot use a std::vector<rvec> since it does not compile.
948 : // we thus use a std::unique_ptr<rvec[]>
949 : auto pos=Tools::make_unique<xdrfile::rvec[]>(natoms);
950 : float prec,lambda;
951 : int ret=xdrfile::exdrOK;
952 105 : if(trajectory_fmt=="xdr-xtc") ret=xdrfile::read_xtc(xd,natoms,&localstep,&time,box,pos.get(),&prec);
953 105 : if(trajectory_fmt=="xdr-trr") ret=xdrfile::read_trr(xd,natoms,&localstep,&time,&lambda,box,pos.get(),NULL,NULL);
954 105 : if(stride==0) step=localstep;
955 105 : if(ret==xdrfile::exdrENDOFFILE) break;
956 103 : if(ret!=xdrfile::exdrOK) break;
957 1222 : for(unsigned i=0; i<3; i++) for(unsigned j=0; j<3; j++) cell[3*i+j]=box[i][j];
958 26486 : for(int i=0; i<natoms; i++) for(unsigned j=0; j<3; j++)
959 19794 : coordinates[3*i+j]=real(pos[i][j]);
960 : } else {
961 38889 : if(trajectory_fmt=="xyz") {
962 26526 : if(!Tools::getline(fp,line)) error("premature end of trajectory file");
963 :
964 26526 : std::vector<double> celld(9,0.0);
965 26526 : if(pbc_cli_given==false) {
966 : std::vector<std::string> words;
967 52316 : words=Tools::getWords(line);
968 26158 : if(words.size()==3) {
969 25581 : Tools::convert(words[0],celld[0]);
970 25581 : Tools::convert(words[1],celld[4]);
971 25581 : Tools::convert(words[2],celld[8]);
972 577 : } else if(words.size()==9) {
973 577 : Tools::convert(words[0],celld[0]);
974 577 : Tools::convert(words[1],celld[1]);
975 577 : Tools::convert(words[2],celld[2]);
976 577 : Tools::convert(words[3],celld[3]);
977 577 : Tools::convert(words[4],celld[4]);
978 577 : Tools::convert(words[5],celld[5]);
979 577 : Tools::convert(words[6],celld[6]);
980 577 : Tools::convert(words[7],celld[7]);
981 577 : Tools::convert(words[8],celld[8]);
982 0 : } else error("needed box in second line of xyz file");
983 26158 : } else { // from command line
984 368 : celld=pbc_cli_box;
985 : }
986 265260 : for(unsigned i=0; i<9; i++)cell[i]=real(celld[i]);
987 : }
988 38889 : if(trajectory_fmt=="dlp4") {
989 20 : std::vector<double> celld(9,0.0);
990 20 : if(pbc_cli_given==false) {
991 20 : if(!Tools::getline(fp,line)) error("error reading vector a of cell");
992 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[0],&celld[1],&celld[2]);
993 20 : if(!Tools::getline(fp,line)) error("error reading vector b of cell");
994 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[3],&celld[4],&celld[5]);
995 20 : if(!Tools::getline(fp,line)) error("error reading vector c of cell");
996 20 : std::sscanf(line.c_str(),"%lf %lf %lf",&celld[6],&celld[7],&celld[8]);
997 : } else {
998 0 : celld=pbc_cli_box;
999 : }
1000 200 : for(auto i=0; i<9; i++)cell[i]=real(celld[i])*0.1;
1001 : }
1002 : int ddist=0;
1003 : // Read coordinates
1004 2808703 : for(int i=0; i<natoms; i++) {
1005 2769814 : bool ok=Tools::getline(fp,line);
1006 2769814 : if(!ok) error("premature end of trajectory file");
1007 : double cc[3];
1008 2769814 : if(trajectory_fmt=="xyz") {
1009 : char dummy[1000];
1010 1882537 : int ret=std::sscanf(line.c_str(),"%999s %100lf %100lf %100lf",dummy,&cc[0],&cc[1],&cc[2]);
1011 1882537 : if(ret!=4) error("cannot read line"+line);
1012 887277 : } else if(trajectory_fmt=="gro") {
1013 : // do the gromacs way
1014 886637 : if(!i) {
1015 : //
1016 : // calculate the distance between dots (as in gromacs gmxlib/confio.c, routine get_w_conf )
1017 : //
1018 : const char *p1, *p2, *p3;
1019 : p1 = std::strchr(line.c_str(), '.');
1020 12343 : if (p1 == NULL) error("seems there are no coordinates in the gro file");
1021 12343 : p2 = std::strchr(&p1[1], '.');
1022 12343 : if (p2 == NULL) error("seems there is only one coordinates in the gro file");
1023 12343 : ddist = p2 - p1;
1024 12343 : p3 = std::strchr(&p2[1], '.');
1025 12343 : if (p3 == NULL)error("seems there are only two coordinates in the gro file");
1026 12343 : if (p3 - p2 != ddist)error("not uniform spacing in fields in the gro file");
1027 : }
1028 886637 : Tools::convert(line.substr(20,ddist),cc[0]);
1029 886637 : Tools::convert(line.substr(20+ddist,ddist),cc[1]);
1030 1773274 : Tools::convert(line.substr(20+ddist+ddist,ddist),cc[2]);
1031 640 : } else if(trajectory_fmt=="dlp4") {
1032 : char dummy[9];
1033 : int idummy;
1034 : double m,c;
1035 640 : std::sscanf(line.c_str(),"%8s %d %lf %lf",dummy,&idummy,&m,&c);
1036 640 : masses[i]=real(m);
1037 640 : charges[i]=real(c);
1038 640 : if(!Tools::getline(fp,line)) error("error reading coordinates");
1039 640 : std::sscanf(line.c_str(),"%lf %lf %lf",&cc[0],&cc[1],&cc[2]);
1040 640 : cc[0]*=0.1;
1041 640 : cc[1]*=0.1;
1042 640 : cc[2]*=0.1;
1043 640 : if(lvl>0) {
1044 640 : if(!Tools::getline(fp,line)) error("error skipping velocities");
1045 : }
1046 640 : if(lvl>1) {
1047 640 : if(!Tools::getline(fp,line)) error("error skipping forces");
1048 : }
1049 0 : } else plumed_error();
1050 2769814 : if(!debug_pd || ( i>=pd_start && i<pd_start+pd_nlocal) ) {
1051 2750374 : coordinates[3*i]=real(cc[0]);
1052 2750374 : coordinates[3*i+1]=real(cc[1]);
1053 2750374 : coordinates[3*i+2]=real(cc[2]);
1054 : }
1055 : }
1056 38889 : if(trajectory_fmt=="gro") {
1057 12343 : if(!Tools::getline(fp,line)) error("premature end of trajectory file");
1058 12343 : std::vector<std::string> words=Tools::getWords(line);
1059 12343 : if(words.size()<3) error("cannot understand box format");
1060 12343 : Tools::convert(words[0],cell[0]);
1061 12343 : Tools::convert(words[1],cell[4]);
1062 12343 : Tools::convert(words[2],cell[8]);
1063 12343 : if(words.size()>3) Tools::convert(words[3],cell[1]);
1064 12343 : if(words.size()>4) Tools::convert(words[4],cell[2]);
1065 12343 : if(words.size()>5) Tools::convert(words[5],cell[3]);
1066 12343 : if(words.size()>6) Tools::convert(words[6],cell[5]);
1067 12343 : if(words.size()>7) Tools::convert(words[7],cell[6]);
1068 12343 : if(words.size()>8) Tools::convert(words[8],cell[7]);
1069 12343 : }
1070 :
1071 : }
1072 :
1073 57629 : p.cmd("setStepLongLong",step);
1074 57629 : p.cmd("setStopFlag",&plumedStopCondition);
1075 :
1076 57629 : if(debug_dd) {
1077 38944 : for(int i=0; i<dd_nlocal; ++i) {
1078 37156 : int kk=dd_gatindex[i];
1079 37156 : dd_coordinates[3*i+0]=coordinates[3*kk+0];
1080 37156 : dd_coordinates[3*i+1]=coordinates[3*kk+1];
1081 37156 : dd_coordinates[3*i+2]=coordinates[3*kk+2];
1082 : }
1083 1788 : p.cmd("setForces",&dd_forces[0], {dd_nlocal,3});
1084 1788 : p.cmd("setPositions",&dd_coordinates[0], {dd_nlocal,3});
1085 1788 : p.cmd("setMasses",&dd_masses[0], {dd_nlocal});
1086 1788 : p.cmd("setCharges",&dd_charges[0], {dd_nlocal});
1087 : } else {
1088 : // this is required to avoid troubles when the last domain
1089 : // contains zero atoms
1090 : // Basically, for empty domains we pass null pointers
1091 : #define fix_pd(xx) (pd_nlocal!=0?&xx:NULL)
1092 111682 : p.cmd("setForces",fix_pd(forces[3*pd_start]), {pd_nlocal,3});
1093 111682 : p.cmd("setPositions",fix_pd(coordinates[3*pd_start]), {pd_nlocal,3});
1094 111682 : p.cmd("setMasses",fix_pd(masses[pd_start]), {pd_nlocal});
1095 111682 : p.cmd("setCharges",fix_pd(charges[pd_start]), {pd_nlocal});
1096 : }
1097 57629 : p.cmd("setBox",cell.data(), {3,3});
1098 57629 : p.cmd("setVirial",virial.data(), {3,3});
1099 : } else {
1100 200582 : p.cmd("setStepLongLong",step);
1101 200582 : p.cmd("setStopFlag",&plumedStopCondition);
1102 : }
1103 258211 : p.cmd("calc");
1104 258211 : if(debugforces.length()>0) {
1105 96 : virial.assign(9,real(0.0));
1106 96 : forces.assign(3*natoms,real(0.0));
1107 96 : p.cmd("prepareCalc");
1108 96 : p.cmd("performCalcNoUpdate");
1109 : }
1110 :
1111 : // this is necessary as only processor zero is adding to the virial:
1112 258211 : intracomm.Bcast(virial,0);
1113 258211 : if(debug_pd) intracomm.Sum(forces);
1114 258211 : if(debug_dd) {
1115 38944 : for(int i=0; i<dd_nlocal; i++) {
1116 37156 : forces[3*dd_gatindex[i]+0]=dd_forces[3*i+0];
1117 37156 : forces[3*dd_gatindex[i]+1]=dd_forces[3*i+1];
1118 37156 : forces[3*dd_gatindex[i]+2]=dd_forces[3*i+2];
1119 : }
1120 1788 : dd_forces.assign(3*natoms,0.0);
1121 1788 : intracomm.Sum(forces);
1122 : }
1123 258211 : if(debug_grex &&step%grex_stride==0) {
1124 114 : p.cmd("GREX savePositions");
1125 114 : if(intracomm.Get_rank()>0) {
1126 57 : p.cmd("GREX prepare");
1127 : } else {
1128 57 : int r=intercomm.Get_rank();
1129 57 : int n=intercomm.Get_size();
1130 57 : int partner=r+(2*((r+step/grex_stride)%2))-1;
1131 : if(partner<0)partner=0;
1132 57 : if(partner>=n) partner=n-1;
1133 57 : p.cmd("GREX setPartner",partner);
1134 57 : p.cmd("GREX calculate");
1135 57 : p.cmd("GREX shareAllDeltaBias");
1136 228 : for(int i=0; i<n; i++) {
1137 171 : std::string s; Tools::convert(i,s);
1138 342 : real a=std::numeric_limits<real>::quiet_NaN(); s="GREX getDeltaBias "+s; p.cmd(s,&a);
1139 171 : if(grex_log) std::fprintf(grex_log," %f",a);
1140 : }
1141 57 : if(grex_log) std::fprintf(grex_log,"\n");
1142 : }
1143 : }
1144 :
1145 :
1146 258211 : if(fp_forces) {
1147 35897 : std::fprintf(fp_forces,"%d\n",natoms);
1148 71794 : std::string fmtv=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
1149 71794 : std::string fmt=dumpforcesFmt+" "+dumpforcesFmt+" "+dumpforcesFmt+"\n";
1150 35897 : if(dumpfullvirial) {
1151 351 : 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]);
1152 : } else {
1153 35546 : std::fprintf(fp_forces,fmt.c_str(),virial[0],virial[4],virial[8]);
1154 : }
1155 35897 : fmt="X "+fmt;
1156 4000472 : for(int i=0; i<natoms; i++)
1157 3964575 : std::fprintf(fp_forces,fmt.c_str(),forces[3*i],forces[3*i+1],forces[3*i+2]);
1158 : }
1159 258211 : if(debugforces.length()>0) {
1160 : // Now call the routine to work out the derivatives numerically
1161 96 : numder.assign(3*natoms+9,real(0.0)); real base=0;
1162 96 : p.cmd("getBias",&base);
1163 96 : if( std::abs(base)<epsilon ) printf("WARNING: bias for configuration appears to be zero so debugging forces is trivial");
1164 96 : evaluateNumericalDerivatives( step, p, coordinates, masses, charges, cell, base, numder );
1165 :
1166 : // And output everything to a file
1167 96 : fp_dforces.fmtField(" " + dumpforcesFmt);
1168 6033 : for(int i=0; i<3*natoms; ++i) {
1169 5937 : fp_dforces.printField("parameter",(int)i);
1170 11874 : fp_dforces.printField("analytical",forces[i]);
1171 5937 : fp_dforces.printField("numerical",-numder[i]);
1172 5937 : fp_dforces.printField();
1173 : }
1174 : // And print the virial
1175 960 : for(int i=0; i<9; ++i) {
1176 864 : fp_dforces.printField("parameter",3*natoms+i );
1177 864 : fp_dforces.printField("analytical",virial[i] );
1178 864 : fp_dforces.printField("numerical",-numder[3*natoms+i]);
1179 864 : fp_dforces.printField();
1180 : }
1181 : }
1182 :
1183 258211 : if(plumedStopCondition) break;
1184 :
1185 258151 : step+=stride;
1186 : }
1187 942 : if(!parseOnly) p.cmd("runFinalJobs");
1188 :
1189 : return 0;
1190 3768 : }
1191 :
1192 : template<typename real>
1193 96 : void Driver<real>::evaluateNumericalDerivatives( const long long int& step, PlumedMain& p, const std::vector<real>& coordinates,
1194 : const std::vector<real>& masses, const std::vector<real>& charges,
1195 : std::vector<real>& cell, const double& base, std::vector<real>& numder ) {
1196 :
1197 96 : int natoms = coordinates.size() / 3; real delta = std::sqrt(epsilon);
1198 96 : std::vector<Vector> pos(natoms); real bias=0;
1199 96 : std::vector<real> fake_forces( 3*natoms ), fake_virial(9);
1200 2075 : for(int i=0; i<natoms; ++i) {
1201 7916 : for(unsigned j=0; j<3; ++j) pos[i][j]=coordinates[3*i+j];
1202 : }
1203 :
1204 2075 : for(int i=0; i<natoms; ++i) {
1205 7916 : for(unsigned j=0; j<3; ++j) {
1206 5937 : pos[i][j]=pos[i][j]+delta;
1207 5937 : p.cmd("setStepLongLong",step);
1208 5937 : p.cmd("setPositions",&pos[0][0], {natoms,3});
1209 5937 : p.cmd("setForces",&fake_forces[0], {natoms,3});
1210 5937 : p.cmd("setMasses",&masses[0], {natoms});
1211 5937 : p.cmd("setCharges",&charges[0], {natoms});
1212 5937 : p.cmd("setBox",&cell[0], {3,3});
1213 5937 : p.cmd("setVirial",&fake_virial[0], {3,3});
1214 5937 : p.cmd("prepareCalc");
1215 5937 : p.cmd("performCalcNoUpdate");
1216 5937 : p.cmd("getBias",&bias);
1217 5937 : pos[i][j]=coordinates[3*i+j];
1218 5937 : numder[3*i+j] = (bias - base) / delta;
1219 : }
1220 : }
1221 :
1222 : // Create the cell
1223 96 : Tensor box( cell[0], cell[1], cell[2], cell[3], cell[4], cell[5], cell[6], cell[7], cell[8] );
1224 : // And the virial
1225 96 : Pbc pbc; pbc.setBox( box ); Tensor nvirial;
1226 1248 : for(unsigned i=0; i<3; i++) for(unsigned k=0; k<3; k++) {
1227 864 : double arg0=box(i,k);
1228 18675 : for(int j=0; j<natoms; ++j) pos[j]=pbc.realToScaled( pos[j] );
1229 864 : cell[3*i+k]=box(i,k)=box(i,k)+delta; pbc.setBox(box);
1230 18675 : for(int j=0; j<natoms; j++) pos[j]=pbc.scaledToReal( pos[j] );
1231 864 : p.cmd("setStepLongLong",step);
1232 864 : p.cmd("setPositions",&pos[0][0], {natoms,3});
1233 864 : p.cmd("setForces",&fake_forces[0], {natoms,3});
1234 864 : p.cmd("setMasses",&masses[0], {natoms});
1235 864 : p.cmd("setCharges",&charges[0], {natoms});
1236 864 : p.cmd("setBox",&cell[0], {3,3});
1237 864 : p.cmd("setVirial",&fake_virial[0], {3,3});
1238 864 : p.cmd("prepareCalc");
1239 864 : p.cmd("performCalcNoUpdate");
1240 864 : p.cmd("getBias",&bias);
1241 864 : cell[3*i+k]=box(i,k)=arg0; pbc.setBox(box);
1242 72108 : for(int j=0; j<natoms; j++) for(unsigned n=0; n<3; ++n) pos[j][n]=coordinates[3*j+n];
1243 864 : nvirial(i,k) = ( bias - base ) / delta;
1244 : }
1245 96 : nvirial=-matmul(box.transpose(),nvirial);
1246 1248 : for(unsigned i=0; i<3; i++) for(unsigned k=0; k<3; k++) numder[3*natoms+3*i+k] = nvirial(i,k);
1247 :
1248 96 : }
1249 :
1250 : }
1251 : }
|
