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