Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 : #include "CLTool.h"
23 : #include "core/CLToolRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "tools/Vector.h"
26 : #include "tools/Random.h"
27 : #include "tools/Communicator.h"
28 : #include <string>
29 : #include <cstdio>
30 : #include <vector>
31 : #include <memory>
32 :
33 : //+PLUMEDOC TOOLS pesmd
34 : /*
35 : Pesmd allows one to do (biased) Langevin dynamics on a two-dimensional potential energy surface.
36 :
37 : The energy landscape that you are moving about on is specified using in a plumed input file.
38 : For example the following example input uses [MATHEVAL](MATHEVAL.md) to define a two dimensional potential.
39 :
40 : ```plumed
41 : d1: DISTANCE ATOMS=1,2 COMPONENTS
42 : ff: MATHEVAL ARG=d1.x,d1.y PERIODIC=NO FUNC=x*x+y*y
43 : bb: BIASVALUE ARG=ff
44 : ```
45 :
46 : Atom 1 is placed at the origin. The x and y components on our surface are the
47 : positions of the particle on our two dimensional energy landscape. By calculating the
48 : vector connecting atom 1 (the origin) to atom 2 (the position of our particle) we are thus
49 : getting the position of the atom on the energy landscape. This is then inserted into the function
50 : that is calculated on the second line. The value of this function is then used as a bias.
51 :
52 : We can also specify a potential on a grid and look at the dynamics on this function using pesmd. To do
53 : so you would use the DISTANCE command as in the input above together with an [EXTERNAL](EXTERNAL.md) action
54 : that takes d1.x and d1.y as arguments. In this way you can use pesmd to do a dynamics on a free energy surface calculated using metadynamics
55 : and sum_hills. Notice also, that once we have defined our potential we can use all the biasing functions
56 : within plumed in addition in order to do a biased dynamics on the potential energy landscape of interest.
57 :
58 : You run a Langevin simulation using pesmd with the following command:
59 :
60 : ```plumed
61 : plumed pesmd < input
62 : ```
63 :
64 : The input to pesmd is specified in an input file that contains one directive per line like the one shown below.
65 :
66 : ```plumed
67 : #TOOL=pesmd
68 : temperature 0.722
69 : tstep 0.005
70 : friction 1
71 : dimension 2
72 : nstep 50
73 : ipos 0.0 0.0
74 : ```
75 :
76 : This file instructs pesmd to do 50 steps of Langevin dynamics on a 2D potential energy surface
77 : at a temperature of 0.722
78 :
79 : If you run the following command a description of all the directives that can be used in the input file will be output.
80 :
81 : ```plumed
82 : plumed pesmd --help
83 : ```
84 :
85 : */
86 : //+ENDPLUMEDOC
87 :
88 : namespace PLMD {
89 : namespace cltools {
90 :
91 : class PesMD : public PLMD::CLTool {
92 7 : std::string description() const override {
93 7 : return "Langevin dynamics on PLUMED energy landscape";
94 : }
95 : public:
96 7321 : static void registerKeywords( Keywords& keys ) {
97 7321 : CLTool::registerKeywords( keys );
98 7321 : keys.add("compulsory","nstep","The number of steps of dynamics you want to run");
99 7321 : keys.add("compulsory","temperature","NVE","the temperature at which you wish to run the simulation in LJ units. With NVE there is no thermostat and the particle starts at rest");
100 7321 : keys.add("compulsory","friction","off","The relaxation time (in LJ units) of the Langevin thermostat that is used to keep the temperature constant, i.e. the inverse of the friction coefficient");
101 7321 : keys.add("compulsory","tstep","0.005","the integration timestep in LJ units");
102 7321 : keys.add("compulsory","dimension","the dimension of your energy landscape");
103 7321 : keys.add("compulsory","plumed","plumed.dat","the name of the plumed input file containing the potential");
104 7321 : keys.add("compulsory","ipos","0.0","the initial position of the system");
105 7321 : keys.add("compulsory","idum","0","The random number seed");
106 7321 : keys.addFlag("periodic",false,"are your input coordinates periodic");
107 7321 : keys.add("optional","min","minimum value the coordinates can take for a periodic domain");
108 7321 : keys.add("optional","max","maximum value the coordinates can take for a periodic domain");
109 7321 : }
110 :
111 18 : explicit PesMD( const CLToolOptions& co ) :
112 18 : CLTool(co) {
113 18 : inputdata=inputType::ifile;
114 : }
115 :
116 : private:
117 :
118 7 : void read_input(double& temperature,
119 : double& tstep,
120 : double& friction,
121 : int& dim,
122 : std::string& plumedin,
123 : std::vector<double>& ipos,
124 : int& nstep,
125 : bool& lperiod,
126 : std::vector<double>& periods,
127 : int& idum) {
128 : // Read everything from input file
129 : std::string tempstr;
130 7 : parse("temperature",tempstr);
131 16 : parse("tstep",tstep);
132 : std::string frictionstr;
133 7 : parse("friction",frictionstr);
134 7 : if( tempstr=="NVE" ) {
135 : // no thermostat, the particle starts at rest
136 1 : temperature=0.0;
137 1 : friction=-1.0; // negative: no thermostat
138 : } else {
139 6 : if(!Tools::convertNoexcept(tempstr,temperature) || temperature<0.0) {
140 0 : error("temperature should be NVE or a non negative number");
141 : }
142 6 : if(frictionstr=="off") {
143 0 : error("Specify friction for thermostat");
144 : }
145 6 : if(!Tools::convertNoexcept(frictionstr,friction) || friction<0.0) {
146 1 : error("friction should be a non negative number");
147 : }
148 : }
149 6 : parse("plumed",plumedin);
150 6 : parse("dimension",dim);
151 6 : if( dim<1 ) {
152 1 : error("dimension should be a positive integer");
153 : }
154 5 : parse("nstep",nstep);
155 5 : parse("idum",idum);
156 5 : ipos.resize( dim );
157 5 : parseVector("ipos",ipos);
158 :
159 5 : parseFlag("periodic",lperiod);
160 5 : if( lperiod ) {
161 2 : if( dim>3 ) {
162 0 : error("can only do three dimensional periodic functions");
163 : }
164 4 : std::vector<double> min( dim );
165 2 : parseVector("min",min);
166 2 : std::vector<double> max( dim );
167 2 : parseVector("max",max);
168 2 : periods.resize( dim );
169 7 : for(int i=0; i<dim; ++i) {
170 5 : if( max[i]<min[i] ) {
171 0 : error("invalid periods specified max is less than min");
172 : }
173 5 : periods[i]=max[i]-min[i];
174 : }
175 : }
176 5 : }
177 :
178 :
179 : public:
180 :
181 7 : int main( FILE* in, FILE* out, PLMD::Communicator& pc) override {
182 : std::string plumedin;
183 7 : std::vector<double> ipos;
184 : double temp, tstep, friction;
185 : bool lperiod;
186 : int dim, nsteps, seed;
187 7 : std::vector<double> periods;
188 : int plumedWantsToStop;
189 7 : Random random;
190 :
191 7 : read_input( temp, tstep, friction, dim, plumedin, ipos, nsteps, lperiod, periods, seed );
192 : // Setup random number generator
193 5 : random.setSeed(seed);
194 :
195 : // Setup box if we have periodic domain
196 7 : std::vector<double> box(9, 0.0);
197 5 : if( lperiod && dim==1 ) {
198 0 : box[0]=box[4]=box[8]=periods[0];
199 5 : } else if( lperiod && dim==2 ) {
200 1 : box[0]=periods[0];
201 1 : box[4]=box[8]=periods[1];
202 4 : } else if( lperiod && dim==3 ) {
203 1 : box[0]=periods[0];
204 1 : box[4]=periods[1];
205 1 : box[8]=periods[2];
206 3 : } else if( lperiod ) {
207 0 : error("invalid dimension for periodic potential must be 1, 2 or 3");
208 : }
209 :
210 : // Create plumed object and initialize
211 : auto plumed=Tools::make_unique<PLMD::PlumedMain>();
212 5 : int s=sizeof(double);
213 5 : plumed->cmd("setRealPrecision",&s);
214 5 : if(Communicator::initialized()) {
215 0 : plumed->cmd("setMPIComm",&pc.Get_comm());
216 : }
217 5 : int natoms=( std::floor(dim/3) + 2 );
218 5 : plumed->cmd("setNatoms",&natoms);
219 5 : plumed->cmd("setNoVirial");
220 5 : plumed->cmd("setMDEngine","pesmd");
221 5 : plumed->cmd("setTimestep",&tstep);
222 5 : plumed->cmd("setPlumedDat",plumedin.c_str());
223 5 : plumed->cmd("init");
224 :
225 : // Now create some fake atoms
226 5 : int nat = std::floor( dim/3 ) + 1;
227 5 : std::vector<double> masses( 1+nat, 1 );
228 5 : std::vector<Vector> velocities( nat ), positions( nat+1 ), forces( nat+1 );
229 : // Will set these properly eventually
230 : // Atom zero is fixed at origin
231 : positions[0].zero();
232 : // the for loop initializes two variable:
233 : // k is used only here, so its scope is limited
234 12 : for(int i=0,k=0; i<nat; ++i) {
235 28 : for(unsigned j=0; j<3; ++j) {
236 21 : if( k<dim ) {
237 12 : positions[1+i][j]=ipos[k];
238 : } else {
239 9 : positions[1+i][j]=0;
240 : }
241 21 : ++k;
242 : }
243 : }
244 : // And initialize the velocities
245 12 : for(int i=0; i<nat; ++i)
246 28 : for(int j=0; j<3; ++j) {
247 21 : velocities[i][j]=random.Gaussian() * std::sqrt( temp );
248 : }
249 : // And calculate the kinetic energy
250 : double tke=0;
251 12 : for(int i=0; i<nat; ++i) {
252 19 : for(int j=0; j<3; ++j) {
253 17 : if( 3*i+j>dim-1 ) {
254 : break;
255 : }
256 : tke += 0.5*velocities[i][j]*velocities[i][j];
257 : }
258 : }
259 :
260 : // Now call plumed to get initial forces
261 : {
262 5 : int istep=0;
263 5 : double zero=0;
264 5 : plumed->cmd("setStep",&istep);
265 5 : plumed->cmd("setMasses",&masses[0]);
266 5 : Tools::set_to_zero(forces);
267 5 : plumed->cmd("setForces",&forces[0][0]);
268 5 : plumed->cmd("setEnergy",&zero);
269 : }
270 5 : if( lperiod ) {
271 2 : plumed->cmd("setBox",&box[0]);
272 : }
273 5 : plumed->cmd("setPositions",&positions[0][0]);
274 5 : plumed->cmd("calc");
275 :
276 5 : double therm_eng=0;
277 : // Langevin thermostat on half a time step (friction is the relaxation time).
278 : // With NVE (negative friction) the velocities are left untouched, so that
279 : // therm_eng stays exactly zero (adding and subtracting the kinetic energy
280 : // can leave a tiny residue, whose sign depends e.g. on the use of FMA)
281 5 : const double lscale=(friction<0.0?1.0:std::exp(-0.5*tstep/friction));
282 5 : const double lrand=std::sqrt((1.-lscale*lscale)*temp);
283 900 : auto thermostat=[&]() {
284 900 : if( friction<0.0 ) {
285 : return;
286 : }
287 1900 : for(int j=0; j<nat; ++j) {
288 3400 : for(int k=0; k<3; ++k) {
289 2900 : if( 3*j+k>dim-1 ) {
290 : break;
291 : }
292 2200 : therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
293 2200 : velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
294 2200 : therm_eng=therm_eng-0.5*velocities[j][k]*velocities[j][k];
295 : }
296 : }
297 5 : };
298 : FILE* fp=nullptr;
299 5 : if( pc.Get_rank()==0 ) {
300 5 : fp=std::fopen("stats.out","w");
301 5 : if(!fp) {
302 0 : error("cannot open file stats.out for writing");
303 : }
304 : }
305 : // call fclose when fp goes out of scope
306 : auto deleter=[](auto f) {
307 : if(f) {
308 5 : std::fclose(f);
309 : }
310 5 : };
311 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
312 :
313 455 : for(int istep=0; istep<nsteps; ++istep) {
314 :
315 450 : if( istep%20==0 && pc.Get_rank()==0 ) {
316 : std::fprintf(out,"Doing step %i\n",istep);
317 : }
318 :
319 : // Langevin thermostat (friction is the relaxation time, a negative value means NVE)
320 450 : thermostat();
321 :
322 : // First step of velocity verlet
323 1150 : for(int j=0; j<nat; ++j) {
324 2000 : for(int k=0; k<3; ++k) {
325 1750 : if( 3*j+k>dim-1 ) {
326 : break;
327 : }
328 1300 : velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
329 1300 : positions[1+j][k] = positions[1+j][k] + tstep*velocities[j][k];
330 : }
331 : }
332 :
333 450 : int istepplusone=istep+1;
334 450 : plumedWantsToStop=0;
335 450 : plumed->cmd("setStep",&istepplusone);
336 450 : plumed->cmd("setMasses",&masses[0]);
337 450 : Tools::set_to_zero(forces);
338 450 : plumed->cmd("setForces",&forces[0][0]);
339 450 : double fenergy=0.0;
340 450 : plumed->cmd("setEnergy",&fenergy);
341 450 : plumed->cmd("setPositions",&positions[0][0]);
342 450 : plumed->cmd("setStopFlag",&plumedWantsToStop);
343 450 : plumed->cmd("calc");
344 450 : if(plumedWantsToStop) {
345 0 : nsteps=istep;
346 : }
347 :
348 : // Second step of velocity verlet
349 1150 : for(int j=0; j<nat; ++j) {
350 2000 : for(int k=0; k<3; ++k) {
351 1750 : if( 3*j+k>dim-1 ) {
352 : break;
353 : }
354 1300 : velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
355 : }
356 : }
357 :
358 : // Langevin thermostat
359 450 : thermostat();
360 : // Calculate total kinetic energy
361 : tke=0;
362 1150 : for(int i=0; i<nat; ++i) {
363 2000 : for(int j=0; j<3; ++j) {
364 1750 : if( 3*i+j>dim-1 ) {
365 : break;
366 : }
367 1300 : tke += 0.5*velocities[i][j]*velocities[i][j];
368 : }
369 : }
370 :
371 : // Print everything
372 : // conserved = potential+1.5*ttt+therm_eng;
373 450 : if( pc.Get_rank()==0 ) {
374 450 : std::fprintf(fp,"%i %f %f %f \n", istep, istep*tstep, tke, therm_eng );
375 : }
376 : }
377 :
378 5 : plumed->cmd("runFinalJobs");
379 5 : return 0;
380 5 : }
381 : };
382 :
383 21981 : PLUMED_REGISTER_CLTOOL(PesMD,"pesmd")
384 :
385 : }
386 : }
|
