Simpatico  v1.10
McEnergyAnalyzer.cpp
1 /*
2 * Simpatico - Simulation Package for Polymeric and Molecular Liquids
3 *
4 * Copyright 2010 - 2017, The Regents of the University of Minnesota
5 * Distributed under the terms of the GNU General Public License.
6 */
7 
8 #include "McEnergyAnalyzer.h"
9 #include <mcMd/mcSimulation/McSimulation.h>
10 #include <mcMd/mcSimulation/McSystem.h>
11 
12 #ifndef SIMP_NOPAIR
13 #include <mcMd/potentials/pair/McPairPotential.h>
14 #endif
15 #ifdef SIMP_BOND
16 #include <mcMd/potentials/bond/BondPotential.h>
17 #endif
18 #ifdef SIMP_ANGLE
19 #include <mcMd/potentials/angle/AnglePotential.h>
20 #endif
21 #ifdef SIMP_DIHEDRAL
22 #include <mcMd/potentials/dihedral/DihedralPotential.h>
23 #endif
24 #ifdef SIMP_EXTERNAL
25 #include <mcMd/potentials/external/ExternalPotential.h>
26 #endif
27 
28 #include <util/accumulators/Average.h>
29 #include <util/format/Int.h>
30 #include <util/format/Dbl.h>
31 #include <util/mpi/MpiLoader.h>
32 #include <util/misc/ioUtil.h>
33 
34 #include <sstream>
35 
36 namespace McMd
37 {
38 
39  using namespace Util;
40 
41  /*
42  * Constructor.
43  */
45  : SystemAnalyzer<McSystem>(system),
46  totalAveragePtr_(0),
47  #ifndef SIMP_NOPAIR
48  pairAveragePtr_(0),
49  #endif
50  #ifdef SIMP_BOND
51  bondAveragePtr_(0),
52  #endif
53  #ifdef SIMP_ANGLE
54  angleAveragePtr_(0),
55  #endif
56  #ifdef SIMP_DIHEDRAL
57  dihedralAveragePtr_(0),
58  #endif
59  #ifdef SIMP_EXTERNAL
60  externalAveragePtr_(0),
61  #endif
62  nSamplePerBlock_(0),
63  isInitialized_(false)
64  { setClassName("McEnergyAnalyzer"); }
65 
66  /*
67  * Read interval and outputFileName.
68  */
69  void McEnergyAnalyzer::readParameters(std::istream& in)
70  {
71  readInterval(in);
73  nSamplePerBlock_ = 0;
74  readOptional<int>(in, "nSamplePerBlock", nSamplePerBlock_);
75 
76  McSystem& sys = system();
77 
78  totalAveragePtr_ = new Average;
79  totalAveragePtr_->setNSamplePerBlock(nSamplePerBlock_);
80  #ifndef SIMP_NOPAIR
81  pairAveragePtr_ = new Average;
82  pairAveragePtr_->setNSamplePerBlock(nSamplePerBlock_);
83  #endif
84  #ifdef SIMP_BOND
85  if (sys.hasBondPotential()) {
86  bondAveragePtr_ = new Average;
87  bondAveragePtr_->setNSamplePerBlock(nSamplePerBlock_);
88  }
89  #endif
90  #ifdef SIMP_ANGLE
91  if (sys.hasAnglePotential()) {
92  angleAveragePtr_ = new Average;
93  angleAveragePtr_->setNSamplePerBlock(nSamplePerBlock_);
94  }
95  #endif
96  #ifdef SIMP_DIHEDRAL
97  if (sys.hasDihedralPotential()) {
98  dihedralAveragePtr_ = new Average;
99  dihedralAveragePtr_->setNSamplePerBlock(nSamplePerBlock_);
100  }
101  #endif
102  #ifdef SIMP_EXTERNAL
103  if (sys.hasExternalPotential()) {
104  externalAveragePtr_ = new Average;
105  externalAveragePtr_->setNSamplePerBlock(nSamplePerBlock_);
106  }
107  #endif
108 
109  isInitialized_ = true;
110  }
111 
112  /*
113  * Load internal state from an archive.
114  */
116  {
117  // Load interval and outputFileName
119 
120  nSamplePerBlock_ = 0; // default value
121  bool isRequired = false;
122  loadParameter<int>(ar, "nSamplePerBlock", nSamplePerBlock_, isRequired);
123 
124  McSystem& sys = system();
125 
126  // Load Average accumulators
127  totalAveragePtr_ = new Average;
128  ar >> *totalAveragePtr_;
129  #ifndef SIMP_NOPAIR
130  pairAveragePtr_ = new Average;
131  ar >> *pairAveragePtr_;
132  #endif
133  #ifdef SIMP_BOND
134  if (sys.hasBondPotential()) {
135  bondAveragePtr_ = new Average;
136  ar >> *bondAveragePtr_;
137  }
138  #endif
139  #ifdef SIMP_ANGLE
140  if (sys.hasAnglePotential()) {
141  angleAveragePtr_ = new Average;
142  ar >> *angleAveragePtr_;
143  }
144  #endif
145  #ifdef SIMP_DIHEDRAL
146  if (sys.hasDihedralPotential()) {
147  dihedralAveragePtr_ = new Average;
148  ar >> *dihedralAveragePtr_;
149  }
150  #endif
151  #ifdef SIMP_EXTERNAL
152  if (sys.hasExternalPotential()) {
153  externalAveragePtr_ = new Average;
154  ar >> *externalAveragePtr_;
155  }
156  #endif
157 
158  // If nSamplePerBlock != 0, open an output file for block averages.
159  if (nSamplePerBlock_ != 0) {
160  fileMaster().openOutputFile(Analyzer::outputFileName(".dat"), outputFile_);
161  }
162 
163  isInitialized_ = true;
164  }
165 
166  /*
167  * Save internal state to an archive.
168  */
170  {
171  Analyzer::save(ar);
172  //saveInterval(ar);
173  //saveOutputFileName(ar);
174  bool isActive = bool(nSamplePerBlock_);
175  Parameter::saveOptional(ar, nSamplePerBlock_, isActive);
176 
177  McSystem& sys = system();
178 
179  // Save average accumulators
180  ar << *totalAveragePtr_;
181  #ifndef SIMP_NOPAIR
182  ar << *pairAveragePtr_;
183  #endif
184  #ifdef SIMP_BOND
185  if (sys.hasBondPotential()) {
186  ar << *bondAveragePtr_;
187  }
188  #endif
189  #ifdef SIMP_ANGLE
190  if (sys.hasAnglePotential()) {
191  assert(angleAveragePtr_);
192  ar << *angleAveragePtr_;
193  }
194  #endif
195  #ifdef SIMP_DIHEDRAL
196  if (sys.hasDihedralPotential()) {
197  assert(dihedralAveragePtr_);
198  ar << *dihedralAveragePtr_;
199  }
200  #endif
201  #ifdef SIMP_EXTERNAL
202  if (sys.hasExternalPotential()) {
203  assert(externalAveragePtr_);
204  ar << *externalAveragePtr_;
205  }
206  #endif
207  }
208 
209  /*
210  * Reset nSample.
211  */
213  {
214  McSystem& sys = system();
215 
216  totalAveragePtr_->clear();
217  #ifndef SIMP_NOPAIR
218  pairAveragePtr_->clear();
219  #endif
220  #ifdef SIMP_BOND
221  if (sys.hasBondPotential()) {
222  bondAveragePtr_->clear();
223  }
224  #endif
225  #ifdef SIMP_ANGLE
226  if (sys.hasAnglePotential()) {
227  assert(angleAveragePtr_);
228  angleAveragePtr_->clear();
229  }
230  #endif
231  #ifdef SIMP_DIHEDRAL
232  if (sys.hasDihedralPotential()) {
233  assert(dihedralAveragePtr_);
234  dihedralAveragePtr_->clear();
235  }
236  #endif
237  #ifdef SIMP_EXTERNAL
238  if (sys.hasExternalPotential()) {
239  assert(externalAveragePtr_);
240  externalAveragePtr_->clear();
241  }
242  #endif
243  }
244 
245  /*
246  * Open outputfile
247  */
249  {
250 
251  if (outputFile_.is_open()) {
252  outputFile_.close();
253  }
254 
255  if (nSamplePerBlock_ > 0) {
256  McSystem& sys = system();
257  Simulation& sim = system().simulation();
258 
259  // Open *.fmt file with format of *.dat data file
260  std::string filename;
261  filename = outputFileName(".fmt");
262  sim.fileMaster().openOutputFile(filename, outputFile_);
263 
264  outputFile_ << " iStep";
265  #ifndef SIMP_NOPAIR
266  outputFile_ << " Pair";
267  #endif
268  #ifdef SIMP_BOND
269  if (sys.hasBondPotential()) {
270  outputFile_ << " Bond";
271  }
272  #endif
273  #ifdef SIMP_ANGLE
274  if (sys.hasAnglePotential()) {
275  outputFile_ << " Angle";
276  }
277  #endif
278  #ifdef SIMP_DIHEDRAL
279  if (sys.hasDihedralPotential()) {
280  outputFile_ << " Dihedral";
281  }
282  #endif
283  #if 0
284  if (sys.hasCoulombPotential()) {
285  outputFile_ << " Coulomb";
286  }
287  #endif
288  #ifdef SIMP_EXTERNAL
289  if (sys.hasExternalPotential()) {
290  outputFile_ << " External";
291  }
292  #endif
293  outputFile_ << " Total";
294  outputFile_.close();
295 
296  // Open *.dat data file, leave open for writing
297  filename = outputFileName(".dat");
298  sim.fileMaster().openOutputFile(filename, outputFile_);
299  }
300 
301  //std::string filename;
302  //filename = outputFileName(".dat");
303  //sim.fileMaster().openOutputFile(filename, outputFile_);
304  }
305 
306  /*
307  * Output energy to file
308  */
309  void McEnergyAnalyzer::sample(long iStep)
310  {
311  if (isAtInterval(iStep)) {
312  McSystem& sys = system();
313 
314  //outputFile_ << Int(iStep, 10);
315 
316  double potential = 0.0;
317  #ifndef SIMP_NOPAIR
318  double pair = sys.pairPotential().energy();
319  potential += pair;
320  pairAveragePtr_->sample(pair);
321  //outputFile_ << Dbl(pair, 15);
322  #endif
323  #ifdef SIMP_BOND
324  if (sys.hasBondPotential()) {
325  double bond = sys.bondPotential().energy();
326  potential += bond;
327  bondAveragePtr_->sample(bond);
328  //outputFile_ << Dbl(bond, 15);
329  }
330  #endif
331  #ifdef SIMP_ANGLE
332  if (sys.hasAnglePotential()) {
333  double angle = sys.anglePotential().energy();
334  potential += angle;
335  assert(angleAveragePtr_);
336  angleAveragePtr_->sample(angle);
337  //outputFile_ << Dbl(angle, 15);
338  }
339  #endif
340  #ifdef SIMP_DIHEDRAL
341  if (sys.hasDihedralPotential()) {
342  double dihedral = sys.dihedralPotential().energy();
343  potential += dihedral;
344  assert(dihedralAveragePtr_);
345  dihedralAveragePtr_->sample(dihedral);
346  // outputFile_ << Dbl(dihedral, 15);
347  }
348  #endif
349  #ifdef SIMP_EXTERNAL
350  if (sys.hasExternalPotential()) {
351  double external = sys.externalPotential().energy();
352  potential += external;
353  assert(externalAveragePtr_);
354  externalAveragePtr_->sample(external);
355  // outputFile_ << Dbl(external, 15);
356  }
357  #endif
358  double total = potential;
359  totalAveragePtr_->sample(total);
360  // outputFile_ << Dbl(total, 20)
361  // << std::endl;
362 
363  // Output block averages, if needed
364  if (nSamplePerBlock_ > 0 && totalAveragePtr_->isBlockComplete()) {
365  int beginStep = iStep - (nSamplePerBlock_ - 1)*interval();
366  outputFile_ << Int(beginStep, 12);
367  #ifndef SIMP_NOPAIR
368  outputFile_ << Dbl(pairAveragePtr_->blockAverage());
369  #endif
370  #ifdef SIMP_BOND
371  if (sys.hasBondPotential()) {
372  outputFile_ << Dbl(bondAveragePtr_->blockAverage());
373  }
374  #endif
375  #ifdef SIMP_ANGLE
376  if (sys.hasAnglePotential()) {
377  assert(angleAveragePtr_);
378  outputFile_ << Dbl(angleAveragePtr_->blockAverage());
379  }
380  #endif
381  #ifdef SIMP_DIHEDRAL
382  if (sys.hasDihedralPotential()) {
383  assert(dihedralAveragePtr_);
384  outputFile_ << Dbl(dihedralAveragePtr_->blockAverage());
385  }
386  #endif
387  #ifdef SIMP_EXTERNAL
388  if (sys.hasExternalPotential()) {
389  assert(externalAveragePtr_);
390  outputFile_ << Dbl(externalAveragePtr_->blockAverage());
391  }
392  #endif
393  outputFile_ << Dbl(totalAveragePtr_->blockAverage());
394  outputFile_ << "\n";
395  }
396 
397  }
398  }
399 
400  /*
401  * Output results to file after simulation is completed.
402  */
404  {
405  McSystem& sys = system();
406  Simulation& sim = sys.simulation();
407 
408  // Close data (*.dat) file, if any
409  if (outputFile_.is_open()) {
410  outputFile_.close();
411  }
412 
413  // Write parameter (*.prm) file
414  sim.fileMaster().openOutputFile(outputFileName(".prm"), outputFile_);
415  ParamComposite::writeParam(outputFile_);
416  outputFile_.close();
417 
418  // Write average (*.ave) file
419  sim.fileMaster().openOutputFile(outputFileName(".ave"), outputFile_);
420  double ave, err;
421  #ifndef SIMP_NOPAIR
422  ave = pairAveragePtr_->average();
423  err = pairAveragePtr_->blockingError();
424  outputFile_ << "Pair " << Dbl(ave) << " +- " << Dbl(err, 9, 2) << "\n";
425  #endif
426  #ifdef SIMP_BOND
427  if (sys.hasBondPotential()) {
428  ave = bondAveragePtr_->average();
429  err = bondAveragePtr_->blockingError();
430  outputFile_ << "Bond " << Dbl(ave) << " +- " << Dbl(err, 9, 2) << "\n";
431  }
432  #endif
433  #ifdef SIMP_ANGLE
434  if (sys.hasAnglePotential()) {
435  assert(angleAveragePtr_);
436  ave = angleAveragePtr_->average();
437  err = angleAveragePtr_->blockingError();
438  outputFile_ << "Angle " << Dbl(ave) << " +- " << Dbl(err, 9, 2) << "\n";
439  }
440  #endif
441  #ifdef SIMP_DIHEDRAL
442  if (sys.hasDihedralPotential()) {
443  assert(dihedralAveragePtr_);
444  ave = dihedralAveragePtr_->average();
445  err = dihedralAveragePtr_->blockingError();
446  outputFile_ << "Dihedral " << Dbl(ave) << " +- " << Dbl(err, 9, 2) << "\n";
447  }
448  #endif
449  #ifdef SIMP_EXTERNAL
450  if (sys.hasExternalPotential()) {
451  assert(externalAveragePtr_);
452  ave = externalAveragePtr_->average();
453  err = externalAveragePtr_->blockingError();
454  outputFile_ << "External " << Dbl(ave) << " +- " << Dbl(err, 9, 2) << "\n";
455  }
456  #endif
457  ave = totalAveragePtr_->average();
458  err = totalAveragePtr_->blockingError();
459  outputFile_ << "Total " << Dbl(ave) << " +- " << Dbl(err, 9, 2) << "\n";
460  outputFile_.close();
461 
462  // Write error analysis (*.aer) file
463  sim.fileMaster().openOutputFile(outputFileName(".aer"), outputFile_);
464  outputFile_ <<
465  "---------------------------------------------------------------------------------\n";
466  outputFile_ << "Pair:\n\n";
467  #ifndef SIMP_NOPAIR
468  pairAveragePtr_->output(outputFile_);
469  #endif
470  #ifdef SIMP_BOND
471  if (sys.hasBondPotential()) {
472  outputFile_ <<
473  "---------------------------------------------------------------------------------\n";
474  outputFile_ << "Bond:\n\n";
475  bondAveragePtr_->output(outputFile_);
476  }
477  #endif
478  #ifdef SIMP_ANGLE
479  if (sys.hasAnglePotential()) {
480  outputFile_ <<
481  "---------------------------------------------------------------------------------\n";
482  outputFile_ << "Angle:\n\n";
483  assert(angleAveragePtr_);
484  angleAveragePtr_->output(outputFile_);
485  }
486  #endif
487  #ifdef SIMP_DIHEDRAL
488  if (sys.hasDihedralPotential()) {
489  outputFile_ <<
490  "---------------------------------------------------------------------------------\n";
491  outputFile_ << "Dihedral:\n\n";
492  assert(dihedralAveragePtr_);
493  dihedralAveragePtr_->output(outputFile_);
494  }
495  #endif
496  #ifdef SIMP_EXTERNAL
497  if (sys.hasExternalPotential()) {
498  outputFile_ <<
499  "---------------------------------------------------------------------------------\n";
500  outputFile_ << "External:\n\n";
501  assert(externalAveragePtr_);
502  externalAveragePtr_->output(outputFile_);
503  }
504  #endif
505  outputFile_ <<
506  "---------------------------------------------------------------------------------\n";
507  outputFile_ << "Total:\n\n";
508  totalAveragePtr_->output(outputFile_);
509  outputFile_.close();
510 
511  }
512 
513 
514 }
virtual void save(Serializable::OArchive &ar)
Load parameters to archive.
virtual void clear()
Clear nSample counter and all accumulators.
void clear()
Clear all accumulators, set to empty initial state.
Definition: Average.cpp:42
A System for use in a Markov chain Monte Carlo simulation.
Definition: McSystem.h:52
bool isBlockComplete() const
Is the current block average complete?
Definition: Average.h:232
virtual double energy(double rsq, int iAtomType, int jAtomType) const =0
Return pair energy for a single pair.
bool isRequired() const
Is this ParamComposite required in the input file?
Calculates the average and variance of a sampled property.
Definition: Average.h:43
double average() const
Return the average of all sampled values.
BondPotential & bondPotential() const
Return the BondPotential by reference.
Definition: McSystem.h:405
bool hasDihedralPotential() const
Does a dihedral potential exist?.
Definition: McSystem.h:433
McEnergyAnalyzer(McSystem &system)
Constructor.
virtual void loadParameters(Serializable::IArchive &ar)
Load internal state from an input archive.
virtual void readParameters(std::istream &in)
Read interval and output file name.
void openOutputFile(const std::string &filename, std::ofstream &out, std::ios_base::openmode mode=std::ios_base::out) const
Open an output file.
Definition: FileMaster.cpp:290
virtual void loadParameters(Serializable::IArchive &ar)
Load parameters from archive.
bool hasAnglePotential() const
Does angle potential exist?.
Definition: McSystem.h:416
McSystem & system()
Return reference to parent system.
Wrapper for a double precision number, for formatted ostream output.
Definition: Dbl.h:39
double blockingError() const
Return estimated error on average from blocking analysis.
Definition: Average.cpp:133
void readOutputFileName(std::istream &in)
Read outputFileName from file.
The main object in a simulation, which coordinates others.
bool isActive() const
Is this parameter active?
Saving / output archive for binary ostream.
virtual double energy(const Vector &R1, const Vector &R2, const Vector &R3, int type) const =0
Returns potential energy for one dihedral.
virtual void setup()
Setup before main loop.
virtual double energy(const Vector &position, int i) const =0
Returns external potential energy of a single particle.
Simulation & simulation() const
Get the parent Simulation by reference.
Definition: System.h:1055
bool hasBondPotential() const
Does a bond potential exist?.
Definition: McSystem.h:399
virtual void writeParam(std::ostream &out)
Write all parameters to an output stream.
ExternalPotential & externalPotential() const
Return ExternalPotential by reference.
Definition: McSystem.h:473
void readInterval(std::istream &in)
Read interval from file, with error checking.
Utility classes for scientific computation.
Definition: accumulators.mod:1
Wrapper for an int, for formatted ostream output.
Definition: Int.h:36
virtual void sample(long iStep)
Write energy to file.
virtual void output()
Write file averages and error analysis to file.
Template for Analyzer associated with one System.
void output(std::ostream &out) const
Output final statistical properties to file.
Definition: Average.cpp:178
McPairPotential & pairPotential() const
Return the McPairPotential by reference.
Definition: McSystem.h:388
virtual double energy(double cosTheta, int type) const =0
Returns potential energy for one angle.
virtual double energy(double rSq, int type) const =0
Returns potential energy for one bond.
void setNSamplePerBlock(int nSamplePerBlock)
Set nSamplePerBlock.
Definition: Average.cpp:63
Saving archive for binary istream.
void sample(double value)
Add a sampled value to the ensemble.
Definition: Average.cpp:94
bool hasCoulombPotential() const
Does a Coulomb potential exist?.
Definition: McSystem.h:450
Single-processor Monte Carlo (MC) and molecular dynamics (MD).
bool hasExternalPotential() const
Does an external potential exist?.
Definition: McSystem.h:467
void setClassName(const char *className)
Set class name string.
FileMaster & fileMaster()
Get the FileMaster by reference.
bool isAtInterval(long counter) const
Return true iff counter is a multiple of the interval.
AnglePotential & anglePotential() const
Return AnglePotential by reference.
Definition: McSystem.h:422
const std::string & outputFileName() const
Return outputFileName string.
int interval() const
Get interval value.
virtual void save(Serializable::OArchive &ar)
Save internal state to an output archive.
static void saveOptional(Serializable::OArchive &ar, Type &value, bool isActive)
Save an optional parameter value to an output archive.
Definition: Parameter.h:224
FileMaster & fileMaster()
Get the FileMaster object.
DihedralPotential & dihedralPotential() const
Return the DihedralPotential by reference.
Definition: McSystem.h:439