Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2013-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             : 
      23             : /*
      24             :  This class was originally written by Marco Jacopo Ferrarotti
      25             :  (marco.ferrarotti@gmail.com) and Giovanni Bussi
      26             : */
      27             : 
      28             : #include "core/Action.h"
      29             : #include "core/ActionPilot.h"
      30             : #include "core/ActionWithValue.h"
      31             : #include "core/ActionSet.h"
      32             : #include "core/ActionRegister.h"
      33             : #include "core/DomainDecomposition.h"
      34             : #include "core/ActionToPutData.h"
      35             : #include "core/PbcAction.h"
      36             : #include "core/PlumedMain.h"
      37             : 
      38             : #include "tools/File.h"
      39             : #include "tools/Pbc.h"
      40             : #include "tools/OpenMP.h"
      41             : 
      42             : #include <algorithm>
      43             : 
      44             : namespace PLMD {
      45             : namespace generic {
      46             : 
      47             : //+PLUMEDOC GENERIC EFFECTIVE_ENERGY_DRIFT
      48             : /*
      49             : Print the effective energy drift
      50             : 
      51             : The method that is used to calculate the effective energy drift here is described in the
      52             : paper in the bibliography.
      53             : 
      54             : ## Examples
      55             : 
      56             : This input is to monitor the effective energy drift for a metadynamics
      57             : simulation on the Debye-Huckel energy. Since this variable is very expensive,
      58             : it could be conveniently computed every second step.
      59             : 
      60             : ```plumed
      61             : dh: DHENERGY GROUPA=1-10 GROUPB=11-20 EPSILON=80.0 I=0.1 TEMP=300.0
      62             : METAD ARG=dh HEIGHT=0.5 SIGMA=0.1 PACE=500 STRIDE=2
      63             : EFFECTIVE_ENERGY_DRIFT PRINT_STRIDE=100 FILE=eff
      64             : ```
      65             : 
      66             : This exampls shows how to monitor if a restraint is too stiff
      67             : 
      68             : ```plumed
      69             : d: DISTANCE ATOMS=10,20
      70             : RESTRAINT ARG=d KAPPA=100000 AT=0.6
      71             : EFFECTIVE_ENERGY_DRIFT PRINT_STRIDE=100 FILE=eff
      72             : ```
      73             : 
      74             : */
      75             : //+ENDPLUMEDOC
      76             : 
      77             : 
      78             : class EffectiveEnergyDrift:
      79             :   public ActionPilot {
      80             :   OFile output;
      81             :   long long int printStride;
      82             :   std::string fmt;
      83             : 
      84             :   double eed;
      85             : 
      86             :   std::vector<ActionWithValue*> biases;
      87             : 
      88             :   unsigned nLocalAtoms;
      89             :   unsigned pNLocalAtoms;
      90             :   std::vector<int> pGatindex;
      91             :   std::vector<double> xpositions;
      92             :   std::vector<double> ypositions;
      93             :   std::vector<double> zpositions;
      94             :   std::vector<Vector> positions;
      95             :   std::vector<Vector> pPositions;
      96             :   std::vector<Vector> forces;
      97             :   std::vector<Vector> pForces;
      98             :   Tensor box,pbox;
      99             :   Tensor fbox,pfbox;
     100             : 
     101             :   const int nProc;
     102             :   std::vector<int> indexCnt;
     103             :   std::vector<int> indexDsp;
     104             :   std::vector<int> dataCnt;
     105             :   std::vector<int> dataDsp;
     106             :   std::vector<int> indexS;
     107             :   std::vector<int> indexR;
     108             :   std::vector<double> dataS;
     109             :   std::vector<double> dataR;
     110             :   std::vector<int> backmap;
     111             : 
     112             :   double initialBias;
     113             :   bool isFirstStep;
     114             : 
     115             :   bool ensemble;
     116             :   PbcAction* pbc_action;
     117             :   DomainDecomposition* domains;
     118             :   ActionToPutData* posx;
     119             :   ActionToPutData* posy;
     120             :   ActionToPutData* posz;
     121             : public:
     122             :   explicit EffectiveEnergyDrift(const ActionOptions&);
     123             :   ~EffectiveEnergyDrift();
     124             : 
     125             :   static void registerKeywords( Keywords& keys );
     126             : 
     127         474 :   void calculate() override {};
     128         474 :   void apply() override {};
     129             :   void update() override;
     130             : };
     131             : 
     132             : PLUMED_REGISTER_ACTION(EffectiveEnergyDrift,"EFFECTIVE_ENERGY_DRIFT")
     133             : 
     134          15 : void EffectiveEnergyDrift::registerKeywords( Keywords& keys ) {
     135          15 :   Action::registerKeywords( keys );
     136          15 :   ActionPilot::registerKeywords( keys );
     137             : 
     138          15 :   keys.add("compulsory","STRIDE","1","should be set to 1. Effective energy drift computation has to be active at each step.");
     139          15 :   keys.add("compulsory", "FILE", "file on which to output the effective energy drift.");
     140          15 :   keys.add("compulsory", "PRINT_STRIDE", "frequency to which output the effective energy drift on FILE");
     141          15 :   keys.addFlag("ENSEMBLE",false,"Set to TRUE if you want to average over multiple replicas.");
     142          15 :   keys.add("optional","FMT","the format that should be used to output real numbers");
     143          15 :   keys.use("RESTART");
     144          15 :   keys.use("UPDATE_FROM");
     145          15 :   keys.use("UPDATE_UNTIL");
     146          15 :   keys.addDOI("10.1021/ct5007086");
     147          15 : }
     148             : 
     149          11 : EffectiveEnergyDrift::EffectiveEnergyDrift(const ActionOptions&ao):
     150             :   Action(ao),
     151             :   ActionPilot(ao),
     152          11 :   fmt("%f"),
     153          11 :   eed(0.0),
     154          11 :   nProc(plumed.comm.Get_size()),
     155          11 :   initialBias(0.0),
     156          11 :   isFirstStep(true),
     157          11 :   ensemble(false),
     158          11 :   pbc_action(NULL),
     159          11 :   domains(NULL),
     160          11 :   posx(NULL),
     161          11 :   posy(NULL),
     162          11 :   posz(NULL) {
     163             :   //stride must be == 1
     164          11 :   if(getStride()!=1) {
     165           0 :     error("EFFECTIVE_ENERGY_DRIFT must have STRIDE=1 to work properly");
     166             :   }
     167             : 
     168             :   //parse and open FILE
     169             :   std::string fileName;
     170          22 :   parse("FILE",fileName);
     171          11 :   if(fileName.length()==0) {
     172           0 :     error("name out output file was not specified\n");
     173             :   }
     174          11 :   output.link(*this);
     175          11 :   output.open(fileName);
     176             : 
     177             :   //parse PRINT_STRIDE
     178          11 :   parse("PRINT_STRIDE",printStride);
     179             : 
     180             : 
     181             :   //parse FMT
     182          11 :   parse("FMT",fmt);
     183          11 :   fmt=" "+fmt;
     184          11 :   log.printf("  with format %s\n",fmt.c_str());
     185             : 
     186             :   //parse ENSEMBLE
     187          11 :   ensemble=false;
     188          11 :   parseFlag("ENSEMBLE",ensemble);
     189          11 :   if(ensemble && getNumberOfReplicas()<2) {
     190           0 :     error("You CANNOT run Replica-Averaged simulations without running multiple replicas!\n");
     191             :   }
     192             : 
     193          11 :   log<<"Bibliography "<<cite("Ferrarotti, Bottaro, Perez-Villa, and Bussi, J. Chem. Theory Comput. 11, 139 (2015)")<<"\n";
     194             : 
     195             :   //construct biases from ActionWithValue with a component named bias
     196          11 :   std::vector<ActionWithValue*> tmpActions=plumed.getActionSet().select<ActionWithValue*>();
     197         128 :   for(unsigned i=0; i<tmpActions.size(); i++)
     198         234 :     if(tmpActions[i]->exists(tmpActions[i]->getLabel()+".bias")) {
     199          13 :       biases.push_back(tmpActions[i]);
     200             :     }
     201             : 
     202             :   //resize counters and displacements useful to communicate with MPI_Allgatherv
     203          11 :   indexCnt.resize(nProc);
     204          11 :   indexDsp.resize(nProc);
     205          11 :   dataCnt.resize(nProc);
     206          11 :   dataDsp.resize(nProc);
     207             :   // Retrieve the box
     208          11 :   pbc_action=plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
     209             :   // Get the domain decomposition object
     210          11 :   std::vector<DomainDecomposition*> ddact=plumed.getActionSet().select<DomainDecomposition*>();
     211          11 :   if( ddact.empty() ) {
     212           0 :     error("cannot find the atoms passed by the MD code");
     213             :   }
     214          11 :   if( ddact.size()>1 ) {
     215           0 :     warning("found more than one interface so don't know get gatindex");
     216             :   }
     217          11 :   domains = ddact[0];
     218          11 :   std::vector<ActionToPutData*> inputs=plumed.getActionSet().select<ActionToPutData*>();
     219          89 :   for(const auto & pp : inputs ) {
     220          78 :     if( pp->getRole()=="x" ) {
     221          11 :       posx = pp;
     222             :     }
     223          78 :     if( pp->getRole()=="y" ) {
     224          11 :       posy = pp;
     225             :     }
     226          78 :     if( pp->getRole()=="z" ) {
     227          11 :       posz = pp;
     228             :     }
     229             :   }
     230          11 :   plumed_assert( posx && posy && posz );
     231             :   //resize the received buffers
     232          11 :   indexR.resize((posx->copyOutput(0))->getShape()[0]);
     233          11 :   dataR.resize((posx->copyOutput(0))->getShape()[0]*6);
     234          11 :   backmap.resize((posx->copyOutput(0))->getShape()[0]);
     235          22 : }
     236             : 
     237          22 : EffectiveEnergyDrift::~EffectiveEnergyDrift() {
     238             : 
     239          22 : }
     240             : 
     241         474 : void EffectiveEnergyDrift::update() {
     242         474 :   Pbc & tpbc(pbc_action->getPbc());
     243             :   bool pbc=tpbc.isSet();
     244             : 
     245             :   //retrieve data of local atoms
     246         474 :   const std::vector<int>& gatindex = domains->getGatindex();
     247         474 :   nLocalAtoms = gatindex.size();
     248         474 :   xpositions.resize( gatindex.size() );
     249         474 :   posx->getLocalValues( xpositions );
     250         474 :   ypositions.resize( gatindex.size() );
     251         474 :   posy->getLocalValues( ypositions );
     252         474 :   zpositions.resize( gatindex.size() );
     253         474 :   posz->getLocalValues( zpositions );
     254         474 :   positions.resize( gatindex.size() );
     255         474 :   forces.resize( gatindex.size() );
     256         474 :   const Value* fx=posx->copyOutput(0);
     257         474 :   const Value* fy=posy->copyOutput(0);
     258         474 :   const Value* fz=posz->copyOutput(0);
     259       17154 :   for(unsigned i=0; i<gatindex.size(); ++i ) {
     260       16680 :     positions[i][0] = xpositions[i];
     261       16680 :     positions[i][1] = ypositions[i];
     262       16680 :     positions[i][2] = zpositions[i];
     263       16680 :     forces[i][0] = fx->getForce( gatindex[i] );
     264       16680 :     forces[i][1] = fy->getForce( gatindex[i] );
     265       16680 :     forces[i][2] = fz->getForce( gatindex[i] );
     266             :   }
     267         474 :   if(pbc) {
     268         474 :     Tensor B=tpbc.getBox();
     269         474 :     Tensor IB=tpbc.getInvBox();
     270         474 :     #pragma omp parallel for num_threads(OpenMP::getGoodNumThreads(positions))
     271             :     for(unsigned i=0; i<positions.size(); ++i) {
     272             :       positions[i]=matmul(positions[i],IB);
     273             :       forces[i]=matmul(B,forces[i]);
     274             :     }
     275         474 :     box=B;
     276         474 :     Tensor virial;
     277         474 :     Value* boxValue = pbc_action->copyOutput(0);
     278        1896 :     for(unsigned i=0; i<3; ++i)
     279        5688 :       for(unsigned j=0; j<3; ++j) {
     280        4266 :         virial[i][j]=boxValue->getForce(3*i+j);
     281             :       }
     282         474 :     fbox=matmul(transpose(inverse(box)),virial);
     283             :   }
     284             : 
     285             :   //init stored data at the first step
     286         474 :   if(isFirstStep) {
     287          11 :     pGatindex = domains->getGatindex();
     288          11 :     pNLocalAtoms = pGatindex.size();
     289          11 :     pPositions=positions;
     290          11 :     pForces=forces;
     291          11 :     pbox=box;
     292          11 :     pfbox=fbox;
     293          11 :     initialBias=plumed.getBias();
     294             : 
     295          11 :     isFirstStep=false;
     296             :   }
     297             : 
     298             :   // The local atoms, or their order, may have changed since the previous step:
     299             :   // domain decomposition, or a serial MD code that reorders its atoms (as LAMMPS does).
     300             :   // In that case the stored data is remapped. The local lists are compared at each step:
     301             :   // DomainDecomposition::getDdStep() cannot be used, since most MD codes pass the
     302             :   // local atoms at every step.
     303         474 :   int changed=(pGatindex.size()!=gatindex.size() || !std::equal(gatindex.begin(),gatindex.end(),pGatindex.begin()));
     304         474 :   if(nProc>1) {
     305         400 :     plumed.comm.Sum(changed);
     306             :   }
     307         474 :   if(changed) {
     308             :     unsigned nrecv=0;
     309         215 :     if(nProc>1) {
     310             :       //prepare the data to be sent
     311         204 :       indexS.resize(pNLocalAtoms);
     312         204 :       dataS.resize(pNLocalAtoms*6);
     313        5712 :       for(unsigned i=0; i<pNLocalAtoms; i++) {
     314        5508 :         indexS[i] = pGatindex[i];
     315        5508 :         dataS[i*6] = pPositions[i][0];
     316        5508 :         dataS[i*6+1] = pPositions[i][1];
     317        5508 :         dataS[i*6+2] = pPositions[i][2];
     318        5508 :         dataS[i*6+3] = pForces[i][0];
     319        5508 :         dataS[i*6+4] = pForces[i][1];
     320        5508 :         dataS[i*6+5] = pForces[i][2];
     321             :       }
     322             :       //setup the counters and displacements for the communication
     323         204 :       plumed.comm.Allgather(&pNLocalAtoms,1,&indexCnt[0],1);
     324         204 :       indexDsp[0] = 0;
     325        1020 :       for(int i=0; i<nProc; i++) {
     326         816 :         dataCnt[i] = indexCnt[i]*6;
     327         816 :         if(i+1<nProc) {
     328         612 :           indexDsp[i+1] = indexDsp[i]+indexCnt[i];
     329             :         }
     330         816 :         dataDsp[i] = indexDsp[i]*6;
     331         816 :         nrecv += indexCnt[i];
     332             :       }
     333         204 :       plumed_assert(nrecv<=indexR.size());
     334             :       //share stored data
     335         402 :       plumed.comm.Allgatherv((!indexS.empty()?&indexS[0]:NULL), pNLocalAtoms, &indexR[0], &indexCnt[0], &indexDsp[0]);
     336         402 :       plumed.comm.Allgatherv((!dataS.empty()?&dataS[0]:NULL), pNLocalAtoms*6, &dataR[0], &dataCnt[0], &dataDsp[0]);
     337             :     } else {
     338             :       // a single process: the stored data is only reordered
     339          11 :       nrecv=pNLocalAtoms;
     340          11 :       plumed_assert(nrecv<=indexR.size());
     341         231 :       for(unsigned i=0; i<pNLocalAtoms; i++) {
     342         220 :         indexR[i] = pGatindex[i];
     343         220 :         dataR[i*6] = pPositions[i][0];
     344         220 :         dataR[i*6+1] = pPositions[i][1];
     345         220 :         dataR[i*6+2] = pPositions[i][2];
     346         220 :         dataR[i*6+3] = pForces[i][0];
     347         220 :         dataR[i*6+4] = pForces[i][1];
     348         220 :         dataR[i*6+5] = pForces[i][2];
     349             :       }
     350             :     }
     351             :     //resize vectors to store the proper amount of data
     352         215 :     pPositions.resize(nLocalAtoms);
     353         215 :     pForces.resize(nLocalAtoms);
     354             :     //compute backmap (atoms that were not stored anywhere keep -1)
     355             :     std::fill(backmap.begin(),backmap.end(),-1);
     356       22467 :     for(unsigned j=0; j<nrecv; j++) {
     357       22252 :       backmap[indexR[j]]=j;
     358             :     }
     359             :     //fill the vectors pPositions and pForces
     360        5943 :     for(unsigned i=0; i<nLocalAtoms; i++) {
     361        5728 :       int glb=backmap[gatindex[i]];
     362        5728 :       if(glb<0) {
     363             :         // no data from the previous step: this atom does not contribute to this step
     364           0 :         pPositions[i] = positions[i];
     365           0 :         pForces[i] = forces[i];
     366           0 :         continue;
     367             :       }
     368        5728 :       pPositions[i][0] = dataR[glb*6];
     369        5728 :       pPositions[i][1] = dataR[glb*6+1];
     370        5728 :       pPositions[i][2] = dataR[glb*6+2];
     371        5728 :       pForces[i][0] = dataR[glb*6+3];
     372        5728 :       pForces[i][1] = dataR[glb*6+4];
     373        5728 :       pForces[i][2] = dataR[glb*6+5];
     374             :     }
     375         215 :     pGatindex=gatindex;
     376             :   }
     377             : 
     378             :   //compute the effective energy drift on local atoms
     379             : 
     380         474 :   double eed_tmp=eed;
     381         474 :   #pragma omp parallel for reduction(+:eed_tmp) num_threads(OpenMP::getGoodNumThreads(positions))
     382             :   for(unsigned i=0; i<nLocalAtoms; i++) {
     383             :     Vector dst=delta(pPositions[i],positions[i]);
     384             :     if(pbc)
     385             :       for(unsigned k=0; k<3; k++) {
     386             :         dst[k]=Tools::pbc(dst[k]);
     387             :       }
     388             :     eed_tmp += dotProduct(dst, forces[i]+pForces[i])*0.5;
     389             :   }
     390             : 
     391         474 :   eed=eed_tmp;
     392             : 
     393         474 :   if(plumed.comm.Get_rank()==0) {
     394         696 :     for(unsigned i=0; i<3; i++)
     395        2088 :       for(unsigned j=0; j<3; j++) {
     396        1566 :         eed-=0.5*(pfbox(i,j)+fbox(i,j))*(box(i,j)-pbox(i,j));
     397             :       }
     398             :   }
     399             : 
     400             : 
     401             :   //print the effective energy drift on FILE with frequency PRINT_STRIDE
     402         474 :   if(plumed.getStep()%printStride==0) {
     403         474 :     double eedSum = eed;
     404             :     double bias = 0.0;
     405             : 
     406             :     //we cannot just use plumed.getBias() because it will be ==0.0 if PRINT_STRIDE
     407             :     //is not a multiple of the bias actions stride
     408         972 :     for(unsigned i=0; i<biases.size(); i++) {
     409         498 :       bias+=biases[i]->getOutputQuantity("bias");
     410             :     }
     411             : 
     412         474 :     plumed.comm.Sum(&eedSum,1);
     413             : 
     414         474 :     double effective = eedSum+bias-initialBias-plumed.getWork();
     415             :     // this is to take into account ensemble averaging
     416         474 :     if(ensemble) {
     417           0 :       if(plumed.comm.Get_rank()==0) {
     418           0 :         plumed.multi_sim_comm.Sum(&effective,1);
     419             :       } else {
     420           0 :         effective=0.;
     421             :       }
     422           0 :       plumed.comm.Sum(&effective,1);
     423             :     }
     424         474 :     output.fmtField(" %f");
     425         474 :     output.printField("time",getTime());
     426         474 :     output.fmtField(fmt);
     427         474 :     output.printField("effective-energy",effective);
     428         474 :     output.printField();
     429             :   }
     430             : 
     431             :   //store the data of the current step
     432         474 :   pNLocalAtoms = nLocalAtoms;
     433             :   pPositions.swap(positions);
     434             :   pForces.swap(forces);
     435         474 :   pbox=box;
     436         474 :   pfbox=fbox;
     437         474 : }
     438             : 
     439             : }
     440             : }