SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
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 (unsigned int frequency, const MPI_Comm &world, const MPI_Comm &comm)
 Constructor. More...
 
void SetCVMask (const std::vector< unsigned int > &mask)
 Sets the collective variable mask.
 
virtual ~Method ()
 Destructor.
 
- Public Member Functions inherited from SSAGES::EventListener
 EventListener (unsigned int frequency)
 Constructor. More...
 
unsigned int 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)
 Build a derived method from JSON node. More...
 
- 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

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< 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...
 
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.
 
unsigned int 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< unsigned int > 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 40 of file ABF.h.

Constructor & Destructor Documentation

◆ ABF()

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 213 of file ABF.h.

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

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

Member Function Documentation

◆ boundsCheck()

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 286 of file ABF.cpp.

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  }
const std::vector< double > GetLower() const
Return the lower edges of the Grid.
Definition: GridBase.h:226
Grid< int > * N_
To store number of local hits at a given CV bin.
Definition: ABF.h:48
const std::vector< double > GetUpper() const
Return the upper edges of the Grid.
Definition: GridBase.h:257
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:116

◆ Build()

ABF * SSAGES::ABF::Build ( const Json::Value &  json,
const MPI_Comm &  world,
const MPI_Comm &  comm,
const std::string &  path 
)
static

Build a derived method from JSON node.

Parameters
jsonJSON Value containing all input information.
worldMPI global communicator.
commMPI local communicator.
pathPath for JSON path specification.
Returns
Pointer to the Method built. nullptr if an unknown error occurred.

This function builds a registered method from a JSON node. The difference between this function and "Build" is that this automatically determines the appropriate derived type based on the JSON node information.

Note
Object lifetime is the caller's responsibility.

Definition at line 299 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().

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  }
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.
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:213
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:

◆ CalcBiasForce()

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 176 of file ABF.cpp.

References SSAGES::Snapshot::GetNumAtoms().

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  }
const T & at(const std::vector< int > &indices) const
Access Grid element read-only.
Definition: GridBase.h:541
std::vector< bool > isperiodic_
For each CV, holds whether that CV has periodic boundaries or not.
Definition: ABF.h:84
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:113
std::vector< std::vector< double > > restraint_
Information for a harmonic restraint to keep CV in the region of interest.
Definition: ABF.h:81
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:116
std::vector< Grid< double > * > Fworld_
Will hold the global total, synced across walkers at every time step.
Definition: ABF.h:70
bool boundsCheck(const std::vector< double > &CVs)
Checks whether the local walker is within CV bounds.
Definition: ABF.cpp:286
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
int min_
Definition: ABF.h:98
Grid< int > * Nworld_
To store number of global hits at a given CV bin.
Definition: ABF.h:55
std::vector< std::vector< double > > periodicboundaries_
Holds periodic boundaries of CVs.
Definition: ABF.h:94
Here is the call graph for this function:

◆ PostIntegration()

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 72 of file ABF.cpp.

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

Referenced by ABF().

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  }
Grid< int > * N_
To store number of local hits at a given CV bin.
Definition: ABF.h:48
size_t size() const
Get the size of the internal storage vector.
Definition: GridBase.h:321
const T & at(const std::vector< int > &indices) const
Access Grid element read-only.
Definition: GridBase.h:541
mxx::comm comm_
Local MPI communicator.
Definition: Method.h:47
unsigned int iteration_
Iteration counter.
Definition: ABF.h:172
mxx::comm world_
Global MPI communicator.
Definition: Method.h:46
void WriteData()
Writes out data to file.
Definition: ABF.cpp:232
double unitconv_
Unit Conversion Constant from W dot P to Force.
Definition: ABF.h:155
std::vector< Grid< double > * > F_
To store running total of the local walker.
Definition: ABF.h:62
bool massweigh_
Mass weighing of bias enabled/disabled.
Definition: ABF.h:110
Eigen::VectorXd Fold_
To hold last iterations F_ value for removing bias.
Definition: ABF.h:107
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:113
void CalcBiasForce(const Snapshot *snapshot, const CVList &cvs, const std::vector< double > &cvVals)
Computes the bias force.
Definition: ABF.cpp:176
std::vector< unsigned int > cvmask_
Mask which identifies which CVs to act on.
Definition: Method.h:50
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:116
std::vector< Grid< double > * > Fworld_
Will hold the global total, synced across walkers at every time step.
Definition: ABF.h:70
bool boundsCheck(const std::vector< double > &CVs)
Checks whether the local walker is within CV bounds.
Definition: ABF.cpp:286
T * data()
Get pointer to the internal data storage vector.
Definition: GridBase.h:333
int min_
Definition: ABF.h:98
int FBackupInterv_
Integer to hold F estimate backup information.
Definition: ABF.h:146
Grid< int > * Nworld_
To store number of global hits at a given CV bin.
Definition: ABF.h:55
Eigen::VectorXd wdotp2_
To hold second to last iteration wdotp value for numerical derivative.
Definition: ABF.h:104
double timestep_
Timestep of integration.
Definition: ABF.h:169
Eigen::VectorXd wdotp1_
To hold last iteration wdotp value for numerical derivative.
Definition: ABF.h:101
Here is the call graph for this function:
Here is the caller graph for this function:

◆ PostSimulation()

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 171 of file ABF.cpp.

Referenced by ABF().

172  {
173  WriteData();
174  }
void WriteData()
Writes out data to file.
Definition: ABF.cpp:232
Here is the caller graph for this function:

◆ PreSimulation()

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  // 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  }
std::string filename_
Definition: ABF.h:123
std::ofstream worldout_
Output stream for F/N world data.
Definition: ABF.h:119
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:107
std::vector< Vector3 > biases_
Biases applied to atoms each timestep.
Definition: ABF.h:113
std::vector< unsigned int > cvmask_
Mask which identifies which CVs to act on.
Definition: Method.h:50
unsigned int dim_
Number of CVs in system.
Definition: ABF.h:116
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:104
Eigen::VectorXd wdotp1_
To hold last iteration wdotp value for numerical derivative.
Definition: ABF.h:101
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SetIteration()

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

Set iteration of the method.

Parameters
iterNew value for the iteration counter.

Definition at line 266 of file ABF.h.

References Build().

267  {
268  iteration_ = iter;
269  }
unsigned int iteration_
Iteration counter.
Definition: ABF.h:172
Here is the call graph for this function:

Member Data Documentation

◆ F_

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 62 of file ABF.h.

◆ FBackupInterv_

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 146 of file ABF.h.

◆ filename_

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

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

Definition at line 123 of file ABF.h.

◆ Fworld_

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 70 of file ABF.h.

◆ histdetails_

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 140 of file ABF.h.

◆ min_

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 98 of file ABF.h.

◆ N_

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 48 of file ABF.h.

◆ Nworld_

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 55 of file ABF.h.

◆ periodicboundaries_

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 94 of file ABF.h.

◆ restraint_

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 81 of file ABF.h.

◆ unitconv_

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 155 of file ABF.h.


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