LCOV - code coverage report
Current view: top level - generic - Plumed.cpp (source / functions) Hit Total Coverage
Test: plumed test coverage Lines: 213 229 93.0 %
Date: 2025-04-08 18:07:56 Functions: 10 12 83.3 %

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

Generated by: LCOV version 1.16