Collective variable on components of the gyration tensor. More...

#include <GyrationTensorCV.h>

Inheritance diagram for SSAGES::GyrationTensorCV:

Public Member Functions
	GyrationTensorCV (const Label &atomids, GyrationTensor component)
	Constructor. 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.

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.

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

virtual double	GetDifference (double location) const

Static Public Member Functions
static GyrationTensorCV *	Build (const Json::Value &json, const std::string &path)

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
Label	atomids_
	IDs of the atoms used for calculation.

GyrationTensor	component_
	Component of gyration tensor to compute.

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 on components of the gyration tensor.

Collective variable on components of gyration tensor. Depending on the user selection, this will specify a principal moment, radius of gyration, or another shape descriptor.

Definition at line 51 of file GyrationTensorCV.h.

Constructor & Destructor Documentation

SSAGES::GyrationTensorCV::GyrationTensorCV	(	const Label &	atomids,
		GyrationTensor	component
	)

inline

Constructor.

Parameters

atomids	IDs of the atoms defining gyration tensor.
component	Specification of component to compute.

Construct a GyrationTensorCV.

Definition at line 66 of file GyrationTensorCV.h.

                                                                          :
         atomids_(atomids), component_(component)
         {
         }

Member Function Documentation

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

inlineoverride

Evaluate the CV.

Parameters

snapshot Current simulation snapshot.

Definition at line 104 of file GyrationTensorCV.h.

References SSAGES::Snapshot::ApplyMinimumImage(), atomids_, SSAGES::Snapshot::CenterOfMass(), component_, SSAGES::Snapshot::GetCommunicator(), SSAGES::Snapshot::GetLocalIndices(), SSAGES::Snapshot::GetMasses(), SSAGES::Snapshot::GetNumAtoms(), SSAGES::Snapshot::GetPositions(), SSAGES::CollectiveVariable::grad_, SSAGES::Snapshot::TotalMass(), and SSAGES::CollectiveVariable::val_.

         {
             using namespace Eigen; 
 
             // Get local atom indices and compute COM. 
             std::vector<int> idx;
             snapshot.GetLocalIndices(atomids_, &idx);
 
             // Get data from snapshot. 
             auto n = snapshot.GetNumAtoms();
             const auto& masses = snapshot.GetMasses();
             const auto& pos = snapshot.GetPositions();
 
             // Initialize gradient.
             std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
             grad_.resize(n, Vector3{0,0,0});
             
             // Compute total and center of mass.
             auto masstot = snapshot.TotalMass(idx);
             Vector3 com = snapshot.CenterOfMass(idx, masstot);
 
             // Gyration tensor and temporary vector to store positions in inertial frame. 
             Matrix3 S = Matrix3::Zero();
             std::vector<Vector3> ris; 
             ris.reserve(idx.size());
             for(auto& i : idx)
             {
                 ris.emplace_back(snapshot.ApplyMinimumImage(pos[i] - com));
                 S.noalias() += masses[i]*ris.back()*ris.back().transpose();
             }
 
             // Reduce gyration tensor across processors and normalize.
             MPI_Allreduce(MPI_IN_PLACE, S.data(), S.size(), MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());
             S /= masstot;
 
             // Perform EVD. The columns are the eigenvectors. 
             // SelfAdjoint solver sorts in ascending order. 
             SelfAdjointEigenSolver<Matrix3> solver;
             solver.computeDirect(S);
             const auto& eigvals = solver.eigenvalues();
             const auto& eigvecs = solver.eigenvectors();
             
             // Assign variables for clarity. l1 is largest.
             auto l1 = eigvals[2], 
                  l2 = eigvals[1], 
                  l3 = eigvals[0];
             const auto& n1 = eigvecs.col(2), 
                         n2 = eigvecs.col(1), 
                         n3 = eigvecs.col(0);
 
             // Compute gradient.
             size_t j = 0;
             val_ = 0;
             for(auto& i : idx)
             {
                 // Compute derivative of each eigenvalue and use combos in components. 
                 auto dl1 = 2.*masses[i]/masstot*(1.-masses[i]/masstot)*ris[j].dot(n1)*n1;
                 auto dl2 = 2.*masses[i]/masstot*(1.-masses[i]/masstot)*ris[j].dot(n2)*n2;
                 auto dl3 = 2.*masses[i]/masstot*(1.-masses[i]/masstot)*ris[j].dot(n3)*n3;
 
                 // It is inefficient to keep reassigning val, but it's better than having two 
                 // switches, and the compiler will probably optimize it out anyways.
                 switch(component_)
                 {
                     case Rg:
                         grad_[i] = dl1 + dl2 + dl3;
                         val_ = l1 + l2 + l3; 
                         break;
                     case principal1:
                         grad_[i] = dl1;
                         val_ = l1;
                         break;
                     case principal2:
                         grad_[i] = dl2;
                         val_ = l2;
                         break;
                     case principal3:
                         grad_[i] = dl3;
                         val_ = l3;
                         break;
                     case asphericity:
                         grad_[i] = dl1 - 0.5*(dl2 + dl3);
                         val_ = l1 - 0.5*(l2 + l3);
                         break;
                     case acylindricity:
                         grad_[i] = dl2 - dl3;
                         val_ = l2 - l3; 
                         break;
                     case shapeaniso:
                         auto l1_2 = l1*l1, l2_2 = l2*l2, l3_2 = l3*l3; 
                         auto sum = l1 + l2 + l3;
                         auto sqsum = l1_2 + l2_2 + l3_2;
                         grad_[i] = 3.*(l1*dl1+l2*dl2+l3*dl3)/(sum*sum) - 
                                    3.*sqsum*(dl1+dl2+dl3)/(sum*sum*sum);
                         val_ = 1.5*sqsum/(sum*sum) - 0.5;
                         break;
                 }
 
                 ++j;
             }
         }

Here is the call graph for this function:

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

inlineoverride

Initialize necessary variables.

Parameters

snapshot Current simulation snapshot.

Definition at line 75 of file GyrationTensorCV.h.

References atomids_, SSAGES::Snapshot::GetCommunicator(), and SSAGES::Snapshot::GetLocalIndex().

         {
             using std::to_string;
 
             auto n = atomids_.size();
 
             // Make sure atom ID's are on at least one processor. 
             std::vector<int> found(n, 0);
             for(size_t i = 0; i < n; ++i)
             {
                 if(snapshot.GetLocalIndex(atomids_[i]) != -1)
                     found[i] = 1;
             }
 
             MPI_Allreduce(MPI_IN_PLACE, found.data(), n, MPI_INT, MPI_SUM, snapshot.GetCommunicator());
             unsigned ntot = std::accumulate(found.begin(), found.end(), 0, std::plus<int>());
             if(ntot != n)
                 throw BuildException({
                     "GyrationTensorCV: Expected to find " + 
                     to_string(n) + 
                     " atoms, but only found " + 
                     to_string(ntot) + "."
                 });     
         }       

Here is the call graph for this function:

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

/home/cody/SSAGES/src/CVs/GyrationTensorCV.h

Public Member Functions

Static Public Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation