SSAGES  0.1
A MetaDynamics Package
Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
SSAGES::ABF Class Reference

Adaptive Biasing Force Algorithm. More...

#include <ABF.h>

Inheritance diagram for SSAGES::ABF:
Inheritance graph
[legend]

Public Member Functions

 ABF (const MPI_Comm &world, const MPI_Comm &comm, Grid< int > *N, Grid< int > *Nworld, std::vector< Grid< double > * > F, std::vector< Grid< double > * > Fworld, std::vector< std::vector< double >> restraint, std::vector< bool > isperiodic, std::vector< std::vector< double >> periodicboundaries, double min, bool massweigh, std::string filename, std::string Nworld_filename, std::string Fworld_filename, const std::vector< std::vector< double >> &histdetails, int FBackupInterv, double unitconv, double timestep, unsigned int frequency)
 Constructor. More...
 
void PreSimulation (Snapshot *snapshot, const class CVManager &cvmanager) override
 Pre-simulation hook. More...
 
void PostIntegration (Snapshot *snapshot, const class CVManager &cvmanager) override
 Post-integration hook. More...
 
void PostSimulation (Snapshot *snapshot, const class CVManager &cvmanager) override
 Post-simulation hook. More...
 
void SetIteration (const int iter)
 Set iteration of the method. More...
 
 ~ABF ()
 Destructor.
 
- Public Member Functions inherited from SSAGES::Method
 Method (uint frequency, const MPI_Comm &world, const MPI_Comm &comm)
 Constructor. More...
 
void SetCVMask (const std::vector< uint > &mask)
 Sets the collective variable mask.
 
virtual ~Method ()
 Destructor.
 
- Public Member Functions inherited from SSAGES::EventListener
 EventListener (uint frequency)
 Constructor. More...
 
uint GetFrequency () const
 Get frequency of event listener. More...
 
virtual ~EventListener ()
 Destructor.
 

Static Public Member Functions

static ABFBuild (const Json::Value &json, const MPI_Comm &world, const MPI_Comm &comm, const std::string &path)
 
- Static Public Member Functions inherited from SSAGES::Method
static MethodBuildMethod (const Json::Value &json, const MPI_Comm &world, const MPI_Comm &comm, const std::string &path)
 Build a derived method from JSON node. More...
 

Private Member Functions

void CalcBiasForce (const Snapshot *snapshot, const CVList &cvs, const std::vector< double > &cvVals)
 Computes the bias force. More...
 
void WriteData ()
 Writes out data to file.
 
bool boundsCheck (const std::vector< double > &CVs)
 Checks whether the local walker is within CV bounds. More...
 

Private Attributes

std::vector< Grid< double > * > F_
 To store running total of the local walker. More...
 
std::vector< Grid< double > * > Fworld_
 Will hold the global total, synced across walkers at every time step. More...
 
Grid< int > * N_
 To store number of local hits at a given CV bin. More...
 
Grid< int > * Nworld_
 To store number of global hits at a given CV bin. More...
 
std::vector< std::vector< double > > restraint_
 Information for a harmonic restraint to keep CV in the region of interest. More...
 
std::vector< bool > isperiodic_
 For each CV, holds whether that CV has periodic boundaries or not.
 
std::vector< std::vector< double > > periodicboundaries_
 Holds periodic boundaries of CVs. More...
 
int min_
 
Eigen::VectorXd wdotp1_
 To hold last iteration wdotp value for numerical derivative.
 
Eigen::VectorXd wdotp2_
 To hold second to last iteration wdotp value for numerical derivative.
 
Eigen::VectorXd Fold_
 To hold last iterations F_ value for removing bias.
 
bool massweigh_
 Mass weighing of bias enabled/disabled.
 
std::vector< Vector3biases_
 Biases applied to atoms each timestep.
 
unsigned int dim_
 Number of CVs in system.
 
std::ofstream worldout_
 Output stream for F/N world data.
 
std::string filename_
 
std::string Nworld_filename_
 Nworld print out filename, for restarts.
 
std::string Fworld_filename_
 Fworld print out filename, for restarts.
 
std::vector< std::vector< double > > histdetails_
 Histogram details. More...
 
int FBackupInterv_
 Integer to hold F estimate backup information. More...
 
double unitconv_
 Unit Conversion Constant from W dot P to Force. More...
 
double timestep_
 Timestep of integration.
 
uint iteration_
 Iteration counter.
 
Eigen::VectorXd mass_
 Mass vector. Empty unless required.
 

Additional Inherited Members

- Protected Attributes inherited from SSAGES::Method
mxx::comm world_
 Global MPI communicator.
 
mxx::comm comm_
 Local MPI communicator.
 
std::vector< uint > cvmask_
 Mask which identifies which CVs to act on.
 

Detailed Description

Adaptive Biasing Force Algorithm.

Implementation of the Adaptive Biasing Force algorithm based on [2]

Definition at line 42 of file ABF.h.

Constructor & Destructor Documentation

SSAGES::ABF::ABF ( const MPI_Comm &  world,
const MPI_Comm &  comm,
Grid< int > *  N,
Grid< int > *  Nworld,
std::vector< Grid< double > * >  F,
std::vector< Grid< double > * >  Fworld,
std::vector< std::vector< double >>  restraint,
std::vector< bool >  isperiodic,
std::vector< std::vector< double >>  periodicboundaries,
double  min,
bool  massweigh,
std::string  filename,
std::string  Nworld_filename,
std::string  Fworld_filename,
const std::vector< std::vector< double >> &  histdetails,
int  FBackupInterv,
double  unitconv,
double  timestep,
unsigned int  frequency 
)
inline

Constructor.

Parameters
worldMPI global communicator.
commMPI local communicator.
NLocal CV histogram.
NworldGlobal CV histogram.
FVector of grids holding local raw generalized force totals per bin per CV.
FworldVector of grids holding global raw generalized force totals per bin per CV.
restraintMinimum, maximum and spring constant for CV restraints.
isperiodicVector of bools that holds whether each CV is periodic or not.
periodicboundariesPeriodic boundaries of each CV.
minMinimum number of hist in a histogram bin before biasing is applied.
massweighTurn mass weighing on/off.
filenameName for output file.
Nworld_filenameName to read/write CV histogram.
Fworld_filenameName to read/write raw generalized force histograms.
histdetailsMinimum, maximum and number of bins for each CV.
FBackupIntervSet how often the adaptive force histograms are backed up.
unitconvUnit conversion from d(momentum)/d(time) to force.
timestepSimulation time step.
frequencyFrequency with which this method is invoked.

Constructs an instance of Adaptive Biasing Force method.

Note
The restraints should be outside of the range defined in histdetails_ by at least one bin size on each side.

Definition at line 216 of file ABF.h.

References PostIntegration(), PostSimulation(), and PreSimulation().

234  :
235  Method(frequency, world, comm), N_(N), Nworld_(Nworld), F_(F), Fworld_(Fworld), restraint_(restraint), isperiodic_(isperiodic), periodicboundaries_(periodicboundaries),
236  min_(min), wdotp1_(), wdotp2_(), Fold_(), massweigh_(massweigh),
237  biases_(), dim_(0), filename_(filename), Nworld_filename_(Nworld_filename), Fworld_filename_(Fworld_filename),
238  histdetails_(histdetails),
239  FBackupInterv_(FBackupInterv), unitconv_(unitconv), timestep_(timestep),
240  iteration_(0)
241  {
242  }
std::string filename_
Definition: ABF.h:126
Grid< int > * N_
To store number of local hits at a given CV bin.
Definition: ABF.h:66
std::string Fworld_filename_
Fworld print out filename, for restarts.
Definition: ABF.h:132
double unitconv_
Unit Conversion Constant from W dot P to Force.
Definition: ABF.h:158
std::vector< std::vector< double > > histdetails_
Histogram details.
Definition: ABF.h:143
Method(uint frequency, const MPI_Comm &world, const MPI_Comm &comm)
Constructor.
Definition: Method.h:61
std::vector< Grid< double > * > F_
To store running total of the local walker.
Definition: ABF.h:51
bool massweigh_
Mass weighing of bias enabled/disabled.
Definition: ABF.h:113
Eigen::VectorXd Fold_
To hold last iterations F_ value for removing bias.
Definition: ABF.h:110
std::vector< bool > isperiodic_
For each CV, holds whether that CV has periodic boundaries or not.
Definition: ABF.h:87
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:116
std::vector< std::vector< double > > restraint_
Information for a harmonic restraint to keep CV in the region of interest.
Definition: ABF.h:84
std::string Nworld_filename_
Nworld print out filename, for restarts.
Definition: ABF.h:129
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:119
std::vector< Grid< double > * > Fworld_
Will hold the global total, synced across walkers at every time step.
Definition: ABF.h:59
int min_
Definition: ABF.h:101
int FBackupInterv_
Integer to hold F estimate backup information.
Definition: ABF.h:149
Grid< int > * Nworld_
To store number of global hits at a given CV bin.
Definition: ABF.h:73
Eigen::VectorXd wdotp2_
To hold second to last iteration wdotp value for numerical derivative.
Definition: ABF.h:107
double timestep_
Timestep of integration.
Definition: ABF.h:172
uint iteration_
Iteration counter.
Definition: ABF.h:175
std::vector< std::vector< double > > periodicboundaries_
Holds periodic boundaries of CVs.
Definition: ABF.h:97
Eigen::VectorXd wdotp1_
To hold last iteration wdotp value for numerical derivative.
Definition: ABF.h:104

Here is the call graph for this function:

Member Function Documentation

bool SSAGES::ABF::boundsCheck ( const std::vector< double > &  CVs)
private

Checks whether the local walker is within CV bounds.

Parameters
CVsCurrent values of CVs.
Returns
Whether the walker is in bounds or not.

Definition at line 289 of file ABF.cpp.

290  {
291  for(size_t i = 0; i < dim_ ; ++i)
292  {
293  if((CVs[i] < N_->GetLower(i)) || (CVs[i] > N_->GetUpper(i)))
294  {
295  return false;
296  }
297  }
298  return true;
299  }
Grid< int > * N_
To store number of local hits at a given CV bin.
Definition: ABF.h:66
const std::vector< double > GetLower() const
Return the lower edges of the Grid.
Definition: GridBase.h:227
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:119
const std::vector< double > GetUpper() const
Return the upper edges of the Grid.
Definition: GridBase.h:258
ABF * SSAGES::ABF::Build ( const Json::Value &  json,
const MPI_Comm &  world,
const MPI_Comm &  comm,
const std::string &  path 
)
static

Definition at line 302 of file ABF.cpp.

References Json::Requirement::GetErrors(), Json::Requirement::HasErrors(), SSAGES::Grid< T >::LoadFromFile(), Json::ObjectRequirement::Parse(), and Json::ObjectRequirement::Validate().

Referenced by SetIteration().

306  {
307  ObjectRequirement validator;
308  Value schema;
309  Reader reader;
310 
311  reader.parse(JsonSchema::ABFMethod, schema);
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  uint wid = mxx::comm(world).rank()/mxx::comm(comm).size();
320 
321  //bool multiwalkerinput = false;
322 
323  // Obtain lower bound for each CV.
324  std::vector<double> minsCV;
325  for(auto& mins : json["CV_lower_bounds"])
326  if(mins.isArray())
327  {
328  throw BuildException({"Separate inputs for multi-walker not fully implemented. Please use global entries for CV ranges"});
329  //multiwalkerinput = true;
330  //for(auto& bound : mins)
331  // minsCV.push_back(bound.asDouble());
332  }
333  else
334  minsCV.push_back(mins.asDouble());
335 
336  // Obtain upper bound for each CV.
337  std::vector<double> maxsCV;
338  for(auto& maxs : json["CV_upper_bounds"])
339  /*if(maxs.isArray())
340  {
341  for(auto& bound : maxs)
342  maxsCV.push_back(bound.asDouble());
343  }*/
344  //else
345  maxsCV.push_back(maxs.asDouble());
346 
347  // Obtain number of bins for each CV dimension.
348  std::vector<int> binsCV;
349  for(auto& bins : json["CV_bins"])
350  binsCV.push_back(bins.asInt());
351 
352  // Obtain lower bounds for restraints for each CV.
353  std::vector<double> minsrestCV;
354  for(auto& mins : json["CV_restraint_minimums"])
355  /*if(mins.isArray())
356  {
357  for(auto& bound : mins)
358  minsrestCV.push_back(bound.asDouble());
359  }*/
360  //else
361  minsrestCV.push_back(mins.asDouble());
362 
363  // Obtain upper bounds for restraints for each CV.
364  std::vector<double> maxsrestCV;
365  for(auto& maxs : json["CV_restraint_maximums"])
366  /*if(maxs.isArray())
367  {
368  for(auto& bound : maxs)
369  maxsrestCV.push_back(bound.asDouble());
370  }*/
371  //else
372  maxsrestCV.push_back(maxs.asDouble());
373 
374  // Obtain harmonic spring constant for restraints for each CV.
375  std::vector<double> springkrestCV;
376  for(auto& bins : json["CV_restraint_spring_constants"])
377  springkrestCV.push_back(bins.asDouble());
378 
379  // Obtain periodicity information for restraints for each CV for the purpose of correctly applying restraints through periodic boundaries.
380  std::vector<bool> isperiodic;
381  for(auto& isperCV : json["CV_isperiodic"])
382  isperiodic.push_back(isperCV.asBool());
383 
384  // Obtain lower periodic boundary for each CV.
385  std::vector<double> minperboundaryCV;
386  for(auto& minsperCV : json["CV_periodic_boundary_lower_bounds"])
387  minperboundaryCV.push_back(minsperCV.asDouble());
388 
389  // Obtain upper periodic boundary for each CV.
390  std::vector<double> maxperboundaryCV;
391  for(auto& maxsperCV : json["CV_periodic_boundary_upper_bounds"])
392  maxperboundaryCV.push_back(maxsperCV.asDouble());
393 
394  // Verify inputs are all correct lengths.
395  if(!(( minsCV.size() == maxsCV.size() &&
396  maxsCV.size() == minsrestCV.size() &&
397  minsrestCV.size() == maxsrestCV.size()) &&
398  (binsCV.size() == springkrestCV.size() &&
399  springkrestCV.size() == isperiodic.size())))
400  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."});
401 
402  // Verify that all periodicity information is provided if CVs are periodic.
403  bool anyperiodic=false;
404  for(size_t i = 0; i<isperiodic.size(); ++i)
405  if(isperiodic[i])
406  {
407  anyperiodic = true;
408  if(!(isperiodic.size() == minperboundaryCV.size() &&
409  minperboundaryCV.size() == maxperboundaryCV.size()))
410  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."});
411  }
412 
413  int dim = binsCV.size();
414 
415  // Read in JSON info.
416  auto FBackupInterv = json.get("output_frequency", 1000).asInt();
417  auto unitconv = json.get("unit_conversion", 0).asDouble();
418  auto timestep = json.get("timestep", 2).asDouble();
419  auto min = json.get("minimum_count", 200).asDouble();
420  auto massweigh = json.get("mass_weighing",false).asBool();
421 
422  std::vector<std::vector<double>> histdetails;
423  std::vector<std::vector<double>> restraint;
424  std::vector<std::vector<double>> periodicboundaries;
425  std::vector<double> temp1(3);
426  std::vector<double> temp2(3);
427  std::vector<double> temp3(2);
428 
429  auto freq = json.get("frequency", 1).asInt();
430  auto filename = json.get("output_file", "F_out").asString();
431  auto Nworld_filename = json.get("Nworld_output_file", "Nworld").asString();
432  auto Fworld_filename = json.get("Fworld_output_file", "Fworld_cv").asString();
433 
434  // Generate the grids based on JSON.
435  Grid<int> *N;
436  Grid<int> *Nworld;
437  std::vector<Grid<double>*> F(dim);
438  std::vector<Grid<double>*> Fworld(dim);
439 
440  // This feature is disabled for now.
441  /*if(multiwalkerinput)
442  {
443  for(size_t i=0; i<dim; ++i)
444  {
445  temp1 = {minsCV[i+wid*dim], maxsCV[i+wid*dim], binsCV[i]};
446  temp2 = {minsrestCV[i+wid*dim], maxsrestCV[i+wid*dim], springkrestCV[i]};
447  histdetails.push_back(temp1);
448  restraint.push_back(temp2);
449  if(anyperiodic)
450  {
451  temp3 = {minperboundaryCV[i], maxperboundaryCV[i]};
452  periodicboundaries.push_back(temp3);
453  }
454  }
455  for(size_t i=0; i<dim; ++i)
456  {
457  minsCVperwalker[i] = histdetails[i][0];
458  maxsCVperwalker[i] = histdetails[i][1];
459  }
460 
461  N= new Grid<int>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
462  Nworld= new Grid<int>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
463 
464  for(auto& grid : F)
465  {
466  grid= new Grid<double>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
467  }
468  for(auto& grid : Fworld)
469  {
470  grid= new Grid<double>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
471  }
472 
473  }
474  */
475  //else
476  //{
477 
478  // Appropriately size the grids.
479 
480  N= new Grid<int>(binsCV, minsCV, maxsCV, isperiodic);
481  Nworld= new Grid<int>(binsCV, minsCV, maxsCV, isperiodic);
482 
483  for(auto& grid : F)
484  {
485  grid= new Grid<double>(binsCV, minsCV, maxsCV, isperiodic);
486  }
487  for(auto& grid : Fworld)
488  {
489  grid= new Grid<double>(binsCV, minsCV, maxsCV, isperiodic);
490  }
491 
492  for(size_t i=0; i<dim; ++i)
493  {
494  temp1 = {minsCV[i], maxsCV[i], double(binsCV[i])};
495  temp2 = {minsrestCV[i], maxsrestCV[i], springkrestCV[i]};
496  histdetails.push_back(temp1);
497  restraint.push_back(temp2);
498  if(anyperiodic)
499  {
500  temp3 = {minperboundaryCV[i], maxperboundaryCV[i]};
501  periodicboundaries.push_back(temp3);
502  }
503  }
504  //}
505 
506  // Check if previously saved grids exist. If so, check that data match and load grids.
507  if(std::ifstream(Nworld_filename) && std::ifstream(Fworld_filename+std::to_string(0)) && wid == 0)
508  {
509  std::cout << "Attempting to load data from a previous run of ABF." << std::endl;
510  N->LoadFromFile(Nworld_filename);
511  for(size_t i=0; i<dim; ++i)
512  if(std::ifstream(Fworld_filename+std::to_string(i)))
513  F[i]->LoadFromFile(Fworld_filename+std::to_string(i));
514  else
515  throw BuildException({"Some, but not all Fworld outputs were found. Please check that these are appropriate inputs, or clean the working directory of other Fworld and Nworld inputs."});
516  }
517 
518 
519  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);
520 
521 
522  if(json.isMember("iteration"))
523  m->SetIteration(json.get("iteration",0).asInt());
524 
525  return m;
526 
527  }
ABF(const MPI_Comm &world, const MPI_Comm &comm, Grid< int > *N, Grid< int > *Nworld, std::vector< Grid< double > * > F, std::vector< Grid< double > * > Fworld, std::vector< std::vector< double >> restraint, std::vector< bool > isperiodic, std::vector< std::vector< double >> periodicboundaries, double min, bool massweigh, std::string filename, std::string Nworld_filename, std::string Fworld_filename, const std::vector< std::vector< double >> &histdetails, int FBackupInterv, double unitconv, double timestep, unsigned int frequency)
Constructor.
Definition: ABF.h:216
bool HasErrors()
Check if errors have occured.
Definition: Requirement.h:86
virtual void Parse(Value json, const std::string &path) override
Parse JSON value to generate Requirement(s).
std::vector< std::string > GetErrors()
Get list of error messages.
Definition: Requirement.h:92
Requirements on an object.
virtual void Validate(const Value &json, const std::string &path) override
Validate JSON value.

Here is the call graph for this function:

Here is the caller graph for this function:

void SSAGES::ABF::CalcBiasForce ( const Snapshot snapshot,
const CVList cvs,
const std::vector< double > &  cvVals 
)
private

Computes the bias force.

Parameters
snapshotCurrent simulation snapshot.
cvsDetails of CVs.
cvValsCurrent values of CVs.

Definition at line 179 of file ABF.cpp.

References SSAGES::Snapshot::GetNumAtoms().

180  {
181  // Reset the bias.
182  biases_.resize(snapshot->GetNumAtoms(), Vector3{0,0,0});
183  for(auto& b : biases_)
184  b.setZero();
185 
186  // Compute bias if within bounds.
187  if(boundsCheck(cvVals))
188  {
189  for(int i = 0; i < dim_; ++i)
190  {
191  auto& grad = cvs[i]->GetGradient();
192  for(size_t j = 0; j < biases_.size(); ++j)
193  biases_[j] -= (Fworld_[i]->at(cvVals))*grad[j]/std::max(min_,Nworld_->at(cvVals));
194  }
195  }
196  // Otherwise, apply harmonic restraint if enabled.
197  else
198  {
199  for(int i = 0; i < dim_; ++i)
200  {
201  double cvVal = cvs[i]->GetValue();
202  auto k = 0.;
203  auto x0 = 0.;
204 
205  if(isperiodic_[i])
206  {
207  double periodsize = periodicboundaries_[i][1]-periodicboundaries_[i][0];
208  double cvRestrMidpoint = (restraint_[i][1]+restraint_[i][0])/2;
209  while((cvVal-cvRestrMidpoint) > periodsize/2)
210  cvVal -= periodsize;
211  while((cvVal-cvRestrMidpoint) < -periodsize/2)
212  cvVal += periodsize;
213  }
214 
215  if(cvVal < restraint_[i][0] && restraint_[i][2] > 0)
216  {
217  k = restraint_[i][2];
218  x0 = restraint_[i][0];
219  }
220  else if (cvVal > restraint_[i][1] && restraint_[i][2] > 0)
221  {
222  k = restraint_[i][2];
223  x0 = restraint_[i][1];
224  }
225 
226  auto& grad = cvs[i]->GetGradient();
227  for(size_t j = 0; j < biases_.size(); ++j)
228  biases_[j] -= (cvVal - x0)*k*grad[j];
229  }
230  }
231  }
std::vector< bool > isperiodic_
For each CV, holds whether that CV has periodic boundaries or not.
Definition: ABF.h:87
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:116
std::vector< std::vector< double > > restraint_
Information for a harmonic restraint to keep CV in the region of interest.
Definition: ABF.h:84
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:119
std::vector< Grid< double > * > Fworld_
Will hold the global total, synced across walkers at every time step.
Definition: ABF.h:59
bool boundsCheck(const std::vector< double > &CVs)
Checks whether the local walker is within CV bounds.
Definition: ABF.cpp:289
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
const T & at(const std::vector< int > &indices) const
Access Grid element read-only.
Definition: GridBase.h:536
int min_
Definition: ABF.h:101
Grid< int > * Nworld_
To store number of global hits at a given CV bin.
Definition: ABF.h:73
std::vector< std::vector< double > > periodicboundaries_
Holds periodic boundaries of CVs.
Definition: ABF.h:97

Here is the call graph for this function:

void SSAGES::ABF::PostIntegration ( Snapshot snapshot,
const class CVManager cvmanager 
)
overridevirtual

Post-integration hook.

Parameters
snapshotCurrent simulation snapshot.
cvmanagerCollective variable manager.

Post-integration hook is where processes carried out every timestep happen. First, information from current snapshot is retrieved and stored in variables as necessary. Then, coordinates in CV space are determined. Then, for each CV, the time derivative of w.p is calculated. Then, information is printed out if its enabled. Finally, bias is applied from current estimate of generalized force.

Coord holds where we are in CV space in current timestep.

Eigen::MatrixXd to hold the CV gradient.

std::vector<double> to hold CV values to address grid.

Implements SSAGES::Method.

Definition at line 75 of file ABF.cpp.

References SSAGES::CVManager::GetCVs(), SSAGES::Snapshot::GetForces(), SSAGES::Snapshot::GetMasses(), SSAGES::Snapshot::GetNumAtoms(), and SSAGES::Snapshot::GetVelocities().

Referenced by ABF().

76  {
77 
78  ++iteration_;
79 
80  // Gather information.
81  auto cvs = cvmanager.GetCVs(cvmask_);
82  auto& vels = snapshot->GetVelocities();
83  auto& mass = snapshot->GetMasses();
84  auto& forces = snapshot->GetForces();
85  auto n = snapshot->GetNumAtoms();
86 
87 
89  //int coord = histCoords(cvs);
90 
92  Eigen::MatrixXd J(dim_, 3*n);
93 
95  std::vector<double> cvVals(dim_);
96 
97  // Fill J and CV. Each column represents grad(CV) with flattened Cartesian elements.
98  for(int i = 0; i < dim_; ++i)
99  {
100  auto& grad = cvs[i]->GetGradient();
101  for(size_t j = 0; j < n; ++j)
102  J.block<1, 3>(i,3*j) = grad[j];
103  cvVals[i] = cvs[i]->GetValue();
104  }
105 
106  //* Calculate W using Darve's approach (http://mc.stanford.edu/cgi-bin/images/0/06/Darve_2008.pdf).
107  Eigen::MatrixXd Jmass = J.transpose();
108  if(massweigh_)
109  {
110  for(size_t i = 0; i < forces.size(); ++i)
111  Jmass.block(3*i, 0, 3, dim_) = Jmass.block(3*i, 0, 3, dim_)/mass[i];
112  }
113 
114  Eigen::MatrixXd Minv = J*Jmass;
115  MPI_Allreduce(MPI_IN_PLACE, Minv.data(), Minv.size(), MPI_DOUBLE, MPI_SUM, comm_);
116  Eigen::MatrixXd Wt = Minv.inverse()*Jmass.transpose();
117  // Fill momenta.
118  Eigen::VectorXd momenta(3*vels.size());
119  for(size_t i = 0; i < vels.size(); ++i)
120  momenta.segment<3>(3*i) = mass[i]*vels[i];
121 
122  // Compute dot(w,p)
123  Eigen::VectorXd wdotp = Wt*momenta;
124 
125 
126  // Reduce dot product across processors.
127  MPI_Allreduce(MPI_IN_PLACE, wdotp.data(), wdotp.size(), MPI_DOUBLE, MPI_SUM, comm_);
128 
129 
130  // Compute d(wdotp)/dt second order backwards finite difference.
131  // Adding old force removes bias.
132  Eigen::VectorXd dwdotpdt = unitconv_*(1.5*wdotp - 2.0*wdotp1_ + 0.5*wdotp2_)/timestep_ + Fold_;
133 
134  // If we are in bounds, sum force into running total.
135  if(boundsCheck(cvVals))
136  {
137  for(size_t i=0; i<dim_; ++i)
138  F_[i]->at(cvVals) += dwdotpdt[i];
139  N_->at(cvVals)++;
140  }
141 
142  // Reduce data across processors.
143  for(size_t i=0; i<dim_; ++i)
144  MPI_Allreduce(F_[i]->data(), Fworld_[i]->data(), (F_[i]->size()), MPI_DOUBLE, MPI_SUM, world_);
145  MPI_Allreduce(N_->data(), Nworld_->data(), N_->size(), MPI_INT, MPI_SUM, world_);
146 
147  // If we are in bounds, store the old summed force.
148  if(boundsCheck(cvVals))
149  {
150  for(size_t i=0; i < dim_; ++i)
151  Fold_[i] = Fworld_[i]->at(cvVals)/std::max(min_, Nworld_->at(cvVals));
152  }
153 
154  // Update finite difference time derivatives.
155  wdotp2_ = wdotp1_;
156  wdotp1_ = wdotp;
157 
158  // Write out data to file.
159  if(iteration_ % FBackupInterv_ == 0){
160  WriteData();
161  }
162 
163  // Calculate the bias from averaged F at current CV coordinates.
164  // Or apply harmonic restraints to return CVs back in bounds.
165  CalcBiasForce(snapshot, cvs, cvVals);
166 
167  // Update the forces in snapshot by adding in the force bias from each
168  // CV to each atom based on the gradient of the CV.
169  for (size_t j = 0; j < forces.size(); ++j)
170  forces[j] += biases_[j];
171  }
Grid< int > * N_
To store number of local hits at a given CV bin.
Definition: ABF.h:66
mxx::comm comm_
Local MPI communicator.
Definition: Method.h:47
mxx::comm world_
Global MPI communicator.
Definition: Method.h:46
void WriteData()
Writes out data to file.
Definition: ABF.cpp:235
double unitconv_
Unit Conversion Constant from W dot P to Force.
Definition: ABF.h:158
size_t size() const
Get the size of the internal storage vector.
Definition: GridBase.h:322
std::vector< Grid< double > * > F_
To store running total of the local walker.
Definition: ABF.h:51
bool massweigh_
Mass weighing of bias enabled/disabled.
Definition: ABF.h:113
Eigen::VectorXd Fold_
To hold last iterations F_ value for removing bias.
Definition: ABF.h:110
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:116
void CalcBiasForce(const Snapshot *snapshot, const CVList &cvs, const std::vector< double > &cvVals)
Computes the bias force.
Definition: ABF.cpp:179
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:119
std::vector< Grid< double > * > Fworld_
Will hold the global total, synced across walkers at every time step.
Definition: ABF.h:59
bool boundsCheck(const std::vector< double > &CVs)
Checks whether the local walker is within CV bounds.
Definition: ABF.cpp:289
std::vector< uint > cvmask_
Mask which identifies which CVs to act on.
Definition: Method.h:50
const T & at(const std::vector< int > &indices) const
Access Grid element read-only.
Definition: GridBase.h:536
T * data()
Get pointer to the internal data storage vector.
Definition: GridBase.h:334
int min_
Definition: ABF.h:101
int FBackupInterv_
Integer to hold F estimate backup information.
Definition: ABF.h:149
Grid< int > * Nworld_
To store number of global hits at a given CV bin.
Definition: ABF.h:73
Eigen::VectorXd wdotp2_
To hold second to last iteration wdotp value for numerical derivative.
Definition: ABF.h:107
double timestep_
Timestep of integration.
Definition: ABF.h:172
uint iteration_
Iteration counter.
Definition: ABF.h:175
Eigen::VectorXd wdotp1_
To hold last iteration wdotp value for numerical derivative.
Definition: ABF.h:104

Here is the call graph for this function:

Here is the caller graph for this function:

void SSAGES::ABF::PostSimulation ( Snapshot snapshot,
const class CVManager cvmanager 
)
overridevirtual

Post-simulation hook.

Parameters
snapshotCurrent simulation snapshot.
cvmanagerCollective variable manager.

Implements SSAGES::Method.

Definition at line 174 of file ABF.cpp.

Referenced by ABF().

175  {
176  WriteData();
177  }
void WriteData()
Writes out data to file.
Definition: ABF.cpp:235

Here is the caller graph for this function:

void SSAGES::ABF::PreSimulation ( Snapshot snapshot,
const class CVManager cvmanager 
)
overridevirtual

Pre-simulation hook.

Parameters
snapshotCurrent simulation snapshot.
cvmanagerCollective variable manager.

Implements SSAGES::Method.

Definition at line 39 of file ABF.cpp.

References SSAGES::CVManager::GetCVs(), and SSAGES::Snapshot::GetPositions().

Referenced by ABF().

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  // Size and initialize Fworld_ and Nworld_
56  auto nel = 1;
57 
58  // Initialize biases.
59  biases_.resize(snapshot->GetPositions().size(), Vector3{0, 0, 0});
60 
61  // Initialize w \dot p's for finite difference.
62  wdotp1_.setZero(dim_);
63  wdotp2_.setZero(dim_);
64  }
std::string filename_
Definition: ABF.h:126
std::ofstream worldout_
Output stream for F/N world data.
Definition: ABF.h:122
mxx::comm world_
Global MPI communicator.
Definition: Method.h:46
Eigen::VectorXd Fold_
To hold last iterations F_ value for removing bias.
Definition: ABF.h:110
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:116
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:119
std::vector< uint > cvmask_
Mask which identifies which CVs to act on.
Definition: Method.h:50
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
Eigen::VectorXd wdotp2_
To hold second to last iteration wdotp value for numerical derivative.
Definition: ABF.h:107
Eigen::VectorXd wdotp1_
To hold last iteration wdotp value for numerical derivative.
Definition: ABF.h:104

Here is the call graph for this function:

Here is the caller graph for this function:

void SSAGES::ABF::SetIteration ( const int  iter)
inline

Set iteration of the method.

Parameters
iterNew value for the iteration counter.

Definition at line 269 of file ABF.h.

References Build().

270  {
271  iteration_ = iter;
272  }
uint iteration_
Iteration counter.
Definition: ABF.h:175

Here is the call graph for this function:

Member Data Documentation

std::vector<Grid<double>*> SSAGES::ABF::F_
private

To store running total of the local walker.

This is a grid that holds dF/dCVx where x is one of the CVs. For N CVs, there will be N grids, and thus this vector will be N long.

Definition at line 51 of file ABF.h.

int SSAGES::ABF::FBackupInterv_
private

Integer to hold F estimate backup information.

Print every how many timesteps? ; -1: Do not backup during simulation

Definition at line 149 of file ABF.h.

std::string SSAGES::ABF::filename_
private

File name for world data. Principal output of the method.

Definition at line 126 of file ABF.h.

std::vector<Grid<double>*> SSAGES::ABF::Fworld_
private

Will hold the global total, synced across walkers at every time step.

This is a grid that holds dF/dCVx where x is one of the CVs. For N CVs, there will be N grids, and thus this vector will be N long. Global version.

Definition at line 59 of file ABF.h.

std::vector<std::vector<double> > SSAGES::ABF::histdetails_
private

Histogram details.

This is a 2 Dimensional object set up as the following: Histdetails is a vector of (nr of CV) vectors, ordered in CV order. Each of those vectors are 3 long. First entry of each vector holds the lower bound for that CV. Second entry of each vector holds the upper bound for that CV. Third entry of each vector holds the nr of bins for that CV dimension.

Definition at line 143 of file ABF.h.

int SSAGES::ABF::min_
private

The minimum number of hits required before full biasing, bias is F_[i]/max(N_[i],min_).

Definition at line 101 of file ABF.h.

Grid<int>* SSAGES::ABF::N_
private

To store number of local hits at a given CV bin.

Stores the number of times each bin was visited in CV space.

Definition at line 66 of file ABF.h.

Grid<int>* SSAGES::ABF::Nworld_
private

To store number of global hits at a given CV bin.

Stores the number of times each bin was visited in CV space. Global version.

Definition at line 73 of file ABF.h.

std::vector<std::vector<double> > SSAGES::ABF::periodicboundaries_
private

Holds periodic boundaries of CVs.

This is a 2 Dimensional object set up as the following: periodicboundaries_ is a vector of (nr of CV) vectors, ordered in CV order. Each of those vectors are 2 long. First entry of each vector holds the lower bound for that CV periodic boundary. Second entry of each vector holds the upper bound for that CV periodic boundary.

Definition at line 97 of file ABF.h.

std::vector<std::vector<double> > SSAGES::ABF::restraint_
private

Information for a harmonic restraint to keep CV in the region of interest.

This is a 2 Dimensional object set up as the following: restraint_ is a vector of (nr of CV) vectors, ordered in CV order. Each of those vectors are 3 long. First entry of each vector holds the lower bound for that CV restraint. Second entry of each vector holds the upper bound for that CV restraint. Third entry of each vector holds the spring constant for that CV restraint.

Definition at line 84 of file ABF.h.

double SSAGES::ABF::unitconv_
private

Unit Conversion Constant from W dot P to Force.

It is crucial that unit conversion is entered correctly. Unit conversion from d(momentum)/d(time) to force for the simulation. For LAMMPS using units real, this is 2390.06 (gram.angstrom/mole.femtosecond^2 -> kcal/mole.angstrom)

Definition at line 158 of file ABF.h.


The documentation for this class was generated from the following files: