Simpatico  v1.10
IntraBondStressAutoCorr.tpp
1 #ifndef MCMD_INTRA_BOND_STRESS_AUTO_CORR_INC_H
2 #define MCMD_INTRA_BOND_STRESS_AUTO_CORR_INC_H
3 
4 /*
5 * Simpatico - Simulation Package for Polymeric and Molecular Liquids
6 *
7 * Copyright 2010 - 2017, The Regents of the University of Minnesota
8 * Distributed under the terms of the GNU General Public License.
9 */
10 
11 #include "IntraBondStressAutoCorr.h"
12 #include <mcMd/simulation/Simulation.h>
13 #include <mcMd/simulation/stress.h>
14 #include <mcMd/chemistry/Molecule.h>
15 #include <mcMd/chemistry/Bond.h>
16 #include <mcMd/chemistry/Atom.h>
17 #include <simp/species/Species.h>
18 #include <simp/boundary/Boundary.h>
19 #include <util/misc/FileMaster.h>
20 
21 #include <util/global.h>
22 
23 namespace McMd
24 {
25 
26  using namespace Util;
27  using namespace Simp;
28 
29  /*
30  * Constructor.
31  */
32  template <class SystemType>
34  : SystemAnalyzer<SystemType>(system),
35  outputFile_(),
36  accumulator_(),
37  speciesId_(-1),
38  nMolecule_(-1),
39  nBond_(-1),
40  capacity_(-1),
41  isInitialized_(false)
42  {}
43 
44  /*
45  * Destructor.
46  */
47  template <class SystemType>
49  {}
50 
51  /*
52  * Read parameters from file.
53  */
54  template <class SystemType>
56  {
57  readInterval(in);
59  read(in, "speciesId", speciesId_);
60  read(in, "capacity", capacity_);
61 
62  // Validate input
63  if (speciesId_ < 0)
64  UTIL_THROW("Negative speciesId");
65  if (speciesId_ >= system().simulation().nSpecies())
66  UTIL_THROW("speciesId >= nSpecies");
67  if (capacity_ <= 0)
68  UTIL_THROW("Negative capacity");
69 
70  speciesPtr_ = &system().simulation().species(speciesId_);
71  nBond_ = speciesPtr_->nBond();
72  if (nBond_ <= 0)
73  UTIL_THROW("Number of bonds per molecule <= 0");
74 
75  // Allocate memory
76  int speciesCapacity = speciesPtr_->capacity();
77  data_.allocate(speciesCapacity);
78  accumulator_.setParam(speciesCapacity, capacity_);
79 
80  // Delay initialization of AutoCorrArray until first call
81  // of sample(), when the number of molecules is known.
82 
83  isInitialized_ = true;
84  }
85 
86  /*
87  * Load internal state from file.
88  */
89  template <class SystemType> void
91  {
93  loadParameter(ar, "speciesId", speciesId_);
94  loadParameter(ar, "capacity", capacity_);
95  ar & nBond_;
96  ar & nMolecule_;
97  ar & accumulator_;
98 
99  // Validate parameters and set speciesPtr_
100  if (speciesId_ < 0) {
101  UTIL_THROW("Negative speciesId");
102  }
103  if (speciesId_ >= system().simulation().nSpecies()) {
104  UTIL_THROW("speciesId >= nSpecies");
105  }
106  if (capacity_ <= 0) {
107  UTIL_THROW("Negative capacity");
108  }
109  if (nBond_ <= 0) {
110  UTIL_THROW("Number of bonds per molecule <= 0");
111  }
112  speciesPtr_ = &system().simulation().species(speciesId_);
113  if (nBond_ != speciesPtr_->nBond()) {
114  UTIL_THROW("Inconsistent values for nBond");
115  }
116 
117  // Allocate data_ array for internal usage.
118  int speciesCapacity = speciesPtr_->capacity();
119  if (speciesCapacity <= 0) {
120  UTIL_THROW("Species capacity <= 0");
121  }
122  data_.allocate(speciesCapacity);
123 
124  isInitialized_ = true;
125  }
126 
127  /*
128  * Save internal state to an archive.
129  */
130  template <class SystemType> void
132  { ar & *this; }
133 
134  /*
135  * Evaluate end-to-end vectors of all chains, add to ensemble.
136  */
137  template <class SystemType>
139  {
140 
141  if (!isInitialized_) {
142  UTIL_THROW("Object is not intitialized");
143  }
144 
145  // Get the actual number of molecules.
146  nMolecule_ = system().nMolecule(speciesId_);
147 
148  // Initialize the AutoCorrArray object
149  accumulator_.setNEnsemble(nMolecule_);
150  accumulator_.clear();
151  }
152 
153  /*
154  * Evaluate end-to-end vectors of all chains, add to ensemble.
155  */
156  template <class SystemType>
158  {
159  if (isAtInterval(iStep)) {
160 
161  Tensor stress;
162  Vector dr, f;
163  double rsq, pressure;
164  int i, j;
165 
166  // Confirm that nMolecule has remained constant
167  if (nMolecule_ != system().nMolecule(speciesId_)) {
168  UTIL_THROW("Number of molecules has changed.");
169  }
170 
171  // Loop over molecules
174  i = 0;
175  system().begin(speciesId_, molIter);
176  for ( ; molIter.notEnd(); ++ molIter) {
177  data_[i].zero();
178 
179  molIter->begin(bondIter);
180  for ( ; bondIter.notEnd(); ++bondIter) {
181  rsq = system().boundary().distanceSq(
182  bondIter->atom(0).position(),
183  bondIter->atom(1).position(), dr);
184  f = dr;
185  f *= system().bondPotential()
186  .forceOverR(rsq, bondIter->typeId());
187  incrementPairStress(f, dr, data_[i]);
188  }
189 
190  // Remove trace
191  pressure = data_[i].trace() / double(Dimension);
192  for (j = 0; j < Dimension; ++j) {
193  data_[i](j, j) -= pressure;
194  }
195 
196  ++i;
197  }
198 
199  accumulator_.sample(data_);
200 
201  } // if isAtInterval
202 
203  }
204 
206  template <class SystemType>
208  {
209 
210  // Echo parameters to analyzer log file
211  fileMaster().openOutputFile(outputFileName(), outputFile_);
212  writeParam(outputFile_);
213  outputFile_.close();
214 
215  // Output statistical analysis to separate data file
216  fileMaster().openOutputFile(outputFileName(".dat"), outputFile_);
217  accumulator_.output(outputFile_);
218  outputFile_.close();
219 
220  }
221 
222 }
223 #endif
const int Dimension
Dimensionality of space.
Definition: Dimension.h:19
A Vector is a Cartesian vector.
Definition: Vector.h:75
virtual void readParameters(std::istream &in)
Read parameters from file.
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.
IntraBondStressAutoCorr(SystemType &system)
Constructor.
SystemType & system()
Return reference to parent system.
Forward iterator for a PArray.
File containing preprocessor macros for error handling.
Classes used by all simpatico molecular simulations.
void readOutputFileName(std::istream &in)
Read outputFileName from file.
A Tensor represents a Cartesian tensor.
Definition: Tensor.h:32
Forward const iterator for an Array or a C array.
Saving / output archive for binary ostream.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
Definition: global.h:51
virtual void writeParam(std::ostream &out)
Write all parameters to an output stream.
virtual void save(Serializable::OArchive &ar)
Save internal state to archive.
void readInterval(std::istream &in)
Read interval from file, with error checking.
Utility classes for scientific computation.
Definition: accumulators.mod:1
virtual void output()
Output results after simulation is completed.
ScalarParam< Type > & loadParameter(Serializable::IArchive &ar, const char *label, Type &value, bool isRequired)
Add and load a new ScalarParam < Type > object.
virtual void loadParameters(Serializable::IArchive &ar)
Load internal state from archive.
Template for Analyzer associated with one System.
bool notEnd() const
Is the current pointer not at the end of the array?
Saving archive for binary istream.
Single-processor Monte Carlo (MC) and molecular dynamics (MD).
int capacity() const
Maximum allowed number of molecules for this Species.
FileMaster & fileMaster()
Get the FileMaster by reference.
bool isAtInterval(long counter) const
Return true iff counter is a multiple of the interval.
virtual void setup()
Set number of molecules and clear accumulator.
const std::string & outputFileName() const
Return outputFileName string.
bool notEnd() const
Is this not the end of the array?
virtual ~IntraBondStressAutoCorr()
Destructor.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
int nBond() const
Get number of bonds per molecule for this Species.
virtual void sample(long iStep)
Evaluate end-to-end vectors of all chains, add to ensemble.