Collective variable is a Rouse mode for a polymer chain comprised of N particle groups. More...

#include <RouseModeCV.h>

Inheritance diagram for SSAGES::RouseModeCV:

Public Member Functions
	RouseModeCV (const std::vector< Label > &groups, int p)
	Basic Constructor for Rouse Mode CV. More...

void	setMasses (const std::vector< Label > &groups, const Snapshot &snapshot)
	Helper function to determine masses of each group. More...

void	Initialize (const Snapshot &snapshot) override
	Initialize necessary variables. More...

void	Evaluate (const Snapshot &snapshot) override
	Evaluate the CV. More...

Public Member Functions inherited from SSAGES::CollectiveVariable
	CollectiveVariable ()
	Constructor.

virtual	~CollectiveVariable ()
	Destructor.

virtual void	Initialize (const class Snapshot &)
	Initialize CV. More...

virtual void	Evaluate (const class Snapshot &)=0
	Evaluate CV. More...

double	GetValue () const
	Get current value of the CV. More...

virtual double	GetMinimumImage (double) const
	Returns the minimum image of a CV based on the input location. More...

virtual double	GetPeriodicValue (double location) const
	Apply periodic boundaries to a given value. More...

const std::vector< Vector3 > &	GetGradient () const
	Get current gradient of the CV. More...

const Matrix3 &	GetBoxGradient () const
	Get gradient contribution to box. More...

const std::array< double, 2 > &	GetBoundaries ()
	Get CV boundaries. More...

virtual double	GetDifference (double location) const

Static Public Member Functions
static RouseModeCV *	Build (const Json::Value &json, const std::string &path)
	Set up collective variable. More...

Static Public Member Functions inherited from SSAGES::CollectiveVariable
static CollectiveVariable *	BuildCV (const Json::Value &json, const std::string &path)
	Set up collective variable. More...

Private Attributes
std::vector< Label >	groups_
	vector of groups of indices to define the particle groups

std::vector< double >	massg_
	vector of the total mass for each particle group

size_t	p_
	index of mode of interest as CV

size_t	N_
	number of Rouse beads in polymer chain

Vector3	xp_
	3d vector for containing vectorial rouse amplitude

std::vector< Vector3 >	r_
	vector of coordinate positions for each bead

Additional Inherited Members
Protected Attributes inherited from SSAGES::CollectiveVariable
std::vector< Vector3 >	grad_
	Gradient vector dCv/dxi.

Matrix3	boxgrad_
	Gradient w.r.t box vectors dCv/dHij.

double	val_
	Current value of CV.

std::array< double, 2 >	bounds_
	Bounds on CV.

Detailed Description

Collective variable is a Rouse mode for a polymer chain comprised of N particle groups.

This CV returns the value for the p'th Rouse mode, computed as Xp(t) ~ (1/N) sum_{i=1}^{N} ri(t)*cos[pi(n-0.5)p/N], where N is the number of particle groups, p is the mode index, ri is the center-of-mass position of a collection of atoms comprising the i'th bead in the N-bead polymer chain

Definition at line 39 of file RouseModeCV.h.

Constructor & Destructor Documentation

◆ RouseModeCV()

SSAGES::RouseModeCV::RouseModeCV	(	const std::vector< Label > &	groups,
		int	p
	)

inline

Basic Constructor for Rouse Mode CV.

Parameters

groups	- vector of vector of atom IDs, each group comprising a bead in the Rouse chain
p	- index for relevant Rouse mode

Definition at line 55 of file RouseModeCV.h.

                                                            : 
         groups_(groups), p_(p), N_( groups_.size())
         { massg_.resize( groups_.size(),0.0 ); }

References groups_, and massg_.

Referenced by Build().

Here is the caller graph for this function:

Member Function Documentation

◆ Build()

static RouseModeCV* SSAGES::RouseModeCV::Build	(	const Json::Value &	json,
		const std::string &	path
	)

inlinestatic

Set up collective variable.

Parameters

json	JSON input value.
path	Path for JSON path specification.

Returns: Pointer to the CV built by this function. nullptr in case of unknown error.

Builds a CV from a JSON node. Returns a pointer to the built cv. If an unknown error is encountered, this function will return a nullptr, but generally it will throw a BuildException on failure.

Warning: Object lifetime is the caller's responsibility.

Definition at line 192 of file RouseModeCV.h.

         {
             Json::ObjectRequirement validator;
             Json::Value schema;
             Json::CharReaderBuilder rbuilder;
             Json::CharReader* reader = rbuilder.newCharReader();
  
             reader->parse(JsonSchema::RouseModeCV.c_str(),
                           JsonSchema::RouseModeCV.c_str() + JsonSchema::RouseModeCV.size(),
                           &schema, nullptr);
             validator.Parse(schema, path);
  
             // Validate inputs.
             validator.Validate(json, path);
             if(validator.HasErrors())
                 throw BuildException(validator.GetErrors());
             
             std::vector<Label> groups;
             for (auto& group : json["groups"]) 
             {
                 groups.push_back({});
                 for(auto& id : group) 
                     groups.back().push_back(id.asInt());
             }
  
             auto mode = json.get("mode",0).asInt();
             return new RouseModeCV( groups, mode);  
         }

References Json::Requirement::GetErrors(), Json::Requirement::HasErrors(), Json::ObjectRequirement::Parse(), RouseModeCV(), and Json::ObjectRequirement::Validate().

Referenced by SSAGES::CollectiveVariable::BuildCV().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Evaluate()

void SSAGES::RouseModeCV::Evaluate ( const Snapshot & snapshot )

inlineoverride

Evaluate the CV.

Parameters

snapshot Current simulation snapshot.

Definition at line 121 of file RouseModeCV.h.

         {
  
             // Get data from snapshot.
             auto          ntot = snapshot.GetNumAtoms(); // total number of atoms
             const auto& masses = snapshot.GetMasses();   // mass of each atom
  
             // Initialize working variables
             double ppi_n = p_ * M_PI / N_;  // constant
             xp_.fill(0.0);      // vectorial Rouse mode
             r_.resize(N_);      // position vector for beads in Rouse chain (unwrapped)
             grad_.resize(ntot, Vector3{0,0,0});  // gradient set to 0.0 for all atoms
             std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
  
             // Iterate over all N_ atom groups and compute the center of mass for each
             std::vector<Vector3>     rcom; // vector of COM positions 
             for (size_t i = 0; i < N_; ++i) {
                 Label idi;  // list of indices
                 snapshot.GetLocalIndices(groups_[i], &idi);
                 rcom.push_back(snapshot.CenterOfMass(idi)); // center of mass for group
             }
  
             // Now compute differences vectors between the neighboring beads
             // accumulate displacements to reconstruct unwrapped polymer chain
             // for simplicity, we consider the first bead to be the reference position
             // in all snapshots
             r_[0] = rcom[0];
             for (size_t i = 1; i < N_; ++i) {
                 Vector3 dri = snapshot.ApplyMinimumImage(rcom[i] - rcom[i-1]);
                 r_[i] = r_[i-1] + dri;  // r_i = r_{i-1} + (r_i - r_{i-1})
             }
  
             // Determine the value of the Rouse coordinate
             // Xp(t) = 1/sqrt(N) * sum_{i=1}^{N} ri, p = 0
             // Xp(t) = sqrt(2/N) * sum_{i=1}^{N} ri * cos[p*pi/N*(i-0.5)], p = 1,...,N-1
             // Note: this solution is valid for homogeneous friction
             xp_.fill(0.0);
             for (size_t i = 0; i < N_; ++i) {
                 xp_ += r_[i]*cos ( ppi_n * (i+0.5) );
             }
             xp_ /= sqrt(N_);
             if ( p_ != 0 ) xp_ *= sqrt(2.0) ;
  
             // Compute Rouse vector norm as the CV
             // CV = sqrt(Xp*Xp), Xp = (Xp1,Xp2,Xp3)
             val_ = xp_.norm();
  
             // Now perform gradient operation
             // dCV/dxjd = (Xpd/CV)*(c/N)**0.5*sum_i=1^N cos[p*pi(i-0.5)/N] mj/Mi *delta_j({i})
             // delta_j({i}) = 1, if j in {i}, 0 otherwise
             Vector3   gradpcon  = xp_ / sqrt(N_) / val_;
             if ( p_ != 0) gradpcon *= sqrt(2.0);         // (Xpd/CV)*(c/N)**0.5
             for (size_t i = 0; i < N_; ++i) {
                 Label idi;  // list of indices
                 snapshot.GetLocalIndices(groups_[i], &idi);
                 // go over each atom in the group and add to its gradient
                 // Note: performance tradeoff here. All gradient elements have a common factor of
                 // Xpd/CV*sqrt(c/N) = prefactor, with c = 2 if p != 0
                 // this could be factored out, but if ntot >> number of atoms in groups
                 // then it won't be worth it to post multiply all gradient terms... 
                 // Could also go over loop again after to do the multiplication, but 
                 // that is troublesome if atom ids appear in multiple groups for some reason
                 double cosval = cos(ppi_n*(i+0.5)) / massg_[i]; // cos[p*pi(i-0.5)/N] / Mi
                 for (auto& id : idi) {
                     grad_[id] +=  gradpcon* cosval * masses[id];
                 }
             }
  
         }

References SSAGES::Snapshot::ApplyMinimumImage(), SSAGES::Snapshot::CenterOfMass(), SSAGES::Snapshot::GetLocalIndices(), SSAGES::Snapshot::GetMasses(), SSAGES::Snapshot::GetNumAtoms(), SSAGES::CollectiveVariable::grad_, groups_, massg_, N_, p_, r_, SSAGES::CollectiveVariable::val_, and xp_.

Here is the call graph for this function:

◆ Initialize()

void SSAGES::RouseModeCV::Initialize ( const Snapshot & snapshot )

inlineoverride

Initialize necessary variables.

Parameters

snapshot Current simulation snapshot.

Definition at line 80 of file RouseModeCV.h.

         {
             using std::to_string;
  
             // Check for valid p_
             if (p_ > groups_.size())
                 throw std::invalid_argument(
                     "RouseModeCV: Expected to find p to be less than " + 
                     to_string(groups_.size()) +" but found p = " + 
                     to_string(p_) 
                 );          
             // Check for valid groups
             for (size_t j = 0; j < groups_.size(); ++j) {
                 auto nj = groups_[j].size();
                 // Make sure atom IDs in the group are somewhere
                 std::vector<int> found(nj,0);
                 for (size_t i = 0; i < nj; ++i) {
                     if(snapshot.GetLocalIndex(groups_[j][i]) != -1)
                         found[i] = 1;
                 }
  
                 MPI_Allreduce(MPI_IN_PLACE, found.data(), nj, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
                 unsigned njtot = std::accumulate(found.begin(), found.end(), 0, std::plus<int>());
                 
                 if(njtot != nj)
                     throw std::invalid_argument(
                         "RouseModeCV: Expected to find " + 
                         to_string(nj) + 
                         " atoms in group " + to_string(j) +", but only found " + 
                         to_string(njtot) + "."
                     );          
  
             }
             // Set the masses of each particle group in massg_
             this->setMasses(groups_, snapshot);
         }

References SSAGES::Snapshot::GetCommunicator(), SSAGES::Snapshot::GetLocalIndex(), groups_, p_, and setMasses().

Here is the call graph for this function:

◆ setMasses()

void SSAGES::RouseModeCV::setMasses	(	const std::vector< Label > &	groups,
		const Snapshot &	snapshot
	)

inline

Helper function to determine masses of each group.

Note: this here assumes that the masses of each group are not changing during the simulation, which is likely typical...

Parameters

snapshot	Current simulation snapshot.
groups	- vector of vector of atom IDs, each group comprising a bead in the Rouse chain

Definition at line 65 of file RouseModeCV.h.

                                                                                   {
             // Compute mass of each group and store to massg_
             Label listi;
             for (size_t i = 0; i < N_; ++i) {
                 listi     = groups[i];     // MW: can this be used directly in Total Mass
                 Label idi;
                 snapshot.GetLocalIndices(listi, &idi);
                 massg_[i] = snapshot.TotalMass(idi);
             } 
         }