8 #include "MdEnergyAnalyzer.h" 9 #include <mcMd/mdSimulation/MdSimulation.h> 10 #include <mcMd/mdSimulation/MdSystem.h> 13 #include <mcMd/potentials/pair/MdPairPotential.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/coulomb/MdCoulombPotential.h> 28 #include <mcMd/potentials/external/ExternalPotential.h> 31 #include <util/accumulators/Average.h> 32 #include <util/format/Int.h> 33 #include <util/format/Dbl.h> 34 #include <util/mpi/MpiLoader.h> 35 #include <util/misc/ioUtil.h> 50 kineticAveragePtr_(0),
51 potentialAveragePtr_(0),
62 dihedralAveragePtr_(0),
65 coulombRSpaceAveragePtr_(0),
66 coulombKSpaceAveragePtr_(0),
67 coulombAveragePtr_(0),
68 coulombComponents_(false),
71 externalAveragePtr_(0),
87 readOptional<int>(in,
"nSamplePerBlock", nSamplePerBlock_);
90 coulombComponents_ =
false;
91 readOptional<bool>(in,
"coulombComponents", coulombComponents_);
97 kineticAveragePtr_ =
new Average;
99 potentialAveragePtr_ =
new Average;
113 angleAveragePtr_ =
new Average;
119 dihedralAveragePtr_ =
new Average;
125 if (coulombComponents_) {
126 coulombRSpaceAveragePtr_ =
new Average;
128 coulombKSpaceAveragePtr_ =
new Average;
131 coulombAveragePtr_ =
new Average;
137 externalAveragePtr_ =
new Average;
142 isInitialized_ =
true;
155 nSamplePerBlock_ = 0;
157 loadParameter<int>(ar,
"nSamplePerBlock", nSamplePerBlock_,
isRequired);
160 loadParameter<bool>(ar,
"coulombComponents", coulombComponents_,
isRequired);
165 totalAveragePtr_ =
new Average;
166 ar >> *totalAveragePtr_;
167 kineticAveragePtr_ =
new Average;
168 ar >> *kineticAveragePtr_;
169 potentialAveragePtr_ =
new Average;
170 ar >> *potentialAveragePtr_;
173 ar >> *pairAveragePtr_;
178 ar >> *bondAveragePtr_;
183 angleAveragePtr_ =
new Average;
184 ar >> *angleAveragePtr_;
189 dihedralAveragePtr_ =
new Average;
190 ar >> *dihedralAveragePtr_;
195 if (coulombComponents_) {
196 coulombRSpaceAveragePtr_ =
new Average;
197 ar >> *coulombRSpaceAveragePtr_;
198 coulombKSpaceAveragePtr_ =
new Average;
199 ar >> *coulombKSpaceAveragePtr_;
201 coulombAveragePtr_ =
new Average;
202 ar >> *coulombAveragePtr_;
207 externalAveragePtr_ =
new Average;
208 ar >> *externalAveragePtr_;
213 if (nSamplePerBlock_ != 0) {
217 isInitialized_ =
true;
228 bool isActive = bool(nSamplePerBlock_);
232 isActive = coulombComponents_;
238 ar << *totalAveragePtr_;
239 ar << *kineticAveragePtr_;
240 ar << *potentialAveragePtr_;
242 ar << *pairAveragePtr_;
246 ar << *bondAveragePtr_;
251 assert(angleAveragePtr_);
252 ar << *angleAveragePtr_;
257 assert(dihedralAveragePtr_);
258 ar << *dihedralAveragePtr_;
263 if (coulombComponents_) {
264 assert(coulombRSpaceAveragePtr_);
265 ar << *coulombRSpaceAveragePtr_;
266 assert(coulombKSpaceAveragePtr_);
267 ar << *coulombKSpaceAveragePtr_;
269 assert(coulombAveragePtr_);
270 ar << *coulombAveragePtr_;
275 assert(externalAveragePtr_);
276 ar << *externalAveragePtr_;
288 totalAveragePtr_->
clear();
289 kineticAveragePtr_->
clear();
290 potentialAveragePtr_->
clear();
292 pairAveragePtr_->
clear();
296 bondAveragePtr_->
clear();
301 assert(angleAveragePtr_);
302 angleAveragePtr_->
clear();
307 assert(dihedralAveragePtr_);
308 dihedralAveragePtr_->
clear();
313 if (coulombComponents_) {
314 assert(coulombRSpaceAveragePtr_);
315 coulombRSpaceAveragePtr_->
clear();
316 assert(coulombKSpaceAveragePtr_);
317 coulombKSpaceAveragePtr_->
clear();
319 assert(coulombAveragePtr_);
320 coulombAveragePtr_->
clear();
325 assert(externalAveragePtr_);
326 externalAveragePtr_->
clear();
339 if (outputFile_.is_open()) {
343 if (nSamplePerBlock_ > 0) {
346 std::string filename;
350 outputFile_ <<
" iStep";
351 outputFile_ <<
" Kinetic";
353 outputFile_ <<
" Pair";
357 outputFile_ <<
" Bond";
362 outputFile_ <<
" Angle";
367 outputFile_ <<
" Dihedral";
372 outputFile_ <<
" Coulomb";
377 outputFile_ <<
" External";
380 outputFile_ <<
" Total";
401 kineticAveragePtr_->
sample(kinetic);
404 double potential = 0.0;
408 pairAveragePtr_->
sample(pair);
415 bondAveragePtr_->
sample(bond);
423 assert(angleAveragePtr_);
424 angleAveragePtr_->
sample(angle);
431 potential += dihedral;
432 assert(dihedralAveragePtr_);
433 dihedralAveragePtr_->
sample(dihedral);
439 if (coulombComponents_) {
442 assert(coulombRSpaceAveragePtr_);
443 coulombRSpaceAveragePtr_->
sample(coulombRSpace);
446 assert(coulombKSpaceAveragePtr_);
447 coulombKSpaceAveragePtr_->
sample(coulombKSpace);
451 assert(coulombAveragePtr_);
452 coulombAveragePtr_->
sample(coulomb);
453 potential += coulomb;
460 potential += external;
461 assert(externalAveragePtr_);
462 externalAveragePtr_->
sample(external);
467 assert(potentialAveragePtr_);
468 potentialAveragePtr_->
sample(potential);
469 double total = kinetic + potential;
470 totalAveragePtr_->
sample(total);
476 int beginStep = iStep - (nSamplePerBlock_ - 1)*
interval();
477 outputFile_ <<
Int(beginStep, 12);
478 outputFile_ <<
Dbl(kineticAveragePtr_->blockAverage());
480 outputFile_ <<
Dbl(pairAveragePtr_->blockAverage());
484 outputFile_ <<
Dbl(bondAveragePtr_->blockAverage());
489 assert(angleAveragePtr_);
490 outputFile_ <<
Dbl(angleAveragePtr_->blockAverage());
495 assert(dihedralAveragePtr_);
496 outputFile_ <<
Dbl(dihedralAveragePtr_->blockAverage());
501 assert(coulombAveragePtr_);
502 outputFile_ <<
Dbl(coulombAveragePtr_->blockAverage());
507 assert(externalAveragePtr_);
508 outputFile_ <<
Dbl(externalAveragePtr_->blockAverage());
511 outputFile_ <<
Dbl(totalAveragePtr_->blockAverage());
527 if (outputFile_.is_open()) {
539 ave = kineticAveragePtr_->
average();
541 outputFile_ <<
"Kinetic " <<
Dbl(ave)
542 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
544 ave = pairAveragePtr_->
average();
546 outputFile_ <<
"Pair " <<
Dbl(ave)
547 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
551 ave = bondAveragePtr_->
average();
553 outputFile_ <<
"Bond " <<
Dbl(ave)
554 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
559 assert(angleAveragePtr_);
560 ave = angleAveragePtr_->
average();
562 outputFile_ <<
"Angle " <<
Dbl(ave)
563 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
568 assert(dihedralAveragePtr_);
569 ave = dihedralAveragePtr_->
average();
571 outputFile_ <<
"Dihedral " <<
Dbl(ave)
572 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
578 assert(coulombAveragePtr_);
579 ave = coulombAveragePtr_->
average();
581 outputFile_ <<
"Coulomb " <<
Dbl(ave)
582 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
584 if (coulombComponents_) {
585 assert(coulombRSpaceAveragePtr_);
586 ave = coulombRSpaceAveragePtr_->
average();
588 outputFile_ <<
"Coulomb(R) " <<
Dbl(ave)
589 <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
590 assert(coulombKSpaceAveragePtr_);
591 ave = coulombKSpaceAveragePtr_->
average();
593 outputFile_ <<
"Coulomb(K) " <<
Dbl(ave) <<
" +- " 594 <<
Dbl(err, 9, 2) <<
"\n";
600 assert(externalAveragePtr_);
601 ave = externalAveragePtr_->
average();
603 outputFile_ <<
"External " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
607 assert(potentialAveragePtr_);
608 ave = potentialAveragePtr_->
average();
610 outputFile_ <<
"Potential " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
612 ave = totalAveragePtr_->
average();
614 outputFile_ <<
"Total " <<
Dbl(ave) <<
" +- " <<
Dbl(err, 9, 2) <<
"\n";
621 "---------------------------------------------------------------------------------\n";
622 outputFile_ <<
"Kinetic:\n\n";
623 kineticAveragePtr_->
output(outputFile_);
625 "---------------------------------------------------------------------------------\n";
626 outputFile_ <<
"Pair:\n\n";
628 pairAveragePtr_->
output(outputFile_);
633 "---------------------------------------------------------------------------------\n";
634 outputFile_ <<
"Bond:\n\n";
635 bondAveragePtr_->
output(outputFile_);
641 "---------------------------------------------------------------------------------\n";
642 outputFile_ <<
"Angle:\n\n";
643 assert(angleAveragePtr_);
644 angleAveragePtr_->
output(outputFile_);
650 "---------------------------------------------------------------------------------\n";
651 outputFile_ <<
"Dihedral:\n\n";
652 assert(dihedralAveragePtr_);
653 dihedralAveragePtr_->
output(outputFile_);
659 "---------------------------------------------------------------------------------\n";
660 outputFile_ <<
"Coulomb:\n\n";
661 assert(coulombAveragePtr_);
662 coulombAveragePtr_->
output(outputFile_);
668 "---------------------------------------------------------------------------------\n";
669 outputFile_ <<
"External:\n\n";
670 assert(externalAveragePtr_);
671 externalAveragePtr_->
output(outputFile_);
676 "---------------------------------------------------------------------------------\n";
677 outputFile_ <<
"Potential:\n\n";
678 assert(potentialAveragePtr_);
679 potentialAveragePtr_->
output(outputFile_);
682 "---------------------------------------------------------------------------------\n";
683 outputFile_ <<
"Total:\n\n";
684 assert(totalAveragePtr_);
685 totalAveragePtr_->
output(outputFile_);
virtual void save(Serializable::OArchive &ar)
Load parameters to archive.
virtual void readParameters(std::istream &in)
Read interval and output file name.
void clear()
Clear all accumulators, set to empty initial state.
DihedralPotential & dihedralPotential() const
Return DihedralPotential by reference.
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.
double rSpaceEnergy()
Return short-range r-space part of Coulomb energy.
bool isRequired() const
Is this ParamComposite required in the input file?
virtual void save(Serializable::OArchive &ar)
Save internal state to an output archive.
Calculates the average and variance of a sampled property.
double average() const
Return the average of all sampled values.
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.
virtual void output()
Write file averages and error analysis to file.
MdSystem & 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.
bool hasBondPotential() const
Does a bond potential exist?.
ExternalPotential & externalPotential() const
Return ExternalPotential by reference.
The main object in a simulation, which coordinates others.
bool isActive() const
Is this parameter active?
bool hasExternalPotential() const
Does an external potential exist?.
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.
MdPairPotential & pairPotential() const
Return MdPairPotential by reference.
virtual double energy(const Vector &position, int i) const =0
Returns external potential energy of a single particle.
double kSpaceEnergy()
Get long-range k-space part of Coulomb energy.
Simulation & simulation() const
Get the parent Simulation by reference.
virtual void writeParam(std::ostream &out)
Write all parameters to an output stream.
void readInterval(std::istream &in)
Read interval from file, with error checking.
virtual void clear()
Clear nSample counter and all accumulators.
Utility classes for scientific computation.
virtual void loadParameters(Serializable::IArchive &ar)
Load internal state from an input archive.
virtual void sample(long iStep)
Write energy to file.
bool hasAnglePotential() const
Does angle potential exist?.
AnglePotential & anglePotential() const
Return AnglePotential by reference.
Wrapper for an int, for formatted ostream output.
bool hasCoulombPotential() const
Does a Coulomb potential exist?.
double kineticEnergy() const
Compute and return total kinetic energy.
Template for Analyzer associated with one System.
void output(std::ostream &out) const
Output final statistical properties to file.
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.
bool hasDihedralPotential() const
Does a dihedral potential exist?.
Saving archive for binary istream.
void sample(double value)
Add a sampled value to the ensemble.
Single-processor Monte Carlo (MC) and molecular dynamics (MD).
MdEnergyAnalyzer(MdSystem &system)
Constructor.
void setClassName(const char *className)
Set class name string.
virtual void setup()
Setup before main loop.
FileMaster & fileMaster()
Get the FileMaster by reference.
A System for Molecular Dynamics simulation.
bool isAtInterval(long counter) const
Return true iff counter is a multiple of the interval.
const std::string & outputFileName() const
Return outputFileName string.
int interval() const
Get interval value.
static void saveOptional(Serializable::OArchive &ar, Type &value, bool isActive)
Save an optional parameter value to an output archive.
double energy()
Get total Coulomb energy.
BondPotential & bondPotential() const
Return BondPotential by reference.
MdCoulombPotential & coulombPotential() const
Return CoulombPotential by reference.
FileMaster & fileMaster()
Get the FileMaster object.