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