SSAGES  0.1
A MetaDynamics Package
Meta.cpp
1 
22 #include "Meta.h"
23 #include <math.h>
24 #include <iostream>
25 #include "Drivers/DriverException.h"
26 #include "CVs/CVManager.h"
27 #include "Validator/ObjectRequirement.h"
28 #include "Drivers/DriverException.h"
29 #include "Snapshot.h"
30 #include "schema.h"
31 
32 using namespace Json;
33 
34 namespace SSAGES
35 {
37 
42  double gaussian(double dx, double sigma)
43  {
44  double arg = (dx * dx) / (2. * sigma * sigma);
45  return exp(-arg);
46  }
47 
49 
54  double gaussianDerv(double dx, double sigma)
55  {
56  double arg = (dx * dx) / (2. * sigma * sigma);
57  double pre = - dx / (sigma * sigma);
58  return pre * exp(-arg);
59  }
60 
61  // Pre-simulation hook.
62  void Meta::PreSimulation(Snapshot* snapshot, const CVManager& cvmanager)
63  {
64  auto cvs = cvmanager.GetCVs(cvmask_);
65  // Write ouput file header.
66  if(world_.rank() == 0)
67  {
68  hillsout_.open("hills.out");
69  hillsout_ << "#Iteration ";
70 
71  for(size_t i = 0; i < cvs.size(); ++i)
72  hillsout_ << "center." << i << " ";
73 
74  for(size_t i = 0; i < cvs.size(); ++i)
75  hillsout_ << "sigma." << i << " ";
76 
77  hillsout_ << "height" << std::endl;
78 
79  hillsout_.close();
80  }
81 
82  // Initialize grid to zero.
83  if(grid_ != nullptr)
84  {
85  Vector vec = Vector::Zero(cvs.size());
86  std::fill(grid_->begin(), grid_->end(), vec);
87  }
88 
89  auto n = snapshot->GetTargetIterations();
90  n = n ? n : 1e5; // Pre-allocate at least something.
91 
92  hills_.reserve(n+1);
93  widths_.reserve(n+1);
94  derivatives_.resize(cvs.size());
95  tder_.resize(cvs.size());
96  dx_.resize(cvs.size());
97  }
98 
99  // Post-integration hook.
100  void Meta::PostIntegration(Snapshot* snapshot, const CVManager& cvmanager)
101  {
102  auto cvs = cvmanager.GetCVs(cvmask_);
103  // Add hills when needed.
104  if(snapshot->GetIteration() % hillfreq_ == 0)
105  AddHill(cvs, snapshot->GetIteration());
106 
107  // Always calculate the current bias.
108  CalcBiasForce(cvs);
109 
110  // Take each CV and add its biased forces to the atoms
111  // using the chain rule.
112  auto& forces = snapshot->GetForces();
113  auto& virial = snapshot->GetVirial();
114 
115  for(size_t i = 0; i < cvs.size(); ++i)
116  {
117  auto& grad = cvs[i]->GetGradient();
118  auto& boxgrad = cvs[i]->GetBoxGradient();
119 
120  // Update the forces in snapshot by adding in the force bias from each
121  // CV to each atom based on the gradient of the CV.
122  for (size_t j = 0; j < forces.size(); ++j)
123  for(size_t k = 0; k < 3; ++k)
124  forces[j][k] -= derivatives_[i]*grad[j][k];
125 
126  virial += derivatives_[i]*boxgrad;
127  }
128  }
129 
130  // Post-simulation hook.
131  void Meta::PostSimulation(Snapshot*, const CVManager&)
132  {
133  }
134 
135  // Load hills from file.
136  void Meta::LoadHills(const std::string& filename)
137  {
138  std::ifstream file(filename);
139  std::string line;
140 
141  auto dim = widths_.size();
142  double iteration = 0, height = 0;
143  std::vector<double> width(dim, 0.), center(dim, 0);
144 
145  // First line is a comment.
146  std::getline(file, line);
147  while(std::getline(file, line))
148  {
149  std::istringstream iss(line);
150  iss >> iteration; // Not really using this.
151 
152  // Load centers.
153  for(size_t i = 0; i < dim; ++i)
154  iss >> center[i];
155 
156  // Load widths.
157  for(size_t i = 0; i < dim; ++i)
158  iss >> width[i];
159 
160  // Load height.
161  iss >> height;
162 
163  hills_.emplace_back(center, width, height);
164  }
165  }
166 
167  // Drop a new hill.
168  void Meta::AddHill(const CVList& cvs, int iteration)
169  {
170  int n = cvs.size();
171 
172  // Assume we have the same number of procs per walker.
173  int nwalkers = world_.size()/comm_.size();
174 
175  // We need to exchange CV values across the walkers
176  // and to each proc on a walker.
177  std::vector<double> cvals(n*nwalkers, 0);
178 
179  if(comm_.rank() == 0)
180  {
181  for(auto i = 0, j = world_.rank()/comm_.size()*n; i < n; ++i,++j)
182  cvals[j] = cvs[i]->GetValue();
183  }
184 
185  // Reduce across all processors and add hills.
186  MPI_Allreduce(MPI_IN_PLACE, cvals.data(), n*nwalkers, MPI_DOUBLE, MPI_SUM, world_);
187 
188  for(int i = 0; i < n*nwalkers; i += n)
189  {
190  std::vector<double> cval(cvals.begin() + i, cvals.begin() + i + n);
191  hills_.emplace_back(cval, widths_, height_);
192 
193  // Write hill to file.
194  if(world_.rank() == 0)
195  PrintHill(hills_.back(), iteration);
196  }
197 
198  // If grid is defined, add bias onto grid.
199  if(grid_ != nullptr)
200  {
201  std::vector<double> dx(n, 0.0), df(n, 1.0);
202  auto& hill = hills_.back();
203  for(auto it = grid_->begin(); it != grid_->end(); ++it)
204  {
205  auto& val = *it;
206  auto coord = it.coordinates();
207 
208  // Compute difference between grid point and current val.
209  for(size_t i = 0; i < n; ++i)
210  {
211  dx[i] = -cvs[i]->GetDifference(coord[i]);
212  df[i] = 1.;
213  }
214 
215  // Compute derivative.
216  for(size_t i = 0; i < n; ++i)
217  {
218  for(size_t j = 0; j < n; ++j)
219  {
220  if(j != i)
221  df[i] *= gaussian(dx[j], hill.width[j]);
222  else
223  df[i] *= gaussianDerv(dx[j], hill.width[j]);
224  }
225  }
226 
227  // Add to grid.
228  for(size_t i = 0; i < n; ++i)
229  val[i] += height_*df[i];
230  }
231  }
232  }
233 
234  // Writes hill to output file. This should only be called by the
235  // world master node.
236  void Meta::PrintHill(const Hill& hill, int iteration)
237  {
238  hillsout_.open("hills.out", std::fstream::app);
239 
240  hillsout_ << iteration << " ";
241  hillsout_.precision(8);
242 
243  for(auto& cv : hill.center)
244  hillsout_ << cv << " ";
245 
246  for(auto& w : hill.width)
247  hillsout_ << w << " ";
248 
249  hillsout_ << height_ << std::endl;
250  hillsout_.close();
251  }
252 
253  void Meta::CalcBiasForce(const CVList& cvs)
254  {
255  // Reset bias and derivatives.
256  double bias = 0.;
257  auto n = cvs.size();
258 
259  // Reset vectors.
260  std::fill(derivatives_.begin(), derivatives_.end(), 0);
261 
262  // Look up and apply grid bias.
263  if(grid_ != nullptr)
264  {
265  bool inbounds = true;
266  std::vector<double> val(n, 0.);
267  for(size_t i = 0; i < n; ++i)
268  {
269  val[i] = cvs[i]->GetValue();
270  if(val[i] < grid_->GetLower(i) || val[i] > grid_->GetUpper(i))
271  inbounds = false;
272  }
273 
274  if(inbounds)
275  {
276  auto frc = (*grid_)[val];
277  for(size_t i = 0; i < n; ++i)
278  derivatives_[i] = frc[i];
279  }
280  else
281  {
282  if(comm_.rank() == 0)
283  {
284  std::cerr << "Metadynamics: out of bounds ( ";
285  for(auto& v : val)
286  std::cerr << v << " ";
287  std::cerr << ")" << std::endl;
288  }
289  }
290  }
291  else
292  {
293  // Loop through hills and calculate the bias force.
294  for(auto& hill : hills_)
295  {
296  auto tbias = 1.;
297  std::fill(tder_.begin(), tder_.end(), 1.0);
298  std::fill(dx_.begin(), dx_.end(), 1.0);
299 
300  for(size_t i = 0; i < n; ++i)
301  {
302  dx_[i] = cvs[i]->GetDifference(hill.center[i]);
303  tbias *= gaussian(dx_[i], hill.width[i]);
304  }
305 
306  for(size_t i = 0; i < n; ++i)
307  for(size_t j = 0; j < n; ++j)
308  {
309  if(j != i)
310  tder_[i] *= gaussian(dx_[j], hill.width[j]);
311  else
312  tder_[i] *= gaussianDerv(dx_[j], hill.width[j]);
313  }
314 
315  bias += height_ * tbias;
316  for(size_t i = 0; i < n; ++i)
317  derivatives_[i] += height_*tder_[i];
318  }
319  }
320 
321  // Restraints.
322  for(size_t i = 0; i < n; ++i)
323  {
324  auto cval = cvs[i]->GetValue();
325  if(cval < lowerb_[i])
326  derivatives_[i] += lowerk_[i]*(cval - lowerb_[i]);
327  else if(cval > upperb_[i])
328  derivatives_[i] += upperk_[i]*(cval - upperb_[i]);
329  }
330  }
331 
332  Meta* Meta::Build(const Json::Value& json,
333  const MPI_Comm& world,
334  const MPI_Comm& comm,
335  const std::string& path)
336  {
337  ObjectRequirement validator;
338  Value schema;
339  Reader reader;
340 
341  reader.parse(JsonSchema::MetadynamicsMethod, schema);
342  validator.Parse(schema, path);
343 
344  // Validate inputs.
345  validator.Validate(json, path);
346  if(validator.HasErrors())
347  throw BuildException(validator.GetErrors());
348 
349  std::vector<double> widths;
350  for(auto& s : json["widths"])
351  widths.push_back(s.asDouble());
352 
353  std::vector<double> lowerb, upperb, lowerk, upperk;
354  Grid<Vector>* grid = nullptr;
355  if(json.isMember("grid"))
356  grid = Grid<Vector>::BuildGrid(json.get("grid", Json::Value()));
357  else if(!json.isMember("lower_bounds") || !json.isMember("upper_bounds"))
358  throw BuildException({
359  "#/Method/Metadynamics: Both upper_bounds and lower_bounds "
360  "must be defined if grid is not being used."});
361 
362  // Assume all vectors are the same size.
363  for(int i = 0; i < json["lower_bound_restraints"].size(); ++i)
364  {
365  lowerk.push_back(json["lower_bound_restraints"][i].asDouble());
366  upperk.push_back(json["upper_bound_restraints"][i].asDouble());
367  lowerb.push_back(json["lower_bounds"][i].asDouble());
368  upperb.push_back(json["upper_bounds"][i].asDouble());
369  }
370 
371  auto height = json.get("height", 1.0).asDouble();
372  auto hillfreq = json.get("hill_frequency", 1).asInt();
373  auto freq = json.get("frequency", 1).asInt();
374 
375  auto* m = new Meta(
376  world, comm, height, widths,
377  lowerb, upperb, lowerk, upperk,
378  grid, hillfreq, freq
379  );
380 
381  if(json.isMember("load_hills"))
382  m->LoadHills(json["load_hills"].asString());
383 
384  return m;
385  }
386 }
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
Collective variable manager.
Definition: CVManager.h:40
int GetIteration() const
Get the current iteration.
Definition: Snapshot.h:103
const Matrix3 & GetVirial() const
Get box virial.
Definition: Snapshot.h:133
std::vector< CollectiveVariable * > CVList
List of Collective Variables.
Definition: types.h:51
Eigen::VectorXd Vector
Arbitrary length vector.
Definition: types.h:30
double gaussian(double dx, double sigma)
Evlauate Gaussian. Helper function.
Definition: Meta.cpp:42
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:43
double gaussianDerv(double dx, double sigma)
Evaluate Gaussian derivative. Helper function.
Definition: Meta.cpp:54
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
Exception to be thrown when building the Driver fails.
int GetTargetIterations() const
Get target iterations.
Definition: Snapshot.h:109
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
std::vector< CollectiveVariable * > GetCVs(const std::vector< uint > &mask=std::vector< uint >()) const
Get CV iterator.
Definition: CVManager.h:80
"Vanilla" multi-dimensional Metadynamics.
Definition: Meta.h:68
Multidimensional hill.
Definition: Meta.h:37
Requirements on an object.
std::vector< double > center
Hill center.
Definition: Meta.h:40
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:362
std::vector< double > width
Hill width.
Definition: Meta.h:43
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.