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