SSAGES  0.8.5
Software Suite for Advanced General Ensemble Simulations
TorsionalCV.h
1 
22 #pragma once
23 
24 #include "CollectiveVariable.h"
25 #include "Validator/ObjectRequirement.h"
26 #include "Drivers/DriverException.h"
27 #include "Snapshot.h"
28 #include "schema.h"
29 #include <array>
30 #include <cmath>
31 
32 namespace SSAGES
33 {
35 
43  {
44  private:
47 
49  bool periodic_;
50 
51  public:
53 
64  TorsionalCV(int atomid1, int atomid2, int atomid3, int atomid4, bool periodic) :
65  atomids_({atomid1, atomid2, atomid3, atomid4}), periodic_(periodic)
66  {
67  }
68 
70 
73  void Initialize(const Snapshot& snapshot) override
74  {
75  using std::to_string;
76 
77  std::vector<int> found;
78  snapshot.GetLocalIndices(atomids_, &found);
79  int nfound = found.size();
80  MPI_Allreduce(MPI_IN_PLACE, &nfound, 1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
81 
82  if(nfound != 4)
83  throw BuildException({
84  "TorsionalCV: Expected to find " +
85  to_string(4) +
86  " atoms, but only found " +
87  to_string(nfound) + "."
88  });
89  }
90 
92 
95  void Evaluate(const Snapshot& snapshot) override
96  {
97  // Get data from snapshot.
98  auto n = snapshot.GetNumAtoms();
99  const auto& pos = snapshot.GetPositions();
100  auto& comm = snapshot.GetCommunicator();
101 
102  // Initialize gradient.
103  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
104  grad_.resize(n, Vector3{0,0,0});
105 
106  Vector3 xi{0, 0, 0}, xj{0, 0, 0}, xk{0, 0, 0}, xl{0, 0, 0};
107 
108  auto iindex = snapshot.GetLocalIndex(atomids_[0]);
109  auto jindex = snapshot.GetLocalIndex(atomids_[1]);
110  auto kindex = snapshot.GetLocalIndex(atomids_[2]);
111  auto lindex = snapshot.GetLocalIndex(atomids_[3]);
112 
113  if(iindex != -1) xi = pos[iindex];
114  if(jindex != -1) xj = pos[jindex];
115  if(kindex != -1) xk = pos[kindex];
116  if(lindex != -1) xl = pos[lindex];
117 
118  // By performing a reduce, we actually collect all. This can
119  // be converted to a more intelligent allgater on rank then bcast.
120  MPI_Allreduce(MPI_IN_PLACE, xi.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
121  MPI_Allreduce(MPI_IN_PLACE, xj.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
122  MPI_Allreduce(MPI_IN_PLACE, xk.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
123  MPI_Allreduce(MPI_IN_PLACE, xl.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
124 
125  //Calculate pertinent vectors
126  auto F = snapshot.ApplyMinimumImage(xi - xj);
127  auto G = snapshot.ApplyMinimumImage(xj - xk);
128  auto H = snapshot.ApplyMinimumImage(xl - xk);
129  auto A = F.cross(G);
130  auto B = H.cross(G);
131 
132  auto y = B.cross(A).dot(G.normalized());
133  auto x = A.dot(B);
134 
135  val_ = atan2(y, x);
136 
137  auto Zed = F.dot(G.normalized())/A.dot(A);
138  auto Ned = H.dot(G.normalized())/B.dot(B);
139 
140  Vector3 gradi{0,0,0}, gradl{0,0,0};
141  if(iindex != -1) gradi = -G.norm()*A/A.dot(A);
142  if(lindex != -1) gradl = G.norm()*B/B.dot(B);
143  MPI_Allreduce(MPI_IN_PLACE, gradi.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
144  MPI_Allreduce(MPI_IN_PLACE, gradl.data(), 3, MPI_DOUBLE, MPI_SUM, comm);
145  if(iindex != -1) grad_[iindex] = gradi;
146  if(lindex != -1) grad_[lindex] = gradl;
147  if(jindex != -1) grad_[jindex] = Zed*A - Ned*B - gradi;
148  if(kindex != -1) grad_[kindex] = Ned*B - Zed*A - gradl;
149  }
150 
152 
161  double GetPeriodicValue(double Location) const override
162  {
163  if(!periodic_)
164  return Location;
165 
166  int n = (int)(Location/(2.0*M_PI));
167  double PeriodicLocation;
168 
169  PeriodicLocation = Location - n*M_PI;
170  if(PeriodicLocation < -M_PI)
171  PeriodicLocation += 2.0*M_PI;
172  else if (PeriodicLocation > M_PI)
173  PeriodicLocation -= 2.0*M_PI;
174 
175  return PeriodicLocation;
176  }
177 
179 
188  double GetDifference(const double Location) const override
189  {
190  double PeriodicDiff = val_ - Location;
191 
192  if(!periodic_)
193  return PeriodicDiff;
194 
195  PeriodicDiff = GetPeriodicValue(PeriodicDiff);
196 
197  if(PeriodicDiff > M_PI)
198  PeriodicDiff -= 2.0*M_PI;
199  else if(PeriodicDiff < -M_PI)
200  PeriodicDiff += 2.0*M_PI;
201 
202  return PeriodicDiff;
203  }
204 
206 
213  double GetMinimumImage(const double Location) const override
214  {
215  double PeriodicDiff = val_ - Location;
216 
217  if(!periodic_)
218  return val_;
219 
220  if(PeriodicDiff > M_PI)
221  return (val_ - 2.0*M_PI);
222  else if(PeriodicDiff < -M_PI)
223  return (val_ + 2.0*M_PI);
224 
225  return val_;
226  }
227 
229  static TorsionalCV* Build(const Json::Value& json, const std::string& path)
230  {
231  Json::ObjectRequirement validator;
232  Json::Value schema;
233  Json::CharReaderBuilder rbuilder;
234  Json::CharReader* reader = rbuilder.newCharReader();
235 
236  reader->parse(JsonSchema::TorsionalCV.c_str(),
237  JsonSchema::TorsionalCV.c_str() + JsonSchema::TorsionalCV.size(),
238  &schema, nullptr);
239  validator.Parse(schema, path);
240 
241  // Validate inputs.
242  validator.Validate(json, path);
243  if(validator.HasErrors())
244  throw BuildException(validator.GetErrors());
245 
246  std::vector<int> atomids;
247  for(auto& s : json["atom_ids"])
248  atomids.push_back(s.asInt());
249 
250  auto periodic = json.get("periodic", true).asBool();
251  return new TorsionalCV(atomids[0], atomids[1], atomids[2], atomids[3], periodic);
252  }
253  };
254 }
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:204
bool periodic_
If True, use periodic boundary conditions.
Definition: TorsionalCV.h:49
double GetMinimumImage(const double Location) const override
Returns the minimum image of a CV based on the input location.
Definition: TorsionalCV.h:213
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:47
Label atomids_
Vector of 4 atom ID&#39;s of interest.
Definition: TorsionalCV.h:46
double GetPeriodicValue(double Location) const override
Return value taking periodic boundary conditions into account.
Definition: TorsionalCV.h:161
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
int GetLocalIndex(int id) const
Gets the local atom index corresponding to an atom ID.
Definition: Snapshot.h:616
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
const mxx::comm & GetCommunicator() const
Get communicator for walker.
Definition: Snapshot.h:188
static TorsionalCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
Definition: TorsionalCV.h:229
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:634
Exception to be thrown when building the Driver fails.
Collective variable on the torsion angles.
Definition: TorsionalCV.h:42
TorsionalCV(int atomid1, int atomid2, int atomid3, int atomid4, bool periodic)
Constructor.
Definition: TorsionalCV.h:64
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Requirements on an object.
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Definition: TorsionalCV.h:73
double GetDifference(const double Location) const override
Get difference taking periodic boundary conditions into account.
Definition: TorsionalCV.h:188
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
Definition: TorsionalCV.h:95
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
void ApplyMinimumImage(Vector3 *v) const
Apply minimum image to a vector.
Definition: Snapshot.h:422
double val_
Current value of CV.
Abstract class for a collective variable.
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:327
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.