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