8 #include "McEnergyAnalyzer.h" 9 #include <mcMd/mcSimulation/McSimulation.h> 10 #include <mcMd/mcSimulation/McSystem.h> 13 #include <mcMd/potentials/pair/McPairPotential.h> 16 #include <mcMd/potentials/bond/BondPotential.h> 19 #include <mcMd/potentials/angle/AnglePotential.h> 22 #include <mcMd/potentials/dihedral/DihedralPotential.h> 25 #include <mcMd/potentials/external/ExternalPotential.h> 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> 57 dihedralAveragePtr_(0),
60 externalAveragePtr_(0),
74 readOptional<int>(in,
"nSamplePerBlock", nSamplePerBlock_);
85 if (sys.hasBondPotential()) {
91 if (sys.hasAnglePotential()) {
97 if (sys.hasDihedralPotential()) {
98 dihedralAveragePtr_ =
new Average;
103 if (sys.hasExternalPotential()) {
104 externalAveragePtr_ =
new Average;
109 isInitialized_ =
true;
120 nSamplePerBlock_ = 0;
122 loadParameter<int>(ar,
"nSamplePerBlock", nSamplePerBlock_,
isRequired);
127 totalAveragePtr_ =
new Average;
128 ar >> *totalAveragePtr_;
131 ar >> *pairAveragePtr_;
134 if (sys.hasBondPotential()) {
136 ar >> *bondAveragePtr_;
140 if (sys.hasAnglePotential()) {
141 angleAveragePtr_ =
new Average;
142 ar >> *angleAveragePtr_;
146 if (sys.hasDihedralPotential()) {
147 dihedralAveragePtr_ =
new Average;
148 ar >> *dihedralAveragePtr_;
152 if (sys.hasExternalPotential()) {
153 externalAveragePtr_ =
new Average;
154 ar >> *externalAveragePtr_;
159 if (nSamplePerBlock_ != 0) {
163 isInitialized_ =
true;
174 bool isActive = bool(nSamplePerBlock_);
180 ar << *totalAveragePtr_;
182 ar << *pairAveragePtr_;
186 ar << *bondAveragePtr_;
191 assert(angleAveragePtr_);
192 ar << *angleAveragePtr_;
197 assert(dihedralAveragePtr_);
198 ar << *dihedralAveragePtr_;
203 assert(externalAveragePtr_);
204 ar << *externalAveragePtr_;
216 totalAveragePtr_->
clear();
218 pairAveragePtr_->
clear();
222 bondAveragePtr_->
clear();
227 assert(angleAveragePtr_);
228 angleAveragePtr_->
clear();
233 assert(dihedralAveragePtr_);
234 dihedralAveragePtr_->
clear();
239 assert(externalAveragePtr_);
240 externalAveragePtr_->
clear();
251 if (outputFile_.is_open()) {
255 if (nSamplePerBlock_ > 0) {
260 std::string filename;
264 outputFile_ <<
" iStep";
266 outputFile_ <<
" Pair";
270 outputFile_ <<
" Bond";
275 outputFile_ <<
" Angle";
280 outputFile_ <<
" Dihedral";
285 outputFile_ <<
" Coulomb";
290 outputFile_ <<
" External";
293 outputFile_ <<
" Total";
316 double potential = 0.0;
320 pairAveragePtr_->
sample(pair);
327 bondAveragePtr_->
sample(bond);
335 assert(angleAveragePtr_);
336 angleAveragePtr_->
sample(angle);
343 potential += dihedral;
344 assert(dihedralAveragePtr_);
345 dihedralAveragePtr_->
sample(dihedral);
352 potential += external;
353 assert(externalAveragePtr_);
354 externalAveragePtr_->
sample(external);
358 double total = potential;
359 totalAveragePtr_->
sample(total);
365 int beginStep = iStep - (nSamplePerBlock_ - 1)*
interval();
366 outputFile_ <<
Int(beginStep, 12);
368 outputFile_ <<
Dbl(pairAveragePtr_->blockAverage());
372 outputFile_ <<
Dbl(bondAveragePtr_->blockAverage());
377 assert(angleAveragePtr_);
378 outputFile_ <<
Dbl(angleAveragePtr_->blockAverage());
383 assert(dihedralAveragePtr_);
384 outputFile_ <<
Dbl(dihedralAveragePtr_->blockAverage());
389 assert(externalAveragePtr_);
390 outputFile_ <<
Dbl(externalAveragePtr_->blockAverage());
393 outputFile_ <<
Dbl(totalAveragePtr_->blockAverage());
409 if (outputFile_.is_open()) {
422 ave = pairAveragePtr_->
average();
424 outputFile_ <<
"Pair " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
428 ave = bondAveragePtr_->
average();
430 outputFile_ <<
"Bond " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
435 assert(angleAveragePtr_);
436 ave = angleAveragePtr_->
average();
438 outputFile_ <<
"Angle " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
443 assert(dihedralAveragePtr_);
444 ave = dihedralAveragePtr_->
average();
446 outputFile_ <<
"Dihedral " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
451 assert(externalAveragePtr_);
452 ave = externalAveragePtr_->
average();
454 outputFile_ <<
"External " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
457 ave = totalAveragePtr_->
average();
459 outputFile_ <<
"Total " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
465 "---------------------------------------------------------------------------------\n";
466 outputFile_ <<
"Pair:\n\n";
468 pairAveragePtr_->
output(outputFile_);
473 "---------------------------------------------------------------------------------\n";
474 outputFile_ <<
"Bond:\n\n";
475 bondAveragePtr_->
output(outputFile_);
481 "---------------------------------------------------------------------------------\n";
482 outputFile_ <<
"Angle:\n\n";
483 assert(angleAveragePtr_);
484 angleAveragePtr_->
output(outputFile_);
490 "---------------------------------------------------------------------------------\n";
491 outputFile_ <<
"Dihedral:\n\n";
492 assert(dihedralAveragePtr_);
493 dihedralAveragePtr_->
output(outputFile_);
499 "---------------------------------------------------------------------------------\n";
500 outputFile_ <<
"External:\n\n";
501 assert(externalAveragePtr_);
502 externalAveragePtr_->
output(outputFile_);
506 "---------------------------------------------------------------------------------\n";
507 outputFile_ <<
"Total:\n\n";
508 totalAveragePtr_->
output(outputFile_);
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.
A System for use in a Markov chain Monte Carlo simulation.
bool isBlockComplete() const
Is the current block average complete?
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.
double average() const
Return the average of all sampled values.
BondPotential & bondPotential() const
Return the BondPotential by reference.
bool hasDihedralPotential() const
Does a dihedral potential exist?.
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.
virtual void loadParameters(Serializable::IArchive &ar)
Load parameters from archive.
bool hasAnglePotential() const
Does angle potential exist?.
McSystem & system()
Return reference to parent system.
Wrapper for a double precision number, for formatted ostream output.
double blockingError() const
Return estimated error on average from blocking analysis.
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.
bool hasBondPotential() const
Does a bond potential exist?.
virtual void writeParam(std::ostream &out)
Write all parameters to an output stream.
ExternalPotential & externalPotential() const
Return ExternalPotential by reference.
void readInterval(std::istream &in)
Read interval from file, with error checking.
Utility classes for scientific computation.
Wrapper for an int, for formatted ostream output.
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.
McPairPotential & pairPotential() const
Return the McPairPotential by reference.
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.
Saving archive for binary istream.
void sample(double value)
Add a sampled value to the ensemble.
bool hasCoulombPotential() const
Does a Coulomb potential exist?.
Single-processor Monte Carlo (MC) and molecular dynamics (MD).
bool hasExternalPotential() const
Does an external potential exist?.
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.
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.
FileMaster & fileMaster()
Get the FileMaster object.
DihedralPotential & dihedralPotential() const
Return the DihedralPotential by reference.