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