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