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