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