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             : 
      41             : #include <algorithm>
      42             : 
      43             : namespace PLMD {
      44             : namespace generic {
      45             : 
      46             : //+PLUMEDOC GENERIC EFFECTIVE_ENERGY_DRIFT
      47             : /*
      48             : Print the effective energy drift
      49             : 
      50             : The method that is used to calculate the effective energy drift here is described in the
      51             : paper in the bibliography.
      52             : 
      53             : ## Examples
      54             : 
      55             : This input is to monitor the effective energy drift for a metadynamics
      56             : simulation on the Debye-Huckel energy. Since this variable is very expensive,
      57             : it could be conveniently computed every second step.
      58             : 
      59             : ```plumed
      60             : dh: DHENERGY GROUPA=1-10 GROUPB=11-20 EPSILON=80.0 I=0.1 TEMP=300.0
      61             : METAD ARG=dh HEIGHT=0.5 SIGMA=0.1 PACE=500 STRIDE=2
      62             : EFFECTIVE_ENERGY_DRIFT PRINT_STRIDE=100 FILE=eff
      63             : ```
      64             : 
      65             : This exampls shows how to monitor if a restraint is too stiff
      66             : 
      67             : ```plumed
      68             : d: DISTANCE ATOMS=10,20
      69             : RESTRAINT ARG=d KAPPA=100000 AT=0.6
      70             : EFFECTIVE_ENERGY_DRIFT PRINT_STRIDE=100 FILE=eff
      71             : ```
      72             : 
      73             : */
      74             : //+ENDPLUMEDOC
      75             : 
      76             : 
      77             : class EffectiveEnergyDrift:
      78             :   public ActionPilot {
      79             :   OFile output;
      80             :   long long int printStride;
      81             :   std::string fmt;
      82             : 
      83             :   double eed;
      84             : 
      85             :   std::vector<ActionWithValue*> biases;
      86             : 
      87             :   long long int pDdStep;
      88             :   int nLocalAtoms;
      89             :   int 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         450 :   void calculate() override {};
     128         450 :   void apply() override {};
     129             :   void update() override;
     130             : };
     131             : 
     132             : PLUMED_REGISTER_ACTION(EffectiveEnergyDrift,"EFFECTIVE_ENERGY_DRIFT")
     133             : 
     134          11 : void EffectiveEnergyDrift::registerKeywords( Keywords& keys ) {
     135          11 :   Action::registerKeywords( keys );
     136          11 :   ActionPilot::registerKeywords( keys );
     137             : 
     138          11 :   keys.add("compulsory","STRIDE","1","should be set to 1. Effective energy drift computation has to be active at each step.");
     139          11 :   keys.add("compulsory", "FILE", "file on which to output the effective energy drift.");
     140          11 :   keys.add("compulsory", "PRINT_STRIDE", "frequency to which output the effective energy drift on FILE");
     141          11 :   keys.addFlag("ENSEMBLE",false,"Set to TRUE if you want to average over multiple replicas.");
     142          11 :   keys.add("optional","FMT","the format that should be used to output real numbers");
     143          11 :   keys.use("RESTART");
     144          11 :   keys.use("UPDATE_FROM");
     145          11 :   keys.use("UPDATE_UNTIL");
     146          11 :   keys.addDOI("10.1021/ct5007086");
     147          11 : }
     148             : 
     149           9 : EffectiveEnergyDrift::EffectiveEnergyDrift(const ActionOptions&ao):
     150             :   Action(ao),
     151             :   ActionPilot(ao),
     152           9 :   fmt("%f"),
     153           9 :   eed(0.0),
     154           9 :   nProc(plumed.comm.Get_size()),
     155           9 :   initialBias(0.0),
     156           9 :   isFirstStep(true),
     157           9 :   ensemble(false),
     158           9 :   pbc_action(NULL),
     159           9 :   domains(NULL),
     160           9 :   posx(NULL),
     161           9 :   posy(NULL),
     162          27 :   posz(NULL) {
     163             :   //stride must be == 1
     164           9 :   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          18 :   parse("FILE",fileName);
     171           9 :   if(fileName.length()==0) {
     172           0 :     error("name out output file was not specified\n");
     173             :   }
     174           9 :   output.link(*this);
     175           9 :   output.open(fileName);
     176             : 
     177             :   //parse PRINT_STRIDE
     178           9 :   parse("PRINT_STRIDE",printStride);
     179             : 
     180             : 
     181             :   //parse FMT
     182           9 :   parse("FMT",fmt);
     183           9 :   fmt=" "+fmt;
     184           9 :   log.printf("  with format %s\n",fmt.c_str());
     185             : 
     186             :   //parse ENSEMBLE
     187           9 :   ensemble=false;
     188           9 :   parseFlag("ENSEMBLE",ensemble);
     189           9 :   if(ensemble&&comm.Get_rank()==0) {
     190           0 :     if(multi_sim_comm.Get_size()<2) {
     191           0 :       error("You CANNOT run Replica-Averaged simulations without running multiple replicas!\n");
     192             :     }
     193             :   }
     194             : 
     195          18 :   log<<"Bibliography "<<cite("Ferrarotti, Bottaro, Perez-Villa, and Bussi, J. Chem. Theory Comput. 11, 139 (2015)")<<"\n";
     196             : 
     197             :   //construct biases from ActionWithValue with a component named bias
     198           9 :   std::vector<ActionWithValue*> tmpActions=plumed.getActionSet().select<ActionWithValue*>();
     199         100 :   for(unsigned i=0; i<tmpActions.size(); i++)
     200         182 :     if(tmpActions[i]->exists(tmpActions[i]->getLabel()+".bias")) {
     201           9 :       biases.push_back(tmpActions[i]);
     202             :     }
     203             : 
     204             :   //resize counters and displacements useful to communicate with MPI_Allgatherv
     205           9 :   indexCnt.resize(nProc);
     206           9 :   indexDsp.resize(nProc);
     207           9 :   dataCnt.resize(nProc);
     208           9 :   dataDsp.resize(nProc);
     209             :   // Retrieve the box
     210           9 :   pbc_action=plumed.getActionSet().selectWithLabel<PbcAction*>("Box");
     211             :   // Get the domain decomposition object
     212           9 :   std::vector<DomainDecomposition*> ddact=plumed.getActionSet().select<DomainDecomposition*>();
     213           9 :   if( ddact.size()>1 ) {
     214           0 :     warning("found more than one interface so don't know get gatindex");
     215             :   }
     216           9 :   domains = ddact[0];
     217           9 :   std::vector<ActionToPutData*> inputs=plumed.getActionSet().select<ActionToPutData*>();
     218          73 :   for(const auto & pp : inputs ) {
     219         128 :     if( pp->getRole()=="x" ) {
     220           9 :       posx = pp;
     221             :     }
     222         128 :     if( pp->getRole()=="y" ) {
     223           9 :       posy = pp;
     224             :     }
     225         128 :     if( pp->getRole()=="z" ) {
     226           9 :       posz = pp;
     227             :     }
     228             :   }
     229           9 :   plumed_assert( posx && posy && posz );
     230             :   //resize the received buffers
     231           9 :   indexR.resize((posx->copyOutput(0))->getShape()[0]);
     232           9 :   dataR.resize((posx->copyOutput(0))->getShape()[0]*6);
     233           9 :   backmap.resize((posx->copyOutput(0))->getShape()[0]);
     234           9 : }
     235             : 
     236          18 : EffectiveEnergyDrift::~EffectiveEnergyDrift() {
     237             : 
     238          18 : }
     239             : 
     240         450 : void EffectiveEnergyDrift::update() {
     241         450 :   Pbc & tpbc(pbc_action->getPbc());
     242             :   bool pbc=tpbc.isSet();
     243             : 
     244             :   //retrieve data of local atoms
     245         450 :   const std::vector<int>& gatindex = domains->getGatindex();
     246         450 :   nLocalAtoms = gatindex.size();
     247         450 :   xpositions.resize( gatindex.size() );
     248         450 :   posx->getLocalValues( xpositions );
     249         450 :   ypositions.resize( gatindex.size() );
     250         450 :   posy->getLocalValues( ypositions );
     251         450 :   zpositions.resize( gatindex.size() );
     252         450 :   posz->getLocalValues( zpositions );
     253         450 :   positions.resize( gatindex.size() );
     254         450 :   forces.resize( gatindex.size() );
     255       16650 :   for(unsigned i=0; i<gatindex.size(); ++i ) {
     256       16200 :     positions[i][0] = xpositions[i];
     257       16200 :     positions[i][1] = ypositions[i];
     258       16200 :     positions[i][2] = zpositions[i];
     259       16200 :     forces[i][0] = (posx->copyOutput(0))->getForce( gatindex[i] );
     260       16200 :     forces[i][1] = (posy->copyOutput(0))->getForce( gatindex[i] );
     261       16200 :     forces[i][2] = (posz->copyOutput(0))->getForce( gatindex[i] );
     262             :   }
     263         450 :   if(pbc) {
     264         450 :     Tensor B=tpbc.getBox();
     265         450 :     Tensor IB=tpbc.getInvBox();
     266         450 :     #pragma omp parallel for
     267             :     for(unsigned i=0; i<positions.size(); ++i) {
     268             :       positions[i]=matmul(positions[i],IB);
     269             :       forces[i]=matmul(B,forces[i]);
     270             :     }
     271         450 :     box=B;
     272         450 :     Tensor virial;
     273         450 :     Value* boxValue = pbc_action->copyOutput(0);
     274        1800 :     for(unsigned i=0; i<3; ++i)
     275        5400 :       for(unsigned j=0; j<3; ++j) {
     276        4050 :         virial[i][j]=boxValue->getForce(3*i+j);
     277             :       }
     278         450 :     fbox=matmul(transpose(inverse(box)),virial);
     279             :   }
     280             : 
     281             :   //init stored data at the first step
     282         450 :   if(isFirstStep) {
     283           9 :     pDdStep=0;
     284           9 :     pGatindex = domains->getGatindex();
     285           9 :     pNLocalAtoms = pGatindex.size();
     286           9 :     pPositions=positions;
     287           9 :     pForces=forces;
     288           9 :     pbox=box;
     289           9 :     pfbox=fbox;
     290           9 :     initialBias=plumed.getBias();
     291             : 
     292           9 :     isFirstStep=false;
     293             :   }
     294             : 
     295             :   //if the dd has changed we have to reshare the stored data
     296         450 :   if(pDdStep<domains->getDdStep() && nLocalAtoms<(posx->copyOutput(0))->getShape()[0]) {
     297             :     //prepare the data to be sent
     298         204 :     indexS.resize(pNLocalAtoms);
     299         204 :     dataS.resize(pNLocalAtoms*6);
     300             : 
     301        5712 :     for(int i=0; i<pNLocalAtoms; i++) {
     302        5508 :       indexS[i] = pGatindex[i];
     303        5508 :       dataS[i*6] = pPositions[i][0];
     304        5508 :       dataS[i*6+1] = pPositions[i][1];
     305        5508 :       dataS[i*6+2] = pPositions[i][2];
     306        5508 :       dataS[i*6+3] = pForces[i][0];
     307        5508 :       dataS[i*6+4] = pForces[i][1];
     308        5508 :       dataS[i*6+5] = pForces[i][2];
     309             :     }
     310             : 
     311             :     //setup the counters and displacements for the communication
     312         204 :     plumed.comm.Allgather(&pNLocalAtoms,1,&indexCnt[0],1);
     313         204 :     indexDsp[0] = 0;
     314        1020 :     for(int i=0; i<nProc; i++) {
     315         816 :       dataCnt[i] = indexCnt[i]*6;
     316             : 
     317         816 :       if(i+1<nProc) {
     318         612 :         indexDsp[i+1] = indexDsp[i]+indexCnt[i];
     319             :       }
     320         816 :       dataDsp[i] = indexDsp[i]*6;
     321             :     }
     322             : 
     323             :     //share stored data
     324         402 :     plumed.comm.Allgatherv((!indexS.empty()?&indexS[0]:NULL), pNLocalAtoms, &indexR[0], &indexCnt[0], &indexDsp[0]);
     325         402 :     plumed.comm.Allgatherv((!dataS.empty()?&dataS[0]:NULL), pNLocalAtoms*6, &dataR[0], &dataCnt[0], &dataDsp[0]);
     326             : 
     327             :     //resize vectors to store the proper amount of data
     328         204 :     pGatindex.resize(nLocalAtoms);
     329         204 :     pPositions.resize(nLocalAtoms);
     330         204 :     pForces.resize(nLocalAtoms);
     331             : 
     332             :     //compute backmap
     333       22236 :     for(unsigned j=0; j<indexR.size(); j++) {
     334       22032 :       backmap[indexR[j]]=j;
     335             :     }
     336             : 
     337             :     //fill the vectors pGatindex, pPositions and pForces
     338        5712 :     for(int i=0; i<nLocalAtoms; i++) {
     339        5508 :       int glb=backmap[gatindex[i]];
     340        5508 :       pGatindex[i] = indexR[glb];
     341        5508 :       pPositions[i][0] = dataR[glb*6];
     342        5508 :       pPositions[i][1] = dataR[glb*6+1];
     343        5508 :       pPositions[i][2] = dataR[glb*6+2];
     344        5508 :       pForces[i][0] = dataR[glb*6+3];
     345        5508 :       pForces[i][1] = dataR[glb*6+4];
     346        5508 :       pForces[i][2] = dataR[glb*6+5];
     347             :     }
     348             :   }
     349             : 
     350             :   //compute the effective energy drift on local atoms
     351             : 
     352         450 :   double eed_tmp=eed;
     353         450 :   #pragma omp parallel for reduction(+:eed_tmp)
     354             :   for(int i=0; i<nLocalAtoms; i++) {
     355             :     Vector dst=delta(pPositions[i],positions[i]);
     356             :     if(pbc)
     357             :       for(unsigned k=0; k<3; k++) {
     358             :         dst[k]=Tools::pbc(dst[k]);
     359             :       }
     360             :     eed_tmp += dotProduct(dst, forces[i]+pForces[i])*0.5;
     361             :   }
     362             : 
     363         450 :   eed=eed_tmp;
     364             : 
     365         450 :   if(plumed.comm.Get_rank()==0) {
     366         600 :     for(unsigned i=0; i<3; i++)
     367        1800 :       for(unsigned j=0; j<3; j++) {
     368        1350 :         eed-=0.5*(pfbox(i,j)+fbox(i,j))*(box(i,j)-pbox(i,j));
     369             :       }
     370             :   }
     371             : 
     372             : 
     373             :   //print the effective energy drift on FILE with frequency PRINT_STRIDE
     374         450 :   if(plumed.getStep()%printStride==0) {
     375         450 :     double eedSum = eed;
     376             :     double bias = 0.0;
     377             : 
     378             :     //we cannot just use plumed.getBias() because it will be ==0.0 if PRINT_STRIDE
     379             :     //is not a multiple of the bias actions stride
     380         900 :     for(unsigned i=0; i<biases.size(); i++) {
     381         450 :       bias+=biases[i]->getOutputQuantity("bias");
     382             :     }
     383             : 
     384         450 :     plumed.comm.Sum(&eedSum,1);
     385             : 
     386         450 :     double effective = eedSum+bias-initialBias-plumed.getWork();
     387             :     // this is to take into account ensemble averaging
     388         450 :     if(ensemble) {
     389           0 :       if(plumed.comm.Get_rank()==0) {
     390           0 :         plumed.multi_sim_comm.Sum(&effective,1);
     391             :       } else {
     392           0 :         effective=0.;
     393             :       }
     394           0 :       plumed.comm.Sum(&effective,1);
     395             :     }
     396         450 :     output.fmtField(" %f");
     397         450 :     output.printField("time",getTime());
     398         450 :     output.fmtField(fmt);
     399         450 :     output.printField("effective-energy",effective);
     400         450 :     output.printField();
     401             :   }
     402             : 
     403             :   //store the data of the current step
     404         450 :   pDdStep = domains->getDdStep();
     405         450 :   pNLocalAtoms = nLocalAtoms;
     406             :   pPositions.swap(positions);
     407             :   pForces.swap(forces);
     408         450 :   pbox=box;
     409         450 :   pfbox=fbox;
     410         450 : }
     411             : 
     412             : }
     413             : }