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