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