Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 <string>
28 : #include <cstdio>
29 : #include <cmath>
30 : #include <vector>
31 :
32 : using namespace std;
33 :
34 : namespace PLMD {
35 : namespace cltools {
36 :
37 : //+PLUMEDOC TOOLS simplemd
38 : /*
39 : simplemd allows one to do molecular dynamics on systems of Lennard-Jones atoms.
40 :
41 : The input to simplemd is spcified in an input file. Configurations are input and
42 : output in xyz format. The input file should contain one directive per line.
43 : The directives available are as follows:
44 :
45 : \par Examples
46 :
47 : You run an MD simulation using simplemd with the following command:
48 : \verbatim
49 : plumed simplemd < in
50 : \endverbatim
51 :
52 : The following is an example of an input file for a simplemd calculation. This file
53 : instructs simplemd to do 50 steps of MD at a temperature of 0.722
54 : \verbatim
55 : nputfile input.xyz
56 : outputfile output.xyz
57 : temperature 0.722
58 : tstep 0.005
59 : friction 1
60 : forcecutoff 2.5
61 : listcutoff 3.0
62 : nstep 50
63 : nconfig 10 trajectory.xyz
64 : nstat 10 energies.dat
65 : \endverbatim
66 :
67 : If you run the following a description of all the directives that can be used in the
68 : input file will be output.
69 : \verbatim
70 : plumed simplemd --help
71 : \endverbatim
72 :
73 : */
74 : //+ENDPLUMEDOC
75 :
76 5 : class SimpleMD:
77 : public PLMD::CLTool
78 : {
79 0 : string description()const {
80 0 : return "run lj code";
81 : }
82 :
83 : bool write_positions_first;
84 : bool write_statistics_first;
85 : int write_statistics_last_time_reopened;
86 : FILE* write_statistics_fp;
87 :
88 :
89 : public:
90 1613 : static void registerKeywords( Keywords& keys ) {
91 6452 : keys.add("compulsory","nstep","The number of steps of dynamics you want to run");
92 8065 : keys.add("compulsory","temperature","NVE","the temperature at which you wish to run the simulation in LJ units");
93 8065 : keys.add("compulsory","friction","off","The friction (in LJ units) for the langevin thermostat that is used to keep the temperature constant");
94 8065 : keys.add("compulsory","tstep","0.005","the integration timestep in LJ units");
95 6452 : keys.add("compulsory","inputfile","An xyz file containing the initial configuration of the system");
96 8065 : keys.add("compulsory","forcecutoff","2.5","");
97 8065 : keys.add("compulsory","listcutoff","3.0","");
98 6452 : keys.add("compulsory","outputfile","An output xyz file containing the final configuration of the system");
99 8065 : keys.add("compulsory","nconfig","10","The frequency with which to write configurations to the trajectory file followed by the name of the trajectory file");
100 8065 : keys.add("compulsory","nstat","1","The frequency with which to write the statistics to the statistics file followed by the name of the statistics file");
101 8065 : keys.add("compulsory","maxneighbours","10000","The maximum number of neighbours an atom can have");
102 8065 : keys.add("compulsory","idum","0","The random number seed");
103 8065 : keys.add("compulsory","ndim","3","The dimensionality of the system (some interesting LJ clusters are two dimensional)");
104 8065 : keys.add("compulsory","wrapatoms","false","If true, atomic coordinates are written wrapped in minimal cell");
105 1613 : }
106 :
107 5 : explicit SimpleMD( const CLToolOptions& co ) :
108 : CLTool(co),
109 : write_positions_first(true),
110 : write_statistics_first(true),
111 : write_statistics_last_time_reopened(0),
112 5 : write_statistics_fp(NULL)
113 : {
114 5 : inputdata=ifile;
115 : }
116 :
117 : private:
118 :
119 : void
120 5 : read_input(double& temperature,
121 : double& tstep,
122 : double& friction,
123 : double& forcecutoff,
124 : double& listcutoff,
125 : int& nstep,
126 : int& nconfig,
127 : int& nstat,
128 : bool& wrapatoms,
129 : string& inputfile,
130 : string& outputfile,
131 : string& trajfile,
132 : string& statfile,
133 : int& maxneighbours,
134 : int& ndim,
135 : int& idum)
136 : {
137 :
138 : // Read everything from input file
139 : char buffer1[256];
140 10 : std::string tempstr; parse("temperature",tempstr);
141 5 : if( tempstr!="NVE" ) Tools::convert(tempstr,temperature);
142 10 : parse("tstep",tstep);
143 10 : std::string frictionstr; parse("friction",frictionstr);
144 5 : if( tempstr!="NVE" ) {
145 5 : if(frictionstr=="off") { fprintf(stderr,"Specify friction for thermostat\n"); exit(1); }
146 5 : Tools::convert(frictionstr,friction);
147 : }
148 10 : parse("forcecutoff",forcecutoff);
149 10 : parse("listcutoff",listcutoff);
150 10 : parse("nstep",nstep);
151 10 : parse("maxneighbours",maxneighbours);
152 10 : parse("idum",idum);
153 :
154 : // Read in stuff with sanity checks
155 10 : parse("inputfile",inputfile);
156 5 : if(inputfile.length()==0) {
157 0 : fprintf(stderr,"Specify input file\n");
158 0 : exit(1);
159 : }
160 10 : parse("outputfile",outputfile);
161 5 : if(outputfile.length()==0) {
162 0 : fprintf(stderr,"Specify output file\n");
163 0 : exit(1);
164 : }
165 10 : std::string nconfstr; parse("nconfig",nconfstr);
166 5 : sscanf(nconfstr.c_str(),"%100d %255s",&nconfig,buffer1);
167 : trajfile=buffer1;
168 5 : if(trajfile.length()==0) {
169 0 : fprintf(stderr,"Specify traj file\n");
170 0 : exit(1);
171 : }
172 10 : std::string nstatstr; parse("nstat",nstatstr);
173 5 : sscanf(nstatstr.c_str(),"%100d %255s",&nstat,buffer1);
174 : statfile=buffer1;
175 5 : if(statfile.length()==0) {
176 0 : fprintf(stderr,"Specify stat file\n");
177 0 : exit(1);
178 : }
179 10 : parse("ndim",ndim);
180 5 : if(ndim<1 || ndim>3) {
181 0 : fprintf(stderr,"ndim should be 1,2 or 3\n");
182 0 : exit(1);
183 : }
184 : std::string w;
185 10 : parse("wrapatoms",w);
186 5 : wrapatoms=false;
187 10 : if(w.length()>0 && (w[0]=='T' || w[0]=='t')) wrapatoms=true;
188 5 : }
189 :
190 5 : void read_natoms(const string & inputfile,int & natoms) {
191 : // read the number of atoms in file "input.xyz"
192 5 : FILE* fp=fopen(inputfile.c_str(),"r");
193 5 : if(!fp) {
194 0 : fprintf(stderr,"ERROR: file %s not found\n",inputfile.c_str());
195 0 : exit(1);
196 : }
197 5 : fscanf(fp,"%1000d",&natoms);
198 5 : fclose(fp);
199 5 : }
200 :
201 5 : void read_positions(const string& inputfile,int natoms,vector<Vector>& positions,double cell[3]) {
202 : // read positions and cell from a file called inputfile
203 : // natoms (input variable) and number of atoms in the file should be consistent
204 5 : FILE* fp=fopen(inputfile.c_str(),"r");
205 5 : if(!fp) {
206 0 : fprintf(stderr,"ERROR: file %s not found\n",inputfile.c_str());
207 0 : exit(1);
208 : }
209 : char buffer[256];
210 : char atomname[256];
211 : fgets(buffer,256,fp);
212 5 : fscanf(fp,"%1000lf %1000lf %1000lf",&cell[0],&cell[1],&cell[2]);
213 883 : for(int i=0; i<natoms; i++) {
214 878 : fscanf(fp,"%255s %1000lf %1000lf %1000lf",atomname,&positions[i][0],&positions[i][1],&positions[i][2]);
215 : // note: atomname is read but not used
216 : }
217 5 : fclose(fp);
218 5 : }
219 :
220 5 : void randomize_velocities(const int natoms,const int ndim,const double temperature,const vector<double>&masses,vector<Vector>& velocities,Random&random) {
221 : // randomize the velocities according to the temperature
222 1754 : for(int iatom=0; iatom<natoms; iatom++) for(int i=0; i<ndim; i++)
223 3930 : velocities[iatom][i]=sqrt(temperature/masses[iatom])*random.Gaussian();
224 5 : }
225 :
226 957956 : void pbc(const double cell[3],const Vector & vin,Vector & vout) {
227 : // apply periodic boundary condition to a vector
228 6705692 : for(int i=0; i<3; i++) {
229 2873868 : vout[i]=vin[i]-floor(vin[i]/cell[i]+0.5)*cell[i];
230 : }
231 957956 : }
232 :
233 2200 : void check_list(const int natoms,const vector<Vector>& positions,const vector<Vector>&positions0,const double listcutoff,
234 : const double forcecutoff,bool & recompute)
235 : {
236 : // check if the neighbour list have to be recomputed
237 2200 : Vector displacement; // displacement from positions0 to positions
238 : double delta2; // square of the 'skin' thickness
239 2200 : recompute=false;
240 2200 : delta2=(0.5*(listcutoff-forcecutoff))*(0.5*(listcutoff-forcecutoff));
241 : // if ANY atom moved more than half of the skin thickness, recompute is set to .true.
242 73400 : for(int iatom=0; iatom<natoms; iatom++) {
243 356000 : for(int k=0; k<3; k++) displacement[k]=positions[iatom][k]-positions0[iatom][k];
244 : double s=0.0;
245 142400 : for(int k=0; k<3; k++) s+=displacement[k]*displacement[k];
246 35600 : if(s>delta2) recompute=true;
247 : }
248 2200 : }
249 :
250 :
251 37 : void compute_list(const int natoms,const int listsize,const vector<Vector>& positions,const double cell[3],const double listcutoff,
252 : vector<int>& point,vector<int>& list) {
253 : // see Allen-Tildesey for a definition of point and list
254 37 : Vector distance; // distance of the two atoms
255 37 : Vector distance_pbc; // minimum-image distance of the two atoms
256 : double listcutoff2; // squared list cutoff
257 37 : listcutoff2=listcutoff*listcutoff;
258 37 : point[0]=0;
259 1572 : for(int iatom=0; iatom<natoms-1; iatom++) {
260 4605 : point[iatom+1]=point[iatom];
261 152771 : for(int jatom=iatom+1; jatom<natoms; jatom++) {
262 756180 : for(int k=0; k<3; k++) distance[k]=positions[iatom][k]-positions[jatom][k];
263 75618 : pbc(cell,distance,distance_pbc);
264 : // if the interparticle distance is larger than the cutoff, skip
265 302472 : double d2=0; for(int k=0; k<3; k++) d2+=distance_pbc[k]*distance_pbc[k];
266 75618 : if(d2>listcutoff2)continue;
267 54532 : if(point[iatom+1]>listsize) {
268 : // too many neighbours
269 0 : fprintf(stderr,"%s","Verlet list size exceeded\n");
270 0 : fprintf(stderr,"%s","Increase maxneighbours\n");
271 0 : exit(1);
272 : }
273 109064 : list[point[iatom+1]]=jatom;
274 54532 : point[iatom+1]++;
275 : }
276 : }
277 37 : }
278 :
279 2205 : void compute_forces(const int natoms,const int listsize,const vector<Vector>& positions,const double cell[3],
280 : double forcecutoff,const vector<int>& point,const vector<int>& list,vector<Vector>& forces,double & engconf)
281 : {
282 2205 : Vector distance; // distance of the two atoms
283 2205 : Vector distance_pbc; // minimum-image distance of the two atoms
284 : double distance_pbc2; // squared minimum-image distance
285 : double forcecutoff2; // squared force cutoff
286 2205 : Vector f; // force
287 : double engcorrection; // energy necessary shift the potential avoiding discontinuities
288 :
289 2205 : forcecutoff2=forcecutoff*forcecutoff;
290 2205 : engconf=0.0;
291 146361 : for(int i=0; i<natoms; i++)for(int k=0; k<3; k++) forces[i][k]=0.0;
292 4410 : engcorrection=4.0*(1.0/pow(forcecutoff2,6.0)-1.0/pow(forcecutoff2,3));
293 36039 : for(int iatom=0; iatom<natoms-1; iatom++) {
294 1863586 : for(int jlist=point[iatom]; jlist<point[iatom+1]; jlist++) {
295 1762084 : int jatom=list[jlist];
296 8810420 : for(int k=0; k<3; k++) distance[k]=positions[iatom][k]-positions[jatom][k];
297 881042 : pbc(cell,distance,distance_pbc);
298 3524168 : distance_pbc2=0.0; for(int k=0; k<3; k++) distance_pbc2+=distance_pbc[k]*distance_pbc[k];
299 : // if the interparticle distance is larger than the cutoff, skip
300 881042 : if(distance_pbc2>forcecutoff2) continue;
301 641298 : double distance_pbc6=distance_pbc2*distance_pbc2*distance_pbc2;
302 641298 : double distance_pbc8=distance_pbc6*distance_pbc2;
303 641298 : double distance_pbc12=distance_pbc6*distance_pbc6;
304 641298 : double distance_pbc14=distance_pbc12*distance_pbc2;
305 641298 : engconf+=4.0*(1.0/distance_pbc12 - 1.0/distance_pbc6) - engcorrection;
306 2565192 : for(int k=0; k<3; k++) f[k]=2.0*distance_pbc[k]*4.0*(6.0/distance_pbc14-3.0/distance_pbc8);
307 : // same force on the two atoms, with opposite sign:
308 4489086 : for(int k=0; k<3; k++) forces[iatom][k]+=f[k];
309 4489086 : for(int k=0; k<3; k++) forces[jatom][k]-=f[k];
310 : }
311 : }
312 2205 : }
313 :
314 2200 : void compute_engkin(const int natoms,const vector<double>& masses,const vector<Vector>& velocities,double & engkin)
315 : {
316 : // calculate the kinetic energy from the velocities
317 2200 : engkin=0.0;
318 144600 : for(int iatom=0; iatom<natoms; iatom++)for(int k=0; k<3; k++) {
319 320400 : engkin+=0.5*masses[iatom]*velocities[iatom][k]*velocities[iatom][k];
320 : }
321 2200 : }
322 :
323 :
324 4400 : void thermostat(const int natoms,const int ndim,const vector<double>& masses,const double dt,const double friction,
325 : const double temperature,vector<Vector>& velocities,double & engint,Random & random) {
326 : // Langevin thermostat, implemented as decribed in Bussi and Parrinello, Phys. Rev. E (2007)
327 : // it is a linear combination of old velocities and new, randomly chosen, velocity,
328 : // with proper coefficients
329 4400 : double c1=exp(-friction*dt);
330 146800 : for(int iatom=0; iatom<natoms; iatom++) {
331 142400 : double c2=sqrt((1.0-c1*c1)*temperature/masses[iatom]);
332 442400 : for(int i=0; i<ndim; i++) {
333 371200 : engint+=0.5*masses[iatom]*velocities[iatom][i]*velocities[iatom][i];
334 371200 : velocities[iatom][i]=c1*velocities[iatom][i]+c2*random.Gaussian();
335 371200 : engint-=0.5*masses[iatom]*velocities[iatom][i]*velocities[iatom][i];
336 : }
337 : }
338 4400 : }
339 :
340 22 : void write_positions(const string& trajfile,int natoms,const vector<Vector>& positions,const double cell[3],const bool wrapatoms)
341 : {
342 : // write positions on file trajfile
343 : // positions are appended at the end of the file
344 22 : Vector pos;
345 : FILE*fp;
346 22 : if(write_positions_first) {
347 5 : fp=fopen(trajfile.c_str(),"w");
348 5 : write_positions_first=false;
349 : } else {
350 17 : fp=fopen(trajfile.c_str(),"a");
351 : }
352 : fprintf(fp,"%d\n",natoms);
353 22 : fprintf(fp,"%f %f %f\n",cell[0],cell[1],cell[2]);
354 4370 : for(int iatom=0; iatom<natoms; iatom++) {
355 : // usually, it is better not to apply pbc here, so that diffusion
356 : // is more easily calculated from a trajectory file:
357 3254 : if(wrapatoms) pbc(cell,positions[iatom],pos);
358 7658 : else for(int k=0; k<3; k++) pos[k]=positions[iatom][k];
359 2174 : fprintf(fp,"Ar %10.7f %10.7f %10.7f\n",pos[0],pos[1],pos[2]);
360 : }
361 22 : fclose(fp);
362 22 : }
363 :
364 5 : void write_final_positions(const string& outputfile,int natoms,const vector<Vector>& positions,const double cell[3],const bool wrapatoms)
365 : {
366 : // write positions on file outputfile
367 5 : Vector pos;
368 : FILE*fp;
369 5 : fp=fopen(outputfile.c_str(),"w");
370 : fprintf(fp,"%d\n",natoms);
371 5 : fprintf(fp,"%f %f %f\n",cell[0],cell[1],cell[2]);
372 883 : for(int iatom=0; iatom<natoms; iatom++) {
373 : // usually, it is better not to apply pbc here, so that diffusion
374 : // is more easily calculated from a trajectory file:
375 655 : if(wrapatoms) pbc(cell,positions[iatom],pos);
376 1561 : else for(int k=0; k<3; k++) pos[k]=positions[iatom][k];
377 439 : fprintf(fp,"Ar %10.7f %10.7f %10.7f\n",pos[0],pos[1],pos[2]);
378 : }
379 5 : fclose(fp);
380 5 : }
381 :
382 :
383 22 : void write_statistics(const string & statfile,const int istep,const double tstep,
384 : const int natoms,const int ndim,const double engkin,const double engconf,const double engint) {
385 : // write statistics on file statfile
386 22 : if(write_statistics_first) {
387 : // first time this routine is called, open the file
388 5 : write_statistics_fp=fopen(statfile.c_str(),"w");
389 5 : write_statistics_first=false;
390 : }
391 22 : if(istep-write_statistics_last_time_reopened>100) {
392 : // every 100 steps, reopen the file to flush the buffer
393 2 : fclose(write_statistics_fp);
394 2 : write_statistics_fp=fopen(statfile.c_str(),"a");
395 2 : write_statistics_last_time_reopened=istep;
396 : }
397 22 : fprintf(write_statistics_fp,"%d %f %f %f %f %f\n",istep,istep*tstep,2.0*engkin/double(ndim*natoms),engconf,engkin+engconf,engkin+engconf+engint);
398 22 : }
399 :
400 :
401 :
402 5 : virtual int main(FILE* in,FILE*out,PLMD::Communicator& pc) {
403 : int natoms; // number of atoms
404 : vector<Vector> positions; // atomic positions
405 : vector<Vector> velocities; // velocities
406 : vector<double> masses; // masses
407 : vector<Vector> forces; // forces
408 : double cell[3]; // cell size
409 : double cell9[3][3]; // cell size
410 :
411 : // neighbour list variables
412 : // see Allen and Tildesey book for details
413 : int listsize; // size of the list array
414 : vector<int> list; // neighbour list
415 : vector<int> point; // pointer to neighbour list
416 : vector<Vector> positions0; // reference atomic positions, i.e. positions when the neighbour list
417 :
418 : // input parameters
419 : // all of them have a reasonable default value, set in read_input()
420 : double tstep; // simulation timestep
421 : double temperature; // temperature
422 : double friction; // friction for Langevin dynamics (for NVE, use 0)
423 : double listcutoff; // cutoff for neighbour list
424 : double forcecutoff; // cutoff for forces
425 : int nstep; // number of steps
426 : int nconfig; // stride for output of configurations
427 : int nstat; // stride for output of statistics
428 : int maxneighbour; // maximum average number of neighbours per atom
429 : int ndim; // dimensionality of the system (1, 2, or 3)
430 : int idum; // seed
431 : int plumedWantsToStop; // stop flag
432 : bool wrapatoms; // if true, atomic coordinates are written wrapped in minimal cell
433 : string inputfile; // name of file with starting configuration (xyz)
434 : string outputfile; // name of file with final configuration (xyz)
435 : string trajfile; // name of the trajectory file (xyz)
436 : string statfile; // name of the file with statistics
437 :
438 : double engkin; // kinetic energy
439 : double engconf; // configurational energy
440 : double engint; // integral for conserved energy in Langevin dynamics
441 :
442 : bool recompute_list; // control if the neighbour list have to be recomputed
443 :
444 5 : Random random; // random numbers stream
445 :
446 : PLMD::Plumed* plumed=NULL;
447 :
448 : // Commenting the next line it is possible to switch-off plumed
449 5 : plumed=new PLMD::Plumed;
450 :
451 5 : if(plumed) {
452 5 : int s=sizeof(double);
453 : plumed->cmd("setRealPrecision",&s);
454 : }
455 :
456 5 : read_input(temperature,tstep,friction,forcecutoff,
457 : listcutoff,nstep,nconfig,nstat,
458 : wrapatoms,inputfile,outputfile,trajfile,statfile,
459 : maxneighbour,ndim,idum);
460 :
461 : // number of atoms is read from file inputfile
462 5 : read_natoms(inputfile,natoms);
463 :
464 : // write the parameters in output so they can be checked
465 : fprintf(out,"%s %s\n","Starting configuration :",inputfile.c_str());
466 : fprintf(out,"%s %s\n","Final configuration :",outputfile.c_str());
467 5 : fprintf(out,"%s %d\n","Number of atoms :",natoms);
468 5 : fprintf(out,"%s %f\n","Temperature :",temperature);
469 5 : fprintf(out,"%s %f\n","Time step :",tstep);
470 5 : fprintf(out,"%s %f\n","Friction :",friction);
471 5 : fprintf(out,"%s %f\n","Cutoff for forces :",forcecutoff);
472 5 : fprintf(out,"%s %f\n","Cutoff for neighbour list :",listcutoff);
473 5 : fprintf(out,"%s %d\n","Number of steps :",nstep);
474 5 : fprintf(out,"%s %d\n","Stride for trajectory :",nconfig);
475 : fprintf(out,"%s %s\n","Trajectory file :",trajfile.c_str());
476 5 : fprintf(out,"%s %d\n","Stride for statistics :",nstat);
477 : fprintf(out,"%s %s\n","Statistics file :",statfile.c_str());
478 5 : fprintf(out,"%s %d\n","Max average number of neighbours :",maxneighbour);
479 5 : fprintf(out,"%s %d\n","Dimensionality :",ndim);
480 5 : fprintf(out,"%s %d\n","Seed :",idum);
481 5 : fprintf(out,"%s %s\n","Are atoms wrapped on output? :",(wrapatoms?"T":"F"));
482 :
483 : // Setting the seed
484 5 : random.setSeed(idum);
485 :
486 : // Since each atom pair is counted once, the total number of pairs
487 : // will be half of the number of neighbours times the number of atoms
488 5 : listsize=maxneighbour*natoms/2;
489 :
490 : // allocation of dynamical arrays
491 5 : positions.resize(natoms);
492 5 : positions0.resize(natoms);
493 5 : velocities.resize(natoms);
494 5 : forces.resize(natoms);
495 5 : masses.resize(natoms);
496 5 : point.resize(natoms);
497 5 : list.resize(listsize);
498 :
499 : // masses are hard-coded to 1
500 883 : for(int i=0; i<natoms; ++i) masses[i]=1.0;
501 :
502 : // energy integral initialized to 0
503 5 : engint=0.0;
504 :
505 : // positions are read from file inputfile
506 5 : read_positions(inputfile,natoms,positions,cell);
507 :
508 : // velocities are randomized according to temperature
509 5 : randomize_velocities(natoms,ndim,temperature,masses,velocities,random);
510 :
511 5 : if(plumed) {
512 : plumed->cmd("setNoVirial");
513 : plumed->cmd("setNatoms",&natoms);
514 : plumed->cmd("setMDEngine","simpleMD");
515 : plumed->cmd("setTimestep",&tstep);
516 : plumed->cmd("setPlumedDat","plumed.dat");
517 5 : int pversion=0;
518 : plumed->cmd("getApiVersion",&pversion);
519 : // setting kbT is only implemented with api>1
520 : // even if not necessary in principle in SimpleMD (which is part of plumed)
521 : // we leave the check here as a reference
522 5 : if(pversion>1) {
523 : plumed->cmd("setKbT",&temperature);
524 : }
525 : plumed->cmd("init");
526 : }
527 :
528 : // neighbour list are computed, and reference positions are saved
529 5 : compute_list(natoms,listsize,positions,cell,listcutoff,point,list);
530 :
531 10 : fprintf(out,"List size: %d\n",point[natoms-1]);
532 3078 : for(int iatom=0; iatom<natoms; ++iatom) for(int k=0; k<3; ++k) positions0[iatom][k]=positions[iatom][k];
533 :
534 : // forces are computed before starting md
535 5 : compute_forces(natoms,listsize,positions,cell,forcecutoff,point,list,forces,engconf);
536 :
537 : // remove forces if ndim<3
538 5 : if(ndim<3)
539 15 : for(int iatom=0; iatom<natoms; ++iatom) for(int k=ndim; k<3; ++k) forces[iatom][k]=0.0;
540 :
541 : // here is the main md loop
542 : // Langevin thermostat is applied before and after a velocity-Verlet integrator
543 : // the overall structure is:
544 : // thermostat
545 : // update velocities
546 : // update positions
547 : // (eventually recompute neighbour list)
548 : // compute forces
549 : // update velocities
550 : // thermostat
551 : // (eventually dump output informations)
552 2205 : for(int istep=0; istep<nstep; istep++) {
553 2200 : thermostat(natoms,ndim,masses,0.5*tstep,friction,temperature,velocities,engint,random);
554 :
555 144600 : for(int iatom=0; iatom<natoms; iatom++) for(int k=0; k<3; k++)
556 427200 : velocities[iatom][k]+=forces[iatom][k]*0.5*tstep/masses[iatom];
557 :
558 144600 : for(int iatom=0; iatom<natoms; iatom++) for(int k=0; k<3; k++)
559 320400 : positions[iatom][k]+=velocities[iatom][k]*tstep;
560 :
561 : // a check is performed to decide whether to recalculate the neighbour list
562 2200 : check_list(natoms,positions,positions0,listcutoff,forcecutoff,recompute_list);
563 2200 : if(recompute_list) {
564 32 : compute_list(natoms,listsize,positions,cell,listcutoff,point,list);
565 7963 : for(int iatom=0; iatom<natoms; ++iatom) for(int k=0; k<3; ++k) positions0[iatom][k]=positions[iatom][k];
566 : fprintf(out,"Neighbour list recomputed at step %d\n",istep);
567 64 : fprintf(out,"List size: %d\n",point[natoms-1]);
568 : }
569 :
570 2200 : compute_forces(natoms,listsize,positions,cell,forcecutoff,point,list,forces,engconf);
571 :
572 2200 : if(plumed) {
573 2200 : int istepplusone=istep+1;
574 2200 : plumedWantsToStop=0;
575 8800 : for(int i=0; i<3; i++)for(int k=0; k<3; k++) cell9[i][k]=0.0;
576 8800 : for(int i=0; i<3; i++) cell9[i][i]=cell[i];
577 : plumed->cmd("setStep",&istepplusone);
578 : plumed->cmd("setMasses",&masses[0]);
579 : plumed->cmd("setForces",&forces[0]);
580 : plumed->cmd("setEnergy",&engconf);
581 : plumed->cmd("setPositions",&positions[0]);
582 : plumed->cmd("setBox",cell9);
583 : plumed->cmd("setStopFlag",&plumedWantsToStop);
584 : plumed->cmd("calc");
585 2200 : if(plumedWantsToStop) nstep=istep;
586 : }
587 : // remove forces if ndim<3
588 2200 : if(ndim<3)
589 30000 : for(int iatom=0; iatom<natoms; ++iatom) for(int k=ndim; k<3; ++k) forces[iatom][k]=0.0;
590 :
591 144600 : for(int iatom=0; iatom<natoms; iatom++) for(int k=0; k<3; k++)
592 427200 : velocities[iatom][k]+=forces[iatom][k]*0.5*tstep/masses[iatom];
593 :
594 2200 : thermostat(natoms,ndim,masses,0.5*tstep,friction,temperature,velocities,engint,random);
595 :
596 : // kinetic energy is calculated
597 2200 : compute_engkin(natoms,masses,velocities,engkin);
598 :
599 : // eventually, write positions and statistics
600 2200 : if((istep+1)%nconfig==0) write_positions(trajfile,natoms,positions,cell,wrapatoms);
601 2200 : if((istep+1)%nstat==0) write_statistics(statfile,istep+1,tstep,natoms,ndim,engkin,engconf,engint);
602 :
603 : }
604 :
605 : // call final plumed jobs
606 : plumed->cmd("runFinalJobs");
607 :
608 : // write final positions
609 5 : write_final_positions(outputfile,natoms,positions,cell,wrapatoms);
610 :
611 : // close the statistic file if it was open:
612 5 : if(write_statistics_fp) fclose(write_statistics_fp);
613 :
614 5 : if(plumed) delete plumed;
615 :
616 5 : return 0;
617 : }
618 :
619 :
620 : };
621 :
622 6462 : PLUMED_REGISTER_CLTOOL(SimpleMD,"simplemd")
623 :
624 : }
625 4839 : }
626 :
627 :
628 :
629 :
|
