SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
ABF.cpp
1 
21 #include "ABF.h"
22 #include "CVs/CVManager.h"
23 #include "Drivers/DriverException.h"
24 #include "Validator/ObjectRequirement.h"
25 #include "Snapshot.h"
26 #include "schema.h"
27 #include <algorithm>
28 #include <fstream>
29 #include <iostream>
30 
31 using namespace Json;
32 // "Adaptive biasing force method for scalar and vector free energy calculations"
33 // Darve, Rodriguez-Gomez, Pohorille
34 // J. Chem. Phys. (2008)
35 namespace SSAGES
36 {
37 
38  // Pre-simulation hook.
39  void ABF::PreSimulation(Snapshot* snapshot, const CVManager& cvmanager)
40  {
41  // Open/close outfile to create it fresh.
42  if(world_.rank() == 0)
43  {
44  worldout_.open(filename_);
45  worldout_.close();
46  }
47 
48  // Convenience. Number of CVs.
49  auto cvs = cvmanager.GetCVs(cvmask_);
50  dim_ = cvs.size();
51 
52  // Size and initialize Fold_
53  Fold_.setZero(dim_);
54 
55  // Initialize biases.
56  biases_.resize(snapshot->GetPositions().size(), Vector3{0, 0, 0});
57 
58  // Initialize w \dot p's for finite difference.
59  wdotp1_.setZero(dim_);
60  wdotp2_.setZero(dim_);
61  }
62 
64 
72  void ABF::PostIntegration(Snapshot* snapshot, const CVManager& cvmanager)
73  {
74 
75  ++iteration_;
76 
77  // Gather information.
78  auto cvs = cvmanager.GetCVs(cvmask_);
79  auto& vels = snapshot->GetVelocities();
80  auto& mass = snapshot->GetMasses();
81  auto& forces = snapshot->GetForces();
82  auto n = snapshot->GetNumAtoms();
83 
84 
86  //int coord = histCoords(cvs);
87 
89  Eigen::MatrixXd J(dim_, 3*n);
90 
92  std::vector<double> cvVals(dim_);
93 
94  // Fill J and CV. Each column represents grad(CV) with flattened Cartesian elements.
95  for(size_t i = 0; i < dim_; ++i)
96  {
97  auto& grad = cvs[i]->GetGradient();
98  for(size_t j = 0; j < n; ++j)
99  J.block<1, 3>(i,3*j) = grad[j];
100  cvVals[i] = cvs[i]->GetValue();
101  }
102 
103  //* Calculate W using Darve's approach (http://mc.stanford.edu/cgi-bin/images/0/06/Darve_2008.pdf).
104  Eigen::MatrixXd Jmass = J.transpose();
105  if(massweigh_)
106  {
107  for(size_t i = 0; i < forces.size(); ++i)
108  Jmass.block(3*i, 0, 3, dim_) = Jmass.block(3*i, 0, 3, dim_)/mass[i];
109  }
110 
111  Eigen::MatrixXd Minv = J*Jmass;
112  MPI_Allreduce(MPI_IN_PLACE, Minv.data(), Minv.size(), MPI_DOUBLE, MPI_SUM, comm_);
113  Eigen::MatrixXd Wt = Minv.inverse()*Jmass.transpose();
114  // Fill momenta.
115  Eigen::VectorXd momenta(3*vels.size());
116  for(size_t i = 0; i < vels.size(); ++i)
117  momenta.segment<3>(3*i) = mass[i]*vels[i];
118 
119  // Compute dot(w,p)
120  Eigen::VectorXd wdotp = Wt*momenta;
121 
122 
123  // Reduce dot product across processors.
124  MPI_Allreduce(MPI_IN_PLACE, wdotp.data(), wdotp.size(), MPI_DOUBLE, MPI_SUM, comm_);
125 
126 
127  // Compute d(wdotp)/dt second order backwards finite difference.
128  // Adding old force removes bias.
129  Eigen::VectorXd dwdotpdt = unitconv_*(1.5*wdotp - 2.0*wdotp1_ + 0.5*wdotp2_)/timestep_ + Fold_;
130 
131  // If we are in bounds, sum force into running total.
132  if(boundsCheck(cvVals) && comm_.rank() == 0)
133  {
134  for(size_t i = 0; i < dim_; ++i)
135  F_[i]->at(cvVals) += dwdotpdt[i];
136  N_->at(cvVals)++;
137  }
138 
139  // Reduce data across processors.
140  for(size_t i = 0; i < dim_; ++i)
141  MPI_Allreduce(F_[i]->data(), Fworld_[i]->data(), (F_[i]->size()), MPI_DOUBLE, MPI_SUM, world_);
142  MPI_Allreduce(N_->data(), Nworld_->data(), N_->size(), MPI_INT, MPI_SUM, world_);
143 
144  // If we are in bounds, store the old summed force.
145  if(boundsCheck(cvVals))
146  {
147  for(size_t i=0; i < dim_; ++i)
148  Fold_[i] = Fworld_[i]->at(cvVals)/std::max(min_, Nworld_->at(cvVals));
149  }
150 
151  // Update finite difference time derivatives.
152  wdotp2_ = wdotp1_;
153  wdotp1_ = wdotp;
154 
155  // Write out data to file.
156  if(iteration_ % FBackupInterv_ == 0){
157  WriteData();
158  }
159 
160  // Calculate the bias from averaged F at current CV coordinates.
161  // Or apply harmonic restraints to return CVs back in bounds.
162  CalcBiasForce(snapshot, cvs, cvVals);
163 
164  // Update the forces in snapshot by adding in the force bias from each
165  // CV to each atom based on the gradient of the CV.
166  for (size_t j = 0; j < forces.size(); ++j)
167  forces[j] += biases_[j];
168  }
169 
170  // Post-simulation hook.
171  void ABF::PostSimulation(Snapshot*, const CVManager&)
172  {
173  WriteData();
174  }
175 
176  void ABF::CalcBiasForce(const Snapshot* snapshot, const CVList& cvs, const std::vector<double> &cvVals)
177  {
178  // Reset the bias.
179  biases_.resize(snapshot->GetNumAtoms(), Vector3{0,0,0});
180  for(auto& b : biases_)
181  b.setZero();
182 
183  // Compute bias if within bounds.
184  if(boundsCheck(cvVals))
185  {
186  for(size_t i = 0; i < dim_; ++i)
187  {
188  auto& grad = cvs[i]->GetGradient();
189  for(size_t j = 0; j < biases_.size(); ++j)
190  biases_[j] -= (Fworld_[i]->at(cvVals))*grad[j]/std::max(min_,Nworld_->at(cvVals));
191  }
192  }
193  // Otherwise, apply harmonic restraint if enabled.
194  else
195  {
196  for(size_t i = 0; i < dim_; ++i)
197  {
198  double cvVal = cvs[i]->GetValue();
199  auto k = 0.;
200  auto x0 = 0.;
201 
202  if(isperiodic_[i])
203  {
204  double periodsize = periodicboundaries_[i][1]-periodicboundaries_[i][0];
205  double cvRestrMidpoint = (restraint_[i][1]+restraint_[i][0])/2;
206  while((cvVal-cvRestrMidpoint) > periodsize/2)
207  cvVal -= periodsize;
208  while((cvVal-cvRestrMidpoint) < -periodsize/2)
209  cvVal += periodsize;
210  }
211 
212  if(cvVal < restraint_[i][0] && restraint_[i][2] > 0)
213  {
214  k = restraint_[i][2];
215  x0 = restraint_[i][0];
216  }
217  else if (cvVal > restraint_[i][1] && restraint_[i][2] > 0)
218  {
219  k = restraint_[i][2];
220  x0 = restraint_[i][1];
221  }
222 
223  auto& grad = cvs[i]->GetGradient();
224  for(size_t j = 0; j < biases_.size(); ++j)
225  biases_[j] -= (cvVal - x0)*k*grad[j];
226  }
227  }
228  }
229 
230  // Write out the average generalized force.
231  // Also write out Fworld and Nworld backups for restarts.
232  void ABF::WriteData()
233  {
234 
235  // Only one processor should be performing I/O.
236  if(world_.rank() != 0)
237  return;
238 
239  Nworld_->syncGrid();
240 
241 
242  // Backup Fworld and Nworld.
243  Nworld_->WriteToFile(Nworld_filename_);
244  for(size_t i = 0 ; i < dim_; ++i)
245  {
246  Fworld_[i]->syncGrid();
247  Fworld_[i]->WriteToFile(Fworld_filename_+std::to_string(i));
248  }
249 
250 
251  // Average the generalized force for each bin and print it out.
252  int gridPoints = 1;
253  for(size_t i = 0 ; i < dim_; ++i)
254  gridPoints = gridPoints * N_->GetNumPoints(i);
255 
256  worldout_.open(filename_, std::ios::out);
257  worldout_ << std::endl;
258  worldout_ << "Iteration: " << iteration_ << std::endl;
259  worldout_ << "Printing out the current Adaptive Biasing Vector Field." << std::endl;
260  worldout_ << "First (Nr of CVs) columns are the coordinates, the next (Nr of CVs) columns are components of the Adaptive Force vector at that point." << std::endl;
261  worldout_ << "The columns are " << gridPoints << " long, mapping out a surface of ";
262 
263  for(size_t i = 0 ; i < dim_-1; ++i)
264  worldout_ << (N_->GetNumPoints())[i] << " by " ;
265  worldout_ << (N_->GetNumPoints(dim_-1)) << " points in " << dim_ << " dimensions." << std::endl;
266  worldout_ << std::endl;
267 
268 
269  for(auto it = Nworld_->begin(); it != Nworld_->end(); ++it)
270  {
271  auto& val = *it;
272  auto coord = it.coordinates();
273  for(auto& c : coord)
274  worldout_ << std::setw(14) << std::setprecision(8) << std::fixed << c << " ";
275  for(size_t i = 0 ; i < dim_; ++i)
276  worldout_ << std::setw(14) << std::setprecision(8) << std::fixed << Fworld_[i]->at(coord)/std::max(val,min_) << " ";
277  worldout_ << std::endl;
278  }
279 
280  worldout_ << std::endl;
281  worldout_ << std::endl;
282  worldout_.close();
283  }
284 
285  // Checks whether walker is within CV bounds.
286  bool ABF::boundsCheck(const std::vector<double> &CVs)
287  {
288  for(size_t i = 0; i < dim_ ; ++i)
289  {
290  if((CVs[i] < N_->GetLower(i)) || (CVs[i] > N_->GetUpper(i)))
291  {
292  return false;
293  }
294  }
295  return true;
296  }
297 
298 
299  ABF* ABF::Build(const Value& json,
300  const MPI_Comm& world,
301  const MPI_Comm& comm,
302  const std::string& path)
303  {
304  ObjectRequirement validator;
305  Value schema;
306  CharReaderBuilder rbuilder;
307  CharReader* reader = rbuilder.newCharReader();
308 
309  reader->parse(JsonSchema::ABFMethod.c_str(),
310  JsonSchema::ABFMethod.c_str() + JsonSchema::ABFMethod.size(),
311  &schema, NULL);
312  validator.Parse(schema, path);
313 
314  // Validate inputs.
315  validator.Validate(json, path);
316  if(validator.HasErrors())
317  throw BuildException(validator.GetErrors());
318 
319  unsigned int wid = mxx::comm(world).rank()/mxx::comm(comm).size(); // walker ID
320  unsigned int wrank = mxx::comm(world).rank()%mxx::comm(comm).size(); // walker rank
321 
322  //bool multiwalkerinput = false;
323 
324  // Obtain lower bound for each CV.
325  std::vector<double> minsCV;
326  for(auto& mins : json["CV_lower_bounds"])
327  if(mins.isArray())
328  {
329  throw BuildException({"Separate inputs for multi-walker not fully implemented. Please use global entries for CV ranges"});
330  //multiwalkerinput = true;
331  //for(auto& bound : mins)
332  // minsCV.push_back(bound.asDouble());
333  }
334  else
335  minsCV.push_back(mins.asDouble());
336 
337  // Obtain upper bound for each CV.
338  std::vector<double> maxsCV;
339  for(auto& maxs : json["CV_upper_bounds"])
340  /*if(maxs.isArray())
341  {
342  for(auto& bound : maxs)
343  maxsCV.push_back(bound.asDouble());
344  }*/
345  //else
346  maxsCV.push_back(maxs.asDouble());
347 
348  // Obtain number of bins for each CV dimension.
349  std::vector<int> binsCV;
350  for(auto& bins : json["CV_bins"])
351  binsCV.push_back(bins.asInt());
352 
353  // Obtain lower bounds for restraints for each CV.
354  std::vector<double> minsrestCV;
355  for(auto& mins : json["CV_restraint_minimums"])
356  /*if(mins.isArray())
357  {
358  for(auto& bound : mins)
359  minsrestCV.push_back(bound.asDouble());
360  }*/
361  //else
362  minsrestCV.push_back(mins.asDouble());
363 
364  // Obtain upper bounds for restraints for each CV.
365  std::vector<double> maxsrestCV;
366  for(auto& maxs : json["CV_restraint_maximums"])
367  /*if(maxs.isArray())
368  {
369  for(auto& bound : maxs)
370  maxsrestCV.push_back(bound.asDouble());
371  }*/
372  //else
373  maxsrestCV.push_back(maxs.asDouble());
374 
375  // Obtain harmonic spring constant for restraints for each CV.
376  std::vector<double> springkrestCV;
377  for(auto& bins : json["CV_restraint_spring_constants"])
378  springkrestCV.push_back(bins.asDouble());
379 
380  // Obtain periodicity information for restraints for each CV for the purpose of correctly applying restraints through periodic boundaries.
381  std::vector<bool> isperiodic;
382  for(auto& isperCV : json["CV_isperiodic"])
383  isperiodic.push_back(isperCV.asBool());
384 
385  // Obtain lower periodic boundary for each CV.
386  std::vector<double> minperboundaryCV;
387  for(auto& minsperCV : json["CV_periodic_boundary_lower_bounds"])
388  minperboundaryCV.push_back(minsperCV.asDouble());
389 
390  // Obtain upper periodic boundary for each CV.
391  std::vector<double> maxperboundaryCV;
392  for(auto& maxsperCV : json["CV_periodic_boundary_upper_bounds"])
393  maxperboundaryCV.push_back(maxsperCV.asDouble());
394 
395  // Verify inputs are all correct lengths.
396  if(!(( minsCV.size() == maxsCV.size() &&
397  maxsCV.size() == minsrestCV.size() &&
398  minsrestCV.size() == maxsrestCV.size()) &&
399  (binsCV.size() == springkrestCV.size() &&
400  springkrestCV.size() == isperiodic.size())))
401  throw BuildException({"CV lower bounds, upper bounds, restrain minimums, restrains maximums must match in size. Bins, spring constants and periodicity info must match in size."});
402 
403  // Verify that all periodicity information is provided if CVs are periodic.
404  bool anyperiodic=false;
405  for(size_t i = 0; i<isperiodic.size(); ++i)
406  if(isperiodic[i])
407  {
408  anyperiodic = true;
409  if(!(isperiodic.size() == minperboundaryCV.size() &&
410  minperboundaryCV.size() == maxperboundaryCV.size()))
411  throw BuildException({"If any CV is defined as periodic, please define the full upper and lower bound vectors. They should both have the same number of entries as CV lower bounds, upper bounds... Entries corresponding to non-periodic CVs will not be used."});
412  }
413 
414  size_t dim = binsCV.size();
415 
416  // Read in JSON info.
417  auto FBackupInterv = json.get("output_frequency", 1000).asInt();
418  auto unitconv = json.get("unit_conversion", 0).asDouble();
419  auto timestep = json.get("timestep", 2).asDouble();
420  auto min = json.get("minimum_count", 200).asDouble();
421  auto massweigh = json.get("mass_weighing",false).asBool();
422 
423  std::vector<std::vector<double>> histdetails;
424  std::vector<std::vector<double>> restraint;
425  std::vector<std::vector<double>> periodicboundaries;
426  std::vector<double> temp1(3);
427  std::vector<double> temp2(3);
428  std::vector<double> temp3(2);
429 
430  auto freq = json.get("frequency", 1).asInt();
431  auto filename = json.get("output_file", "F_out").asString();
432  auto Nworld_filename = json.get("Nworld_output_file", "Nworld").asString();
433  auto Fworld_filename = json.get("Fworld_output_file", "Fworld_cv").asString();
434 
435  // Generate the grids based on JSON.
436  Grid<int> *N;
437  Grid<int> *Nworld;
438  std::vector<Grid<double>*> F(dim);
439  std::vector<Grid<double>*> Fworld(dim);
440 
441  // This feature is disabled for now.
442  /*if(multiwalkerinput)
443  {
444  for(size_t i=0; i<dim; ++i)
445  {
446  temp1 = {minsCV[i+wid*dim], maxsCV[i+wid*dim], binsCV[i]};
447  temp2 = {minsrestCV[i+wid*dim], maxsrestCV[i+wid*dim], springkrestCV[i]};
448  histdetails.push_back(temp1);
449  restraint.push_back(temp2);
450  if(anyperiodic)
451  {
452  temp3 = {minperboundaryCV[i], maxperboundaryCV[i]};
453  periodicboundaries.push_back(temp3);
454  }
455  }
456  for(size_t i=0; i<dim; ++i)
457  {
458  minsCVperwalker[i] = histdetails[i][0];
459  maxsCVperwalker[i] = histdetails[i][1];
460  }
461 
462  N= new Grid<int>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
463  Nworld= new Grid<int>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
464 
465  for(auto& grid : F)
466  {
467  grid= new Grid<double>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
468  }
469  for(auto& grid : Fworld)
470  {
471  grid= new Grid<double>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
472  }
473 
474  }
475  */
476  //else
477  //{
478 
479  // Appropriately size the grids.
480 
481  N = new Grid<int>(binsCV, minsCV, maxsCV, isperiodic);
482  Nworld = new Grid<int>(binsCV, minsCV, maxsCV, isperiodic);
483 
484  for(auto& grid : F)
485  {
486  grid = new Grid<double>(binsCV, minsCV, maxsCV, isperiodic);
487  }
488  for(auto& grid : Fworld)
489  {
490  grid = new Grid<double>(binsCV, minsCV, maxsCV, isperiodic);
491  }
492 
493  for(size_t i = 0; i < dim; ++i)
494  {
495  temp1 = {minsCV[i], maxsCV[i], double(binsCV[i])};
496  temp2 = {minsrestCV[i], maxsrestCV[i], springkrestCV[i]};
497  histdetails.push_back(temp1);
498  restraint.push_back(temp2);
499  if(anyperiodic)
500  {
501  temp3 = {minperboundaryCV[i], maxperboundaryCV[i]};
502  periodicboundaries.push_back(temp3);
503  }
504  }
505  //}
506 
507  // Check if a restart is requested.
508  bool restart = json.get("restart",false).asBool();
509 
510  // Either load previous grid, or back it up.
511  if (wrank == 0 && restart)
512  {
513 
514  std::cout << "Attempting to load data from a previous run of ABF..." << std::endl;
515  N->LoadFromFile(Nworld_filename);
516  for(size_t i = 0; i < dim; ++i)
517  {
518  F[i]->LoadFromFile(Fworld_filename+std::to_string(i));
519  }
520 
521  }
522  else if (wrank == 0)
523  {
524  std::ifstream Nworldbackupsource(Nworld_filename, std::ios::binary);
525  if(Nworldbackupsource)
526  {
527  std::cout << "Backing up previous copy of Nworld." << std::endl;
528  std::ofstream Nworldbackuptarget(Nworld_filename+"_backup", std::ios::binary);
529  Nworldbackuptarget << Nworldbackupsource.rdbuf();
530  }
531  for(size_t i = 0; i < dim; ++i)
532  {
533  std::ifstream Fworldbackupsource(Fworld_filename+std::to_string(i), std::ios::binary);
534  if(Fworldbackupsource)
535  {
536  std::cout << "Backing up previous copy of Fworld"+std::to_string(i)+"." << std::endl;
537  std::ofstream Fworldbackuptarget(Fworld_filename+std::to_string(i)+"_backup", std::ios::binary);
538  Fworldbackuptarget << Fworldbackupsource.rdbuf();
539  }
540  }
541  }
542 
543 
544  auto* m = new ABF(world, comm, N, Nworld, F, Fworld, restraint, isperiodic, periodicboundaries, min, massweigh, filename, Nworld_filename, Fworld_filename, histdetails, FBackupInterv, unitconv, timestep, freq);
545 
546 
547  if(json.isMember("iteration"))
548  m->SetIteration(json.get("iteration",0).asInt());
549 
550  return m;
551 
552  }
553 
554 }
555 
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
unsigned GetNumAtoms() const
Get number of atoms in this snapshot.
Definition: Snapshot.h:200
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
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.
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Adaptive Biasing Force Algorithm.
Definition: ABF.h:40
Requirements on an object.
const std::vector< Vector3 > & GetVelocities() const
Access the particle velocities.
Definition: Snapshot.h:349
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
void LoadFromFile(const std::string &filename)
Builds a grid from file.
Definition: Grid.h:629
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:323
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.