Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 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 "wrapper/Plumed.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 :
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 : using namespace std;
96 :
97 : namespace PLMD {
98 : namespace cltools {
99 :
100 1 : class PesMD : public PLMD::CLTool {
101 0 : string description() const {
102 0 : return "langevin dynamics on PLUMED energy landscape";
103 : }
104 : public:
105 1613 : static void registerKeywords( Keywords& keys ) {
106 6452 : keys.add("compulsory","nstep","The number of steps of dynamics you want to run");
107 8065 : keys.add("compulsory","temperature","NVE","the temperature at which you wish to run the simulation in LJ units");
108 8065 : keys.add("compulsory","friction","off","The friction (in LJ units) for the langevin thermostat that is used to keep the temperature constant");
109 8065 : keys.add("compulsory","tstep","0.005","the integration timestep in LJ units");
110 6452 : keys.add("compulsory","dimension","the dimension of your energy landscape");
111 8065 : keys.add("compulsory","plumed","plumed.dat","the name of the plumed input file containing the potential");
112 8065 : keys.add("compulsory","ipos","0.0","the initial position of the system");
113 8065 : keys.add("compulsory","idum","0","The random number seed");
114 4839 : keys.addFlag("periodic","false","are your input coordinates periodic");
115 6452 : keys.add("optional","min","minimum value the coordinates can take for a periodic domain");
116 6452 : keys.add("optional","max","maximum value the coordinates can take for a periodic domain");
117 1613 : }
118 :
119 1 : explicit PesMD( const CLToolOptions& co ) :
120 1 : CLTool(co)
121 : {
122 1 : inputdata=ifile;
123 : }
124 :
125 : private:
126 :
127 1 : void read_input(double& temperature,
128 : double& tstep,
129 : double& friction,
130 : int& dim,
131 : std::string& plumedin,
132 : std::vector<double>& ipos,
133 : int& nstep,
134 : bool& lperiod,
135 : std::vector<double>& periods,
136 : int& idum)
137 : {
138 : // Read everything from input file
139 2 : std::string tempstr; parse("temperature",tempstr);
140 1 : if( tempstr!="NVE" ) Tools::convert(tempstr,temperature);
141 2 : parse("tstep",tstep);
142 2 : std::string frictionstr; parse("friction",frictionstr);
143 1 : if( tempstr!="NVE" ) {
144 1 : if(frictionstr=="off") { fprintf(stderr,"Specify friction for thermostat\n"); exit(1); }
145 1 : Tools::convert(frictionstr,friction);
146 : }
147 3 : parse("plumed",plumedin); parse("dimension",dim);
148 3 : parse("nstep",nstep); parse("idum",idum);
149 2 : ipos.resize( dim ); parseVector("ipos",ipos);
150 :
151 2 : parseFlag("periodic",lperiod);
152 1 : if( lperiod ) {
153 0 : if( dim>3 ) error("can only do three dimensional periodic functions");
154 0 : std::vector<double> min( dim ); parseVector("min",min);
155 0 : std::vector<double> max( dim ); parseVector("max",max);
156 0 : periods.resize( dim );
157 0 : for(unsigned i=0; i<dim; ++i) {
158 0 : if( max[i]<min[i] ) error("invalid periods specified max is less than min");
159 0 : periods[i]=max[i]-min[i];
160 : }
161 : }
162 1 : }
163 :
164 :
165 : public:
166 :
167 1 : int main( FILE* in, FILE* out, PLMD::Communicator& pc) {
168 : std::string plumedin; std::vector<double> ipos;
169 : double temp, tstep, friction; bool lperiod;
170 : int dim, nsteps, seed; std::vector<double> periods;
171 : int plumedWantsToStop;
172 1 : Random random;
173 :
174 1 : read_input( temp, tstep, friction, dim, plumedin, ipos, nsteps, lperiod, periods, seed );
175 : // Setup random number generator
176 1 : random.setSeed(seed);
177 :
178 : // Setup box if we have periodic domain
179 1 : std::vector<double> box(9, 0.0);
180 1 : if( lperiod && dim==1 ) { box[0]=box[5]=box[9]=periods[0]; }
181 1 : else if( lperiod && dim==2 ) { box[0]=periods[0]; box[5]=box[9]=periods[1]; }
182 1 : else if( lperiod && dim==3 ) { box[0]=periods[0]; box[5]=periods[1]; box[9]=periods[2]; }
183 1 : else if( lperiod ) error("invalid dimension for periodic potential must be 1, 2 or 3");
184 :
185 : // Create plumed object and initialize
186 2 : PLMD::Plumed* plumed=new PLMD::Plumed; int s=sizeof(double);
187 : plumed->cmd("setRealPrecision",&s);
188 1 : if(Communicator::initialized()) plumed->cmd("setMPIComm",&pc.Get_comm());
189 : plumed->cmd("setNoVirial");
190 2 : int natoms=( std::floor(dim/3) + 2 );
191 : plumed->cmd("setNatoms",&natoms);
192 : plumed->cmd("setMDEngine","pesmd");
193 : plumed->cmd("setTimestep",&tstep);
194 : plumed->cmd("setPlumedDat",plumedin.c_str());
195 : plumed->cmd("init");
196 :
197 : // Now create some fake atoms
198 2 : unsigned nat = std::floor( dim/3 ) + 1;
199 1 : std::vector<double> masses( 1+nat, 1 );
200 1 : std::vector<Vector> velocities( nat ), positions( nat+1 ), forces( nat+1 );
201 : // Will set these properly eventually
202 1 : unsigned k=0; positions[0].zero(); // Atom zero is fixed at origin
203 5 : for(unsigned i=0; i<nat; ++i) for(unsigned j=0; j<3; ++j) {
204 7 : if( k<dim ) { positions[1+i][j]=ipos[k]; } else { positions[1+i][j]=0;}
205 3 : k++;
206 : }
207 : // And initialize the velocities
208 5 : for(unsigned i=0; i<nat; ++i) for(unsigned j=0; j<3; ++j) velocities[i][j]=random.Gaussian() * sqrt( temp );
209 : // And calcualte the kinetic energy
210 : double tke=0;
211 3 : for(unsigned i=0; i<nat; ++i) {
212 5 : for(unsigned j=0; j<3; ++j) {
213 3 : if( 3*i+j>dim ) break;
214 : tke += 0.5*velocities[i][j]*velocities[i][j];
215 : }
216 : }
217 :
218 : // Now call plumed to get initial forces
219 1 : int istep=0; double zero=0;
220 : plumed->cmd("setStep",&istep);
221 : plumed->cmd("setMasses",&masses[0]);
222 8 : for(unsigned i=0; i<forces.size(); ++i) forces[i].zero();
223 : plumed->cmd("setForces",&forces[0]);
224 : plumed->cmd("setEnergy",&zero);
225 1 : if( lperiod ) plumed->cmd("setBox",&box[0]);
226 : plumed->cmd("setPositions",&positions[0]);
227 : plumed->cmd("calc");
228 :
229 :
230 : // potential=calc_energy(positions,forces);
231 : double therm_eng=0;
232 :
233 1 : FILE* fp=fopen("stats.out","w+");
234 : // double conserved = potential+1.5*ttt+therm_eng; FILE* fp=fopen("stats.out","w+");
235 : // if( pc.Get_rank()==0 ) fprintf(fp,"%d %f %f \n", 0, 0., tke, therm_eng );
236 :
237 101 : for(unsigned istep=0; istep<nsteps; ++istep) {
238 :
239 50 : if( istep%20==0 && pc.Get_rank()==0 ) printf("Doing step %u\n",istep);
240 :
241 : // Langevin thermostat
242 50 : double lscale=exp(-0.5*tstep/friction);
243 50 : double lrand=sqrt((1.-lscale*lscale)*temp);
244 150 : for(unsigned j=0; j<nat; ++j) {
245 150 : for(unsigned k=0; k<3; ++k) {
246 100 : if( 3*j+k>dim-1 ) break;
247 100 : therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
248 100 : velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
249 50 : therm_eng=therm_eng-0.5*velocities[j][k]*velocities[0][k];
250 : }
251 : }
252 :
253 : // First step of velocity verlet
254 150 : for(unsigned j=0; j<nat; ++j) {
255 150 : for(unsigned k=0; k<3; ++k) {
256 100 : if( 3*j+k>dim-1 ) break;
257 150 : velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
258 100 : positions[1+j][k] = positions[1+j][k] + tstep*velocities[j][k];
259 : // Apply pbc
260 : // if( positions[0][k]>pi ) positions[0][k]-=2*pi;
261 : // if( positions[0][k]<=-pi ) positions[0][k]+=2*pi;
262 : }
263 : }
264 :
265 50 : int istepplusone=istep+1;
266 50 : plumedWantsToStop=0;
267 : plumed->cmd("setStep",&istepplusone);
268 : plumed->cmd("setMasses",&masses[0]);
269 400 : for(unsigned i=0; i<forces.size(); ++i) forces[i].zero();
270 : plumed->cmd("setForces",&forces[0]);
271 50 : double fenergy=0.0;
272 : plumed->cmd("setEnergy",&fenergy);
273 : plumed->cmd("setPositions",&positions[0]);
274 : plumed->cmd("setStopFlag",&plumedWantsToStop);
275 : plumed->cmd("calc");
276 : // if(istep%2000==0) plumed->cmd("writeCheckPointFile");
277 50 : if(plumedWantsToStop) nsteps=istep;
278 :
279 : // Second step of velocity verlet
280 150 : for(unsigned j=0; j<nat; ++j) {
281 150 : for(unsigned k=0; k<3; ++k) {
282 100 : if( 3*j+k>dim-1 ) break;
283 150 : velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
284 : }
285 : }
286 :
287 : // Langevin thermostat
288 50 : lscale=exp(-0.5*tstep/friction);
289 50 : lrand=sqrt((1.-lscale*lscale)*temp);
290 150 : for(unsigned j=0; j<nat; ++j) {
291 150 : for(unsigned k=0; k<3; ++k) {
292 100 : if( 3*j+k>dim-1 ) break;
293 100 : therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
294 100 : velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
295 50 : therm_eng=therm_eng-0.5*velocities[j][k]*velocities[j][k];
296 : }
297 : }
298 : // Calculate total kinetic energy
299 : tke=0;
300 150 : for(unsigned i=0; i<nat; ++i) {
301 150 : for(unsigned j=0; j<3; ++j) {
302 100 : if( 3*i+j>dim-1 ) break;
303 100 : tke += 0.5*velocities[i][j]*velocities[i][j];
304 : }
305 : }
306 :
307 : // Print everything
308 : // conserved = potential+1.5*ttt+therm_eng;
309 50 : if( pc.Get_rank()==0 ) fprintf(fp,"%u %f %f %f \n", istep, istep*tstep, tke, therm_eng );
310 : }
311 :
312 1 : delete plumed;
313 1 : fclose(fp);
314 :
315 1 : return 0;
316 : }
317 : };
318 :
319 6454 : PLUMED_REGISTER_CLTOOL(PesMD,"pesmd")
320 :
321 : }
322 4839 : }
|
