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