SSAGES::RMSDCV Class Reference

Collective variable to calculate root mean square displacement. More...

#include <RMSDCV.h>

Inheritance diagram for SSAGES::RMSDCV:

Public Member Functions
	RMSDCV (const Label &atomids, std::string xyzfile, bool use_range=false)
	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.
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 RMSDCV *	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
Label	atomids_
	IDs of the atoms used for RMSD calculation.
std::string	xyzfile_
	Name of model structure.
std::vector< Vector3 >	refcoord_
	Store reference structure coordinates.
Vector3	COMref_
	Center of mass of reference.

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 to calculate root mean square displacement.

RMSD calculation performed as reported in Coutsias, E. A., Seok, C., and Dill, K. A., "Using Quaternions to Calculate RMSD", J. Comput. Chem. 25: 1849-1857, 2004

Definition at line 42 of file RMSDCV.h.

Constructor & Destructor Documentation

RMSDCV()

SSAGES::RMSDCV::RMSDCV ( const Label &  atomids, std::string  xyzfile, bool  use_range = false  )

inline

Constructor.

Parameters

atomids	IDs of the atoms defining RMSD.
xyzfile	String determining the reference file.
use_range	If True Use range of atoms defined by the two atoms in atomids.

Construct a RMSD CV.

Todo:

Bounds needs to be an input and periodic boundary conditions

Definition at line 68 of file RMSDCV.h.

                                                                                :
        atomids_(atomids), xyzfile_(xyzfile)
        {
            if(use_range)
            {
                if(atomids.size() != 2)
                {
                    throw BuildException({
                        "RMSDCV: With use_range, expected 2 atoms, but found " +
                        std::to_string(atomids.size()) + "."
                    });
                }
 
                atomids_.clear();
                if(atomids[0] >= atomids[1])
                {
                    throw BuildException({
                        "RMSDCV: Atom range must be strictly increasing."
                    });
                }
                for(int i = atomids[0]; i <= atomids[1]; ++i)
                {
                    atomids_.push_back(i);
                }
            }
        }

References atomids_.

Referenced by Build().

Here is the caller graph for this function:

Member Function Documentation

Build()

static RMSDCV* SSAGES::RMSDCV::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 244 of file RMSDCV.h.

        {
            Json::ObjectRequirement validator;
            Json::Value schema;
            Json::CharReaderBuilder rbuilder;
            Json::CharReader* reader = rbuilder.newCharReader();
 
            reader->parse(JsonSchema::RMSDCV.c_str(),
                          JsonSchema::RMSDCV.c_str() + JsonSchema::RMSDCV.size(),
                          &schema, NULL);
            validator.Parse(schema, path);
 
            //Validate inputs
            validator.Validate(json, path);
            if(validator.HasErrors())
                    throw BuildException(validator.GetErrors());
 
            std::vector<int> atom_ids;
            for(auto& s : json["atom_ids"])
                atom_ids.push_back(s.asInt());
            auto reference = json["reference"].asString();
            auto use_range = json.get("use_range",false).asBool();
 
            return new RMSDCV(atom_ids, reference, use_range);
        }

References Json::Requirement::GetErrors(), Json::Requirement::HasErrors(), Json::ObjectRequirement::Parse(), RMSDCV(), 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::RMSDCV::Evaluate ( const Snapshot &  snapshot )

inlineoverride

Evaluate the CV.

Parameters

snapshot

Current simulation snapshot.

Definition at line 162 of file RMSDCV.h.

        {
            // Get local atom indices.
            Label idx;
            snapshot.GetLocalIndices(atomids_, &idx);
 
            // Get data from snapshot.
            const auto& masses = snapshot.GetMasses();
            const auto& pos = snapshot.GetPositions();
            double masstot = snapshot.TotalMass(idx);
 
            // Initialize gradient.
            std::fill(grad_.begin(), grad_.end(), Vector3{0,0,0});
            grad_.resize(snapshot.GetNumAtoms(), Vector3{0,0,0});
 
            // Calculate center of mass
            Vector3 COM_ = snapshot.CenterOfMass(idx);
 
            // Build correlation matrix
            Matrix3 R = Matrix3::Zero();
 
            Vector3 diff, diff_ref;
            double part_RMSD = 0;
 
            for(auto& i : idx)
            {
                diff = snapshot.ApplyMinimumImage(pos[i] - COM_);
                diff_ref = refcoord_[i] - COMref_; // Reference has no "box" or minimum image
 
                R.noalias() += masses[i]*diff*diff_ref.transpose();
 
                part_RMSD += masses[i]*(diff.squaredNorm() + diff_ref.squaredNorm());
            }
            MPI_Allreduce(MPI_IN_PLACE, R.data(), R.size(), MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());
            MPI_Allreduce(MPI_IN_PLACE, &part_RMSD, 1, MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());
            R /= masstot;
            part_RMSD /= masstot;
 
            Eigen::Matrix4d F;
            F(0,0) = R(0,0) + R(1,1) + R(2,2);
            F(1,0) = R(1,2) - R(2,1);
            F(2,0) = R(2,0) - R(0,2);
            F(3,0) = R(0,1) - R(1,0);
 
            F(0,1) = F(1,0);
            F(1,1) = R(0,0) - R(1,1) - R(2,2);
            F(2,1) = R(0,1) + R(1,0);
            F(3,1) = R(0,2) + R(2,0);
 
            F(0,2) = F(2,0);
            F(1,2) = F(2,1);
            F(2,2) = -R(0,0) + R(1,1) - R(2,2);
            F(3,2) = R(1,2) + R(2,1);
 
            F(0,3) = F(3,0);
            F(1,3) = F(3,1);
            F(2,3) = F(3,2);
            F(3,3) = -R(0,0) - R(1,1) + R(2,2);
 
            //Find eigenvalues
            Eigen::SelfAdjointEigenSolver<Eigen::Matrix4d> es(F);
            auto max_lambda = es.eigenvalues().real()[3];
 
            // Calculate RMSD
            val_ = sqrt(part_RMSD - (2*max_lambda));
 
            // If RMSD is zero, we have nothing to do.
            if(val_ == 0) return;
 
            // Calculate gradient
            auto ev = es.eigenvectors().real().col(3);
            auto q = Eigen::Quaterniond(ev(0),ev(1),ev(2),ev(3));
            Matrix3 RotMatrix = q.toRotationMatrix();
 
            for(auto& i : idx)
            {
                auto rotref = RotMatrix.transpose()*(refcoord_[i] - COMref_);
                grad_[i] = masses[i]/masstot*(1-masses[i]/masstot)*(snapshot.ApplyMinimumImage(pos[i] - COM_) - rotref)/val_;
            }
        }

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

Here is the call graph for this function:

Initialize()

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

inlineoverride

Initialize necessary variables.

Parameters

snapshot

Current simulation snapshot.

Definition at line 99 of file RMSDCV.h.

        {
            const auto& masses = snapshot.GetMasses();
            auto n = atomids_.size();
 
            // Make sure atom IDs 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({
                    "RMSDCV: Expected to find " + std::to_string(n) +
                    " atoms, but only found " + std::to_string(ntot) + "."
                });
            }
 
            Label idx;
            snapshot.GetLocalIndices(atomids_, &idx);
 
            std::vector<std::array<double,4>> xyzinfo = ReadFile::ReadXYZ(xyzfile_);
 
            refcoord_.resize(snapshot.GetNumAtoms(), Vector3{0,0,0});
 
            // Loop through atom positions
            for(size_t i = 0; i < xyzinfo.size(); ++i)
            {
                int id = snapshot.GetLocalIndex(xyzinfo[i][0]);
                if(id == -1) continue; // Atom not found
                for(size_t j = 0; j < atomids_.size(); ++j)
                {
                    if(atomids_[j] == xyzinfo[i][0])
                    {
                        refcoord_[id][0] = xyzinfo[i][1];
                        refcoord_[id][1] = xyzinfo[i][2];
                        refcoord_[id][2] = xyzinfo[i][3];
                    }
                }
            }
 
            Vector3 mass_pos_prod_ref = Vector3::Zero();
            double total_mass = 0;
 
            for(auto& i : idx)
            {
                mass_pos_prod_ref += masses[i]*refcoord_[i];
                total_mass += masses[i];
            }
            MPI_Allreduce(MPI_IN_PLACE, mass_pos_prod_ref.data(), mass_pos_prod_ref.size(), MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());
            MPI_Allreduce(MPI_IN_PLACE, &total_mass, 1, MPI_DOUBLE, MPI_SUM, snapshot.GetCommunicator());
 
            COMref_ = mass_pos_prod_ref/total_mass;
        }

References atomids_, COMref_, SSAGES::Snapshot::GetCommunicator(), SSAGES::Snapshot::GetLocalIndex(), SSAGES::Snapshot::GetLocalIndices(), SSAGES::Snapshot::GetMasses(), SSAGES::Snapshot::GetNumAtoms(), SSAGES::ReadFile::ReadXYZ(), refcoord_, and xyzfile_.

Here is the call graph for this function:

---

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

• /home/michael/development/SSAGES/src/CVs/RMSDCV.h