Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "core/CLToolRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "tools/Vector.h"
26 : #include "tools/Random.h"
27 : #include "tools/OpenMP.h"
28 : #include <string>
29 : #include <cstdio>
30 : #include <cmath>
31 : #include <vector>
32 : #include <memory>
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 specified 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 : inputfile 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 : // simple static function to close a file
77 : // defined once here since it's used in many places in this file
78 : // in addition, this seems the only way to use it in the write_statistics_fp_deleter member
79 : static void (*deleter)(FILE* f) = [](auto f) { if(f) std::fclose(f); };
80 :
81 : class SimpleMD:
82 : public PLMD::CLTool
83 : {
84 4 : std::string description()const override {
85 4 : return "run lj code";
86 : }
87 :
88 : bool write_positions_first;
89 : bool write_statistics_first;
90 : int write_statistics_last_time_reopened;
91 : FILE* write_statistics_fp;
92 : std::unique_ptr<FILE,decltype(deleter)> write_statistics_fp_deleter{nullptr,deleter};
93 :
94 :
95 : public:
96 5316 : static void registerKeywords( Keywords& keys ) {
97 10632 : keys.add("compulsory","nstep","The number of steps of dynamics you want to run");
98 10632 : keys.add("compulsory","temperature","NVE","the temperature at which you wish to run the simulation in LJ units");
99 10632 : keys.add("compulsory","friction","off","The friction (in LJ units) for the Langevin thermostat that is used to keep the temperature constant");
100 10632 : keys.add("compulsory","tstep","0.005","the integration timestep in LJ units");
101 10632 : keys.add("compulsory","epsilon","1.0","LJ parameter");
102 10632 : keys.add("compulsory","sigma","1.0","LJ parameter");
103 10632 : keys.add("compulsory","inputfile","An xyz file containing the initial configuration of the system");
104 10632 : keys.add("compulsory","forcecutoff","2.5","");
105 10632 : keys.add("compulsory","listcutoff","3.0","");
106 10632 : keys.add("compulsory","outputfile","An output xyz file containing the final configuration of the system");
107 10632 : keys.add("compulsory","nconfig","10","The frequency with which to write configurations to the trajectory file followed by the name of the trajectory file");
108 10632 : 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");
109 10632 : keys.add("compulsory","idum","0","The random number seed");
110 10632 : keys.add("compulsory","ndim","3","The dimensionality of the system (some interesting LJ clusters are two dimensional)");
111 10632 : keys.add("compulsory","wrapatoms","false","If true, atomic coordinates are written wrapped in minimal cell");
112 5316 : }
113 :
114 13 : explicit SimpleMD( const CLToolOptions& co ) :
115 : CLTool(co),
116 13 : write_positions_first(true),
117 13 : write_statistics_first(true),
118 13 : write_statistics_last_time_reopened(0),
119 13 : write_statistics_fp(NULL)
120 : {
121 13 : inputdata=ifile;
122 13 : }
123 :
124 : private:
125 :
126 : void
127 9 : read_input(double& temperature,
128 : double& tstep,
129 : double& friction,
130 : double& forcecutoff,
131 : double& listcutoff,
132 : int& nstep,
133 : int& nconfig,
134 : int& nstat,
135 : bool& wrapatoms,
136 : std::string& inputfile,
137 : std::string& outputfile,
138 : std::string& trajfile,
139 : std::string& statfile,
140 : int& ndim,
141 : int& idum,
142 : double& epsilon,
143 : double& sigma)
144 : {
145 :
146 : // Read everything from input file
147 : char buffer1[256];
148 18 : std::string tempstr; parse("temperature",tempstr);
149 9 : if( tempstr!="NVE" ) Tools::convert(tempstr,temperature);
150 18 : parse("tstep",tstep);
151 18 : std::string frictionstr; parse("friction",frictionstr);
152 9 : if( tempstr!="NVE" ) {
153 9 : if(frictionstr=="off") error("Specify friction for thermostat");
154 9 : Tools::convert(frictionstr,friction);
155 : }
156 9 : parse("forcecutoff",forcecutoff);
157 9 : parse("listcutoff",listcutoff);
158 9 : parse("nstep",nstep);
159 9 : parse("idum",idum);
160 9 : parse("epsilon",epsilon);
161 9 : parse("sigma",sigma);
162 :
163 : // Read in stuff with sanity checks
164 18 : parse("inputfile",inputfile);
165 9 : if(inputfile.length()==0) error("Specify input file");
166 18 : parse("outputfile",outputfile);
167 9 : if(outputfile.length()==0) error("Specify output file");
168 18 : std::string nconfstr; parse("nconfig",nconfstr);
169 9 : std::sscanf(nconfstr.c_str(),"%100d %255s",&nconfig,buffer1);
170 : trajfile=buffer1;
171 9 : if(trajfile.length()==0) error("Specify traj file");
172 18 : std::string nstatstr; parse("nstat",nstatstr);
173 9 : std::sscanf(nstatstr.c_str(),"%100d %255s",&nstat,buffer1);
174 : statfile=buffer1;
175 9 : if(statfile.length()==0) error("Specify stat file");
176 9 : parse("ndim",ndim);
177 9 : if(ndim<1 || ndim>3) error("ndim should be 1,2 or 3");
178 : std::string w;
179 9 : parse("wrapatoms",w);
180 9 : wrapatoms=false;
181 9 : if(w.length()>0 && (w[0]=='T' || w[0]=='t')) wrapatoms=true;
182 9 : }
183 :
184 9 : void read_natoms(const std::string & inputfile,int & natoms) {
185 : // read the number of atoms in file "input.xyz"
186 9 : FILE* fp=std::fopen(inputfile.c_str(),"r");
187 9 : if(!fp) error(std::string("file ") + inputfile + std::string(" not found"));
188 :
189 : // call fclose when fp goes out of scope
190 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
191 :
192 9 : int ret=std::fscanf(fp,"%1000d",&natoms);
193 9 : if(ret==0) plumed_error() <<"Error reading number of atoms from file "<<inputfile;
194 9 : }
195 :
196 9 : void read_positions(const std::string& inputfile,int natoms,std::vector<Vector>& positions,double cell[3]) {
197 : // read positions and cell from a file called inputfile
198 : // natoms (input variable) and number of atoms in the file should be consistent
199 9 : FILE* fp=std::fopen(inputfile.c_str(),"r");
200 9 : if(!fp) error(std::string("file ") + inputfile + std::string(" not found"));
201 : // call fclose when fp goes out of scope
202 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
203 :
204 : char buffer[256];
205 : char atomname[256];
206 : char* cret=fgets(buffer,256,fp);
207 9 : if(cret==nullptr) plumed_error() <<"Error reading buffer from file "<<inputfile;
208 9 : int ret=std::fscanf(fp,"%1000lf %1000lf %1000lf",&cell[0],&cell[1],&cell[2]);
209 9 : if(ret==0) plumed_error() <<"Error reading cell line from file "<<inputfile;
210 880 : for(int i=0; i<natoms; i++) {
211 871 : ret=std::fscanf(fp,"%255s %1000lf %1000lf %1000lf",atomname,&positions[i][0],&positions[i][1],&positions[i][2]);
212 : // note: atomname is read but not used
213 871 : if(ret==0) plumed_error() <<"Error reading atom line from file "<<inputfile;
214 : }
215 9 : }
216 :
217 9 : void randomize_velocities(const int natoms,const int ndim,const double temperature,const std::vector<double>&masses,std::vector<Vector>& velocities,Random&random) {
218 : // randomize the velocities according to the temperature
219 3486 : for(int iatom=0; iatom<natoms; iatom++) for(int i=0; i<ndim; i++)
220 2606 : velocities[iatom][i]=std::sqrt(temperature/masses[iatom])*random.Gaussian();
221 9 : }
222 :
223 10165978 : void pbc(const double cell[3],const Vector & vin,Vector & vout) {
224 : // apply periodic boundary condition to a vector
225 40663912 : for(int i=0; i<3; i++) {
226 30497934 : vout[i]=vin[i]-std::floor(vin[i]/cell[i]+0.5)*cell[i];
227 : }
228 10165978 : }
229 :
230 3750 : void check_list(const int natoms,const std::vector<Vector>& positions,const std::vector<Vector>&positions0,const double listcutoff,
231 : const double forcecutoff,bool & recompute)
232 : {
233 : // check if the neighbour list have to be recomputed
234 3750 : recompute=false;
235 3750 : auto delta2=(0.5*(listcutoff-forcecutoff))*(0.5*(listcutoff-forcecutoff));
236 : // if ANY atom moved more than half of the skin thickness, recompute is set to .true.
237 206750 : for(int iatom=0; iatom<natoms; iatom++) {
238 203000 : if(modulo2(positions[iatom]-positions0[iatom])>delta2) recompute=true;
239 : }
240 3750 : }
241 :
242 :
243 437 : void compute_list(const int natoms,const std::vector<Vector>& positions,const double cell[3],const double listcutoff,
244 : std::vector<std::vector<int> >& list) {
245 437 : double listcutoff2=listcutoff*listcutoff; // squared list cutoff
246 437 : list.assign(natoms,std::vector<int>());
247 437 : # pragma omp parallel for num_threads(OpenMP::getNumThreads()) schedule(static,1)
248 : for(int iatom=0; iatom<natoms-1; iatom++) {
249 : for(int jatom=iatom+1; jatom<natoms; jatom++) {
250 : auto distance=positions[iatom]-positions[jatom];
251 : Vector distance_pbc; // minimum-image distance of the two atoms
252 : pbc(cell,distance,distance_pbc);
253 : // if the interparticle distance is larger than the cutoff, skip
254 : if(modulo2(distance_pbc)>listcutoff2)continue;
255 : list[iatom].push_back(jatom);
256 : }
257 : }
258 437 : }
259 :
260 3759 : void compute_forces(const int natoms,double epsilon, double sigma,const std::vector<Vector>& positions,const double cell[3],
261 : double forcecutoff,const std::vector<std::vector<int> >& list,std::vector<Vector>& forces,double & engconf)
262 : {
263 3759 : double forcecutoff2=forcecutoff*forcecutoff; // squared force cutoff
264 3759 : engconf=0.0;
265 : double engconf_tmp=0.0; // separate variable for reduction (gcc 4.8 does not make reductions on references)
266 819243 : for(int i=0; i<natoms; i++)for(int k=0; k<3; k++) forces[i][k]=0.0;
267 3759 : double engcorrection=4.0*epsilon*(1.0/std::pow(forcecutoff2/(sigma*sigma),6.0)-1.0/std::pow(forcecutoff2/(sigma*sigma),3)); // energy necessary shift the potential avoiding discontinuities
268 3759 : # pragma omp parallel num_threads(OpenMP::getNumThreads())
269 : {
270 : std::vector<Vector> omp_forces(forces.size());
271 : #pragma omp for reduction(+:engconf_tmp) schedule(static,1) nowait
272 : for(int iatom=0; iatom<natoms-1; iatom++) {
273 : for(int jlist=0; jlist<list[iatom].size(); jlist++) {
274 : const int jatom=list[iatom][jlist];
275 : auto distance=positions[iatom]-positions[jatom];
276 : Vector distance_pbc; // minimum-image distance of the two atoms
277 : pbc(cell,distance,distance_pbc);
278 : auto distance_pbc2=modulo2(distance_pbc); // squared minimum-image distance
279 : // if the interparticle distance is larger than the cutoff, skip
280 : if(distance_pbc2>forcecutoff2) continue;
281 : auto distance_pbcm2=sigma*sigma/distance_pbc2;
282 : auto distance_pbcm6=distance_pbcm2*distance_pbcm2*distance_pbcm2;
283 : auto distance_pbcm8=distance_pbcm6*distance_pbcm2;
284 : auto distance_pbcm12=distance_pbcm6*distance_pbcm6;
285 : auto distance_pbcm14=distance_pbcm12*distance_pbcm2;
286 : engconf_tmp+=4.0*epsilon*(distance_pbcm12 - distance_pbcm6) - engcorrection;
287 : auto f=24.0*distance_pbc*(2.0*distance_pbcm14-distance_pbcm8)*epsilon/sigma;
288 : omp_forces[iatom]+=f;
289 : omp_forces[jatom]-=f;
290 : }
291 : }
292 : # pragma omp critical
293 : for(unsigned i=0; i<omp_forces.size(); i++) forces[i]+=omp_forces[i];
294 : }
295 :
296 3759 : engconf=engconf_tmp;
297 :
298 3759 : }
299 :
300 3750 : void compute_engkin(const int natoms,const std::vector<double>& masses,const std::vector<Vector>& velocities,double & engkin)
301 : {
302 : // calculate the kinetic energy from the velocities
303 3750 : engkin=0.0;
304 206750 : for(int iatom=0; iatom<natoms; iatom++) {
305 203000 : engkin+=0.5*masses[iatom]*modulo2(velocities[iatom]);
306 : }
307 3750 : }
308 :
309 :
310 7500 : void thermostat(const int natoms,const int ndim,const std::vector<double>& masses,const double dt,const double friction,
311 : const double temperature,std::vector<Vector>& velocities,double & engint,Random & random) {
312 : // Langevin thermostat, implemented as described in Bussi and Parrinello, Phys. Rev. E (2007)
313 : // it is a linear combination of old velocities and new, randomly chosen, velocity,
314 : // with proper coefficients
315 7500 : double c1=std::exp(-friction*dt);
316 413500 : for(int iatom=0; iatom<natoms; iatom++) {
317 406000 : double c2=std::sqrt((1.0-c1*c1)*temperature/masses[iatom]);
318 1596000 : for(int i=0; i<ndim; i++) {
319 1190000 : engint+=0.5*masses[iatom]*velocities[iatom][i]*velocities[iatom][i];
320 1190000 : velocities[iatom][i]=c1*velocities[iatom][i]+c2*random.Gaussian();
321 1190000 : engint-=0.5*masses[iatom]*velocities[iatom][i]*velocities[iatom][i];
322 : }
323 : }
324 7500 : }
325 :
326 42 : void write_positions(const std::string& trajfile,int natoms,const std::vector<Vector>& positions,const double cell[3],const bool wrapatoms)
327 : {
328 : // write positions on file trajfile
329 : // positions are appended at the end of the file
330 42 : Vector pos;
331 : FILE*fp;
332 42 : if(write_positions_first) {
333 9 : fp=std::fopen(trajfile.c_str(),"w");
334 9 : write_positions_first=false;
335 : } else {
336 33 : fp=std::fopen(trajfile.c_str(),"a");
337 : }
338 42 : plumed_assert(fp);
339 : // call fclose when fp goes out of scope
340 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
341 : std::fprintf(fp,"%d\n",natoms);
342 42 : std::fprintf(fp,"%f %f %f\n",cell[0],cell[1],cell[2]);
343 4376 : for(int iatom=0; iatom<natoms; iatom++) {
344 : // usually, it is better not to apply pbc here, so that diffusion
345 : // is more easily calculated from a trajectory file:
346 4334 : if(wrapatoms) pbc(cell,positions[iatom],pos);
347 3254 : else pos=positions[iatom];
348 4334 : std::fprintf(fp,"Ar %10.7f %10.7f %10.7f\n",pos[0],pos[1],pos[2]);
349 : }
350 42 : }
351 :
352 9 : void write_final_positions(const std::string& outputfile,int natoms,const std::vector<Vector>& positions,const double cell[3],const bool wrapatoms)
353 : {
354 : // write positions on file outputfile
355 9 : Vector pos;
356 : FILE*fp;
357 9 : fp=std::fopen(outputfile.c_str(),"w");
358 9 : plumed_assert(fp);
359 : // call fclose when fp goes out of scope
360 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
361 : std::fprintf(fp,"%d\n",natoms);
362 9 : std::fprintf(fp,"%f %f %f\n",cell[0],cell[1],cell[2]);
363 880 : for(int iatom=0; iatom<natoms; iatom++) {
364 : // usually, it is better not to apply pbc here, so that diffusion
365 : // is more easily calculated from a trajectory file:
366 871 : if(wrapatoms) pbc(cell,positions[iatom],pos);
367 655 : else pos=positions[iatom];
368 871 : std::fprintf(fp,"Ar %10.7f %10.7f %10.7f\n",pos[0],pos[1],pos[2]);
369 : }
370 9 : }
371 :
372 :
373 177 : void write_statistics(const std::string & statfile,const int istep,const double tstep,
374 : const int natoms,const int ndim,const double engkin,const double engconf,const double engint) {
375 : // write statistics on file statfile
376 177 : if(write_statistics_first) {
377 : // first time this routine is called, open the file
378 9 : write_statistics_fp=std::fopen(statfile.c_str(),"w");
379 : write_statistics_fp_deleter.reset(write_statistics_fp);
380 9 : write_statistics_first=false;
381 : }
382 177 : if(istep-write_statistics_last_time_reopened>100) {
383 : // every 100 steps, reopen the file to flush the buffer
384 : write_statistics_fp_deleter.reset(nullptr); // close file
385 14 : write_statistics_fp=std::fopen(statfile.c_str(),"a");
386 : write_statistics_fp_deleter.reset(write_statistics_fp);
387 14 : write_statistics_last_time_reopened=istep;
388 : }
389 177 : std::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);
390 177 : }
391 :
392 :
393 :
394 9 : int main(FILE* in,FILE*out,PLMD::Communicator& pc) override {
395 : int natoms; // number of atoms
396 : std::vector<Vector> positions; // atomic positions
397 : std::vector<Vector> velocities; // velocities
398 : std::vector<double> masses; // masses
399 : std::vector<Vector> forces; // forces
400 : double cell[3]; // cell size
401 : double cell9[3][3]; // cell size
402 :
403 : // neighbour list variables
404 : // see Allen and Tildesey book for details
405 : std::vector< std::vector<int> > list; // neighbour list
406 : std::vector<Vector> positions0; // reference atomic positions, i.e. positions when the neighbour list
407 :
408 : // input parameters
409 : // all of them have a reasonable default value, set in read_input()
410 : double tstep; // simulation timestep
411 : double temperature; // temperature
412 : double friction; // friction for Langevin dynamics (for NVE, use 0)
413 : double listcutoff; // cutoff for neighbour list
414 : double forcecutoff; // cutoff for forces
415 : int nstep; // number of steps
416 : int nconfig; // stride for output of configurations
417 : int nstat; // stride for output of statistics
418 : int ndim; // dimensionality of the system (1, 2, or 3)
419 : int idum; // seed
420 : int plumedWantsToStop; // stop flag
421 : bool wrapatoms; // if true, atomic coordinates are written wrapped in minimal cell
422 : std::string inputfile; // name of file with starting configuration (xyz)
423 : std::string outputfile; // name of file with final configuration (xyz)
424 : std::string trajfile; // name of the trajectory file (xyz)
425 : std::string statfile; // name of the file with statistics
426 :
427 : double epsilon, sigma; // LJ parameters
428 :
429 : double engkin; // kinetic energy
430 : double engconf; // configurational energy
431 : double engint; // integral for conserved energy in Langevin dynamics
432 :
433 : bool recompute_list; // control if the neighbour list have to be recomputed
434 :
435 9 : Random random; // random numbers stream
436 :
437 9 : std::unique_ptr<PlumedMain> plumed;
438 :
439 : // Commenting the next line it is possible to switch-off plumed
440 9 : plumed=Tools::make_unique<PLMD::PlumedMain>();
441 :
442 9 : if(plumed) plumed->cmd("setRealPrecision",int(sizeof(double)));
443 :
444 9 : read_input(temperature,tstep,friction,forcecutoff,
445 : listcutoff,nstep,nconfig,nstat,
446 : wrapatoms,inputfile,outputfile,trajfile,statfile,
447 : ndim,idum,epsilon,sigma);
448 :
449 : // number of atoms is read from file inputfile
450 9 : read_natoms(inputfile,natoms);
451 :
452 : // write the parameters in output so they can be checked
453 : std::fprintf(out,"%s %s\n","Starting configuration :",inputfile.c_str());
454 : std::fprintf(out,"%s %s\n","Final configuration :",outputfile.c_str());
455 9 : std::fprintf(out,"%s %d\n","Number of atoms :",natoms);
456 9 : std::fprintf(out,"%s %f\n","Temperature :",temperature);
457 9 : std::fprintf(out,"%s %f\n","Time step :",tstep);
458 9 : std::fprintf(out,"%s %f\n","Friction :",friction);
459 9 : std::fprintf(out,"%s %f\n","Cutoff for forces :",forcecutoff);
460 9 : std::fprintf(out,"%s %f\n","Cutoff for neighbour list :",listcutoff);
461 9 : std::fprintf(out,"%s %d\n","Number of steps :",nstep);
462 9 : std::fprintf(out,"%s %d\n","Stride for trajectory :",nconfig);
463 : std::fprintf(out,"%s %s\n","Trajectory file :",trajfile.c_str());
464 9 : std::fprintf(out,"%s %d\n","Stride for statistics :",nstat);
465 : std::fprintf(out,"%s %s\n","Statistics file :",statfile.c_str());
466 9 : std::fprintf(out,"%s %d\n","Dimensionality :",ndim);
467 9 : std::fprintf(out,"%s %d\n","Seed :",idum);
468 9 : std::fprintf(out,"%s %s\n","Are atoms wrapped on output? :",(wrapatoms?"T":"F"));
469 9 : std::fprintf(out,"%s %f\n","Epsilon :",epsilon);
470 9 : std::fprintf(out,"%s %f\n","Sigma :",sigma);
471 :
472 : // Setting the seed
473 9 : random.setSeed(idum);
474 :
475 : // allocation of dynamical arrays
476 9 : positions.resize(natoms);
477 9 : positions0.resize(natoms);
478 9 : velocities.resize(natoms);
479 9 : forces.resize(natoms);
480 9 : masses.resize(natoms);
481 9 : list.resize(natoms);
482 :
483 : // masses are hard-coded to 1
484 880 : for(int i=0; i<natoms; ++i) masses[i]=1.0;
485 :
486 : // energy integral initialized to 0
487 9 : engint=0.0;
488 :
489 : // positions are read from file inputfile
490 9 : read_positions(inputfile,natoms,positions,cell);
491 :
492 : // velocities are randomized according to temperature
493 9 : randomize_velocities(natoms,ndim,temperature,masses,velocities,random);
494 :
495 9 : if(plumed) {
496 9 : plumed->cmd("setNatoms",natoms);
497 9 : plumed->cmd("setNoVirial");
498 9 : plumed->cmd("setMDEngine","simpleMD");
499 9 : plumed->cmd("setTimestep",tstep);
500 9 : plumed->cmd("setPlumedDat","plumed.dat");
501 9 : int pversion=0;
502 9 : plumed->cmd("getApiVersion",&pversion);
503 : // setting kbT is only implemented with api>1
504 : // even if not necessary in principle in SimpleMD (which is part of plumed)
505 : // we leave the check here as a reference
506 9 : if(pversion>1) plumed->cmd("setKbT",temperature);
507 9 : plumed->cmd("init");
508 : }
509 :
510 : // neighbour list are computed, and reference positions are saved
511 9 : compute_list(natoms,positions,cell,listcutoff,list);
512 :
513 : int list_size=0;
514 880 : for(int i=0; i<list.size(); i++) list_size+=list[i].size();
515 : std::fprintf(out,"List size: %d\n",list_size);
516 880 : for(int iatom=0; iatom<natoms; ++iatom) positions0[iatom]=positions[iatom];
517 :
518 : // forces are computed before starting md
519 9 : compute_forces(natoms,epsilon,sigma,positions,cell,forcecutoff,list,forces,engconf);
520 :
521 : // remove forces if ndim<3
522 9 : if(ndim<3)
523 15 : for(int iatom=0; iatom<natoms; ++iatom) for(int k=ndim; k<3; ++k) forces[iatom][k]=0.0;
524 :
525 : // here is the main md loop
526 : // Langevin thermostat is applied before and after a velocity-Verlet integrator
527 : // the overall structure is:
528 : // thermostat
529 : // update velocities
530 : // update positions
531 : // (eventually recompute neighbour list)
532 : // compute forces
533 : // update velocities
534 : // thermostat
535 : // (eventually dump output informations)
536 3759 : for(int istep=0; istep<nstep; istep++) {
537 3750 : thermostat(natoms,ndim,masses,0.5*tstep,friction,temperature,velocities,engint,random);
538 :
539 206750 : for(int iatom=0; iatom<natoms; iatom++)
540 203000 : velocities[iatom]+=forces[iatom]*0.5*tstep/masses[iatom];
541 :
542 206750 : for(int iatom=0; iatom<natoms; iatom++)
543 203000 : positions[iatom]+=velocities[iatom]*tstep;
544 :
545 : // a check is performed to decide whether to recalculate the neighbour list
546 3750 : check_list(natoms,positions,positions0,listcutoff,forcecutoff,recompute_list);
547 3750 : if(recompute_list) {
548 428 : compute_list(natoms,positions,cell,listcutoff,list);
549 44329 : for(int iatom=0; iatom<natoms; ++iatom) positions0[iatom]=positions[iatom];
550 : int list_size=0;
551 44329 : for(int i=0; i<list.size(); i++) list_size+=list[i].size();
552 : std::fprintf(out,"List size: %d\n",list_size);
553 : }
554 :
555 3750 : compute_forces(natoms,epsilon,sigma,positions,cell,forcecutoff,list,forces,engconf);
556 :
557 3750 : if(plumed) {
558 3750 : plumedWantsToStop=0;
559 48750 : for(int i=0; i<3; i++)for(int k=0; k<3; k++) cell9[i][k]=0.0;
560 15000 : for(int i=0; i<3; i++) cell9[i][i]=cell[i];
561 3750 : plumed->cmd("setStep",istep+1);
562 3750 : plumed->cmd("setMasses",&masses[0]);
563 3750 : plumed->cmd("setForces",&forces[0][0]);
564 3750 : plumed->cmd("setEnergy",engconf);
565 3750 : plumed->cmd("setPositions",&positions[0][0]);
566 3750 : plumed->cmd("setBox",&cell9[0][0]);
567 3750 : plumed->cmd("setStopFlag",&plumedWantsToStop);
568 3750 : plumed->cmd("calc");
569 3750 : if(plumedWantsToStop) nstep=istep;
570 : }
571 : // remove forces if ndim<3
572 3750 : if(ndim<3)
573 30000 : for(int iatom=0; iatom<natoms; ++iatom) for(int k=ndim; k<3; ++k) forces[iatom][k]=0.0;
574 :
575 206750 : for(int iatom=0; iatom<natoms; iatom++)
576 203000 : velocities[iatom]+=forces[iatom]*0.5*tstep/masses[iatom];
577 :
578 3750 : thermostat(natoms,ndim,masses,0.5*tstep,friction,temperature,velocities,engint,random);
579 :
580 : // kinetic energy is calculated
581 3750 : compute_engkin(natoms,masses,velocities,engkin);
582 :
583 : // eventually, write positions and statistics
584 3750 : if((istep+1)%nconfig==0) write_positions(trajfile,natoms,positions,cell,wrapatoms);
585 3750 : if((istep+1)%nstat==0) write_statistics(statfile,istep+1,tstep,natoms,ndim,engkin,engconf,engint);
586 :
587 : }
588 :
589 : // call final plumed jobs
590 9 : plumed->cmd("runFinalJobs");
591 :
592 : // write final positions
593 9 : write_final_positions(outputfile,natoms,positions,cell,wrapatoms);
594 :
595 9 : return 0;
596 18 : }
597 :
598 :
599 : };
600 :
601 15961 : PLUMED_REGISTER_CLTOOL(SimpleMD,"simplemd")
602 :
603 : }
604 : }
605 :
606 :
607 :
608 :
|
