SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
RouseModeCV.h
1 
21 #pragma once
22 
23 #include "CollectiveVariable.h"
24 #include "Validator/ObjectRequirement.h"
25 #include "Drivers/DriverException.h"
26 #include "Snapshot.h"
27 #include "schema.h"
28 
29 namespace SSAGES
30 {
32 
40  {
41  private:
42  std::vector<Label> groups_;
43  std::vector<double> massg_;
44  size_t p_;
45  size_t N_;
47  std::vector<Vector3> r_;
48 
49  public:
51 
55  RouseModeCV(const std::vector<Label>& groups, int p) :
56  groups_(groups), p_(p), N_( groups_.size())
57  { massg_.resize( groups_.size(),0.0 ); }
58 
60 
65  void setMasses(const std::vector<Label>& groups, const Snapshot& snapshot ) {
66  // Compute mass of each group and store to massg_
67  Label listi;
68  for (size_t i = 0; i < N_; ++i) {
69  listi = groups[i]; // MW: can this be used directly in Total Mass
70  Label idi;
71  snapshot.GetLocalIndices(listi, &idi);
72  massg_[i] = snapshot.TotalMass(idi);
73  }
74  }
75 
77 
80  void Initialize(const Snapshot& snapshot) override
81  {
82  using std::to_string;
83 
84  // Check for valid p_
85  if (p_ > groups_.size())
86  throw std::invalid_argument(
87  "RouseModeCV: Expected to find p to be less than " +
88  to_string(groups_.size()) +" but found p = " +
89  to_string(p_)
90  );
91  // Check for valid groups
92  for (size_t j = 0; j < groups_.size(); ++j) {
93  auto nj = groups_[j].size();
94  // Make sure atom IDs in the group are somewhere
95  std::vector<int> found(nj,0);
96  for (size_t i = 0; i < nj; ++i) {
97  if(snapshot.GetLocalIndex(groups_[j][i]) != -1)
98  found[i] = 1;
99  }
100 
101  MPI_Allreduce(MPI_IN_PLACE, found.data(), nj, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
102  unsigned njtot = std::accumulate(found.begin(), found.end(), 0, std::plus<int>());
103 
104  if(njtot != nj)
105  throw std::invalid_argument(
106  "RouseModeCV: Expected to find " +
107  to_string(nj) +
108  " atoms in group " + to_string(j) +", but only found " +
109  to_string(njtot) + "."
110  );
111 
112  }
113  // Set the masses of each particle group in massg_
114  this->setMasses(groups_, snapshot);
115  }
116 
118 
121  void Evaluate(const Snapshot& snapshot) override
122  {
123 
124  // Get data from snapshot.
125  auto ntot = snapshot.GetNumAtoms(); // total number of atoms
126  const auto& masses = snapshot.GetMasses(); // mass of each atom
127 
128  // Initialize working variables
129  double ppi_n = p_ * M_PI / N_; // constant
130  xp_.fill(0.0); // vectorial Rouse mode
131  r_.resize(N_); // position vector for beads in Rouse chain (unwrapped)
132  grad_.resize(ntot, Vector3{0,0,0}); // gradient set to 0.0 for all atoms
133  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
134 
135  // Iterate over all N_ atom groups and compute the center of mass for each
136  std::vector<Vector3> rcom; // vector of COM positions
137  for (size_t i = 0; i < N_; ++i) {
138  Label idi; // list of indices
139  snapshot.GetLocalIndices(groups_[i], &idi);
140  rcom.push_back(snapshot.CenterOfMass(idi,massg_[i])); // center of mass for group
141  }
142 
143  // Now compute differences vectors between the neighboring beads
144  // accumulate displacements to reconstruct unwrapped polymer chain
145  // for simplicity, we consider the first bead to be the reference position
146  // in all snapshots
147  r_[0] = rcom[0];
148  for (size_t i = 1; i < N_; ++i) {
149  Vector3 dri = snapshot.ApplyMinimumImage(rcom[i] - rcom[i-1]);
150  r_[i] = r_[i-1] + dri; // r_i = r_{i-1} + (r_i - r_{i-1})
151  }
152 
153  // Determine the value of the Rouse coordinate
154  // Xp(t) = 1/sqrt(N) * sum_{i=1}^{N} ri, p = 0
155  // Xp(t) = sqrt(2/N) * sum_{i=1}^{N} ri * cos[p*pi/N*(i-0.5)], p = 1,...,N-1
156  // Note: this solution is valid for homogeneous friction
157  xp_.fill(0.0);
158  for (size_t i = 0; i < N_; ++i) {
159  xp_ += r_[i]*cos ( ppi_n * (i+0.5) );
160  }
161  xp_ /= sqrt(N_);
162  if ( p_ != 0 ) xp_ *= sqrt(2.0) ;
163 
164  // Compute Rouse vector norm as the CV
165  // CV = sqrt(Xp*Xp), Xp = (Xp1,Xp2,Xp3)
166  val_ = xp_.norm();
167 
168  // Now perform gradient operation
169  // dCV/dxjd = (Xpd/CV)*(c/N)**0.5*sum_i=1^N cos[p*pi(i-0.5)/N] mj/Mi *delta_j({i})
170  // delta_j({i}) = 1, if j in {i}, 0 otherwise
171  Vector3 gradpcon = xp_ / sqrt(N_) / val_;
172  if ( p_ != 0) gradpcon *= sqrt(2.0); // (Xpd/CV)*(c/N)**0.5
173  for (size_t i = 0; i < N_; ++i) {
174  Label idi; // list of indices
175  snapshot.GetLocalIndices(groups_[i], &idi);
176  // go over each atom in the group and add to its gradient
177  // Note: performance tradeoff here. All gradient elements have a common factor of
178  // Xpd/CV*sqrt(c/N) = prefactor, with c = 2 if p != 0
179  // this could be factored out, but if ntot >> number of atoms in groups
180  // then it won't be worth it to post multiply all gradient terms...
181  // Could also go over loop again after to do the multiplication, but
182  // that is troublesome if atom ids appear in multiple groups for some reason
183  double cosval = cos(ppi_n*(i+0.5)) / massg_[i]; // cos[p*pi(i-0.5)/N] / Mi
184  for (auto& id : idi) {
185  grad_[id] += gradpcon* cosval * masses[id];
186  }
187  }
188 
189  }
190 
192  static RouseModeCV* Build(const Json::Value& json, const std::string& path)
193  {
194  Json::ObjectRequirement validator;
195  Json::Value schema;
196  Json::CharReaderBuilder rbuilder;
197  Json::CharReader* reader = rbuilder.newCharReader();
198 
199  reader->parse(JsonSchema::RouseModeCV.c_str(),
200  JsonSchema::RouseModeCV.c_str() + JsonSchema::RouseModeCV.size(),
201  &schema, NULL);
202  validator.Parse(schema, path);
203 
204  // Validate inputs.
205  validator.Validate(json, path);
206  if(validator.HasErrors())
207  throw BuildException(validator.GetErrors());
208 
209  std::vector<Label> groups;
210  for (auto& group : json["groups"])
211  {
212  groups.push_back({});
213  for(auto& id : group)
214  groups.back().push_back(id.asInt());
215  }
216 
217  auto mode = json.get("mode",0).asInt();
218  return new RouseModeCV( groups, mode);
219  }
220  };
221 }
std::vector< Vector3 > r_
vector of coordinate positions for each bead
Definition: RouseModeCV.h:47
std::vector< int > Label
List of integers.
Definition: types.h:48
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:200
void ApplyMinimumImage(Vector3 *v) const
Apply minimum image to a vector.
Definition: Snapshot.h:418
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
Definition: RouseModeCV.h:121
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:184
Collective variable is a Rouse mode for a polymer chain comprised of N particle groups.
Definition: RouseModeCV.h:39
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:43
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:555
std::vector< Label > groups_
vector of groups of indices to define the particle groups
Definition: RouseModeCV.h:42
double TotalMass(const Label &indices) const
Compute the total mass of a group of particles based on index.
Definition: Snapshot.h:451
std::vector< double > massg_
vector of the total mass for each particle group
Definition: RouseModeCV.h:43
size_t N_
number of Rouse beads in polymer chain
Definition: RouseModeCV.h:45
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
RouseModeCV(const std::vector< Label > &groups, int p)
Basic Constructor for Rouse Mode CV.
Definition: RouseModeCV.h:55
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
Exception to be thrown when building the Driver fails.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Requirements on an object.
Vector3 CenterOfMass(const Label &indices) const
Definition: Snapshot.h:466
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Definition: RouseModeCV.h:80
static RouseModeCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
Definition: RouseModeCV.h:192
size_t p_
index of mode of interest as CV
Definition: RouseModeCV.h:44
const std::vector< double > & GetMasses() const
Const access to the particle masses.
Definition: Snapshot.h:378
void setMasses(const std::vector< Label > &groups, const Snapshot &snapshot)
Helper function to determine masses of each group.
Definition: RouseModeCV.h:65
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
Vector3 xp_
3d vector for containing vectorial rouse amplitude
Definition: RouseModeCV.h:46
double val_
Current value of CV.
Abstract class for a collective variable.
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:573
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.