SSAGES::ABF Class Reference

Adaptive Biasing Force Algorithm. More...

#include <ABF.h>

Inheritance diagram for SSAGES::ABF:

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
	Method call prior to simulation initiation. More...
void	PostIntegration (Snapshot *snapshot, const class CVManager &cvmanager) override
	Method call post integration. More...
void	PostSimulation (Snapshot *snapshot, const class CVManager &cvmanager) override
	Method call prior to simulation initiation. 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. More...
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 ABF *	Build (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 Method *	BuildMethod (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::EventListener
static unsigned int	GetWalkerID (const MPI_Comm &world, const MPI_Comm &comm)
	Get walker ID number of specified communicator. More...
static unsigned int	GetNumWalkers (const MPI_Comm &world, const MPI_Comm &comm)
	Get total number of walkers in the simulation. More...
static bool	IsMasterRank (const MPI_Comm &comm)
	Check if current processor is master. 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< Vector3 >	biases_
	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

world	MPI global communicator.
comm	MPI local communicator.
N	Local CV histogram.
Nworld	Global CV histogram.
F	Vector of grids holding local raw generalized force totals per bin per CV.
Fworld	Vector of grids holding global raw generalized force totals per bin per CV.
restraint	Minimum, maximum and spring constant for CV restraints.
isperiodic	Vector of bools that holds whether each CV is periodic or not.
periodicboundaries	Periodic boundaries of each CV.
min	Minimum number of hist in a histogram bin before biasing is applied.
massweigh	Turn mass weighing on/off.
filename	Name for output file.
Nworld_filename	Name to read/write CV histogram.
Fworld_filename	Name to read/write raw generalized force histograms.
histdetails	Minimum, maximum and number of bins for each CV.
FBackupInterv	Set how often the adaptive force histograms are backed up.
unitconv	Unit conversion from d(momentum)/d(time) to force.
timestep	Simulation time step.
frequency	Frequency 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.

                                    :
        Method(frequency, world, comm), N_(N), Nworld_(Nworld), F_(F), Fworld_(Fworld),  restraint_(restraint), isperiodic_(isperiodic), periodicboundaries_(periodicboundaries),
        min_(min), wdotp1_(), wdotp2_(), Fold_(), massweigh_(massweigh),
        biases_(), dim_(0), filename_(filename), Nworld_filename_(Nworld_filename), Fworld_filename_(Fworld_filename),
        histdetails_(histdetails), 
        FBackupInterv_(FBackupInterv), unitconv_(unitconv), timestep_(timestep),
        iteration_(0)
        {
        }

Member Function Documentation

boundsCheck()

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

private

Checks whether the local walker is within CV bounds.

Parameters

CVs	Current values of CVs.

Returns

Whether the walker is in bounds or not.

Definition at line 284 of file ABF.cpp.

    {
        for(size_t i = 0; i < dim_ ; ++i)
        {
            if((CVs[i] < N_->GetLower(i)) || (CVs[i] > N_->GetUpper(i)))
            {
                return false;
            }
        }
        return true;
    }

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

json	JSON Value containing all input information.
world	MPI global communicator.
comm	MPI local communicator.
path	Path 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 297 of file ABF.cpp.

    {
        ObjectRequirement validator;
        Value schema;
        CharReaderBuilder rbuilder;
        CharReader* reader = rbuilder.newCharReader();
 
        reader->parse(JsonSchema::ABFMethod.c_str(),
                      JsonSchema::ABFMethod.c_str() + JsonSchema::ABFMethod.size(),
                      &schema, nullptr);
        validator.Parse(schema, path);
 
        // Validate inputs.
        validator.Validate(json, path);
        if(validator.HasErrors())
            throw BuildException(validator.GetErrors());
 
        //bool multiwalkerinput = false;
 
        // Obtain lower bound for each CV.
        std::vector<double> minsCV;
        for(auto& mins : json["CV_lower_bounds"])
            if(mins.isArray())
                {
                throw BuildException({"Separate inputs for multi-walker not fully implemented. Please use global entries for CV ranges"});
                //multiwalkerinput = true;
                //for(auto& bound : mins)
                //  minsCV.push_back(bound.asDouble());
                }
            else
                minsCV.push_back(mins.asDouble());
            
        // Obtain upper bound for each CV.
        std::vector<double> maxsCV;
        for(auto& maxs : json["CV_upper_bounds"])
            /*if(maxs.isArray())
                {
                for(auto& bound : maxs)
                    maxsCV.push_back(bound.asDouble());
                }*/
            //else
                maxsCV.push_back(maxs.asDouble());
 
        // Obtain number of bins for each CV dimension.
        std::vector<int> binsCV;
        for(auto& bins : json["CV_bins"])
            binsCV.push_back(bins.asInt());
 
        // Obtain lower bounds for restraints for each CV.
        std::vector<double> minsrestCV;
        for(auto& mins : json["CV_restraint_minimums"])
            /*if(mins.isArray())
                {
                for(auto& bound : mins)
                    minsrestCV.push_back(bound.asDouble());
                }*/
            //else
                minsrestCV.push_back(mins.asDouble());
 
        // Obtain upper bounds for restraints for each CV.      
        std::vector<double> maxsrestCV;
        for(auto& maxs : json["CV_restraint_maximums"])
            /*if(maxs.isArray())
                {
                for(auto& bound : maxs)
                    maxsrestCV.push_back(bound.asDouble());
                }*/
            //else
                maxsrestCV.push_back(maxs.asDouble());
        
        // Obtain harmonic spring constant for restraints for each CV.
        std::vector<double> springkrestCV;
        for(auto& springkrest : json["CV_restraint_spring_constants"])
        {
            if(springkrest.asDouble() < 0)
            {
                throw BuildException({"Restraint spring constants must be non-negative."});
            }
            springkrestCV.push_back(springkrest.asDouble());
        }
 
        // Obtain periodicity information for restraints for each CV for the purpose of correctly applying restraints through periodic boundaries.
        std::vector<bool> isperiodic;
        for(auto& isperCV : json["CV_isperiodic"])
            isperiodic.push_back(isperCV.asBool());
        
        // Obtain lower periodic boundary for each CV.
        std::vector<double> minperboundaryCV;
        for(auto& minsperCV : json["CV_periodic_boundary_lower_bounds"])
            minperboundaryCV.push_back(minsperCV.asDouble());
 
        // Obtain upper periodic boundary for each CV.
        std::vector<double> maxperboundaryCV;
        for(auto& maxsperCV : json["CV_periodic_boundary_upper_bounds"])
            maxperboundaryCV.push_back(maxsperCV.asDouble());
 
        // Verify inputs are all correct lengths.
        if(!((   minsCV.size() == maxsCV.size() && 
             maxsCV.size() == minsrestCV.size() &&
             minsrestCV.size() == maxsrestCV.size()) &&
            (binsCV.size() == springkrestCV.size() &&
             springkrestCV.size() == isperiodic.size())))       
            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."});
 
        // Verify that all periodicity information is provided if CVs are periodic.
        bool anyperiodic=false;
        for(size_t i = 0; i<isperiodic.size(); ++i)
            if(isperiodic[i])
            {
                anyperiodic = true;
                if(!(isperiodic.size() == minperboundaryCV.size() &&
                     minperboundaryCV.size() == maxperboundaryCV.size()))
                    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."});
            }
 
        size_t dim = binsCV.size();
        
        // Read in JSON info.
        auto FBackupInterv = json.get("output_frequency", 1000).asInt();
        auto unitconv = json.get("unit_conversion", 0).asDouble();
        auto timestep = json.get("timestep", 2).asDouble();
        auto min = json.get("minimum_count", 200).asDouble();
        auto massweigh = json.get("mass_weighing",false).asBool();          
 
        std::vector<std::vector<double>> histdetails;
        std::vector<std::vector<double>> restraint;
        std::vector<std::vector<double>> periodicboundaries;
        std::vector<double> temp1(3);
        std::vector<double> temp2(3);
        std::vector<double> temp3(2);
 
        auto freq = json.get("frequency", 1).asInt();
        auto filename = json.get("output_file", "F_out").asString();
        auto Nworld_filename = json.get("Nworld_output_file", "Nworld").asString();
        auto Fworld_filename = json.get("Fworld_output_file", "Fworld_cv").asString();
 
        // Generate the grids based on JSON.
        Grid<int> *N;
        Grid<int> *Nworld;
        std::vector<Grid<double>*> F(dim);
        std::vector<Grid<double>*> Fworld(dim);
 
        //unsigned int wid = GetWalkerID(world, comm);
 
        // This feature is disabled for now.
 
        /*if(multiwalkerinput)
        {
            for(size_t i=0; i<dim; ++i)
            {
                temp1 = {minsCV[i+wid*dim], maxsCV[i+wid*dim], binsCV[i]};
                temp2 = {minsrestCV[i+wid*dim], maxsrestCV[i+wid*dim], springkrestCV[i]};
                histdetails.push_back(temp1);
                restraint.push_back(temp2);
                if(anyperiodic)
                {
                    temp3 = {minperboundaryCV[i], maxperboundaryCV[i]};
                    periodicboundaries.push_back(temp3);
                }
            }
            for(size_t i=0; i<dim; ++i)
            {
                minsCVperwalker[i] = histdetails[i][0];
                maxsCVperwalker[i] = histdetails[i][1];
            }
 
            N= new Grid<int>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
            Nworld= new Grid<int>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
            
            for(auto& grid : F)
            {
                grid= new Grid<double>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
            }
            for(auto& grid : Fworld)
            {
                grid= new Grid<double>(binsCV, minsCVperwalker, maxsCVperwalker, isperiodic);
            }
            
        }
        */
        //else
        //{
        
        // Appropriately size the grids.
        
            N = new Grid<int>(binsCV, minsCV, maxsCV, isperiodic);
            Nworld = new Grid<int>(binsCV, minsCV, maxsCV, isperiodic);
            
            for(auto& grid : F)
            {
                grid = new Grid<double>(binsCV, minsCV, maxsCV, isperiodic);
            }
            for(auto& grid : Fworld)
            {
                grid = new Grid<double>(binsCV, minsCV, maxsCV, isperiodic);
            }
 
            for(size_t i = 0; i < dim; ++i)
            {
                temp1 = {minsCV[i], maxsCV[i], double(binsCV[i])};
                temp2 = {minsrestCV[i], maxsrestCV[i], springkrestCV[i]};
                histdetails.push_back(temp1);
                restraint.push_back(temp2);
                if(anyperiodic)
                {
                    temp3 = {minperboundaryCV[i], maxperboundaryCV[i]};
                    periodicboundaries.push_back(temp3);
                }
            }
        //}
 
        // Check if a restart is requested.
        bool restart = json.get("restart",false).asBool();
 
        // Either load previous grid, or back it up.
        if(IsMasterRank(world))
        {
            if(restart)
            {
                std::cout << "Attempting to load data from a previous run of ABF..." << std::endl;
                N->LoadFromFile(Nworld_filename);
                for(size_t i = 0; i < dim; ++i)
                {
                    F[i]->LoadFromFile(Fworld_filename+std::to_string(i));
                }
            }
            else
            {
                std::ifstream Nworldbackupsource(Nworld_filename, std::ios::binary);
                if(Nworldbackupsource)
                {
                    std::cout << "Backing up previous copy of Nworld." << std::endl;
                    std::ofstream Nworldbackuptarget(Nworld_filename+"_backup", std::ios::binary);
                    Nworldbackuptarget << Nworldbackupsource.rdbuf();
                }
                for(size_t i = 0; i < dim; ++i)
                {
                    std::ifstream Fworldbackupsource(Fworld_filename+std::to_string(i), std::ios::binary);
                    if(Fworldbackupsource)
                    {
                        std::cout << "Backing up previous copy of Fworld_cv"+std::to_string(i)+"." << std::endl;
                        std::ofstream Fworldbackuptarget(Fworld_filename+std::to_string(i)+"_backup", std::ios::binary);
                        Fworldbackuptarget << Fworldbackupsource.rdbuf();
                    }
                }
            }
        }
 
        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);
 
        if(json.isMember("iteration"))
            m->SetIteration(json.get("iteration",0).asInt());
 
        return m;
 
    }

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

Here is the call 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

snapshot	Current simulation snapshot.
cvs	Details of CVs.
cvVals	Current values of CVs.

Definition at line 176 of file ABF.cpp.

    {
        // Reset the bias.
        biases_.resize(snapshot->GetNumAtoms(), Vector3{0,0,0});
        for(auto& b : biases_)
            b.setZero();
        
        // Compute bias if within bounds.
        if(boundsCheck(cvVals))
        {
            for(size_t i = 0; i < dim_; ++i)
            {
                auto& grad = cvs[i]->GetGradient();
                for(size_t j = 0; j < biases_.size(); ++j)
                    biases_[j] -= (Fworld_[i]->at(cvVals))*grad[j]/std::max(min_,Nworld_->at(cvVals));
            }
        }
        // Otherwise, apply harmonic restraint if enabled.
        else
        {
            for(size_t i = 0; i < dim_; ++i)
            {
                double cvVal = cvs[i]->GetValue();
                auto k = 0.;
                auto x0 = 0.;
                
                if(isperiodic_[i])
                {
                    double periodsize = periodicboundaries_[i][1]-periodicboundaries_[i][0];
                    double cvRestrMidpoint = (restraint_[i][1]+restraint_[i][0])/2;
                    while((cvVal-cvRestrMidpoint) > periodsize/2)
                        cvVal -= periodsize;
                    while((cvVal-cvRestrMidpoint) < -periodsize/2)
                        cvVal += periodsize;
                }
 
                if(cvVal < restraint_[i][0])
                {
                    k = restraint_[i][2];
                    x0 = restraint_[i][0];
                }
                else if (cvVal > restraint_[i][1])
                {
                    k = restraint_[i][2];
                    x0 = restraint_[i][1];
                }
 
                auto& grad = cvs[i]->GetGradient();
                for(size_t j = 0; j < biases_.size(); ++j)
                    biases_[j] -= (cvVal - x0)*k*grad[j];
            }   
        }
    }

References SSAGES::Snapshot::GetNumAtoms().

Here is the call graph for this function:

PostIntegration()

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

overridevirtual

Method call post integration.

Post-integration hook.

Parameters

snapshot	Pointer to the simulation snapshot.
cvmanager	Collective variable manager.

This function will be called after each integration step.

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.

    {
        
        ++iteration_;
 
        // Gather information.
        auto cvs = cvmanager.GetCVs(cvmask_);
        auto& vels = snapshot->GetVelocities();
        auto& mass = snapshot->GetMasses();
        auto& forces = snapshot->GetForces();
        auto n = snapshot->GetNumAtoms();
        
 
        //int coord = histCoords(cvs);
 
        Eigen::MatrixXd J(dim_, 3*n);
        
        std::vector<double> cvVals(dim_);
 
        // Fill J and CV. Each column represents grad(CV) with flattened Cartesian elements. 
        for(size_t i = 0; i < dim_; ++i)
        {
            auto& grad = cvs[i]->GetGradient();
            for(size_t j = 0; j < n; ++j)
                J.block<1, 3>(i,3*j) = grad[j];
            cvVals[i] = cvs[i]->GetValue();
        }
 
        //* Calculate W using Darve's approach (http://mc.stanford.edu/cgi-bin/images/0/06/Darve_2008.pdf).
        Eigen::MatrixXd Jmass = J.transpose();
        if(massweigh_)
        {
            for(size_t i = 0; i < forces.size(); ++i)
                Jmass.block(3*i, 0, 3, dim_) = Jmass.block(3*i, 0, 3, dim_)/mass[i];
        }
 
        Eigen::MatrixXd Minv = J*Jmass;
        MPI_Allreduce(MPI_IN_PLACE, Minv.data(), Minv.size(), MPI_DOUBLE, MPI_SUM, comm_);
        Eigen::MatrixXd Wt = Minv.inverse()*Jmass.transpose();  
        // Fill momenta.
        Eigen::VectorXd momenta(3*vels.size());
        for(size_t i = 0; i < vels.size(); ++i)
            momenta.segment<3>(3*i) = mass[i]*vels[i];
        
        // Compute dot(w,p)
        Eigen::VectorXd wdotp = Wt*momenta;
 
        
        // Reduce dot product across processors.
        MPI_Allreduce(MPI_IN_PLACE, wdotp.data(), wdotp.size(), MPI_DOUBLE, MPI_SUM, comm_);
        
    
        // Compute d(wdotp)/dt second order backwards finite difference. 
        // Adding old force removes bias. 
        Eigen::VectorXd dwdotpdt = unitconv_*(1.5*wdotp - 2.0*wdotp1_ + 0.5*wdotp2_)/timestep_ + Fold_;
 
        // If we are in bounds, sum force into running total.
        if(boundsCheck(cvVals) && IsMasterRank(comm_))
        {
            for(size_t i = 0; i < dim_; ++i)
                F_[i]->at(cvVals) += dwdotpdt[i];
            N_->at(cvVals)++;               
        }
        
        // Reduce data across processors.
        for(size_t i = 0; i < dim_; ++i)                
            MPI_Allreduce(F_[i]->data(), Fworld_[i]->data(), (F_[i]->size()), MPI_DOUBLE, MPI_SUM, world_); 
        MPI_Allreduce(N_->data(), Nworld_->data(), N_->size(), MPI_INT, MPI_SUM, world_);
 
        // If we are in bounds, store the old summed force.
        if(boundsCheck(cvVals))
        {
            for(size_t i=0; i < dim_; ++i)
                Fold_[i] = Fworld_[i]->at(cvVals)/std::max(min_, Nworld_->at(cvVals));
        }
    
        // Update finite difference time derivatives.
        wdotp2_ = wdotp1_;
        wdotp1_ = wdotp;        
 
        // Write out data to file.
        if(iteration_ % FBackupInterv_ == 0){
            WriteData();
        }
 
        // Calculate the bias from averaged F at current CV coordinates.
        // Or apply harmonic restraints to return CVs back in bounds.
        CalcBiasForce(snapshot, cvs, cvVals);       
 
        // Update the forces in snapshot by adding in the force bias from each
        // CV to each atom based on the gradient of the CV.
        for (size_t j = 0; j < forces.size(); ++j)
            forces[j] += biases_[j];    
    }

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

Here is the call graph for this function:

PostSimulation()

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

overridevirtual

Method call prior to simulation initiation.

Parameters

snapshot	Pointer to the simulation snapshot.
cvmanager	Collective variable manager.

This function will be called before the simulation is started.

Implements SSAGES::Method.

Definition at line 171 of file ABF.cpp.

    {
        WriteData();
    }

PreSimulation()

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

overridevirtual

Method call prior to simulation initiation.

Parameters

snapshot	Pointer to the simulation snapshot.
cvmanager	Collective variable manager.

This function will be called before the simulation is started.

Implements SSAGES::Method.

Definition at line 39 of file ABF.cpp.

    {
        // Open/close outfile to create it fresh. 
        if(IsMasterRank(world_))
        {
            worldout_.open(filename_);
            worldout_.close();
        }
        
        // Convenience. Number of CVs.
        auto cvs = cvmanager.GetCVs(cvmask_);
        dim_ = cvs.size();
 
        // Size and initialize Fold_
        Fold_.setZero(dim_);
 
        // Initialize biases.
        biases_.resize(snapshot->GetPositions().size(), Vector3{0, 0, 0});
        
        // Initialize w \dot p's for finite difference. 
        wdotp1_.setZero(dim_);
        wdotp2_.setZero(dim_);  
    }

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

Here is the call graph for this function:

SetIteration()

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

inline

Set iteration of the method.

Parameters

iter	New value for the iteration counter.

Definition at line 254 of file ABF.h.

        {
            iteration_ = iter;
        }

References iteration_.

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:

• /home/michael/development/SSAGES/src/Methods/ABF.h
• /home/michael/development/SSAGES/src/Methods/ABF.cpp