SSAGES  0.8.5
Software Suite for Advanced General Ensemble Simulations
RMSDCV.h
1 
22 #pragma once
23 
24 #include "Snapshot.h"
25 #include "CollectiveVariable.h"
26 #include "../Utility/ReadFile.h"
27 #include <array>
28 #include <Eigen/Dense>
29 #include <Eigen/Eigenvalues>
30 
31 namespace SSAGES
32 {
34 
39  class RMSDCV : public CollectiveVariable
40  {
41  private:
43  std::vector<int> atomids_;
44 
46  std::vector<int> pertatoms_;
47 
49  std::string molecule_;
50 
52  std::vector<Vector3> refcoord_;
53 
56 
57  public:
59 
68  RMSDCV(std::vector<int> atomids, std::string molxyz, bool use_range = false) :
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  }
88 
90 
93  void Initialize(const Snapshot& snapshot) override
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  }
137 
139 
142  void Evaluate(const Snapshot& snapshot) override
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  }
302  };
303 }
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
const std::vector< Integer3 > & GetImageFlags() const
Access the particles image flags.
Definition: Snapshot.h:340
std::string molecule_
Name of model structure.
Definition: RMSDCV.h:49
Class containing a snapshot of the current simulation in time.
Definition: Snapshot.h:47
RMSDCV(std::vector< int > atomids, std::string molxyz, bool use_range=false)
Constructor.
Definition: RMSDCV.h:68
std::vector< Vector3 > grad_
Gradient vector dCv/dxi.
const Label & GetAtomIDs() const
Access the atom IDs.
Definition: Snapshot.h:602
void Evaluate(const Snapshot &snapshot) override
Evaluate the CV.
Definition: RMSDCV.h:142
const std::vector< double > & GetMasses() const
Const access to the particle masses.
Definition: Snapshot.h:382
void Initialize(const Snapshot &snapshot) override
Initialize necessary variables.
Definition: RMSDCV.h:93
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
Abstract class for a collective variable.
static std::vector< std::array< double, 4 > > ReadXYZ(std::string FileName)
Read xyz file.
Definition: ReadFile.h:70
Vector3 UnwrapVector(const Vector3 &v, const Integer3 &image) const
Unwrap a vector&#39;s real coordinates according to its image replica count.
Definition: Snapshot.h:413
std::vector< Vector3 > refcoord_
Store reference structure coordinates.
Definition: RMSDCV.h:52
const std::vector< Vector3 > & GetPositions() const
Access the particle positions.
Definition: Snapshot.h:327
Collective variable to calculate root mean square displacement.
Definition: RMSDCV.h:39