SSAGES  0.8.3
Software Suite for Advanced General Ensemble Simulations
Public Member Functions | Private Attributes | List of all members
SSAGES::RMSDCV Class Reference

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

#include <RMSDCV.h>

Inheritance diagram for SSAGES::RMSDCV:
Inheritance graph
[legend]

Public Member Functions

 RMSDCV (std::vector< int > atomids, std::string molxyz, 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 Matrix3GetBoxGradient () const
 Get gradient contribution to box.
 
const std::array< double, 2 > & GetBoundaries ()
 Get CV boundaries. More...
 
virtual double GetDifference (double location) const
 

Private Attributes

std::vector< int > atomids_
 IDs of the atoms used for Rg calculation.
 
std::vector< int > pertatoms_
 Array to store indices of atoms of interest.
 
std::string molecule_
 Name of model structure.
 
std::vector< Vector3refcoord_
 Store reference structure coordinates.
 
Vector3 COMref_
 Center of mass of reference.
 

Additional Inherited Members

- Static Public Member Functions inherited from SSAGES::CollectiveVariable
static CollectiveVariableBuildCV (const Json::Value &json, const std::string &path)
 Set up collective variable. More...
 
- Protected Attributes inherited from SSAGES::CollectiveVariable
std::vector< Vector3grad_
 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 39 of file RMSDCV.h.

Constructor & Destructor Documentation

◆ RMSDCV()

SSAGES::RMSDCV::RMSDCV ( std::vector< int >  atomids,
std::string  molxyz,
bool  use_range = false 
)
inline

Constructor.

Parameters
atomidsIDs of the atoms defining Rg.
molxyzString determining the molecule.
use_rangeIf 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.

68  :
69  atomids_(atomids), molecule_(molxyz)
70  {
71  if(use_range)
72  {
73  if(atomids.size() != 2)
74  { std::cout<<"RMSDCV: If using range, must define only two atoms!"<<std::endl;
75  exit(0);
76  }
77 
78  atomids_.clear();
79 
80  if(atomids[0] >= atomids[1])
81  { std::cout<<"RMSDCV: Please reverse atom range or check that atom range is not equal!"<<std::endl;
82  exit(0);
83  }
84  for(int i = atomids[0]; i <= atomids[1];i++)
85  atomids_.push_back(i);
86  }
87  }
std::string molecule_
Name of model structure.
Definition: RMSDCV.h:49
std::vector< int > atomids_
IDs of the atoms used for Rg calculation.
Definition: RMSDCV.h:43

Member Function Documentation

◆ Evaluate()

void SSAGES::RMSDCV::Evaluate ( const Snapshot snapshot)
inlineoverride

Evaluate the CV.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 142 of file RMSDCV.h.

References COMref_, SSAGES::Snapshot::GetAtomIDs(), SSAGES::Snapshot::GetImageFlags(), SSAGES::Snapshot::GetMasses(), SSAGES::Snapshot::GetPositions(), SSAGES::CollectiveVariable::grad_, SSAGES::Snapshot::UnwrapVector(), and SSAGES::CollectiveVariable::val_.

143  {
144  const auto& pos = snapshot.GetPositions();
145  const auto& ids = snapshot.GetAtomIDs();
146  const auto& mass = snapshot.GetMasses();
147  const auto& image_flags = snapshot.GetImageFlags();
148  Vector3 mass_pos_prod = {0,0,0};
149  double total_mass=0;
150  int i;
151  // Loop through atom positions
152  for( size_t i = 0; i < pos.size(); ++i)
153  {
154  grad_[i].setZero();
155  // Loop through pertinent atoms
156  for(size_t j = 0; j < atomids_.size(); j++)
157  {
158  if(ids[i] == atomids_[j])
159  {
160  pertatoms_[j] = i;
161  auto u_coord = snapshot.UnwrapVector(pos[i], image_flags[i]);
162 
163  mass_pos_prod += mass[i]*u_coord;
164  total_mass += mass[i];
165  break;
166  }
167  }
168  }
169 
170 
171  // Calculate center of mass
172  //Vector3 mass_pos_prod;
173  //mass_pos_prod.setZero();
174  total_mass = 0;
175  Vector3 COM_;
176 
177  // Need to unwrap coordinates for appropriate COM calculation.
178 
179  //for( size_t j = 0; j < pertatoms_.size(); ++j)
180  //{
181  // i = pertatoms_[j];
182  // auto u_coord = UnwrapCoordinates(snapshot.GetLatticeConstants(), image_flags[i], pos[i]);
183 
184  // mass_pos_prod[0] += mass[i]*u_coord[0];
185  // mass_pos_prod[1] += mass[i]*u_coord[1];
186  // mass_pos_prod[2] += mass[i]*u_coord[2];
187  // total_mass += mass[i];
188  //}
189 
190  COM_ = mass_pos_prod/total_mass;
191 
192  // Build correlation matrix
193  Matrix3 R;
194 
195  Vector3 diff;
196  Vector3 diff_ref;
197  double part_RMSD = 0;
198 
199  for( size_t j = 0; j < pertatoms_.size(); ++j)
200  {
201  i = pertatoms_[j];
202  auto u_coord = snapshot.UnwrapVector(pos[i], image_flags[i]);
203 
204  diff = u_coord -COM_;
205  diff_ref = refcoord_[j] - COMref_;
206 
207  R(0,0) += diff[0]*diff_ref[0];
208  R(0,1) += diff[0]*diff_ref[1];
209  R(0,2) += diff[0]*diff_ref[2];
210  R(1,0) += diff[1]*diff_ref[0];
211  R(1,1) += diff[1]*diff_ref[1];
212  R(1,2) += diff[1]*diff_ref[2];
213  R(2,0) += diff[2]*diff_ref[0];
214  R(2,1) += diff[2]*diff_ref[1];
215  R(2,2) += diff[2]*diff_ref[2];
216 
217 
218  auto normdiff2 = diff.norm()*diff.norm();
219  auto normref2 = diff_ref.norm()*diff_ref.norm();
220 
221  part_RMSD += (normdiff2 + normref2);
222  }
223 
224  Eigen::Matrix4d F;
225  F(0,0)= R(0,0) + R(1,1) + R(2,2);
226  F(1,0)= R(1,2) - R(2,1);
227  F(2,0)= R(2,0) - R(0,2);
228  F(3,0)= R(0,1) - R(1,0);
229 
230  F(0,1)= R(1,2) - R(2,1);
231  F(1,1)= R(0,0) - R(1,1) - R(2,2);
232  F(2,1)= R(0,1) + R(1,0);
233  F(3,1)= R(0,2) + R(2,0);
234 
235  F(0,2)= R(2,0) - R(0,2);
236  F(1,2)= R(0,1) - R(1,0);
237  F(2,2)= -R(0,0) + R(1,1) - R(2,2);
238  F(3,2)= R(1,2) + R(2,1);
239 
240  F(0,3)= R(0,1) - R(1,0);
241  F(1,3)= R(0,2) - R(2,0);
242  F(2,3)= R(1,2) - R(2,1);
243  F(3,3)= -R(0,0) - R(1,1) + R(2,2);
244 
245  //Find eigenvalues
246  Eigen::EigenSolver<Eigen::Matrix4d> es(F);
247  //EigenSolver<F> es;
248  //es.compute(F);
249  //auto max_lambda = es.eigenvalues().maxCoeff();
250  auto lambda = es.eigenvalues().real();
251  //double max_lambda;
252  auto max_lambda = lambda.maxCoeff();
253  int column;
254 
255 
256  for (i =0; i < 4; ++i)
257  if(es.eigenvalues()[i] == max_lambda)
258  column=i;
259 
260  // Calculate RMSD
261 
262 
263  auto RMSD = sqrt((part_RMSD - (2*max_lambda))/(atomids_.size()));
264 
265  val_ = RMSD;
266 
267  // Calculate gradient
268 
269  auto eigenvector = es.eigenvectors().col(column);
270 
271  auto q0 = eigenvector[0].real();
272  auto q1 = eigenvector[1].real();
273  auto q2 = eigenvector[2].real();
274  auto q3 = eigenvector[3].real();
275 
276  Matrix3 RotMatrix(3,3);
277  RotMatrix(0,0) = q0*q0+q1*q1-q2*q2-q3*q3;
278  RotMatrix(1,0) = 2*(q1*q2+q0*q3);
279  RotMatrix(2,0) = 2*(q1*q3-q0*q2);
280 
281  RotMatrix(0,1) = 2*(q1*q2-q0*q3);
282  RotMatrix(1,1) = q0*q0-q1*q1+q2*q2-q3*q3;
283  RotMatrix(2,1) = 1*(q2*q3+q0*q1);
284 
285  RotMatrix(0,2) = 2*(q1*q3+q0*q2);
286  RotMatrix(1,2) = 2*(q2*q3-q0*q1);
287  RotMatrix(2,2) = q0*q0-q1*q1-q2*q2+q3*q3;
288 
289  //double trans[3];
290  for(size_t j=0; j<pertatoms_.size();++j)
291  {
292  auto trans = RotMatrix.transpose()*refcoord_[j];
293 
294  i = pertatoms_[j];
295  auto u_coord = pos[i]; //UnwrapCoordinates(snapshot.GetLatticeConstants(), image_flags[i], pos[i]);
296 
297  grad_[i][0] = (u_coord[0] - trans[0])/((atomids_.size())*val_);
298  grad_[i][1] = (u_coord[1] - trans[1])/((atomids_.size())*val_);
299  grad_[i][2] = (u_coord[2] - trans[2])/((atomids_.size())*val_);
300  }
301  }
std::vector< int > pertatoms_
Array to store indices of atoms of interest.
Definition: RMSDCV.h:46
Eigen::Matrix3d Matrix3
3x3 matrix.
Definition: types.h:42
Vector3 COMref_
Center of mass of reference.
Definition: RMSDCV.h:55
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
double val_
Current value of CV.
std::vector< int > atomids_
IDs of the atoms used for Rg calculation.
Definition: RMSDCV.h:43
std::vector< Vector3 > refcoord_
Store reference structure coordinates.
Definition: RMSDCV.h:52
Here is the call graph for this function:

◆ Initialize()

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

Initialize necessary variables.

Parameters
snapshotCurrent simulation snapshot.

Definition at line 93 of file RMSDCV.h.

References SSAGES::Snapshot::GetAtomIDs(), SSAGES::Snapshot::GetMasses(), SSAGES::Snapshot::GetPositions(), SSAGES::CollectiveVariable::grad_, and SSAGES::ReadFile::ReadXYZ().

94  {
95 
96  const auto& mass = snapshot.GetMasses();
97  const auto& ids = snapshot.GetAtomIDs();
98 
99  // Initialize gradient
100  auto n = snapshot.GetPositions().size();
101  grad_.resize(n);
102  refcoord_.resize(atomids_.size());
103 
104  std::vector<std::array<double,4>> xyzinfo = ReadFile::ReadXYZ(molecule_);
105 
106  if(refcoord_.size() != xyzinfo.size())
107  throw std::runtime_error("Reference structure and input structure atom size do not match!");
108 
109  for(size_t i = 0; i < xyzinfo.size(); i++)
110  {
111  refcoord_[i][0] = xyzinfo[i][1];
112  refcoord_[i][1] = xyzinfo[i][2];
113  refcoord_[i][2] = xyzinfo[i][3];
114  }
115 
116  Vector3 mass_pos_prod_ref;
117  mass_pos_prod_ref.setZero();
118  double total_mass = 0;
119 
120  // Loop through atom positions
121  for( size_t i = 0; i < mass.size(); ++i)
122  {
123  // Loop through pertinent atoms
124  for(size_t j = 0; j < atomids_.size(); j++)
125  {
126  if(ids[i] == atomids_[j])
127  {
128  mass_pos_prod_ref += mass[i]*refcoord_[j];
129  total_mass += mass[i];
130  break;
131  }
132  }
133  }
134 
135  COMref_ = mass_pos_prod_ref/total_mass;
136  }
Vector3 COMref_
Center of mass of reference.
Definition: RMSDCV.h:55
std::string molecule_
Name of model structure.
Definition: RMSDCV.h:49
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
Eigen::Vector3d Vector3
Three-dimensional vector.
Definition: types.h:33
std::vector< int > atomids_
IDs of the atoms used for Rg calculation.
Definition: RMSDCV.h:43
static std::vector< std::array< double, 4 > > ReadXYZ(std::string FileName)
Read xyz file.
Definition: ReadFile.h:70
std::vector< Vector3 > refcoord_
Store reference structure coordinates.
Definition: RMSDCV.h:52
Here is the call graph for this function:

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