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