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