Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2018-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 "core/ActionAtomistic.h"
23 : #include "core/ActionWithValue.h"
24 : #include "core/ActionPilot.h"
25 : #include "core/ActionRegister.h"
26 : #include "tools/Tools.h"
27 : #include "tools/PlumedHandle.h"
28 : #include "tools/Communicator.h"
29 : #include "core/PlumedMain.h"
30 : #include <cstring>
31 : #ifdef __PLUMED_HAS_DLOPEN
32 : #include <dlfcn.h>
33 : #endif
34 :
35 : namespace PLMD {
36 : namespace generic {
37 :
38 : //+PLUMEDOC GENERIC PLUMED
39 : /*
40 : Embed a separate PLUMED instance.
41 :
42 : This command can be used to embed a separate PLUMED instance.
43 : Only required atoms will be passed to that instance, using an interface
44 : that is similar to the one used when calling PLUMED from the NAMD engine.
45 :
46 : Notice that the two instances are running in the same UNIX process, so that they cannot be perfectly isolated.
47 : However, most of the features are expected to work correctly.
48 :
49 : Notes:
50 : - The [LOAD](LOAD.md) action will not work correctly since registers will be shared among the two instances.
51 : In particular, the loaded actions will be visible to both guest and host irrespective of where they are loaded from.
52 : This can be fixed and will probably be fixed in a later version.
53 : - `CHDIR` is not thread safe.
54 : However, in most implementations there will be a single process running PLUMED, with perhaps multiple OpenMP threads
55 : spawn in order to parallelize the calculation of individual variables. So, this is likely not a problem.
56 : - MPI is working correctly. However, notice that the guest PLUMED will always run with a single process.
57 : Multiple replicas should be handled correctly.
58 :
59 : As an advanced feature, one can use the option `KERNEL` to select the version of the guest PLUMED instance.
60 : In particular, an empty `KERNEL` (default) implies that the guest PLUMED instance is the same as the host one
61 : (no library is loaded).
62 : On the other hand, `KERNEL=/path/to/libplumedKernel.so` will allow specifying a library to be loaded for the
63 : guest instance.
64 : In addition to those mentioned above, this feature has limitations mostly related to
65 : clashes in the symbols defined in the different instances of the PLUMED library.
66 : Clashes might be due by synonymous symbols from the PLUMED library or synonymous symbols
67 : from libraries linked to PLUMED, and would differ depending on the operating system you
68 : are using:
69 : - On OSX:
70 : - Symbols from the PLUMED library would clash. The only way to avoid is to load PLUMED dynamically:
71 : - If you are using PLUMED with an MD code, it should be patched with `--runtime`.
72 : - If you are using PLUMED driver, you should launch the `plumed-runtime` executable (contained in the
73 : `prefix/lib/plumed/` directory) and export `PLUMED_KERNEL` equal to the path of the host kernel library
74 : (as usual in runtime loading). Notice that as of PLUMED 2.10 we load kernels with RTLD_LOCAL by default.
75 : - Symbols from dependent libraries should not clash since, as of version 2.5, we are using
76 : two-level namespace. Problems when loading previous versions should anyway be solved using runtime
77 : mode as explained above.
78 : - On Linux, any `KERNEL` should in principle work correctly. To achieve namespace separation we are loading
79 : the guest kernel with `RTLD_DEEPBIND`. However, this might create difficult to track problems in other linked libraries.
80 : - On Unix systems where `RTLD_DEEPBIND` is not available kernels will not load correctly.
81 : - In general, there might be unexpected crashes. Particularly difficult are situations where different
82 : kernels were compiled with different libraries.
83 :
84 :
85 : The following things need to be done to improve this action:
86 :
87 : - Add support for multiple time stepping (`STRIDE` different from 1).
88 : - Add the possibility to import CVs calculated in the host PLUMED instance into the guest PLUMED instance.
89 : Will be possible after \issue{83} will be closed.
90 : - Add the possibility to export CVs calculated in the guest PLUMED instance into the host PLUMED instance.
91 : Could be implemented using the `DataFetchingObject` class.
92 :
93 : ## Examples
94 :
95 : Here an example plumed file:
96 :
97 : ```plumed
98 : # plumed.dat
99 : p: PLUMED FILE=plumed2.dat
100 : PRINT ARG=p.bias FILE=COLVAR
101 : ```
102 :
103 : `plumed2.dat` can be an arbitrary plumed input file, for instance
104 :
105 : ```plumed
106 : #SETTINGS FILENAME=plumed2.dat
107 : # plumed2.dat
108 : d: DISTANCE ATOMS=1,10
109 : RESTRAINT ARG=d KAPPA=10 AT=2
110 : ```
111 :
112 : Now a more useful example.
113 : Imagine that you ran simulations using two different PLUMED input files.
114 : The files are long and complex and there are some clashes in the name of the variables (that is: same names
115 : are used in both files, same files are written, etc). In addition, files might have been written using different units (see [UNITS](UNITS.md)).
116 : If you want to run a single simulation with a bias potential
117 : that is the sum of the two bias potentials, you can:
118 : - Place the two input files, as well as all the files required by plumed, in separate directories `directory1` and `directory2`.
119 : - Run with the following input file in the parent directory:
120 :
121 : ```plumed
122 : # plumed.dat
123 : PLUMED FILE=plumed.dat CHDIR=directory1
124 : PLUMED FILE=plumed.dat CHDIR=directory2
125 : ```
126 :
127 : */
128 : //+ENDPLUMEDOC
129 :
130 : class Plumed:
131 : public ActionAtomistic,
132 : public ActionWithValue,
133 : public ActionPilot {
134 : /// True on root processor
135 : const bool root;
136 : /// Separate directory.
137 : const std::string directory;
138 : /// Interface to underlying plumed object.
139 : PlumedHandle p;
140 : /// API number.
141 : const int API;
142 : /// Self communicator
143 : Communicator comm_self;
144 : /// Intercommunicator
145 : Communicator intercomm;
146 : /// Detect first usage.
147 : bool first=true;
148 : /// Stop flag, used to stop e.g. in committor analysis
149 : int stop=0;
150 : /// Index of requested atoms.
151 : std::vector<int> index;
152 : /// Masses of requested atoms.
153 : std::vector<double> masses;
154 : /// Charges of requested atoms.
155 : std::vector<double> charges;
156 : /// Forces on requested atoms.
157 : std::vector<double> forces;
158 : /// Requested positions.
159 : std::vector<double> positions;
160 : /// Applied virial.
161 : Tensor virial;
162 : public:
163 : /// Constructor.
164 : explicit Plumed(const ActionOptions&);
165 : /// Documentation.
166 : static void registerKeywords( Keywords& keys );
167 : void prepare() override;
168 : void calculate() override;
169 : void apply() override;
170 : void update() override;
171 0 : unsigned getNumberOfDerivatives() override {
172 0 : return 0;
173 : }
174 92 : bool actionHasForces() override {
175 92 : return true;
176 : }
177 : };
178 :
179 : PLUMED_REGISTER_ACTION(Plumed,"PLUMED")
180 :
181 17 : void Plumed::registerKeywords( Keywords& keys ) {
182 17 : Action::registerKeywords( keys );
183 17 : ActionPilot::registerKeywords( keys );
184 17 : ActionAtomistic::registerKeywords( keys );
185 17 : keys.add("compulsory","STRIDE","1","stride different from 1 are not supported yet");
186 17 : keys.add("optional","FILE","input file for the guest PLUMED instance");
187 17 : keys.add("optional","KERNEL","kernel to be used for the guest PLUMED instance (USE WITH CAUTION!)");
188 17 : keys.add("optional","LOG","log file for the guest PLUMED instance. By default the host log is used");
189 17 : keys.add("optional","CHDIR","run guest in a separate directory");
190 17 : keys.addFlag("NOREPLICAS",false,"run multiple replicas as isolated ones, without letting them know that the host has multiple replicas");
191 34 : keys.addOutputComponent("bias","default","scalar","the instantaneous value of the bias potential");
192 17 : }
193 :
194 13 : Plumed::Plumed(const ActionOptions&ao):
195 : Action(ao),
196 : ActionAtomistic(ao),
197 : ActionWithValue(ao),
198 : ActionPilot(ao),
199 26 : root(comm.Get_rank()==0),
200 27 : directory([&]() {
201 : std::string dir;
202 26 : parse("CHDIR",dir);
203 13 : if(dir.length()>0) {
204 1 : log<<" running on separate directory "<<dir<<"\n";
205 : }
206 13 : return dir;
207 : }()),
208 27 : p([&]() {
209 : std::string kernel;
210 26 : parse("KERNEL",kernel);
211 13 : if(kernel.length()==0) {
212 12 : log<<" using the current kernel\n";
213 12 : return PlumedHandle();
214 : } else {
215 1 : log<<" using the kernel "<<kernel<<"\n";
216 1 : return PlumedHandle::dlopen(kernel.c_str());
217 : }
218 : }()),
219 40 : API([&]() {
220 13 : int api=0;
221 13 : p.cmd("getApiVersion",&api);
222 13 : log<<" reported API version is "<<api<<"\n";
223 : // note: this is required in order to have cmd performCalcNoUpdate and cmd update
224 : // as a matter of fact, any version <2.5 will not even load due to namespace pollution
225 13 : plumed_assert(api>3) << "API>3 is required for the PLUMED action to work correctly\n";
226 13 : return api;
227 26 : }()) {
228 13 : Tools::DirectoryChanger directoryChanger(directory.c_str());
229 :
230 : bool noreplicas;
231 12 : parseFlag("NOREPLICAS",noreplicas);
232 12 : int nreps=getNumberOfReplicas();
233 12 : if(nreps>1) {
234 6 : if(noreplicas) {
235 0 : log<<" running replicas as independent (no suffix used)\n";
236 : } else {
237 6 : log<<" running replicas as standard multi replic (with suffix)\n";
238 6 : if(root) {
239 6 : intercomm.Set_comm(&multi_sim_comm.Get_comm());
240 6 : p.cmd("GREX setMPIIntercomm",&intercomm.Get_comm());
241 6 : p.cmd("GREX setMPIIntracomm",&comm_self.Get_comm());
242 3 : p.cmd("GREX init");
243 : }
244 : }
245 : } else {
246 6 : if(noreplicas) {
247 0 : log<<" WARNING: flag NOREPLICAS ignored since we are running without replicas\n";
248 : }
249 : }
250 :
251 12 : int natoms=getTotAtoms();
252 :
253 12 : plumed_assert(getStride()==1) << "currently only supports STRIDE=1";
254 :
255 12 : double dt=getTimeStep();
256 :
257 : std::string file;
258 24 : parse("FILE",file);
259 12 : if(file.length()>0) {
260 12 : log<<" with input file "<<file<<"\n";
261 : } else {
262 0 : plumed_error() << "you must provide an input file\n";
263 : }
264 :
265 : bool inherited_logfile=false;
266 : std::string logfile;
267 24 : parse("LOG",logfile);
268 12 : if(logfile.length()>0) {
269 0 : log<<" with log file "<<logfile<<"\n";
270 0 : if(root) {
271 0 : p.cmd("setLogFile",logfile.c_str());
272 : }
273 12 : } else if(log.getFILE()) {
274 12 : log<<" with inherited log file\n";
275 12 : if(root) {
276 6 : p.cmd("setLog",log.getFILE());
277 : }
278 : inherited_logfile=true;
279 : } else {
280 0 : log<<" with log on stdout\n";
281 0 : if(root) {
282 0 : p.cmd("setLog",stdout);
283 : }
284 : }
285 :
286 12 : checkRead();
287 :
288 12 : if(root) {
289 6 : p.cmd("setMDEngine","plumed");
290 : }
291 :
292 12 : double engunits=getUnits().getEnergy();
293 12 : if(root) {
294 6 : p.cmd("setMDEnergyUnits",&engunits);
295 : }
296 :
297 12 : double lenunits=getUnits().getLength();
298 12 : if(root) {
299 6 : p.cmd("setMDLengthUnits",&lenunits);
300 : }
301 :
302 12 : double timunits=getUnits().getTime();
303 12 : if(root) {
304 6 : p.cmd("setMDTimeUnits",&timunits);
305 : }
306 :
307 12 : double chaunits=getUnits().getCharge();
308 12 : if(root) {
309 6 : p.cmd("setMDChargeUnits",&chaunits);
310 : }
311 12 : double masunits=getUnits().getMass();
312 12 : if(root) {
313 6 : p.cmd("setMDMassUnits",&masunits);
314 : }
315 :
316 12 : double kbt=getkBT();
317 12 : if(root) {
318 6 : p.cmd("setKbT",&kbt);
319 : }
320 :
321 12 : int res=0;
322 12 : if(getRestart()) {
323 0 : res=1;
324 : }
325 12 : if(root) {
326 6 : p.cmd("setRestart",&res);
327 : }
328 :
329 12 : if(root) {
330 6 : p.cmd("setNatoms",&natoms);
331 : }
332 12 : if(root) {
333 6 : p.cmd("setTimestep",&dt);
334 : }
335 12 : if(root) {
336 6 : p.cmd("setPlumedDat",file.c_str());
337 : }
338 :
339 24 : addComponentWithDerivatives("bias");
340 12 : componentIsNotPeriodic("bias");
341 :
342 12 : if(inherited_logfile) {
343 12 : log<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
344 : }
345 12 : if(root) {
346 6 : p.cmd("init");
347 : }
348 12 : if(inherited_logfile) {
349 12 : log<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
350 : }
351 23 : }
352 :
353 164 : void Plumed::prepare() {
354 164 : Tools::DirectoryChanger directoryChanger(directory.c_str());
355 164 : int step=getStep();
356 164 : if(root) {
357 89 : p.cmd("setStep",&step);
358 : }
359 164 : if(root) {
360 89 : p.cmd("prepareDependencies");
361 : }
362 164 : int ene=0;
363 164 : if(root) {
364 89 : p.cmd("isEnergyNeeded",&ene);
365 : }
366 164 : if(ene) {
367 0 : plumed_error()<<"It is not currently possible to use ENERGY in a guest PLUMED";
368 : }
369 164 : int n=0;
370 164 : if(root) {
371 89 : p.cmd("createFullList",&n);
372 : }
373 164 : const int *pointer=nullptr;
374 164 : if(root) {
375 89 : p.cmd("getFullList",&pointer);
376 : }
377 164 : bool redo=(index.size()!=static_cast<unsigned>(n));
378 164 : if(first) {
379 : redo=true;
380 : }
381 164 : first=false;
382 164 : if(root && !redo)
383 4245 : for(int i=0; i<n; i++)
384 4178 : if(index[i]!=pointer[i]) {
385 : redo=true;
386 : break;
387 : };
388 164 : if(root && redo) {
389 22 : index.resize(n);
390 22 : masses.resize(n);
391 22 : forces.resize(3*n);
392 22 : positions.resize(3*n);
393 22 : charges.resize(n);
394 1043 : for(int i=0; i<n; i++) {
395 1021 : index[i]=pointer[i];
396 : };
397 22 : p.cmd("setAtomsNlocal",&n);
398 22 : p.cmd("setAtomsGatindex",index.data(), {index.size()});
399 : }
400 164 : if(root) {
401 89 : p.cmd("clearFullList");
402 : }
403 164 : int tmp=0;
404 164 : if(root && redo) {
405 22 : tmp=1;
406 : }
407 164 : comm.Bcast(tmp,0);
408 164 : if(tmp) {
409 34 : int s=index.size();
410 34 : comm.Bcast(s,0);
411 34 : if(!root) {
412 12 : index.resize(s);
413 : }
414 34 : comm.Bcast(index,0);
415 34 : std::vector<AtomNumber> numbers;
416 34 : numbers.reserve(index.size());
417 2138 : for(auto i : index) {
418 2104 : numbers.emplace_back(AtomNumber::index(i));
419 : }
420 34 : requestAtoms(numbers);
421 34 : }
422 164 : }
423 :
424 134 : void Plumed::calculate() {
425 134 : Tools::DirectoryChanger directoryChanger(directory.c_str());
426 134 : if(root) {
427 74 : p.cmd("setStopFlag",&stop);
428 : }
429 134 : Tensor box=getPbc().getBox();
430 134 : if(root)
431 74 : p.cmd("setBox",&box[0][0], {3,3});
432 :
433 : virial.zero();
434 : std::fill(forces.begin(),forces.end(),0.0);
435 : std::copy(getMasses().begin(),
436 : getMasses().begin()+masses.size(),
437 : masses.begin());
438 134 : std::copy(getCharges().begin(),
439 134 : getCharges().begin()+charges.size(),
440 : charges.begin());
441 :
442 134 : if(root)
443 74 : p.cmd("setMasses",masses.data(), {masses.size()});
444 134 : if(root)
445 74 : p.cmd("setCharges",charges.data(), {charges.size()});
446 134 : if(root)
447 74 : p.cmd("setPositions",positions.data(), {masses.size(),3});
448 134 : if(root) {
449 74 : p.cmd("setForces",forces.data(),forces.size());
450 : }
451 134 : if(root)
452 74 : p.cmd("setVirial",virial.data(), {3,3});
453 :
454 :
455 134 : if(root)
456 4178 : for(unsigned i=0; i<getNumberOfAtoms(); i++) {
457 4104 : positions[3*i+0]=getPosition(i)[0];
458 4104 : positions[3*i+1]=getPosition(i)[1];
459 4104 : positions[3*i+2]=getPosition(i)[2];
460 : }
461 :
462 134 : if(root) {
463 74 : p.cmd("shareData");
464 : }
465 134 : if(root) {
466 74 : p.cmd("performCalcNoUpdate");
467 : }
468 :
469 134 : int s=forces.size();
470 134 : comm.Bcast(s,0);
471 134 : if(!root) {
472 60 : forces.resize(s);
473 : }
474 134 : comm.Bcast(forces,0);
475 134 : comm.Bcast(virial,0);
476 :
477 134 : double bias=0.0;
478 134 : if(root) {
479 74 : p.cmd("getBias",&bias);
480 : }
481 134 : comm.Bcast(bias,0);
482 134 : getPntrToComponent("bias")->set(bias);
483 134 : }
484 :
485 104 : void Plumed::apply() {
486 104 : Tools::DirectoryChanger directoryChanger(directory.c_str());
487 104 : std::vector<double> fforces( forces.size() + 9, 0 );
488 19580 : for(unsigned i=0; i<forces.size(); i++) {
489 19476 : fforces[i] += forces[i];
490 : }
491 416 : for(unsigned i=0; i<3; ++i)
492 1248 : for(unsigned j=0; j<3; ++j) {
493 936 : fforces[forces.size()+3*i+j] = virial[i][j];
494 : }
495 104 : unsigned ind=0;
496 104 : setForcesOnAtoms( fforces, ind );
497 104 : }
498 :
499 104 : void Plumed::update() {
500 104 : Tools::DirectoryChanger directoryChanger(directory.c_str());
501 104 : if(root) {
502 59 : p.cmd("update");
503 : }
504 104 : comm.Bcast(stop,0);
505 104 : if(stop) {
506 0 : log<<" Action " << getLabel()<<" asked to stop\n";
507 0 : plumed.stop();
508 : }
509 104 : }
510 :
511 : }
512 : }
|
