SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
BasisFunc.cpp
1 
21 #include "BasisFunc.h"
22 #include "Validator/ObjectRequirement.h"
23 #include "Drivers/DriverException.h"
24 #include "CVs/CVManager.h"
25 #include "Snapshot.h"
26 #include <cmath>
27 #include <iostream>
28 #include <iomanip>
29 
30 using namespace Json;
31 
32 namespace SSAGES
33 {
34 
35  // Pre-simulation hook.
36  void BFS::PreSimulation(Snapshot* snapshot, const CVManager& cvmanager)
37  {
38  auto cvs = cvmanager.GetCVs(cvmask_);
39  // For print statements and file I/O, the walker IDs are used
40  mpiid_ = snapshot->GetWalkerID();
41 
42  // Make sure the iteration index is set correctly
43  iteration_ = 0;
44  step_ = 0;
45 
46  // Set the temporary sweep frequency to the same as the input
47  cyclefreqTmp_ = cyclefreq_;
48 
49  // There are a few error messages / checks that are in place with
50  // defining CVs and grids
51  if(h_->GetDimension() != cvs.size())
52  {
53  std::cerr<<"ERROR: Grid dimensions doesn't match number of CVS."<<std::endl;
54  MPI_Abort(world_, EXIT_FAILURE);
55  }
56 
57  // This is to check for non-periodic bounds. It comes into play in the update bias function
58  bounds_ = true;
59  unbias_.resize(h_->size(),0);
60  std::fill(unbias_.begin(),unbias_.end(),1.0);
61 
62  // Initialize the mapping for the hist function
63  for (auto &val : *h_) {
64  val = 0;
65  }
66 
67  // Reset the values in the force grid
68  for (auto &val : *f_) {
69  for (size_t i = 0; i <val.size(); i++) {
70  val[i] = 0;
71  }
72  }
73 
74  // Reset the values in the bias potential grid
75  for (auto &val : *b_) {
76  val = 0;
77  }
78  }
79 
80  // Post-integration hook.
81  void BFS::PostIntegration(Snapshot* snapshot, const CVManager& cvmanager)
82  {
83  auto cvs = cvmanager.GetCVs(cvmask_);
84  std::vector<double> x(cvs.size(),0);
85 
86  /*The binned cv space is updated at every step
87  *After a certain number of steps has been passed, the system updates a
88  *bias projection based on the visited histogram states
89  */
90  for(size_t i = 0; i < cvs.size(); ++i)
91  {
92  x[i] = cvs[i]->GetValue();
93  }
94 
95  // The histogram is updated based on the index
96  InBounds(cvs);
97  if (bounds_) {h_->at(x) += 1;}
98  else {cyclefreq_++;}
99 
100  // Update the basis projection after the requisite number of steps
101  int steps = snapshot->GetIteration()-step_;
102  if(steps % cyclefreq_ == 0) {
103  step_ = snapshot->GetIteration();
104  cyclefreq_ = cyclefreqTmp_;
105  double beta;
106  beta = 1.0 / (temperature_ * snapshot->GetKb());
107 
108  // For systems with poorly defined temperature (ie: 1 particle) the
109  // user needs to define their own temperature. This is a hack that
110  // will be removed in future versions.
111 
112  iteration_+= 1;
113  ProjectBias(cvs,beta);
114  std::cout<<"Node: ["<<mpiid_<<"]"<<std::setw(10)<<"\tSweep: "<<iteration_<<std::endl;
115  }
116 
117  // This gets the bias force from the grid
118  std::vector<double> bias_grad(cvs.size(),0);
119 
120  if (bounds_)
121  {
122  auto& bias = f_->at(x);
123  for (size_t ii = 0; ii<bias.size(); ii++)
124  bias_grad[ii] = bias[ii];
125  }
126  else
127  {
128  // This is where the wall potentials are going to be thrown into the method if the system is not a periodic CV
129  for(size_t j = 0; j < cvs.size(); ++j)
130  {
131  if(!h_->GetPeriodic(j))
132  {
133  if(x[j] > boundUp_[j])
134  bias_grad[j] = -restraint_[j] * (x[j] - boundUp_[j]);
135  else if(x[j] < boundLow_[j])
136  bias_grad[j] = restraint_[j] * (x[j] - boundLow_[j]);
137  }
138  }
139  }
140 
141  // Take each CV and add its biased forces to the atoms using the chain rule
142  auto& forces = snapshot->GetForces();
143  auto& virial = snapshot->GetVirial();
144  for(size_t i = 0; i < cvs.size(); ++i)
145  {
146  auto& grad = cvs[i]->GetGradient();
147  auto& boxgrad = cvs[i]->GetBoxGradient();
148 
149  /* Update the forces in snapshot by adding in the force bias from each
150  *CV to each atom based on the gradient of the CV.
151  */
152  for (size_t j = 0; j < forces.size(); ++j)
153  forces[j] += bias_grad[i]*grad[j];
154 
155  virial -= bias_grad[i]*boxgrad;
156  }
157  }
158 
159  // Post-simulation hook.
160  void BFS::PostSimulation(Snapshot*, const CVManager&)
161  {
162  std::cout<<"Run has finished"<<std::endl;
163  }
164 
165  // Update the coefficients/bias projection
166  void BFS::ProjectBias(const CVList& cvs, const double beta)
167  {
168  double sum = 0.0;
169 
170  // For multiple walkers, the struct is unpacked
171  Grid<unsigned int> histlocal(*h_);
172 
173  // Summed between all walkers
174  MPI_Allreduce(histlocal.data(), h_->data(), h_->size(), MPI_INT, MPI_SUM, world_);
175 
176  h_->syncGrid();
177 
178  // Construct the biased histogram
179  size_t i = 0;
180  for (Grid<unsigned int>::iterator it2 = h_->begin(); it2 != h_->end(); ++it2, ++i)
181  {
182  /* The evaluation of the biased histogram which projects the histogram to the
183  * current bias of CV space.
184  */
185  unbias_[i] += (*it2) * exp(beta * b_->at(it2.coordinates()));
186  }
187 
188  // Reset histogram
189  for (auto &val : *h_) {
190  val = 0;
191  }
192 
193  // Create the log array and send that to the integrator
194  std::vector<double> z (unbias_.size(),0);
195  for (i = 0; i < unbias_.size(); i++)
196  // Make sure the log of the biased histogram is a number
197  z[i] = 1.0/beta*log(unbias_[i] * weight_);
198 
199  //Update the coefficients and determine the difference from the previous iteration
200  sum = evaluator_.UpdateCoeff(z,h_);
201  coeffArr_ = evaluator_.GetCoeff();
202  //Update both the gradient and the bias on the grids
203  evaluator_.UpdateBias(b_,f_);
204 
205  if(world_.rank() == 0) {
206  // Write coeff at this step, but only one walker
207  std::cout<<"Coefficient difference is: " <<sum<<std::endl;
208  PrintBias(cvs,beta);
209  }
210 
211  // The convergence tolerance and whether the user wants to exit are incorporated here
212  if(sum < tol_)
213  {
214  std::cout<<"System has converged"<<std::endl;
215  if(convergeExit_)
216  {
217  std::cout<<"User has elected to exit. System is now exiting"<<std::endl;
218  exit(EXIT_SUCCESS);
219  }
220  }
221  }
222 
223  /*The coefficients are printed out for the purpose of saving the free energy space
224  *Additionally, the current basis projection is printed so that the user can view
225  *the current free energy space
226  */
227  void BFS::PrintBias(const CVList& cvs, const double beta)
228  {
229  // The filenames will have a standard name, with a user-defined suffix
230  std::string filename1 = "basis"+bnme_+".out";
231  std::string filename2 = "restart"+bnme_+".out";
232  std::ofstream basisout;
233  std::ofstream coeffout;
234  basisout.precision(5);
235  coeffout.precision(5);
236  basisout.open(filename1.c_str());
237  coeffout.open(filename2.c_str());
238 
239  // The CV values, PMF projection, PMF, and biased histogram are output for the user
240  basisout << "The information stored in this file is the output of a Basis Function Sampling run" << std::endl;
241  basisout << "CV Values" << std::setw(15*cvs.size()) << "Basis Set Bias" << std::setw(15) << "PMF Estimate" << std::endl;
242  coeffout << "The information stored in this file is for the purpose of restarting simulations with BFS" << std::endl;
243  coeffout << "***COEFFICIENTS***" << std::endl;
244 
245  size_t j = 0;
246  for(Grid<double>::iterator it = b_->begin(); it != b_->end(); ++it, ++j)
247  {
248  for(size_t k = 0; k < cvs.size(); ++k)
249  {
250  // Evaluate the CV values for printing purposes
251  basisout << it.coordinate(k) << std::setw(15);
252  }
253  basisout << -(*it) << std::setw(15) <<
254  -1.0/beta*log(unbias_[j]) <<
255  std::endl;
256  }
257 
258  std::vector<double> coeff = evaluator_.GetCoeff();
259  for (auto& val : coeff) {
260  coeffout << val << std::endl;
261  }
262  coeffout << "***HISTOGRAM***" << std::endl;
263  for (auto& val : unbias_) {
264  coeffout << val << std::endl;
265  }
266 
267  basisout.close();
268  coeffout.close();
269  }
270 
271  // The forces are calculated by chain rule, first the derivatives of the basis set, then in the PostIntegration function, the derivative of the CV is evaluated
272  void BFS::InBounds(const CVList& cvs)
273  {
274  std::vector<double> x (cvs.size(),0);
275  //This is calculating the derivatives for the bias force
276  for (size_t j = 0; j < cvs.size(); ++j)
277  {
278  x[j] = cvs[j]->GetValue();
279  double min = h_->GetLower(j);
280  double max = h_->GetUpper(j);
281 
282  if(!h_->GetPeriodic(j))
283  {
284  // In order to prevent the index for the histogram from going out of bounds a check is in place
285  if(x[j] > max && bounds_)
286  {
287  //std::cout<<"WARNING: CV is above the maximum boundary."<<std::endl;
288  //std::cout<<"Statistics will not be gathered during this interval"<<std::endl;
289  bounds_ = false;
290  break;
291  }
292  else if(x[j] < min && bounds_)
293  {
294  //std::cout<<"WARNING: CV is below the minimum boundary."<<std::endl;
295  //std::cout<<"Statistics will not be gathered during this interval"<<std::endl;
296  bounds_ = false;
297  break;
298  }
299  else if(x[j] < max && x[j] > min && !bounds_)
300  {
301  //std::cout<<"CV has returned in between bounds. Run is resuming"<<std::endl;
302  bounds_ = true;
303  }
304  }
305  }
306  }
307 
308  BFS* BFS::Build(const Json::Value& json,
309  const MPI_Comm& world,
310  const MPI_Comm& comm,
311  const std::string& path)
312  {
313  ObjectRequirement validator;
314  Value schema;
315  CharReaderBuilder rbuilder;
316  CharReader* reader = rbuilder.newCharReader();
317 
318  reader->parse(JsonSchema::BFSMethod.c_str(),
319  JsonSchema::BFSMethod.c_str() + JsonSchema::BFSMethod.size(),
320  &schema, NULL);
321  validator.Parse(schema, path);
322 
323  //Validate Inputs
324  validator.Validate(json, path);
325  if(validator.HasErrors())
326  throw BuildException(validator.GetErrors());
327 
328  std::vector<double> restrCV(0);
329  for(auto& restr : json["CV_restraint_spring_constants"])
330  restrCV.push_back(restr.asDouble());
331 
332  std::vector<double> boundLow(0);
333  for(auto& bndl : json["CV_restraint_minimums"])
334  boundLow.push_back(bndl.asDouble());
335 
336  std::vector<double> boundUp(0);
337  for(auto& bndu : json["CV_restraint_maximums"])
338  boundUp.push_back(bndu.asDouble());
339 
340  auto cyclefreq = json.get("cycle_frequency", 100000).asInt();
341  auto freq = json.get("frequency", 1).asInt();
342  auto wght = json.get("weight", 1.0).asDouble();
343  auto bnme = json.get("basis_filename", "").asString();
344  auto temp = json.get("temperature", 0.0).asDouble();
345  auto tol = json.get("tolerance", 1e-6).asDouble();
346  auto conv = json.get("convergence_exit", false).asBool();
347 
349  json.get("grid", Json::Value()) );
350 
352  json.get("grid", Json::Value()) );
353 
355  json.get("grid", Json::Value()) );
356 
357  size_t ii = 0;
358  std::vector<BasisFunction*> functions;
359  for(auto& m : json["basis_functions"])
360  {
361  auto *bf = BasisFunction::Build(m, path, b->GetNumPoints(ii));
362  functions.push_back(bf);
363  ii++;
364  }
365 
366  // If restraints were not specified then set them equal to zero
367  if (restrCV.size() == 0) {
368  // Choose simplest array that should have equivalent size
369  for(size_t i = 0; i<functions.size(); i++)
370  restrCV.push_back(0); // Push back value of 0
371  }
372 
373  // This also needs to be the same for upper bounds and lower bounds
374  if (boundLow.size() == 0) {
375  // Choose simplest array that should have equivalent size
376  for(size_t i = 0; i<functions.size(); i++)
377  boundLow.push_back(0); // Push back value of 0
378  }
379 
380  // This also needs to be the same for upper bounds and lower bounds
381  if (boundUp.size() == 0) {
382  // Choose simplest array that should have equivalent size
383  for(size_t i = 0; i<functions.size(); i++)
384  boundUp.push_back(0); // Push back value of 0
385  }
386 
387  // Check to make sure that the arrays for boundaries and restraints are of the same length
388  if (boundUp.size() != boundLow.size() && boundLow.size() != restrCV.size()) {
389  std::cerr<<"ERROR: Number of restraint boundaries do not match number of spring constants."<<std::endl;
390  std::cerr<<"Boundary size is " << boundUp.size() << " boundary low size is " << boundLow.size() << " restraints size is" << restrCV.size() << std::endl;
391  MPI_Abort(world, EXIT_FAILURE);
392  }
393 
394  auto* m = new BFS(world, comm, h, f, b, functions, restrCV, boundUp, boundLow,
395  cyclefreq, freq, bnme, temp, tol, wght,
396  conv);
397 
398  if(json.isMember("iteration"))
399  m->SetIteration(json.get("iteration",0).asInt());
400 
401  // Check if previously saved grids exist. If so, check that data match and load grids.
402  std::ifstream restrFile;
403  restrFile.open("restart"+bnme+".out");
404  if(restrFile.is_open())
405  {
406  std::string line;
407  std::vector<double> coeff;
408  std::vector<double> unbias;
409  bool coeffFlag = false;
410  bool basisFlag = false;
411  std::cout << "Attempting to load data from a previous run of BFS." << std::endl;
412 
413  while(!std::getline(restrFile,line).eof())
414  {
415  if(line.find("COEFFICIENTS") != std::string::npos) {
416  std::getline(restrFile,line);
417  coeffFlag = true;
418  }
419  else if (line.find("HISTOGRAM") != std::string::npos) {
420  std::getline(restrFile,line);
421  basisFlag = true;
422  coeffFlag = false;
423  }
424  if(coeffFlag) {
425  coeff.push_back(std::stod(line));
426  }
427  else if (basisFlag) {
428  unbias.push_back(std::stod(line));
429  }
430  }
431  m->SetBasis(coeff, unbias);
432  restrFile.close();
433  }
434 
435  return m;
436  }
437 
438 }
std::vector< CollectiveVariable * > GetCVs(const std::vector< unsigned int > &mask=std::vector< unsigned int >()) const
Get CV iterator.
Definition: CVManager.h:80
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
Collective variable manager.
Definition: CVManager.h:40
std::vector< CollectiveVariable * > CVList
List of Collective Variables.
Definition: types.h:51
const std::vector< Vector3 > & GetForces() const
Access the per-particle forces.
Definition: Snapshot.h:362
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:43
double GetKb() const
Get system Kb.
Definition: Snapshot.h:163
size_t GetIteration() const
Get the current iteration.
Definition: Snapshot.h:103
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
iterator begin()
Return iterator at first grid point.
Definition: Grid.h:527
Exception to be thrown when building the Driver fails.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
const std::vector< int > GetNumPoints() const
Get the number of points for all dimensions.
Definition: GridBase.h:200
Requirements on an object.
const Matrix3 & GetVirial() const
Get box virial.
Definition: Snapshot.h:133
Basis Function Sampling Algorithm.
Definition: BasisFunc.h:38
unsigned GetWalkerID() const
Get walker ID.
Definition: Snapshot.h:193
T * data()
Get pointer to the internal data storage vector.
Definition: GridBase.h:333
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.