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