SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
GyrationTensorCV.h
1 
21 #pragma once
22 
23 #include "CVs/CollectiveVariable.h"
24 #include "Validator/ObjectRequirement.h"
25 #include "Drivers/DriverException.h"
26 #include "Snapshot.h"
27 #include "schema.h"
28 #include <Eigen/Eigenvalues>
29 
30 namespace SSAGES
31 {
33 
38  {
39  Rg = 0, // Radius of gyration (squared)
40  principal1 = 1, // First (largest) principal moment
41  principal2 = 2, // Second (middle) principal moment
42  principal3 = 3, // Third (smallest) principal moment
43  asphericity = 4, // Asphericity
44  acylindricity = 5, // Acylindricity
45  shapeaniso = 6 // Relative shape anisotropy
46  };
47 
49 
57  {
58  private:
61 
64 
65  public:
67 
74  GyrationTensorCV(const Label& atomids, GyrationTensor component) :
75  atomids_(atomids), component_(component), dim_{true, true, true}
76  {
77  }
78 
80 
90  GyrationTensorCV(const Label& atomids, GyrationTensor component, bool dimx, bool dimy, bool dimz) :
91  atomids_(atomids), component_(component), dim_{dimx, dimy, dimz}
92  {
93  }
94 
96 
99  void Initialize(const Snapshot& snapshot) override
100  {
101  using std::to_string;
102 
103  auto n = atomids_.size();
104 
105  // Make sure atom ID's are on at least one processor.
106  std::vector<int> found(n, 0);
107  for(size_t i = 0; i < n; ++i)
108  {
109  if(snapshot.GetLocalIndex(atomids_[i]) != -1)
110  found[i] = 1;
111  }
112 
113  MPI_Allreduce(MPI_IN_PLACE, found.data(), n, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
114  unsigned ntot = std::accumulate(found.begin(), found.end(), 0, std::plus<int>());
115  if(ntot != n)
116  throw BuildException({
117  "GyrationTensorCV: Expected to find " +
118  to_string(n) +
119  " atoms, but only found " +
120  to_string(ntot) + "."
121  });
122  }
123 
125 
128  void Evaluate(const Snapshot& snapshot) override
129  {
130  using namespace Eigen;
131 
132  // Get local atom indices and compute COM.
133  std::vector<int> idx;
134  snapshot.GetLocalIndices(atomids_, &idx);
135 
136  // Get data from snapshot.
137  auto n = snapshot.GetNumAtoms();
138  const auto& masses = snapshot.GetMasses();
139  const auto& pos = snapshot.GetPositions();
140 
141  // Initialize gradient.
142  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
143  grad_.resize(n, Vector3{0,0,0});
144 
145  // Compute total and center of mass.
146  auto masstot = snapshot.TotalMass(idx);
147  Vector3 com = snapshot.CenterOfMass(idx, masstot);
148 
149  // Gyration tensor and temporary vector to store positions in inertial frame.
150  Matrix3 S = Matrix3::Zero();
151  std::vector<Vector3> ris;
152  ris.reserve(idx.size());
153  for(auto& i : idx)
154  {
155  ris.emplace_back(snapshot.ApplyMinimumImage(pos[i] - com).cwiseProduct(dim_.cast<double>()));
156  S.noalias() += masses[i]*ris.back()*ris.back().transpose();
157  }
158 
159  // Reduce gyration tensor across processors and normalize.
160  MPI_Allreduce(MPI_IN_PLACE, S.data(), S.size(), MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());
161  S /= masstot;
162 
163  // Perform EVD. The columns are the eigenvectors.
164  // SelfAdjoint solver sorts in ascending order.
165  SelfAdjointEigenSolver<Matrix3> solver;
166  solver.computeDirect(S);
167  const auto& eigvals = solver.eigenvalues();
168  const auto& eigvecs = solver.eigenvectors();
169 
170  // Assign variables for clarity. l1 is largest.
171  auto l1 = eigvals[2],
172  l2 = eigvals[1],
173  l3 = eigvals[0];
174  const auto& n1 = eigvecs.col(2),
175  n2 = eigvecs.col(1),
176  n3 = eigvecs.col(0);
177 
178  val_ = 0;
179  auto sum = l1 + l2 + l3;
180  auto sqsum = l1*l1 + l2*l2 + l3*l3;
181  switch(component_)
182  {
183  case Rg:
184  val_ = sum;
185  break;
186  case principal1:
187  val_ = l1;
188  break;
189  case principal2:
190  val_ = l2;
191  break;
192  case principal3:
193  val_ = l3;
194  break;
195  case asphericity:
196  val_ = l1 - 0.5*(l2 + l3);
197  break;
198  case acylindricity:
199  val_ = l2 - l3;
200  break;
201  case shapeaniso:
202  val_ = 1.5*sqsum/(sum*sum) - 0.5;
203  break;
204  }
205 
206  // Compute gradient.
207  size_t j = 0;
208  for(auto& i : idx)
209  {
210  // Compute derivative of each eigenvalue and use combos in components.
211  auto dl1 = 2.*masses[i]/masstot*(1.-masses[i]/masstot)*ris[j].dot(n1)*n1;
212  auto dl2 = 2.*masses[i]/masstot*(1.-masses[i]/masstot)*ris[j].dot(n2)*n2;
213  auto dl3 = 2.*masses[i]/masstot*(1.-masses[i]/masstot)*ris[j].dot(n3)*n3;
214 
215  switch(component_)
216  {
217  case Rg:
218  grad_[i] = dl1 + dl2 + dl3;
219  break;
220  case principal1:
221  grad_[i] = dl1;
222  break;
223  case principal2:
224  grad_[i] = dl2;
225  break;
226  case principal3:
227  grad_[i] = dl3;
228  break;
229  case asphericity:
230  grad_[i] = dl1 - 0.5*(dl2 + dl3);
231  break;
232  case acylindricity:
233  grad_[i] = dl2 - dl3;
234  break;
235  case shapeaniso:
236  grad_[i] = 3.*(l1*dl1+l2*dl2+l3*dl3)/(sum*sum) -
237  3.*sqsum*(dl1+dl2+dl3)/(sum*sum*sum);
238  break;
239  }
240 
241  ++j;
242  }
243  }
244 
246  static GyrationTensorCV* Build(const Json::Value& json, const std::string& path)
247  {
248  Json::ObjectRequirement validator;
249  Json::Value schema;
250  Json::CharReaderBuilder rbuilder;
251  Json::CharReader* reader = rbuilder.newCharReader();
252 
253  reader->parse(JsonSchema::GyrationTensorCV.c_str(),
254  JsonSchema::GyrationTensorCV.c_str() + JsonSchema::GyrationTensorCV.size(),
255  &schema, NULL);
256  validator.Parse(schema, path);
257 
258  // Validate inputs.
259  validator.Validate(json, path);
260  if(validator.HasErrors())
261  throw BuildException(validator.GetErrors());
262 
263  std::vector<int> atomids;
264  for(auto& s : json["atom_ids"])
265  atomids.push_back(s.asInt());
266 
267  GyrationTensor component = Rg;
268  auto comp = json["component"].asString();
269 
270  if(comp == "Rg")
271  component = Rg;
272  else if(comp == "principal1")
273  component = principal1;
274  else if(comp == "principal2")
275  component = principal2;
276  else if(comp == "principal3")
277  component = principal3;
278  else if(comp == "asphericity")
279  component = asphericity;
280  else if(comp == "acylindricity")
281  component = acylindricity;
282  else if(comp == "shapeaniso")
283  component = shapeaniso;
284 
285  GyrationTensorCV* c;
286  if(json.isMember("dimension"))
287  {
288  auto dimx = json["dimension"][0].asBool();
289  auto dimy = json["dimension"][1].asBool();
290  auto dimz = json["dimension"][2].asBool();
291  c = new GyrationTensorCV(atomids, component, dimx, dimy, dimz);
292  }
293  else
294  c = new GyrationTensorCV(atomids, component);
295 
296  return c;
297  }
298  };
299 }
std::vector< int > Label
List of integers.
Definition: types.h:48
Eigen::Matrix3d Matrix3
3x3 matrix.
Definition: types.h:42
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
static GyrationTensorCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
GyrationTensor component_
Component of gyration tensor to compute.
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
GyrationTensorCV(const Label &atomids, GyrationTensor component)
Constructor.
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
Collective variable on components of the gyration tensor.
GyrationTensor
Define components of gyration tensor.
double TotalMass(const Label &indices) const
Compute the total mass of a group of particles based on index.
Definition: Snapshot.h:451
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.
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Exception to be thrown when building the Driver fails.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Label atomids_
IDs of the atoms used for calculation.
Requirements on an object.
Vector3 CenterOfMass(const Label &indices) const
Definition: Snapshot.h:466
GyrationTensorCV(const Label &atomids, GyrationTensor component, bool dimx, bool dimy, bool dimz)
Constructor.
const std::vector< double > & GetMasses() const
Const access to the particle masses.
Definition: Snapshot.h:378
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
double val_
Current value of CV.
Abstract class for a collective variable.
Bool3 dim_
Each dimension determines if it is applied by the CV.
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
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.
Eigen::Matrix< bool, 3, 1 > Bool3
Three-dimensional boolean.
Definition: types.h:36