SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
AngleCV.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 #include <array>
30 #include <cmath>
31 
32 namespace SSAGES
33 {
35  class AngleCV : public CollectiveVariable
36  {
37  private:
40 
41  public:
42 
44 
53  AngleCV(int atomid1, int atomid2, int atomid3) :
54  atomids_({atomid1, atomid2, atomid3})
55  {
56  bounds_ = {{0,M_PI}};
57  }
58 
60 
63  void Initialize(const Snapshot& snapshot) override
64  {
65  using std::to_string;
66 
67  std::vector<int> found;
68  snapshot.GetLocalIndices(atomids_, &found);
69  int nfound = found.size();
70  MPI_Allreduce(MPI_IN_PLACE, &nfound, 1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
71 
72  if(nfound != 3)
73  throw BuildException({
74  "TorsionalCV: Expected to find " +
75  to_string(3) +
76  " atoms, but only found " +
77  to_string(nfound) + "."
78  });
79  }
80 
82 
85  void Evaluate(const Snapshot& snapshot) override
86  {
87  // Get data from snapshot.
88  auto n = snapshot.GetNumAtoms();
89  const auto& pos = snapshot.GetPositions();
90  auto& comm = snapshot.GetCommunicator();
91 
92  // Initialize gradient.
93  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
94  grad_.resize(n, Vector3{0,0,0});
95 
96  Vector3 xi{0, 0, 0}, xj{0, 0, 0}, xk{0, 0, 0};
97 
98  auto iindex = snapshot.GetLocalIndex(atomids_[0]);
99  auto jindex = snapshot.GetLocalIndex(atomids_[1]);
100  auto kindex = snapshot.GetLocalIndex(atomids_[2]);
101 
102  if(iindex != -1) xi = pos[iindex];
103  if(jindex != -1) xj = pos[jindex];
104  if(kindex != -1) xk = pos[kindex];
105 
106  // By performing a reduce, we actually collect all. This can
107  // be converted to a more intelligent allgater on rank then bcast.
108  MPI_Allreduce(MPI_IN_PLACE, xi.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
109  MPI_Allreduce(MPI_IN_PLACE, xj.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
110  MPI_Allreduce(MPI_IN_PLACE, xk.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
111 
112  // Two vectors
113  auto rij = snapshot.ApplyMinimumImage(xi - xj);
114  auto rkj = snapshot.ApplyMinimumImage(xk - xj);
115 
116  auto dotP = rij.dot(rkj);
117  auto nrij = rij.norm();
118  auto nrkj = rkj.norm();
119 
120  val_ = acos(dotP/(nrij*nrkj));
121 
122  // Calculate gradients
123  auto prefactor = -1.0/(sqrt(1. - dotP/(nrij*nrkj))*nrij*nrkj);
124 
125  Vector3 gradi{0, 0, 0}, gradk{0, 0, 0};
126  if(iindex != -1) gradi = prefactor*(rkj - dotP*rij/(nrij*nrij));
127  if(kindex != -1) gradk = prefactor*(rij - dotP*rkj/(nrkj*nrkj));
128  MPI_Allreduce(MPI_IN_PLACE, gradi.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
129  MPI_Allreduce(MPI_IN_PLACE, gradk.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
130  if(iindex != -1) grad_[iindex] = gradi;
131  if(kindex != -1) grad_[kindex] = gradk;
132  if(jindex != -1) grad_[jindex] = -gradi - gradk;
133  }
134 
136  static AngleCV* Build(const Json::Value& json, const std::string& path)
137  {
138  Json::ObjectRequirement validator;
139  Json::Value schema;
140  Json::CharReaderBuilder rbuilder;
141  Json::CharReader* reader = rbuilder.newCharReader();
142 
143  reader->parse(JsonSchema::AngleCV.c_str(),
144  JsonSchema::AngleCV.c_str() + JsonSchema::AngleCV.size(),
145  &schema, NULL);
146  validator.Parse(schema, path);
147 
148  // Validate inputs.
149  validator.Validate(json, path);
150  if(validator.HasErrors())
151  throw BuildException(validator.GetErrors());
152 
153  std::vector<int> atomids;
154  for(auto& s : json["atom_ids"])
155  atomids.push_back(s.asInt());
156 
157  return new AngleCV(atomids[0], atomids[1], atomids[2]);
158  }
159  };
160 }
std::vector< int > Label
List of integers.
Definition: types.h:48
AngleCV(int atomid1, int atomid2, int atomid3)
Constructor.
Definition: AngleCV.h:53
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
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:184
std::array< double, 2 > bounds_
Bounds on CV.
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
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Definition: AngleCV.h:63
Collective variable to calculate angle.
Definition: AngleCV.h:35
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
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.
static AngleCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
Definition: AngleCV.h:136
Label atomids_
Vector of 3 atom ID&#39;s of interest.
Definition: AngleCV.h:39
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
double val_
Current value of CV.
Abstract class for a collective variable.
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
Definition: AngleCV.h:85
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
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.