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 5 : std::string description() const override {
93 5 : return "Langevin dynamics on PLUMED energy landscape";
94 : }
95 : public:
96 5683 : static void registerKeywords( Keywords& keys ) {
97 5683 : CLTool::registerKeywords( keys );
98 5683 : keys.add("compulsory","nstep","The number of steps of dynamics you want to run");
99 5683 : keys.add("compulsory","temperature","NVE","the temperature at which you wish to run the simulation in LJ units");
100 5683 : keys.add("compulsory","friction","off","The friction (in LJ units) for the Langevin thermostat that is used to keep the temperature constant");
101 5683 : keys.add("compulsory","tstep","0.005","the integration timestep in LJ units");
102 5683 : keys.add("compulsory","dimension","the dimension of your energy landscape");
103 5683 : keys.add("compulsory","plumed","plumed.dat","the name of the plumed input file containing the potential");
104 5683 : keys.add("compulsory","ipos","0.0","the initial position of the system");
105 5683 : keys.add("compulsory","idum","0","The random number seed");
106 5683 : keys.addFlag("periodic",false,"are your input coordinates periodic");
107 5683 : keys.add("optional","min","minimum value the coordinates can take for a periodic domain");
108 5683 : keys.add("optional","max","maximum value the coordinates can take for a periodic domain");
109 5683 : }
110 :
111 8 : explicit PesMD( const CLToolOptions& co ) :
112 8 : CLTool(co) {
113 8 : inputdata=inputType::ifile;
114 : }
115 :
116 : private:
117 :
118 3 : 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 6 : parse("temperature",tempstr);
131 3 : if( tempstr!="NVE" ) {
132 3 : Tools::convert(tempstr,temperature);
133 : }
134 6 : parse("tstep",tstep);
135 : std::string frictionstr;
136 6 : parse("friction",frictionstr);
137 3 : if( tempstr!="NVE" ) {
138 3 : if(frictionstr=="off") {
139 0 : error("pecify friction for thermostat");
140 : }
141 3 : Tools::convert(frictionstr,friction);
142 : }
143 3 : parse("plumed",plumedin);
144 3 : parse("dimension",dim);
145 3 : parse("nstep",nstep);
146 3 : parse("idum",idum);
147 3 : ipos.resize( dim );
148 3 : parseVector("ipos",ipos);
149 :
150 3 : parseFlag("periodic",lperiod);
151 3 : if( lperiod ) {
152 2 : if( dim>3 ) {
153 0 : error("can only do three dimensional periodic functions");
154 : }
155 2 : std::vector<double> min( dim );
156 2 : parseVector("min",min);
157 2 : std::vector<double> max( dim );
158 2 : parseVector("max",max);
159 2 : periods.resize( dim );
160 7 : for(int i=0; i<dim; ++i) {
161 5 : if( max[i]<min[i] ) {
162 0 : error("invalid periods specified max is less than min");
163 : }
164 5 : periods[i]=max[i]-min[i];
165 : }
166 : }
167 3 : }
168 :
169 :
170 : public:
171 :
172 3 : int main( FILE* in, FILE* out, PLMD::Communicator& pc) override {
173 : std::string plumedin;
174 : std::vector<double> ipos;
175 : double temp, tstep, friction;
176 : bool lperiod;
177 : int dim, nsteps, seed;
178 : std::vector<double> periods;
179 : int plumedWantsToStop;
180 3 : Random random;
181 :
182 3 : read_input( temp, tstep, friction, dim, plumedin, ipos, nsteps, lperiod, periods, seed );
183 : // Setup random number generator
184 3 : random.setSeed(seed);
185 :
186 : // Setup box if we have periodic domain
187 3 : std::vector<double> box(9, 0.0);
188 3 : if( lperiod && dim==1 ) {
189 0 : box[0]=box[4]=box[8]=periods[0];
190 3 : } else if( lperiod && dim==2 ) {
191 1 : box[0]=periods[0];
192 1 : box[4]=box[8]=periods[1];
193 2 : } else if( lperiod && dim==3 ) {
194 1 : box[0]=periods[0];
195 1 : box[4]=periods[1];
196 1 : box[8]=periods[2];
197 1 : } else if( lperiod ) {
198 0 : error("invalid dimension for periodic potential must be 1, 2 or 3");
199 : }
200 :
201 : // Create plumed object and initialize
202 : auto plumed=Tools::make_unique<PLMD::PlumedMain>();
203 3 : int s=sizeof(double);
204 3 : plumed->cmd("setRealPrecision",&s);
205 3 : if(Communicator::initialized()) {
206 0 : plumed->cmd("setMPIComm",&pc.Get_comm());
207 : }
208 3 : int natoms=( std::floor(dim/3) + 2 );
209 3 : plumed->cmd("setNatoms",&natoms);
210 3 : plumed->cmd("setNoVirial");
211 3 : plumed->cmd("setMDEngine","pesmd");
212 3 : plumed->cmd("setTimestep",&tstep);
213 3 : plumed->cmd("setPlumedDat",plumedin.c_str());
214 3 : plumed->cmd("init");
215 :
216 : // Now create some fake atoms
217 3 : int nat = std::floor( dim/3 ) + 1;
218 3 : std::vector<double> masses( 1+nat, 1 );
219 3 : std::vector<Vector> velocities( nat ), positions( nat+1 ), forces( nat+1 );
220 : // Will set these properly eventually
221 : // Atom zero is fixed at origin
222 : positions[0].zero();
223 : // the for loop initializes two variable:
224 : // k is used only here, so its scope is limited
225 7 : for(int i=0,k=0; i<nat; ++i) {
226 16 : for(unsigned j=0; j<3; ++j) {
227 12 : if( k<dim ) {
228 6 : positions[1+i][j]=ipos[k];
229 : } else {
230 6 : positions[1+i][j]=0;
231 : }
232 12 : ++k;
233 : }
234 : }
235 : // And initialize the velocities
236 7 : for(int i=0; i<nat; ++i)
237 16 : for(int j=0; j<3; ++j) {
238 12 : velocities[i][j]=random.Gaussian() * std::sqrt( temp );
239 : }
240 : // And calculate the kinetic energy
241 : double tke=0;
242 7 : for(int i=0; i<nat; ++i) {
243 10 : for(int j=0; j<3; ++j) {
244 9 : if( 3*i+j>dim-1 ) {
245 : break;
246 : }
247 : tke += 0.5*velocities[i][j]*velocities[i][j];
248 : }
249 : }
250 :
251 : // Now call plumed to get initial forces
252 : {
253 3 : int istep=0;
254 3 : double zero=0;
255 3 : plumed->cmd("setStep",&istep);
256 3 : plumed->cmd("setMasses",&masses[0]);
257 3 : Tools::set_to_zero(forces);
258 3 : plumed->cmd("setForces",&forces[0][0]);
259 3 : plumed->cmd("setEnergy",&zero);
260 : }
261 3 : if( lperiod ) {
262 2 : plumed->cmd("setBox",&box[0]);
263 : }
264 3 : plumed->cmd("setPositions",&positions[0][0]);
265 3 : plumed->cmd("calc");
266 :
267 : double therm_eng=0;
268 3 : FILE* fp=fopen("stats.out","w+");
269 :
270 153 : for(int istep=0; istep<nsteps; ++istep) {
271 :
272 150 : if( istep%20==0 && pc.Get_rank()==0 ) {
273 : std::printf("Doing step %i\n",istep);
274 : }
275 :
276 : // Langevin thermostat
277 150 : double lscale=std::exp(-0.5*tstep/friction);
278 150 : double lrand=std::sqrt((1.-lscale*lscale)*temp);
279 350 : for(int j=0; j<nat; ++j) {
280 500 : for(int k=0; k<3; ++k) {
281 450 : if( 3*j+k>dim-1 ) {
282 : break;
283 : }
284 300 : therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
285 300 : velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
286 300 : therm_eng=therm_eng-0.5*velocities[j][k]*velocities[0][k];
287 : }
288 : }
289 :
290 : // First step of velocity verlet
291 350 : for(int j=0; j<nat; ++j) {
292 500 : for(int k=0; k<3; ++k) {
293 450 : if( 3*j+k>dim-1 ) {
294 : break;
295 : }
296 300 : velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
297 300 : positions[1+j][k] = positions[1+j][k] + tstep*velocities[j][k];
298 : }
299 : }
300 :
301 150 : int istepplusone=istep+1;
302 150 : plumedWantsToStop=0;
303 150 : plumed->cmd("setStep",&istepplusone);
304 150 : plumed->cmd("setMasses",&masses[0]);
305 150 : Tools::set_to_zero(forces);
306 150 : plumed->cmd("setForces",&forces[0][0]);
307 150 : double fenergy=0.0;
308 150 : plumed->cmd("setEnergy",&fenergy);
309 150 : plumed->cmd("setPositions",&positions[0][0]);
310 150 : plumed->cmd("setStopFlag",&plumedWantsToStop);
311 150 : plumed->cmd("calc");
312 : // if(istep%2000==0) plumed->cmd("writeCheckPointFile");
313 150 : if(plumedWantsToStop) {
314 0 : nsteps=istep;
315 : }
316 :
317 : // Second step of velocity verlet
318 350 : for(int j=0; j<nat; ++j) {
319 500 : for(int k=0; k<3; ++k) {
320 450 : if( 3*j+k>dim-1 ) {
321 : break;
322 : }
323 300 : velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
324 : }
325 : }
326 :
327 : // Langevin thermostat
328 150 : lscale=std::exp(-0.5*tstep/friction);
329 150 : lrand=std::sqrt((1.-lscale*lscale)*temp);
330 350 : for(int j=0; j<nat; ++j) {
331 500 : for(int k=0; k<3; ++k) {
332 450 : if( 3*j+k>dim-1) {
333 : break;
334 : }
335 300 : therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
336 300 : velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
337 300 : therm_eng=therm_eng-0.5*velocities[j][k]*velocities[j][k];
338 : }
339 : }
340 : // Calculate total kinetic energy
341 : tke=0;
342 350 : for(int i=0; i<nat; ++i) {
343 500 : for(int j=0; j<3; ++j) {
344 450 : if( 3*i+j>dim-1 ) {
345 : break;
346 : }
347 300 : tke += 0.5*velocities[i][j]*velocities[i][j];
348 : }
349 : }
350 :
351 : // Print everything
352 : // conserved = potential+1.5*ttt+therm_eng;
353 150 : if( pc.Get_rank()==0 ) {
354 150 : std::fprintf(fp,"%i %f %f %f \n", istep, istep*tstep, tke, therm_eng );
355 : }
356 : }
357 :
358 3 : fclose(fp);
359 :
360 3 : return 0;
361 3 : }
362 : };
363 :
364 17057 : PLUMED_REGISTER_CLTOOL(PesMD,"pesmd")
365 :
366 : }
367 : }
|
