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 : #include <algorithm>
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 : You can use simplemd to do a molecular dynamics simulation on a system of Lennard-Jones particles by using the following command
43 :
44 : ```plumed
45 : plumed simplemd < in
46 : ```
47 :
48 : The input to simplemd is specified in an input file that contains one directive per line like the one shown below.
49 :
50 : ```plumed
51 : #TOOL=simplemd
52 : inputfile input.xyz
53 : outputfile output.xyz
54 : temperature 0.722
55 : tstep 0.005
56 : friction 1
57 : forcecutoff 2.5
58 : listcutoff 3.0
59 : nstep 50
60 : nconfig 10 trajectory.xyz
61 : nstat 10 energies.dat
62 : ```
63 :
64 : 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).
65 : Notice that the atomic configurations are input and output in xyz format.
66 :
67 : If you run the following command a description of all the directives that can be used in the input file will be output.
68 :
69 : ```plumed
70 : plumed simplemd --help
71 : ```
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) {
80 : if(f) {
81 : std::fclose(f);
82 : }
83 : };
84 :
85 : class SimpleMD:
86 : public PLMD::CLTool {
87 7 : std::string description()const override {
88 7 : return "run lj code";
89 : }
90 :
91 : bool write_positions_first;
92 : bool write_statistics_first;
93 : int write_statistics_last_time_reopened;
94 : FILE* write_statistics_fp;
95 : std::unique_ptr<FILE,decltype(deleter)> write_statistics_fp_deleter{nullptr,deleter};
96 :
97 :
98 : public:
99 7321 : static void registerKeywords( Keywords& keys ) {
100 7321 : CLTool::registerKeywords( keys );
101 7321 : keys.add("compulsory","nstep","The number of steps of dynamics you want to run");
102 7321 : keys.add("compulsory","temperature","NVE","the temperature at which you wish to run the simulation in LJ units. With NVE there is no thermostat and the atoms start at rest");
103 7321 : keys.add("compulsory","friction","off","The friction (in LJ units) for the Langevin thermostat that is used to keep the temperature constant");
104 7321 : keys.add("compulsory","tstep","0.005","the integration timestep in LJ units");
105 7321 : keys.add("compulsory","epsilon","1.0","LJ parameter");
106 7321 : keys.add("compulsory","sigma","1.0","LJ parameter");
107 7321 : keys.add("compulsory","inputfile","An xyz file containing the initial configuration of the system");
108 7321 : keys.add("compulsory","forcecutoff","2.5","the cutoff for the Lennard-Jones interactions in LJ units (the potential is shifted to zero at the cutoff)");
109 7321 : keys.add("compulsory","listcutoff","3.0","the cutoff for the neighbor list in LJ units, should be larger than forcecutoff");
110 7321 : keys.add("compulsory","outputfile","An output xyz file containing the final configuration of the system");
111 7321 : keys.add("compulsory","nconfig","10 trajectory.xyz","The frequency with which to write configurations to the trajectory file followed by the name of the trajectory file");
112 7321 : keys.add("compulsory","nstat","1 energies.dat","The frequency with which to write the statistics to the statistics file followed by the name of the statistics file");
113 7321 : keys.add("compulsory","idum","0","The random number seed");
114 7321 : keys.add("compulsory","ndim","3","The dimensionality of the system (some interesting LJ clusters are two dimensional)");
115 7321 : keys.add("compulsory","wrapatoms","false","If true, atomic coordinates are written wrapped in minimal cell");
116 7321 : }
117 :
118 21 : explicit SimpleMD( const CLToolOptions& co ) :
119 : CLTool(co),
120 21 : write_positions_first(true),
121 21 : write_statistics_first(true),
122 21 : write_statistics_last_time_reopened(0),
123 21 : write_statistics_fp(NULL) {
124 21 : inputdata=inputType::ifile;
125 21 : }
126 :
127 : private:
128 :
129 : void
130 14 : read_input(double& temperature,
131 : double& tstep,
132 : double& friction,
133 : double& forcecutoff,
134 : double& listcutoff,
135 : int& nstep,
136 : int& nconfig,
137 : int& nstat,
138 : bool& wrapatoms,
139 : std::string& inputfile,
140 : std::string& outputfile,
141 : std::string& trajfile,
142 : std::string& statfile,
143 : int& ndim,
144 : int& idum,
145 : double& epsilon,
146 : double& sigma) {
147 :
148 : // Read everything from input file
149 : char buffer1[256];
150 : std::string tempstr;
151 31 : parse("temperature",tempstr);
152 : std::string frictionstr;
153 14 : parse("friction",frictionstr);
154 14 : if( tempstr=="NVE" ) {
155 : // no thermostat, atoms start at rest
156 3 : temperature=0.0;
157 3 : friction=0.0;
158 : } else {
159 11 : if(!Tools::convertNoexcept(tempstr,temperature) || temperature<0.0) {
160 0 : error("temperature should be NVE or a non negative number");
161 : }
162 11 : if(frictionstr=="off") {
163 1 : error("Specify friction for thermostat");
164 : }
165 10 : if(!Tools::convertNoexcept(frictionstr,friction) || friction<0.0) {
166 0 : error("friction should be a non negative number");
167 : }
168 : }
169 13 : parse("tstep",tstep);
170 13 : parse("forcecutoff",forcecutoff);
171 13 : parse("listcutoff",listcutoff);
172 13 : parse("nstep",nstep);
173 13 : parse("idum",idum);
174 13 : parse("epsilon",epsilon);
175 13 : parse("sigma",sigma);
176 :
177 : // Read in stuff with sanity checks
178 26 : parse("inputfile",inputfile);
179 13 : if(inputfile.length()==0) {
180 0 : error("Specify input file");
181 : }
182 13 : if(sigma<=0.0) {
183 0 : error("sigma should be positive");
184 : }
185 13 : if(listcutoff<forcecutoff) {
186 1 : error("listcutoff should not be smaller than forcecutoff");
187 : }
188 24 : parse("outputfile",outputfile);
189 12 : if(outputfile.length()==0) {
190 3 : error("Specify output file");
191 : }
192 : std::string nconfstr;
193 12 : parse("nconfig",nconfstr);
194 12 : buffer1[0]='\0';
195 12 : if(std::sscanf(nconfstr.c_str(),"%100d %255s",&nconfig,buffer1)!=2 || nconfig<=0) {
196 2 : error("nconfig should be followed by a positive stride and by the name of the trajectory file");
197 : }
198 : trajfile=buffer1;
199 : std::string nstatstr;
200 11 : parse("nstat",nstatstr);
201 11 : buffer1[0]='\0';
202 11 : if(std::sscanf(nstatstr.c_str(),"%100d %255s",&nstat,buffer1)!=2 || nstat<=0) {
203 0 : error("nstat should be followed by a positive stride and by the name of the statistics file");
204 : }
205 : statfile=buffer1;
206 11 : parse("ndim",ndim);
207 11 : if(ndim<1 || ndim>3) {
208 0 : error("ndim should be 1,2 or 3");
209 : }
210 : std::string w;
211 11 : parse("wrapatoms",w);
212 11 : wrapatoms=false;
213 11 : if(w.length()>0 && (w[0]=='T' || w[0]=='t')) {
214 2 : wrapatoms=true;
215 : }
216 11 : }
217 :
218 11 : void read_natoms(const std::string & inputfile,int & natoms) {
219 : // read the number of atoms in file "input.xyz"
220 11 : FILE* fp=std::fopen(inputfile.c_str(),"r");
221 11 : if(!fp) {
222 0 : error(std::string("file ") + inputfile + std::string(" not found"));
223 : }
224 :
225 : // call fclose when fp goes out of scope
226 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
227 :
228 11 : int ret=std::fscanf(fp,"%1000d",&natoms);
229 11 : if(ret!=1 || natoms<=0) {
230 0 : plumed_error() <<"Error reading number of atoms from file "<<inputfile;
231 : }
232 11 : }
233 :
234 11 : void read_positions(const std::string& inputfile,int natoms,std::vector<Vector>& positions,double cell[3]) {
235 : // read positions and cell from a file called inputfile
236 : // natoms (input variable) and number of atoms in the file should be consistent
237 11 : FILE* fp=std::fopen(inputfile.c_str(),"r");
238 11 : if(!fp) {
239 0 : error(std::string("file ") + inputfile + std::string(" not found"));
240 : }
241 : // call fclose when fp goes out of scope
242 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
243 :
244 : char buffer[256];
245 : char atomname[256];
246 : char* cret=fgets(buffer,256,fp);
247 11 : if(cret==nullptr) {
248 0 : plumed_error() <<"Error reading buffer from file "<<inputfile;
249 : }
250 11 : int ret=std::fscanf(fp,"%1000lf %1000lf %1000lf",&cell[0],&cell[1],&cell[2]);
251 11 : if(ret!=3) {
252 0 : plumed_error() <<"Error reading cell line from file "<<inputfile;
253 : }
254 1206 : for(int i=0; i<natoms; i++) {
255 1195 : ret=std::fscanf(fp,"%255s %1000lf %1000lf %1000lf",atomname,&positions[i][0],&positions[i][1],&positions[i][2]);
256 : // note: atomname is read but not used
257 1195 : if(ret!=4) {
258 0 : plumed_error() <<"Error reading atom line from file "<<inputfile;
259 : }
260 : }
261 11 : }
262 :
263 11 : void randomize_velocities(const int natoms,const int ndim,const double temperature,const std::vector<double>&masses,std::vector<Vector>& velocities,Random&random) {
264 : // randomize the velocities according to the temperature
265 1206 : for(int iatom=0; iatom<natoms; iatom++)
266 4773 : for(int i=0; i<ndim; i++) {
267 3578 : velocities[iatom][i]=std::sqrt(temperature/masses[iatom])*random.Gaussian();
268 : }
269 11 : }
270 :
271 10690043 : void pbc(const double cell[3],const Vector & vin,Vector & vout) {
272 : // apply periodic boundary condition to a vector
273 42760172 : for(int i=0; i<3; i++) {
274 32070129 : vout[i]=vin[i]-std::floor(vin[i]/cell[i]+0.5)*cell[i];
275 : }
276 10690043 : }
277 :
278 3900 : void check_list(const int natoms,const std::vector<Vector>& positions,const std::vector<Vector>&positions0,const double listcutoff,
279 : const double forcecutoff,bool & recompute) {
280 : // check if the neighbour list have to be recomputed
281 3900 : recompute=false;
282 3900 : auto delta2=(0.5*(listcutoff-forcecutoff))*(0.5*(listcutoff-forcecutoff));
283 : // if ANY atom moved more than half of the skin thickness, recompute is set to .true.
284 192715 : for(int iatom=0; iatom<natoms; iatom++) {
285 189245 : if(modulo2(positions[iatom]-positions0[iatom])>delta2) {
286 430 : recompute=true;
287 430 : break;
288 : }
289 : }
290 3900 : }
291 :
292 :
293 441 : void compute_list(const int natoms,const std::vector<Vector>& positions,const double cell[3],const double listcutoff,
294 : std::vector<std::vector<int> >& list) {
295 441 : double listcutoff2=listcutoff*listcutoff; // squared list cutoff
296 441 : list.assign(natoms,std::vector<int>());
297 : // linked cells: they are used when there are at least 3 cells in every direction
298 : // (so that the 27 neighboring cells are distinct). Cells are slightly larger
299 : // than listcutoff so that rounding cannot miss a pair. The list of each atom is then
300 : // sorted, so that the list (and thus the forces) are identical to the O(N^2) search
301 : int ncells[3];
302 : bool usecells=true;
303 1764 : for(int k=0; k<3; k++) {
304 1323 : const double nc=std::floor(cell[k]/(listcutoff*(1.0+1e-6)));
305 1323 : ncells[k]=(nc<64?int(nc):64); // at most 64^3 cells, larger cells are still correct
306 1323 : if(ncells[k]<3) {
307 : usecells=false;
308 : }
309 : }
310 441 : if(usecells) {
311 25 : const int ntot=ncells[0]*ncells[1]*ncells[2];
312 25 : std::vector<int> atomcell(natoms);
313 25 : std::vector<int> cellstart(ntot+1,0);
314 25 : std::vector<int> cellatoms(natoms);
315 409 : for(int iatom=0; iatom<natoms; iatom++) {
316 : int c[3];
317 1536 : for(int k=0; k<3; k++) {
318 1152 : double sc=positions[iatom][k]/cell[k];
319 1152 : sc-=std::floor(sc);
320 1152 : c[k]=int(sc*ncells[k]);
321 1152 : if(c[k]>=ncells[k]) {
322 0 : c[k]=ncells[k]-1;
323 : }
324 1152 : if(c[k]<0) {
325 0 : c[k]=0;
326 : }
327 : }
328 384 : atomcell[iatom]=(c[0]*ncells[1]+c[1])*ncells[2]+c[2];
329 384 : cellstart[atomcell[iatom]+1]++;
330 : }
331 862540 : for(int ic=0; ic<ntot; ic++) {
332 862515 : cellstart[ic+1]+=cellstart[ic];
333 : }
334 : {
335 25 : std::vector<int> fill(cellstart.begin(),cellstart.end()-1);
336 409 : for(int iatom=0; iatom<natoms; iatom++) {
337 384 : cellatoms[fill[atomcell[iatom]]++]=iatom;
338 : }
339 25 : }
340 25 : # pragma omp parallel for num_threads(OpenMP::getNumThreads()) schedule(static,1)
341 : for(int iatom=0; iatom<natoms-1; iatom++) {
342 : const int ic=atomcell[iatom];
343 : const int c0=ic/(ncells[1]*ncells[2]);
344 : const int c1=(ic/ncells[2])%ncells[1];
345 : const int c2=ic%ncells[2];
346 : for(int d0=-1; d0<=1; d0++)
347 : for(int d1=-1; d1<=1; d1++)
348 : for(int d2=-1; d2<=1; d2++) {
349 : const int jc=(((c0+d0+ncells[0])%ncells[0])*ncells[1]+(c1+d1+ncells[1])%ncells[1])*ncells[2]+(c2+d2+ncells[2])%ncells[2];
350 : for(int jj=cellstart[jc]; jj<cellstart[jc+1]; jj++) {
351 : const int jatom=cellatoms[jj];
352 : if(jatom<=iatom) {
353 : continue;
354 : }
355 : auto distance=positions[iatom]-positions[jatom];
356 : Vector distance_pbc; // minimum-image distance of the two atoms
357 : pbc(cell,distance,distance_pbc);
358 : if(modulo2(distance_pbc)>listcutoff2) {
359 : continue;
360 : }
361 : list[iatom].push_back(jatom);
362 : }
363 : }
364 : std::sort(list[iatom].begin(),list[iatom].end());
365 : }
366 : return;
367 25 : }
368 416 : # pragma omp parallel for num_threads(OpenMP::getNumThreads()) schedule(static,1)
369 : for(int iatom=0; iatom<natoms-1; iatom++) {
370 : for(int jatom=iatom+1; jatom<natoms; jatom++) {
371 : auto distance=positions[iatom]-positions[jatom];
372 : Vector distance_pbc; // minimum-image distance of the two atoms
373 : pbc(cell,distance,distance_pbc);
374 : // if the interparticle distance is larger than the cutoff, skip
375 : if(modulo2(distance_pbc)>listcutoff2) {
376 : continue;
377 : }
378 : list[iatom].push_back(jatom);
379 : }
380 : }
381 : }
382 :
383 3911 : void compute_forces(const int natoms,double epsilon, double sigma,const std::vector<Vector>& positions,const double cell[3],
384 : double forcecutoff,const std::vector<std::vector<int> >& list,std::vector<Vector>& forces,double & engconf) {
385 3911 : double forcecutoff2=forcecutoff*forcecutoff; // squared force cutoff
386 3911 : 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
387 : // one buffer (forces and energy) per thread: they are summed in thread order after the loop,
388 : // so that the result does not depend on the order in which the threads finish
389 3911 : const unsigned nt=OpenMP::getNumThreads();
390 3911 : std::vector<std::vector<Vector>> omp_forces_all(nt);
391 3911 : std::vector<double> omp_engconf_all(nt,0.0);
392 3911 : # pragma omp parallel num_threads(nt)
393 : {
394 : const unsigned tid=OpenMP::getThreadNum();
395 : std::vector<Vector>& omp_forces(omp_forces_all[tid]);
396 : omp_forces.assign(forces.size(),Vector(0.0,0.0,0.0));
397 : double omp_engconf=0.0;
398 : #pragma omp for schedule(static,1)
399 : for(int iatom=0; iatom<natoms-1; iatom++) {
400 : for(unsigned jlist=0; jlist<list[iatom].size(); jlist++) {
401 : const int jatom=list[iatom][jlist];
402 : auto distance=positions[iatom]-positions[jatom];
403 : Vector distance_pbc; // minimum-image distance of the two atoms
404 : pbc(cell,distance,distance_pbc);
405 : auto distance_pbc2=modulo2(distance_pbc); // squared minimum-image distance
406 : // if the interparticle distance is larger than the cutoff, skip
407 : if(distance_pbc2>forcecutoff2) {
408 : continue;
409 : }
410 : auto distance_pbcm2=sigma*sigma/distance_pbc2;
411 : auto distance_pbcm6=distance_pbcm2*distance_pbcm2*distance_pbcm2;
412 : auto distance_pbcm8=distance_pbcm6*distance_pbcm2;
413 : auto distance_pbcm12=distance_pbcm6*distance_pbcm6;
414 : auto distance_pbcm14=distance_pbcm12*distance_pbcm2;
415 : omp_engconf+=4.0*epsilon*(distance_pbcm12 - distance_pbcm6) - engcorrection;
416 : // F = 24 epsilon (2 sigma^12/r^14 - sigma^6/r^8) r = 24 epsilon (2 m14 - m8) r / sigma^2, with m=sigma^2/r^2
417 : auto f=24.0*distance_pbc*(2.0*distance_pbcm14-distance_pbcm8)*epsilon/(sigma*sigma);
418 : omp_forces[iatom]+=f;
419 : omp_forces[jatom]-=f;
420 : }
421 : }
422 : omp_engconf_all[tid]=omp_engconf;
423 : }
424 3911 : engconf=0.0;
425 235106 : for(int i=0; i<natoms; i++) {
426 231195 : forces[i].zero();
427 : }
428 10230 : for(unsigned t=0; t<nt; t++) {
429 6319 : if(omp_forces_all[t].empty()) {
430 0 : continue; // this thread did not run
431 : }
432 6319 : engconf+=omp_engconf_all[t];
433 306385 : for(int i=0; i<natoms; i++) {
434 300066 : forces[i]+=omp_forces_all[t][i];
435 : }
436 : }
437 3911 : }
438 :
439 3900 : void compute_engkin(const int natoms,const std::vector<double>& masses,const std::vector<Vector>& velocities,double & engkin) {
440 : // calculate the kinetic energy from the velocities
441 3900 : engkin=0.0;
442 233900 : for(int iatom=0; iatom<natoms; iatom++) {
443 230000 : engkin+=0.5*masses[iatom]*modulo2(velocities[iatom]);
444 : }
445 3900 : }
446 :
447 :
448 7800 : void thermostat(const int natoms,const int ndim,const std::vector<double>& masses,const double dt,const double friction,
449 : const double temperature,std::vector<Vector>& velocities,double & engint,Random & random) {
450 : // Langevin thermostat, implemented as described in Bussi and Parrinello, Phys. Rev. E (2007)
451 : // it is a linear combination of old velocities and new, randomly chosen, velocity,
452 : // with proper coefficients
453 7800 : double c1=std::exp(-friction*dt);
454 467800 : for(int iatom=0; iatom<natoms; iatom++) {
455 460000 : double c2=std::sqrt((1.0-c1*c1)*temperature/masses[iatom]);
456 1812000 : for(int i=0; i<ndim; i++) {
457 1352000 : engint+=0.5*masses[iatom]*velocities[iatom][i]*velocities[iatom][i];
458 1352000 : velocities[iatom][i]=c1*velocities[iatom][i]+c2*random.Gaussian();
459 1352000 : engint-=0.5*masses[iatom]*velocities[iatom][i]*velocities[iatom][i];
460 : }
461 : }
462 7800 : }
463 :
464 49 : void write_positions(const std::string& trajfile,int natoms,const std::vector<Vector>& positions,const double cell[3],const bool wrapatoms) {
465 : // write positions on file trajfile
466 : // positions are appended at the end of the file
467 49 : Vector pos;
468 : FILE*fp;
469 49 : if(write_positions_first) {
470 11 : fp=std::fopen(trajfile.c_str(),"w");
471 11 : write_positions_first=false;
472 : } else {
473 38 : fp=std::fopen(trajfile.c_str(),"a");
474 : }
475 49 : plumed_assert(fp)<<"cannot open file "<<trajfile;
476 : // call fclose when fp goes out of scope
477 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
478 : std::fprintf(fp,"%d\n",natoms);
479 49 : std::fprintf(fp,"%f %f %f\n",cell[0],cell[1],cell[2]);
480 5355 : for(int iatom=0; iatom<natoms; iatom++) {
481 : // usually, it is better not to apply pbc here, so that diffusion
482 : // is more easily calculated from a trajectory file:
483 5306 : if(wrapatoms) {
484 1080 : pbc(cell,positions[iatom],pos);
485 : } else {
486 4226 : pos=positions[iatom];
487 : }
488 5306 : std::fprintf(fp,"Ar %10.7f %10.7f %10.7f\n",pos[0],pos[1],pos[2]);
489 : }
490 49 : }
491 :
492 11 : void write_final_positions(const std::string& outputfile,int natoms,const std::vector<Vector>& positions,const double cell[3],const bool wrapatoms) {
493 : // write positions on file outputfile
494 11 : Vector pos;
495 : FILE*fp;
496 11 : fp=std::fopen(outputfile.c_str(),"w");
497 11 : plumed_assert(fp)<<"cannot open file "<<outputfile;
498 : // call fclose when fp goes out of scope
499 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
500 : std::fprintf(fp,"%d\n",natoms);
501 11 : std::fprintf(fp,"%f %f %f\n",cell[0],cell[1],cell[2]);
502 1206 : for(int iatom=0; iatom<natoms; iatom++) {
503 : // usually, it is better not to apply pbc here, so that diffusion
504 : // is more easily calculated from a trajectory file:
505 1195 : if(wrapatoms) {
506 216 : pbc(cell,positions[iatom],pos);
507 : } else {
508 979 : pos=positions[iatom];
509 : }
510 1195 : std::fprintf(fp,"Ar %10.7f %10.7f %10.7f\n",pos[0],pos[1],pos[2]);
511 : }
512 11 : }
513 :
514 :
515 192 : void write_statistics(const std::string & statfile,const int istep,const double tstep,
516 : const int natoms,const int ndim,const double engkin,const double engconf,const double engint) {
517 : // write statistics on file statfile
518 192 : if(write_statistics_first) {
519 : // first time this routine is called, open the file
520 11 : write_statistics_fp=std::fopen(statfile.c_str(),"w");
521 11 : plumed_assert(write_statistics_fp)<<"cannot open file "<<statfile;
522 : write_statistics_fp_deleter.reset(write_statistics_fp);
523 11 : write_statistics_first=false;
524 : }
525 192 : if(istep-write_statistics_last_time_reopened>100) {
526 : // every 100 steps, reopen the file to flush the buffer
527 : write_statistics_fp_deleter.reset(nullptr); // close file
528 14 : write_statistics_fp=std::fopen(statfile.c_str(),"a");
529 14 : plumed_assert(write_statistics_fp)<<"cannot open file "<<statfile;
530 : write_statistics_fp_deleter.reset(write_statistics_fp);
531 14 : write_statistics_last_time_reopened=istep;
532 : }
533 192 : 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);
534 192 : }
535 :
536 :
537 :
538 14 : int main(FILE* in,FILE*out,PLMD::Communicator& pc) override {
539 : int natoms; // number of atoms
540 14 : std::vector<Vector> positions; // atomic positions
541 14 : std::vector<Vector> velocities; // velocities
542 14 : std::vector<double> masses; // masses
543 14 : std::vector<Vector> forces; // forces
544 : double cell[3]; // cell size
545 : double cell9[3][3]; // cell size
546 :
547 : // neighbour list variables
548 : // see Allen and Tildesey book for details
549 14 : std::vector< std::vector<int> > list; // neighbour list
550 14 : std::vector<Vector> positions0; // reference atomic positions, i.e. positions when the neighbour list
551 :
552 : // input parameters
553 : // all of them have a reasonable default value, set in read_input()
554 : double tstep; // simulation timestep
555 : double temperature; // temperature
556 : double friction; // friction for Langevin dynamics (for NVE, use 0)
557 : double listcutoff; // cutoff for neighbour list
558 : double forcecutoff; // cutoff for forces
559 : int nstep; // number of steps
560 : int nconfig; // stride for output of configurations
561 : int nstat; // stride for output of statistics
562 : int ndim; // dimensionality of the system (1, 2, or 3)
563 : int idum; // seed
564 : int plumedWantsToStop; // stop flag
565 : bool wrapatoms; // if true, atomic coordinates are written wrapped in minimal cell
566 : std::string inputfile; // name of file with starting configuration (xyz)
567 : std::string outputfile; // name of file with final configuration (xyz)
568 : std::string trajfile; // name of the trajectory file (xyz)
569 : std::string statfile; // name of the file with statistics
570 :
571 : double epsilon, sigma; // LJ parameters
572 :
573 : double engkin; // kinetic energy
574 : double engconf; // configurational energy
575 : double engint; // integral for conserved energy in Langevin dynamics
576 :
577 : bool recompute_list; // control if the neighbour list have to be recomputed
578 :
579 14 : Random random; // random numbers stream
580 :
581 14 : std::unique_ptr<PlumedMain> plumed;
582 :
583 : // Commenting the next line it is possible to switch-off plumed
584 17 : plumed=Tools::make_unique<PLMD::PlumedMain>();
585 :
586 14 : if(plumed) {
587 14 : plumed->cmd("setRealPrecision",int(sizeof(double)));
588 : }
589 :
590 14 : read_input(temperature,tstep,friction,forcecutoff,
591 : listcutoff,nstep,nconfig,nstat,
592 : wrapatoms,inputfile,outputfile,trajfile,statfile,
593 : ndim,idum,epsilon,sigma);
594 :
595 : // number of atoms is read from file inputfile
596 11 : read_natoms(inputfile,natoms);
597 :
598 : // write the parameters in output so they can be checked
599 : std::fprintf(out,"%s %s\n","Starting configuration :",inputfile.c_str());
600 : std::fprintf(out,"%s %s\n","Final configuration :",outputfile.c_str());
601 11 : std::fprintf(out,"%s %d\n","Number of atoms :",natoms);
602 11 : std::fprintf(out,"%s %f\n","Temperature :",temperature);
603 11 : std::fprintf(out,"%s %f\n","Time step :",tstep);
604 11 : std::fprintf(out,"%s %f\n","Friction :",friction);
605 11 : std::fprintf(out,"%s %f\n","Cutoff for forces :",forcecutoff);
606 11 : std::fprintf(out,"%s %f\n","Cutoff for neighbour list :",listcutoff);
607 11 : std::fprintf(out,"%s %d\n","Number of steps :",nstep);
608 11 : std::fprintf(out,"%s %d\n","Stride for trajectory :",nconfig);
609 : std::fprintf(out,"%s %s\n","Trajectory file :",trajfile.c_str());
610 11 : std::fprintf(out,"%s %d\n","Stride for statistics :",nstat);
611 : std::fprintf(out,"%s %s\n","Statistics file :",statfile.c_str());
612 11 : std::fprintf(out,"%s %d\n","Dimensionality :",ndim);
613 11 : std::fprintf(out,"%s %d\n","Seed :",idum);
614 11 : std::fprintf(out,"%s %s\n","Are atoms wrapped on output? :",(wrapatoms?"T":"F"));
615 11 : std::fprintf(out,"%s %f\n","Epsilon :",epsilon);
616 11 : std::fprintf(out,"%s %f\n","Sigma :",sigma);
617 :
618 : // Setting the seed
619 11 : random.setSeed(idum);
620 :
621 : // allocation of dynamical arrays
622 11 : positions.resize(natoms);
623 11 : positions0.resize(natoms);
624 11 : velocities.resize(natoms);
625 11 : forces.resize(natoms);
626 11 : masses.resize(natoms);
627 11 : list.resize(natoms);
628 :
629 : // masses are hard-coded to 1
630 1206 : for(int i=0; i<natoms; ++i) {
631 1195 : masses[i]=1.0;
632 : }
633 :
634 : // energy integral initialized to 0
635 11 : engint=0.0;
636 :
637 : // positions are read from file inputfile
638 11 : read_positions(inputfile,natoms,positions,cell);
639 :
640 : // velocities are randomized according to temperature
641 11 : randomize_velocities(natoms,ndim,temperature,masses,velocities,random);
642 :
643 11 : if(plumed) {
644 11 : plumed->cmd("setNatoms",natoms);
645 11 : plumed->cmd("setNoVirial");
646 11 : plumed->cmd("setMDEngine","simpleMD");
647 11 : plumed->cmd("setTimestep",tstep);
648 11 : plumed->cmd("setPlumedDat","plumed.dat");
649 11 : int pversion=0;
650 11 : plumed->cmd("getApiVersion",&pversion);
651 : // setting kbT is only implemented with api>1
652 : // even if not necessary in principle in SimpleMD (which is part of plumed)
653 : // we leave the check here as a reference
654 11 : if(pversion>1 && temperature>0.0) {
655 10 : plumed->cmd("setKbT",temperature);
656 : }
657 11 : plumed->cmd("init");
658 : }
659 :
660 : // neighbour list are computed, and reference positions are saved
661 11 : compute_list(natoms,positions,cell,listcutoff,list);
662 : {
663 : int list_size=0;
664 1206 : for(const auto &element: list) {
665 1195 : list_size+=element.size();
666 : }
667 : std::fprintf(out,"List size: %d\n",list_size);
668 : }
669 1206 : for(int iatom=0; iatom<natoms; ++iatom) {
670 1195 : positions0[iatom]=positions[iatom];
671 : }
672 :
673 : // forces are computed before starting md
674 11 : compute_forces(natoms,epsilon,sigma,positions,cell,forcecutoff,list,forces,engconf);
675 :
676 : // remove forces if ndim<3
677 11 : if(ndim<3)
678 8 : for(int iatom=0; iatom<natoms; ++iatom)
679 14 : for(int k=ndim; k<3; ++k) {
680 7 : forces[iatom][k]=0.0;
681 : }
682 :
683 : // here is the main md loop
684 : // Langevin thermostat is applied before and after a velocity-Verlet integrator
685 : // the overall structure is:
686 : // thermostat
687 : // update velocities
688 : // update positions
689 : // (eventually recompute neighbour list)
690 : // compute forces
691 : // update velocities
692 : // thermostat
693 : // (eventually dump output information)
694 3911 : for(int istep=0; istep<nstep; istep++) {
695 3900 : thermostat(natoms,ndim,masses,0.5*tstep,friction,temperature,velocities,engint,random);
696 :
697 233900 : for(int iatom=0; iatom<natoms; iatom++) {
698 230000 : velocities[iatom]+=forces[iatom]*0.5*tstep/masses[iatom];
699 : }
700 :
701 233900 : for(int iatom=0; iatom<natoms; iatom++) {
702 230000 : positions[iatom]+=velocities[iatom]*tstep;
703 : }
704 :
705 : // a check is performed to decide whether to recalculate the neighbour list
706 3900 : check_list(natoms,positions,positions0,listcutoff,forcecutoff,recompute_list);
707 3900 : if(recompute_list) {
708 430 : compute_list(natoms,positions,cell,listcutoff,list);
709 44547 : for(int iatom=0; iatom<natoms; ++iatom) {
710 44117 : positions0[iatom]=positions[iatom];
711 : }
712 : int list_size=0;
713 44547 : for(const auto &element: list) {
714 44117 : list_size+=element.size();
715 : }
716 : std::fprintf(out,"List size: %d\n",list_size);
717 : }
718 :
719 3900 : compute_forces(natoms,epsilon,sigma,positions,cell,forcecutoff,list,forces,engconf);
720 :
721 3900 : if(plumed) {
722 3900 : plumedWantsToStop=0;
723 15600 : for(int i=0; i<3; i++)
724 46800 : for(int k=0; k<3; k++) {
725 35100 : cell9[i][k]=0.0;
726 : }
727 15600 : for(int i=0; i<3; i++) {
728 11700 : cell9[i][i]=cell[i];
729 : }
730 3900 : plumed->cmd("setStep",istep+1);
731 3900 : plumed->cmd("setMasses",&masses[0]);
732 3900 : plumed->cmd("setForces",&forces[0][0]);
733 3900 : plumed->cmd("setEnergy",engconf);
734 3900 : plumed->cmd("setPositions",&positions[0][0]);
735 3900 : plumed->cmd("setBox",&cell9[0][0]);
736 3900 : plumed->cmd("setStopFlag",&plumedWantsToStop);
737 3900 : plumed->cmd("calc");
738 3900 : if(plumedWantsToStop) {
739 0 : nstep=istep;
740 : }
741 : }
742 : // remove forces if ndim<3
743 3900 : if(ndim<3)
744 16000 : for(int iatom=0; iatom<natoms; ++iatom)
745 28000 : for(int k=ndim; k<3; ++k) {
746 14000 : forces[iatom][k]=0.0;
747 : }
748 :
749 233900 : for(int iatom=0; iatom<natoms; iatom++) {
750 230000 : velocities[iatom]+=forces[iatom]*0.5*tstep/masses[iatom];
751 : }
752 :
753 3900 : thermostat(natoms,ndim,masses,0.5*tstep,friction,temperature,velocities,engint,random);
754 :
755 : // kinetic energy is calculated
756 3900 : compute_engkin(natoms,masses,velocities,engkin);
757 :
758 : // eventually, write positions and statistics
759 3900 : if((istep+1)%nconfig==0) {
760 49 : write_positions(trajfile,natoms,positions,cell,wrapatoms);
761 : }
762 3900 : if((istep+1)%nstat==0) {
763 192 : write_statistics(statfile,istep+1,tstep,natoms,ndim,engkin,engconf,engint);
764 : }
765 :
766 : }
767 :
768 : // call final plumed jobs
769 11 : if(plumed) {
770 11 : plumed->cmd("runFinalJobs");
771 : }
772 :
773 : // write final positions
774 11 : write_final_positions(outputfile,natoms,positions,cell,wrapatoms);
775 :
776 11 : return 0;
777 51 : }
778 :
779 :
780 : };
781 :
782 21984 : PLUMED_REGISTER_CLTOOL(SimpleMD,"simplemd")
783 :
784 : }
785 : }
786 :
787 :
788 :
789 :
|
