Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module 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 : The VES code module 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 the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "BasisFunctions.h"
24 : #include "LinearBasisSetExpansion.h"
25 : #include "CoeffsVector.h"
26 : #include "GridIntegrationWeights.h"
27 :
28 : #include "cltools/CLTool.h"
29 : #include "cltools/CLToolRegister.h"
30 : #include "tools/Vector.h"
31 : #include "tools/Random.h"
32 : #include "tools/Grid.h"
33 : #include "tools/Communicator.h"
34 : #include "tools/FileBase.h"
35 : #include "core/PlumedMain.h"
36 : #include "core/ActionRegister.h"
37 : #include "core/ActionSet.h"
38 : #include "core/Value.h"
39 :
40 : #include <string>
41 : #include <cstdio>
42 : #include <cmath>
43 : #include <vector>
44 : #include <iostream>
45 :
46 : #ifdef __PLUMED_HAS_MPI
47 : #include <mpi.h>
48 : #endif
49 :
50 :
51 : namespace PLMD {
52 : namespace ves {
53 :
54 : //+PLUMEDOC VES_TOOLS ves_md_linearexpansion
55 : /*
56 : Simple MD code for dynamics on a potential energy surface given by a linear basis set expansion.
57 :
58 : This is simple MD code that allows running dynamics of a single particle on a
59 : potential energy surface given by some linear basis set expansion in one to three
60 : dimensions.
61 :
62 : It is possible to run more than one replica of the system in parallel.
63 :
64 : \par Examples
65 :
66 : In the following example we perform dynamics on the
67 : Wolfe-Quapp potential that is defined as
68 : \f[
69 : U(x,y) = x^4 + y^4 - 2 x^2 - 4 y^2 + xy + 0.3 x + 0.1 y
70 : \f]
71 : To define the potential we employ polynomial power basis
72 : functions (\ref BF_POWERS). The input file is given as
73 : \verbatim
74 : nstep 10000
75 : tstep 0.005
76 : temperature 1.0
77 : friction 10.0
78 : random_seed 4525
79 : plumed_input plumed.dat
80 : dimension 2
81 : replicas 1
82 : basis_functions_1 BF_POWERS ORDER=4 MINIMUM=-3.0 MAXIMUM=+3.0
83 : basis_functions_2 BF_POWERS ORDER=4 MINIMUM=-3.0 MAXIMUM=+3.0
84 : input_coeffs pot_coeffs_input.data
85 : initial_position -1.174,+1.477
86 : output_potential potential.data
87 : output_potential_grid 150
88 : output_histogram histogram.data
89 :
90 : # Wolfe-Quapp potential given by the equation
91 : # U(x,y) = x**4 + y**4 - 2.0*x**2 - 4.0*y**2 + x*y + 0.3*x + 0.1*y
92 : # Minima around (-1.174,1.477); (-0.831,-1.366); (1.124,-1.486)
93 : # Maxima around (0.100,0.050)
94 : # Saddle points around (-1.013,-0.036); (0.093,0.174); (-0.208,-1.407)
95 : \endverbatim
96 :
97 : This input is then run by using the following command.
98 : \verbatim
99 : plumed ves_md_linearexpansion input
100 : \endverbatim
101 :
102 : The corresponding pot_coeffs_input.data file is
103 : \verbatim
104 : #! FIELDS idx_dim1 idx_dim2 pot.coeffs index description
105 : #! SET type LinearBasisSet
106 : #! SET ndimensions 2
107 : #! SET ncoeffs_total 25
108 : #! SET shape_dim1 5
109 : #! SET shape_dim2 5
110 : 0 0 0.0000000000000000e+00 0 1*1
111 : 1 0 0.3000000000000000e+00 1 s^1*1
112 : 2 0 -2.0000000000000000e+00 2 s^2*1
113 : 4 0 1.0000000000000000e+00 4 s^4*1
114 : 0 1 0.1000000000000000e+00 5 1*s^1
115 : 1 1 +1.0000000000000000e+00 6 s^1*s^1
116 : 0 2 -4.0000000000000000e+00 10 1*s^2
117 : 0 4 1.0000000000000000e+00 20 1*s^4
118 : #!-------------------
119 : \endverbatim
120 :
121 : One then uses the (x,y) position of the particle as CVs by using the \ref POSITION
122 : action as shown in the following PLUMED input
123 : \plumedfile
124 : p: POSITION ATOM=1
125 : ene: ENERGY
126 : PRINT ARG=p.x,p.y,ene FILE=colvar.data FMT=%8.4f
127 : \endplumedfile
128 :
129 :
130 :
131 : */
132 : //+ENDPLUMEDOC
133 :
134 38 : class MD_LinearExpansionPES : public PLMD::CLTool {
135 : public:
136 0 : std::string description() const {return "MD of a one particle on a linear expansion PES";}
137 : static void registerKeywords( Keywords& keys );
138 : explicit MD_LinearExpansionPES( const CLToolOptions& co );
139 : int main( FILE* in, FILE* out, PLMD::Communicator& pc);
140 : private:
141 : unsigned int dim;
142 : LinearBasisSetExpansion* potential_expansion_pntr;
143 : //
144 : double calc_energy( const std::vector<Vector>&, std::vector<Vector>& );
145 : double calc_temp( const std::vector<Vector>& );
146 : };
147 :
148 7432 : PLUMED_REGISTER_CLTOOL(MD_LinearExpansionPES,"ves_md_linearexpansion")
149 :
150 1839 : void MD_LinearExpansionPES::registerKeywords( Keywords& keys ) {
151 1839 : CLTool::registerKeywords( keys );
152 9195 : keys.add("compulsory","nstep","10","The number of steps of dynamics you want to run.");
153 9195 : keys.add("compulsory","tstep","0.005","The integration timestep.");
154 9195 : keys.add("compulsory","temperature","1.0","The temperature to perform the simulation at. For multiple replica you can give a separate value for each replica.");
155 9195 : keys.add("compulsory","friction","10.","The friction of the Langevin thermostat. For multiple replica you can give a separate value for each replica.");
156 9195 : keys.add("compulsory","random_seed","5293818","Value of random number seed.");
157 9195 : keys.add("compulsory","plumed_input","plumed.dat","The name of the plumed input file(s). For multiple replica you can give a separate value for each replica.");
158 9195 : keys.add("compulsory","dimension","1","Number of dimensions, supports 1 to 3.");
159 7356 : keys.add("compulsory","initial_position","Initial position of the particle. For multiple replica you can give a separate value for each replica.");
160 9195 : keys.add("compulsory","replicas","1","Number of replicas.");
161 7356 : keys.add("compulsory","basis_functions_1","Basis functions for dimension 1.");
162 7356 : keys.add("optional","basis_functions_2","Basis functions for dimension 2 if needed.");
163 7356 : keys.add("optional","basis_functions_3","Basis functions for dimension 3 if needed.");
164 9195 : keys.add("compulsory","input_coeffs","potential-coeffs.in.data","Filename of the input coefficient file for the potential. For multiple replica you can give a separate value for each replica.");
165 9195 : keys.add("compulsory","output_coeffs","potential-coeffs.out.data","Filename of the output coefficient file for the potential.");
166 9195 : keys.add("compulsory","output_coeffs_fmt","%30.16e","Format of the output coefficient file for the potential. Useful for regtests.");
167 7356 : keys.add("optional","coeffs_prefactor","prefactor for multiplying the coefficients with. For multiple replica you can give a separate value for each replica.");
168 7356 : keys.add("optional","template_coeffs_file","only generate a template coefficient file with the filename given and exit.");
169 9195 : keys.add("compulsory","output_potential_grid","100","The number of grid points used for the potential and histogram output files.");
170 9195 : keys.add("compulsory","output_potential","potential.data","Filename of the potential output file.");
171 9195 : keys.add("compulsory","output_histogram","histogram.data","Filename of the histogram output file.");
172 1839 : }
173 :
174 :
175 38 : MD_LinearExpansionPES::MD_LinearExpansionPES( const CLToolOptions& co ):
176 : CLTool(co),
177 : dim(0),
178 38 : potential_expansion_pntr(NULL)
179 : {
180 38 : inputdata=ifile; //commandline;
181 38 : }
182 :
183 : inline
184 3838 : double MD_LinearExpansionPES::calc_energy( const std::vector<Vector>& pos, std::vector<Vector>& forces) {
185 3838 : std::vector<double> pos_tmp(dim);
186 3838 : std::vector<double> forces_tmp(dim,0.0);
187 12726 : for(unsigned int j=0; j<dim; ++j) {
188 8888 : pos_tmp[j]=pos[0][j];
189 : }
190 3838 : bool all_inside = true;
191 3838 : double potential = potential_expansion_pntr->getBiasAndForces(pos_tmp,all_inside,forces_tmp);
192 12726 : for(unsigned int j=0; j<dim; ++j) {
193 8888 : forces[0][j] = forces_tmp[j];
194 : }
195 3838 : return potential;
196 : }
197 :
198 :
199 : inline
200 3838 : double MD_LinearExpansionPES::calc_temp( const std::vector<Vector>& vel) {
201 : double total_KE=0.0;
202 : //! Double the total kinetic energy of the system
203 12726 : for(unsigned int j=0; j<dim; ++j) {
204 4444 : total_KE+=vel[0][j]*vel[0][j];
205 : }
206 3838 : return total_KE / (double) dim; // total_KE is actually 2*KE
207 : }
208 :
209 38 : int MD_LinearExpansionPES::main( FILE* in, FILE* out, PLMD::Communicator& pc) {
210 : int plumedWantsToStop;
211 38 : Random random;
212 : unsigned int stepWrite=1000;
213 :
214 : PLMD::PlumedMain* plumed=NULL;
215 : PLMD::PlumedMain* plumed_bf=NULL;
216 :
217 : unsigned int replicas;
218 : unsigned int coresPerReplica;
219 76 : parse("replicas",replicas);
220 38 : if(replicas==1) {
221 7 : coresPerReplica = pc.Get_size();
222 : } else {
223 31 : if(pc.Get_size()%replicas!=0) {
224 0 : error("the number of MPI processes is not a multiple of the number of replicas.");
225 : }
226 31 : coresPerReplica = pc.Get_size()/replicas;
227 : }
228 : // create intra and inter communicators
229 76 : Communicator intra, inter;
230 38 : if(Communicator::initialized()) {
231 33 : int iworld=(pc.Get_rank() / coresPerReplica);
232 33 : pc.Split(iworld,0,intra);
233 33 : pc.Split(intra.Get_rank(),0,inter);
234 : }
235 :
236 : unsigned int nsteps;
237 76 : parse("nstep",nsteps);
238 : double tstep;
239 76 : parse("tstep",tstep);
240 : // initialize to solve a cppcheck 1.86 warning
241 38 : double temp=0.0;
242 38 : std::vector<double> temps_vec(0);
243 76 : parseVector("temperature",temps_vec);
244 38 : if(temps_vec.size()==1) {
245 34 : temp = temps_vec[0];
246 : }
247 4 : else if(replicas > 1 && temps_vec.size()==replicas) {
248 12 : temp = temps_vec[inter.Get_rank()];
249 : }
250 : else {
251 0 : error("problem with temperature keyword, you need to give either one value or a value for each replica.");
252 : }
253 : //
254 : double friction;
255 38 : std::vector<double> frictions_vec(0);
256 76 : parseVector("friction",frictions_vec);
257 38 : if(frictions_vec.size()==1) {
258 34 : friction = frictions_vec[0];
259 : }
260 4 : else if(frictions_vec.size()==replicas) {
261 8 : friction = frictions_vec[inter.Get_rank()];
262 : }
263 : else {
264 0 : error("problem with friction keyword, you need to give either one value or a value for each replica.");
265 : }
266 : //
267 : int seed;
268 38 : std::vector<int> seeds_vec(0);
269 76 : parseVector("random_seed",seeds_vec);
270 226 : for(unsigned int i=0; i<seeds_vec.size(); i++) {
271 50 : if(seeds_vec[i]>0) {seeds_vec[i] = -seeds_vec[i];}
272 : }
273 38 : if(replicas==1) {
274 7 : if(seeds_vec.size()>1) {error("problem with random_seed keyword, for a single replica you should only give one value");}
275 7 : seed = seeds_vec[0];
276 : }
277 : else {
278 31 : if(seeds_vec.size()!=1 && seeds_vec.size()!=replicas) {
279 0 : error("problem with random_seed keyword, for multiple replicas you should give either one value or a separate value for each replica");
280 : }
281 31 : if(seeds_vec.size()==1) {
282 27 : seeds_vec.resize(replicas);
283 300 : for(unsigned int i=1; i<seeds_vec.size(); i++) {seeds_vec[i] = seeds_vec[0] + i;}
284 : }
285 62 : seed = seeds_vec[inter.Get_rank()];
286 : }
287 :
288 : //
289 76 : parse("dimension",dim);
290 :
291 38 : std::vector<std::string> plumed_inputfiles;
292 76 : parseVector("plumed_input",plumed_inputfiles);
293 38 : if(plumed_inputfiles.size()!=1 && plumed_inputfiles.size()!=replicas) {
294 0 : error("in plumed_input you should either give one file or separate files for each replica.");
295 : }
296 :
297 38 : std::vector<Vector> initPos(replicas);
298 : std::vector<double> initPosTmp;
299 76 : parseVector("initial_position",initPosTmp);
300 38 : if(initPosTmp.size()==dim) {
301 76 : for(unsigned int i=0; i<replicas; i++) {
302 100 : for(unsigned int k=0; k<dim; k++) {
303 102 : initPos[i][k]=initPosTmp[k];
304 : }
305 : }
306 : }
307 26 : else if(initPosTmp.size()==dim*replicas) {
308 226 : for(unsigned int i=0; i<replicas; i++) {
309 332 : for(unsigned int k=0; k<dim; k++) {
310 348 : initPos[i][k]=initPosTmp[i*dim+k];
311 : }
312 : }
313 : }
314 : else {
315 0 : error("problem with initial_position keyword, you need to give either one value or a value for each replica.");
316 : }
317 :
318 :
319 38 : plumed_bf = new PLMD::PlumedMain;
320 38 : unsigned int nn=1;
321 38 : FILE* file_dummy = fopen("/dev/null","w+");
322 76 : plumed_bf->cmd("setNatoms",&nn);
323 76 : plumed_bf->cmd("setLog",file_dummy);
324 76 : plumed_bf->cmd("init",&nn);
325 38 : std::vector<BasisFunctions*> basisf_pntrs(dim);
326 76 : std::vector<std::string> basisf_keywords(dim);
327 38 : std::vector<Value*> args(dim);
328 76 : std::vector<bool> periodic(dim);
329 38 : std::vector<double> interval_min(dim);
330 38 : std::vector<double> interval_max(dim);
331 38 : std::vector<double> interval_range(dim);
332 126 : for(unsigned int i=0; i<dim; i++) {
333 : std::string bf_keyword;
334 44 : std::string is; Tools::convert(i+1,is);
335 88 : parse("basis_functions_"+is,bf_keyword);
336 44 : if(bf_keyword.size()==0) {
337 0 : error("basis_functions_"+is+" is needed");
338 : }
339 46 : if(bf_keyword.at(0)=='{' && bf_keyword.at(bf_keyword.size()-1)=='}') {
340 4 : bf_keyword = bf_keyword.substr(1,bf_keyword.size()-2);
341 : }
342 44 : basisf_keywords[i] = bf_keyword;
343 132 : plumed_bf->readInputLine(bf_keyword+" LABEL=dim"+is);
344 132 : basisf_pntrs[i] = plumed_bf->getActionSet().selectWithLabel<BasisFunctions*>("dim"+is);
345 132 : args[i] = new Value(NULL,"dim"+is,false);
346 44 : args[i]->setNotPeriodic();
347 44 : periodic[i] = basisf_pntrs[i]->arePeriodic();
348 88 : interval_min[i] = basisf_pntrs[i]->intervalMin();
349 88 : interval_max[i] = basisf_pntrs[i]->intervalMax();
350 88 : interval_range[i] = basisf_pntrs[i]->intervalMax()-basisf_pntrs[i]->intervalMin();
351 : }
352 76 : Communicator comm_dummy;
353 76 : CoeffsVector* coeffs_pntr = new CoeffsVector("pot.coeffs",args,basisf_pntrs,comm_dummy,false);
354 76 : potential_expansion_pntr = new LinearBasisSetExpansion("potential",1.0/temp,comm_dummy,args,basisf_pntrs,coeffs_pntr);
355 :
356 38 : std::string template_coeffs_fname="";
357 76 : parse("template_coeffs_file",template_coeffs_fname);
358 38 : if(template_coeffs_fname.size()>0) {
359 0 : OFile ofile_coeffstmpl;
360 0 : ofile_coeffstmpl.link(pc);
361 0 : ofile_coeffstmpl.open(template_coeffs_fname);
362 0 : coeffs_pntr->writeToFile(ofile_coeffstmpl,true);
363 0 : ofile_coeffstmpl.close();
364 0 : error("Only generating a template coefficient file - Should stop now.");
365 : }
366 :
367 76 : std::vector<std::string> input_coeffs_fnames(0);
368 76 : parseVector("input_coeffs",input_coeffs_fnames);
369 : std::string input_coeffs_fname;
370 : bool diff_input_coeffs = false;
371 38 : if(input_coeffs_fnames.size()==1) {
372 : input_coeffs_fname = input_coeffs_fnames[0];
373 : }
374 9 : else if(replicas > 1 && input_coeffs_fnames.size()==replicas) {
375 : diff_input_coeffs = true;
376 9 : input_coeffs_fname = input_coeffs_fnames[inter.Get_rank()];
377 : }
378 : else {
379 0 : error("problem with coeffs_file keyword, you need to give either one value or a value for each replica.");
380 : }
381 38 : coeffs_pntr->readFromFile(input_coeffs_fname,true,true);
382 38 : std::vector<double> coeffs_prefactors(0);
383 76 : parseVector("coeffs_prefactor",coeffs_prefactors);
384 38 : if(coeffs_prefactors.size()>0) {
385 : double coeffs_prefactor = 1.0;
386 7 : if(coeffs_prefactors.size()==1) {
387 3 : coeffs_prefactor = coeffs_prefactors[0];
388 : }
389 4 : else if(replicas > 1 && coeffs_prefactors.size()==replicas) {
390 : diff_input_coeffs = true;
391 8 : coeffs_prefactor = coeffs_prefactors[inter.Get_rank()];
392 : }
393 : else {
394 0 : error("problem with coeffs_prefactor keyword, you need to give either one value or a value for each replica.");
395 : }
396 7 : coeffs_pntr->scaleAllValues(coeffs_prefactor);
397 : }
398 : unsigned int pot_grid_bins;
399 76 : parse("output_potential_grid",pot_grid_bins);
400 38 : potential_expansion_pntr->setGridBins(pot_grid_bins);
401 38 : potential_expansion_pntr->setupBiasGrid(false);
402 38 : potential_expansion_pntr->updateBiasGrid();
403 38 : potential_expansion_pntr->setBiasMinimumToZero();
404 38 : potential_expansion_pntr->updateBiasGrid();
405 :
406 76 : OFile ofile_potential;
407 38 : ofile_potential.link(pc);
408 : std::string output_potential_fname;
409 76 : parse("output_potential",output_potential_fname);
410 38 : if(diff_input_coeffs) {
411 13 : ofile_potential.link(intra);
412 : std::string suffix;
413 13 : Tools::convert(inter.Get_rank(),suffix);
414 39 : output_potential_fname = FileBase::appendSuffix(output_potential_fname,"."+suffix);
415 : }
416 38 : ofile_potential.open(output_potential_fname);
417 38 : potential_expansion_pntr->writeBiasGridToFile(ofile_potential);
418 38 : ofile_potential.close();
419 :
420 114 : Grid histo_grid(*potential_expansion_pntr->getPntrToBiasGrid());
421 114 : std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(&histo_grid);
422 : double norm=0.0;
423 292914 : for(Grid::index_t i=0; i<histo_grid.getSize(); i++) {
424 146438 : double value = integration_weights[i]*exp(-histo_grid.getValue(i)/temp);
425 146438 : norm += value;
426 146438 : histo_grid.setValue(i,value);
427 : }
428 38 : histo_grid.scaleAllValuesAndDerivatives(1.0/norm);
429 76 : OFile ofile_histogram;
430 38 : ofile_histogram.link(pc);
431 : std::string output_histogram_fname;
432 76 : parse("output_histogram",output_histogram_fname);
433 63 : if(diff_input_coeffs || temps_vec.size()>1) {
434 17 : ofile_histogram.link(intra);
435 : std::string suffix;
436 17 : Tools::convert(inter.Get_rank(),suffix);
437 51 : output_histogram_fname = FileBase::appendSuffix(output_histogram_fname,"."+suffix);
438 : }
439 38 : ofile_histogram.open(output_histogram_fname);
440 38 : histo_grid.writeToFile(ofile_histogram);
441 38 : ofile_histogram.close();
442 :
443 : std::string output_coeffs_fname;
444 76 : parse("output_coeffs",output_coeffs_fname);
445 : std::string output_coeffs_fmt;
446 76 : parse("output_coeffs_fmt",output_coeffs_fmt);
447 : coeffs_pntr->setOutputFmt(output_coeffs_fmt);
448 76 : OFile ofile_coeffsout;
449 38 : ofile_coeffsout.link(pc);
450 38 : if(diff_input_coeffs) {
451 13 : ofile_coeffsout.link(intra);
452 : std::string suffix;
453 13 : Tools::convert(inter.Get_rank(),suffix);
454 39 : output_coeffs_fname = FileBase::appendSuffix(output_coeffs_fname,"."+suffix);
455 : }
456 38 : ofile_coeffsout.open(output_coeffs_fname);
457 38 : coeffs_pntr->writeToFile(ofile_coeffsout,true);
458 38 : ofile_coeffsout.close();
459 :
460 38 : if(pc.Get_rank() == 0) {
461 14 : fprintf(out,"Replicas %u\n",replicas);
462 : fprintf(out,"Cores per replica %u\n",coresPerReplica);
463 14 : fprintf(out,"Number of steps %u\n",nsteps);
464 14 : fprintf(out,"Timestep %f\n",tstep);
465 14 : fprintf(out,"Temperature %f",temps_vec[0]);
466 40 : for(unsigned int i=1; i<temps_vec.size(); i++) {fprintf(out,",%f",temps_vec[i]);}
467 : fprintf(out,"\n");
468 14 : fprintf(out,"Friction %f",frictions_vec[0]);
469 40 : for(unsigned int i=1; i<frictions_vec.size(); i++) {fprintf(out,",%f",frictions_vec[i]);}
470 : fprintf(out,"\n");
471 14 : fprintf(out,"Random seed %d",seeds_vec[0]);
472 120 : for(unsigned int i=1; i<seeds_vec.size(); i++) {fprintf(out,",%d",seeds_vec[i]);}
473 : fprintf(out,"\n");
474 14 : fprintf(out,"Dimensions %u\n",dim);
475 31 : for(unsigned int i=0; i<dim; i++) {
476 34 : fprintf(out,"Basis Function %u %s\n",i+1,basisf_keywords[i].c_str());
477 : }
478 : fprintf(out,"PLUMED input %s",plumed_inputfiles[0].c_str());
479 31 : for(unsigned int i=1; i<plumed_inputfiles.size(); i++) {fprintf(out,",%s",plumed_inputfiles[i].c_str());}
480 : fprintf(out,"\n");
481 : fprintf(out,"kBoltzmann taken as 1, use NATURAL_UNITS in the plumed input\n");
482 14 : if(diff_input_coeffs) {fprintf(out,"using different coefficients for each replica\n");}
483 : }
484 :
485 :
486 38 : plumed=new PLMD::PlumedMain;
487 :
488 :
489 :
490 38 : if(plumed) {
491 38 : int s=sizeof(double);
492 76 : plumed->cmd("setRealPrecision",&s);
493 38 : if(replicas>1) {
494 31 : if (Communicator::initialized()) {
495 62 : plumed->cmd("GREX setMPIIntracomm",&intra.Get_comm());
496 31 : if (intra.Get_rank()==0) {
497 62 : plumed->cmd("GREX setMPIIntercomm",&inter.Get_comm());
498 : }
499 62 : plumed->cmd("GREX init");
500 62 : plumed->cmd("setMPIComm",&intra.Get_comm());
501 : } else {
502 0 : error("More than 1 replica but no MPI");
503 : }
504 : } else {
505 9 : if(Communicator::initialized()) plumed->cmd("setMPIComm",&pc.Get_comm());
506 : }
507 : }
508 :
509 38 : std::string plumed_logfile = "plumed.log";
510 38 : std::string stats_filename = "stats.out";
511 : std::string plumed_input = plumed_inputfiles[0];
512 38 : if(inter.Get_size()>1) {
513 : std::string suffix;
514 31 : Tools::convert(inter.Get_rank(),suffix);
515 93 : plumed_logfile = FileBase::appendSuffix(plumed_logfile,"."+suffix);
516 93 : stats_filename = FileBase::appendSuffix(stats_filename,"."+suffix);
517 31 : if(plumed_inputfiles.size()>1) {
518 2 : plumed_input = plumed_inputfiles[inter.Get_rank()];
519 : }
520 : }
521 :
522 38 : if(plumed) {
523 76 : plumed->cmd("setNoVirial");
524 38 : int natoms=1;
525 76 : plumed->cmd("setNatoms",&natoms);
526 76 : plumed->cmd("setMDEngine","mdrunner_linearexpansion");
527 76 : plumed->cmd("setTimestep",&tstep);
528 76 : plumed->cmd("setPlumedDat",plumed_input.c_str());
529 76 : plumed->cmd("setLogFile",plumed_logfile.c_str());
530 76 : plumed->cmd("setKbT",&temp);
531 38 : double energyunits=1.0;
532 76 : plumed->cmd("setMDEnergyUnits",&energyunits);
533 76 : plumed->cmd("init");
534 : }
535 :
536 : // Setup random number generator
537 38 : random.setSeed(seed);
538 :
539 38 : double potential, therm_eng=0; std::vector<double> masses(1,1);
540 38 : std::vector<Vector> positions(1), velocities(1), forces(1);
541 126 : for(unsigned int k=0; k<dim; k++) {
542 132 : positions[0][k] = initPos[inter.Get_rank()][k];
543 88 : if(periodic[k]) {
544 12 : positions[0][k] = positions[0][k] - floor((positions[0][k]-interval_min[k])/interval_range[k])*interval_range[k];
545 : }
546 : else {
547 80 : if(positions[0][k]>interval_max[k]) {positions[0][k]=interval_max[k];}
548 80 : if(positions[0][k]<interval_min[k]) {positions[0][k]=interval_min[k];}
549 : }
550 : }
551 :
552 :
553 126 : for(unsigned k=0; k<dim; ++k) {
554 88 : velocities[0][k]=random.Gaussian() * sqrt( temp );
555 : }
556 :
557 38 : potential=calc_energy(positions,forces); double ttt=calc_temp(velocities);
558 :
559 38 : FILE* fp=fopen(stats_filename.c_str(),"w+");
560 38 : double conserved = potential+1.5*ttt+therm_eng;
561 : //fprintf(fp,"%d %f %f %f %f %f %f %f %f \n", 0, 0., positions[0][0], positions[0][1], positions[0][2], conserved, ttt, potential, therm_eng );
562 38 : if( intra.Get_rank()==0 ) {
563 74 : fprintf(fp,"%d %f %f %f %f %f %f %f %f \n", 0, 0., positions[0][0], positions[0][1], positions[0][2], conserved, ttt, potential, therm_eng );
564 : }
565 :
566 38 : if(plumed) {
567 38 : int step_tmp = 0;
568 76 : plumed->cmd("setStep",&step_tmp);
569 114 : plumed->cmd("setMasses",&masses[0]);
570 114 : plumed->cmd("setForces",&forces[0]);
571 76 : plumed->cmd("setEnergy",&potential);
572 114 : plumed->cmd("setPositions",&positions[0]);
573 76 : plumed->cmd("calc");
574 : }
575 :
576 7638 : for(unsigned int istep=0; istep<nsteps; ++istep) {
577 : //if( istep%20==0 && pc.Get_rank()==0 ) printf("Doing step %d\n",istep);
578 :
579 : // Langevin thermostat
580 3800 : double lscale=exp(-0.5*tstep*friction); //exp(-0.5*tstep/friction);
581 3800 : double lrand=sqrt((1.-lscale*lscale)*temp);
582 12600 : for(unsigned k=0; k<dim; ++k) {
583 4400 : therm_eng=therm_eng+0.5*velocities[0][k]*velocities[0][k];
584 8800 : velocities[0][k]=lscale*velocities[0][k]+lrand*random.Gaussian();
585 4400 : therm_eng=therm_eng-0.5*velocities[0][k]*velocities[0][k];
586 : }
587 :
588 : // First step of velocity verlet
589 12600 : for(unsigned k=0; k<dim; ++k) {
590 8800 : velocities[0][k] = velocities[0][k] + 0.5*tstep*forces[0][k];
591 8800 : positions[0][k] = positions[0][k] + tstep*velocities[0][k];
592 :
593 8800 : if(periodic[k]) {
594 1200 : positions[0][k] = positions[0][k] - floor((positions[0][k]-interval_min[k])/interval_range[k])*interval_range[k];
595 : }
596 : else {
597 8000 : if(positions[0][k]>interval_max[k]) {
598 7 : positions[0][k]=interval_max[k];
599 14 : velocities[0][k]=-std::abs(velocities[0][k]);
600 : }
601 8000 : if(positions[0][k]<interval_min[k]) {
602 2 : positions[0][k]=interval_min[k];
603 4 : velocities[0][k]=-std::abs(velocities[0][k]);
604 : }
605 : }
606 : }
607 :
608 3800 : potential=calc_energy(positions,forces);
609 :
610 3800 : if(plumed) {
611 3800 : int istepplusone=istep+1;
612 3800 : plumedWantsToStop=0;
613 7600 : plumed->cmd("setStep",&istepplusone);
614 11400 : plumed->cmd("setMasses",&masses[0]);
615 11400 : plumed->cmd("setForces",&forces[0]);
616 7600 : plumed->cmd("setEnergy",&potential);
617 11400 : plumed->cmd("setPositions",&positions[0]);
618 7600 : plumed->cmd("setStopFlag",&plumedWantsToStop);
619 7600 : plumed->cmd("calc");
620 : //if(istep%2000==0) plumed->cmd("writeCheckPointFile");
621 3800 : if(plumedWantsToStop) nsteps=istep;
622 : }
623 :
624 : // Second step of velocity verlet
625 12600 : for(unsigned k=0; k<dim; ++k) {
626 8800 : velocities[0][k] = velocities[0][k] + 0.5*tstep*forces[0][k];
627 : }
628 :
629 : // Langevin thermostat
630 3800 : lscale=exp(-0.5*tstep*friction); //exp(-0.5*tstep/friction);
631 3800 : lrand=sqrt((1.-lscale*lscale)*temp);
632 12600 : for(unsigned k=0; k<dim; ++k) {
633 4400 : therm_eng=therm_eng+0.5*velocities[0][k]*velocities[0][k];
634 8800 : velocities[0][k]=lscale*velocities[0][k]+lrand*random.Gaussian();
635 4400 : therm_eng=therm_eng-0.5*velocities[0][k]*velocities[0][k];
636 : }
637 :
638 : // Print everything
639 3800 : ttt = calc_temp( velocities );
640 3800 : conserved = potential+1.5*ttt+therm_eng;
641 3800 : if( (intra.Get_rank()==0) && ((istep % stepWrite)==0) ) {
642 74 : fprintf(fp,"%u %f %f %f %f %f %f %f %f \n", istep, istep*tstep, positions[0][0], positions[0][1], positions[0][2], conserved, ttt, potential, therm_eng );
643 : }
644 : }
645 :
646 38 : if(plumed) {delete plumed;}
647 38 : if(plumed_bf) {delete plumed_bf;}
648 38 : if(potential_expansion_pntr) {delete potential_expansion_pntr;}
649 38 : if(coeffs_pntr) {delete coeffs_pntr;}
650 208 : for(unsigned int i=0; i<args.size(); i++) {delete args[i];}
651 : args.clear();
652 : //printf("Rank: %d, Size: %d \n", pc.Get_rank(), pc.Get_size() );
653 : //printf("Rank: %d, Size: %d, MultiSimCommRank: %d, MultiSimCommSize: %d \n", pc.Get_rank(), pc.Get_size(), multi_sim_comm.Get_rank(), multi_sim_comm.Get_size() );
654 38 : fclose(fp);
655 38 : fclose(file_dummy);
656 :
657 38 : return 0;
658 : }
659 :
660 : }
661 5517 : }
|
