SSAGES  0.8.5
Software Suite for Advanced General Ensemble Simulations
PairwiseCV.h
1 
21 #pragma once
22 
23 #include "CVs/CollectiveVariable.h"
24 #include "Validator/ObjectRequirement.h"
25 #include "Drivers/DriverException.h"
26 #include "Utility/PairwiseKernel.h"
27 #include "Snapshot.h"
28 #include "schema.h"
29 
30 namespace SSAGES
31 {
33 
41  {
42  private:
46 
47  public:
49 
57  PairwiseCV(const Label& group1, const Label& group2, PairwiseKernel* pk) :
58  group1_(group1), group2_(group2), pk_(pk)
59  {
60  }
61 
63 
66  void Initialize(const Snapshot& snapshot) override
67  {
68  using std::to_string;
69 
70  auto n1 = group1_.size(), n2 = group2_.size();
71 
72  // Make sure atom ID's are on at least one processor.
73  std::vector<int> found1(n1, 0), found2(n2, 0);
74  for(size_t i = 0; i < n1; ++i)
75  {
76  if(snapshot.GetLocalIndex(group1_[i]) != -1)
77  found1[i] = 1;
78  }
79 
80  for(size_t i = 0; i < n2; ++i)
81  {
82  if(snapshot.GetLocalIndex(group2_[i]) != -1)
83  found2[i] = 1;
84  }
85 
86  MPI_Allreduce(MPI_IN_PLACE, found1.data(), n1, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
87  MPI_Allreduce(MPI_IN_PLACE, found2.data(), n2, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
88 
89  unsigned ntot1 = std::accumulate(found1.begin(), found1.end(), 0, std::plus<int>());
90  unsigned ntot2 = std::accumulate(found2.begin(), found2.end(), 0, std::plus<int>());
91  if(ntot1 != n1)
92  throw BuildException({
93  "Pairwise: Expected to find " +
94  to_string(n1) +
95  " atoms in group 1, but only found " +
96  to_string(ntot1) + "."
97  });
98 
99  if(ntot2 != n2)
100  throw BuildException({
101  "Pairwise: Expected to find " +
102  to_string(n2) +
103  " atoms in group 2, but only found " +
104  to_string(ntot2) + "."
105  });
106  }
107 
109 
112  void Evaluate(const Snapshot& snapshot) override
113  {
114  // Get local atom indices.
115  std::vector<int> idx1, idx2;
116  snapshot.GetLocalIndices(group1_, &idx1);
117  snapshot.GetLocalIndices(group2_, &idx2);
118 
119  // Get data from snapshot.
120  auto n = snapshot.GetNumAtoms();
121  auto& atomids = snapshot.GetAtomIDs();
122  auto& positions = snapshot.GetPositions();
123  auto& comm = snapshot.GetCommunicator();
124 
125  // Initialize gradient.
126  std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
127  grad_.resize(n, Vector3{0,0,0});
128  boxgrad_ = Matrix3::Zero();
129 
130  // Fill eigen matrix with coordinate data.
131  std::vector<double> pos1(3*idx1.size()), pos2(3*idx2.size());
132  std::vector<int> id1(idx1.size()), id2(idx2.size());
133  for(size_t i = 0; i < idx1.size(); ++i)
134  {
135  pos1[3*i+0] = positions[idx1[i]][0];
136  pos1[3*i+1] = positions[idx1[i]][1];
137  pos1[3*i+2] = positions[idx1[i]][2];
138 
139  id1[i] = atomids[idx1[i]];
140  }
141 
142  for(size_t i = 0; i < idx2.size(); ++i)
143  {
144  pos2[3*i+0] = positions[idx2[i]][0];
145  pos2[3*i+1] = positions[idx2[i]][1];
146  pos2[3*i+2] = positions[idx2[i]][2];
147 
148  id2[i] = atomids[idx2[i]];
149  }
150 
151  // Gather across all procs.
152  pos1 = mxx::allgatherv(pos1.data(), pos1.size(), comm);
153  pos2 = mxx::allgatherv(pos2.data(), pos2.size(), comm);
154  id1 = mxx::allgatherv(id1.data(), id1.size(), comm);
155  id2 = mxx::allgatherv(id2.data(), id2.size(), comm);
156 
157  // Create Eigen map for ease of use.
158  using Map = Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor>>;
159  Map gpos1(pos1.data(), group1_.size(), 3), gpos2(pos2.data(), group2_.size(), 3);
160 
161  val_ = 0;
162  double df = 0.;
163  // Compute gradient and value.
164  for(size_t i = 0; i < group1_.size(); ++i)
165  {
166  auto t1 = id1[i];
167  const auto& pi = gpos1.row(i);
168  for(size_t j = 0; j < group2_.size(); ++j)
169  {
170  auto t2 = id2[j];
171  const auto& pj = gpos2.row(j);
172 
173  // Skip identical pairs.
174  if(t1 == t2)
175  continue;
176 
177  // Compute distance and pairwise function.
178  Vector3 rij = pi - pj;
179  snapshot.ApplyMinimumImage(&rij);
180  auto r = rij.norm();
181  val_ += pk_->Evaluate(r, df);
182 
183  // Get local indices and sum gradients.
184  auto lid1 = snapshot.GetLocalIndex(t1);
185  if(lid1 != -1)
186  grad_[lid1] += df*rij/r;
187 
188  auto lid2 = snapshot.GetLocalIndex(t2);
189  if(lid2 != -1)
190  grad_[lid2] -= df*rij/r;
191 
192  // Only sum box gradient on a single proc.
193  if(comm.rank() == 0)
194  boxgrad_ += df*rij/r*rij.transpose();
195  }
196  }
197  }
198 
200  static PairwiseCV* Build(const Json::Value& json, const std::string& path)
201  {
202  Json::ObjectRequirement validator;
203  Json::Value schema;
204  Json::CharReaderBuilder rbuilder;
205  Json::CharReader* reader = rbuilder.newCharReader();
206 
207  reader->parse(JsonSchema::PairwiseCV.c_str(),
208  JsonSchema::PairwiseCV.c_str() + JsonSchema::PairwiseCV.size(),
209  &schema, nullptr);
210  validator.Parse(schema, path);
211 
212  // Validate inputs.
213  validator.Validate(json, path);
214  if(validator.HasErrors())
215  throw BuildException(validator.GetErrors());
216 
217  std::vector<int> group1, group2;
218 
219  for(auto& s : json["group1"])
220  group1.push_back(s.asInt());
221 
222  for(auto& s : json["group2"])
223  group2.push_back(s.asInt());
224 
225  return new PairwiseCV(group1, group2, PairwiseKernel::Build(json["kernel"], path));
226  }
227 
230  {
231  delete pk_;
232  }
233  };
234 }
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
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:47
Label group1_
IDs of the first group of atoms.
Definition: PairwiseCV.h:43
~PairwiseCV()
Destructor.
Definition: PairwiseCV.h:229
static PairwiseCV * Build(const Json::Value &json, const std::string &path)
Set up collective variable.
Definition: PairwiseCV.h:200
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
Map for histogram and coefficients.
Definition: Basis.h:41
void GetLocalIndices(const Label &ids, Label *indices) const
Definition: Snapshot.h:634
Exception to be thrown when building the Driver fails.
virtual double Evaluate(double rij, double &df) const =0
Evaluate the pairwise kernel function.
const Label & GetAtomIDs() const
Access the atom IDs.
Definition: Snapshot.h:602
PairwiseCV(const Label &group1, const Label &group2, PairwiseKernel *pk)
Constructor.
Definition: PairwiseCV.h:57
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Definition: PairwiseCV.h:66
Requirements on an object.
Pairwise kernel base class.
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
Definition: PairwiseCV.h:112
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.
Matrix3 boxgrad_
Gradient w.r.t box vectors dCv/dHij.
Generalized collective variable based on pairwise properties of atoms.
Definition: PairwiseCV.h:40
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:327
PairwiseKernel * pk_
Pairwise kernel function used for CV.
Definition: PairwiseCV.h:45
Label group2_
IDs of the second group of atoms.
Definition: PairwiseCV.h:44
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.
static PairwiseKernel * Build(const Json::Value &json, const std::string &path)
Build PairwiseKernel from JSON value.