Public Member Functions | Public Attributes | Private Attributes | List of all members
PLMD::Kearsley Class Reference

#include <Kearsley.h>

Public Member Functions

 Kearsley (const std::vector< Vector > &p0, const std::vector< Vector > &p1, const std::vector< double > &align, Log *&log)
 constructor: need the two structure, the alignment vector and the log reference More...
 
void assignP0 (const std::vector< Vector > &p0)
 switch the assignment of the structure p0 (e.g. at each md step) More...
 
void assignP1 (const std::vector< Vector > &p1)
 derivatives: derivative of the error respect p1 More...
 
void assignAlign (const std::vector< double > &align)
 transfer the alignment vector More...
 
void finiteDifferenceInterface (bool rmsd)
 finite differences of all the relevant quantities: takes a bool which decides if giving back rmsd or not (msd in this case) More...
 
double calculate (bool rmsd)
 

Public Attributes

double err
 error: the distance between two frames (might be rmsd/msd. See below) More...
 
std::vector< Vectordiff0on1
 displacement: the vector that goes from the p0 onto p1 More...
 
std::vector< Vectordiff1on0
 displacement: the vector that goes from the p1 onto p0 (via inverse rotation) More...
 
Vector com0
 center of mass of p0 More...
 
Vector com1
 center of mass of p1 More...
 
std::vector< Vectorp0reset
 position resetted wrt coms p0 More...
 
std::vector< Vectorp1reset
 position resetted wrt coms p1 More...
 
std::vector< Vectorp0rotated
 position rotated: p0 More...
 
std::vector< Vectorp1rotated
 position rotated: p1 More...
 
Tensor rotmat0on1
 rotation matrices p0 on p1 and reverse (p1 over p0) More...
 
Tensor rotmat1on0
 
std::vector< Vectorderrdp0
 derivatives: derivative of the error respect p0 More...
 
std::vector< Vectorderrdp1
 derivatives: derivative of the error respect p1 More...
 
std::vector< double > dmatdp0
 derivative of the rotation matrix note the dimension 3x3 x 3 x N More...
 
std::vector< double > dmatdp1
 

Private Attributes

Loglog
 general log reference that needs to be initialized when constructed More...
 
std::vector< Vectorp0
 position of atoms (first frame. In md is the running frame) More...
 
std::vector< Vectorp1
 position of atoms (second frame. In md is the reference frame) More...
 
std::vector< double > align
 alignment weight: the rmsd/msd that it provides is only based on this scalar More...
 
bool com0_is_removed
 
bool com1_is_removed
 

Detailed Description

A class that implements Kearsley's calculation which is optimal alignment via quaternion and analytical derivatives via perturbation theory

In Kearsley algorithm (see, S. K. Kearsley, Acta Crystallogr., Sect. A: Found. Crystallogr. 45, 208 1989 ), here adapted to included the COM reset, one first calculates the COM reset position respect to frame \( \{x_0,y_0,z_0\} \) (running frame), being \( w_{tot}^{al}=\sum_i w_{al}^i \)

\begin{eqnarray} \tilde x_0^i=x_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} x_0^i\\ \tilde y_0^i=y_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} y_0^i\\ \tilde z_0^i=z_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} z_0^i \end{eqnarray}

and the same is done with the reference one \( \{x_1,y_1,z_1\} \)

\begin{eqnarray} \tilde x_1^i=x_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} x_1^i\\ \tilde y_1^i=y_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} y_1^i\\ \tilde z_1^i=z_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} z_1^i \end{eqnarray}

Then one can build the \( \{x_p,y_p,z_p\} \) and \( \{x_m,y_m,z_m\} \) vectors of weighted summation and differences:

\begin{eqnarray} x_m^i=w_{al}^i(\tilde x_0^i-\tilde x_1^i)\\ y_m^i=w_{al}^i(\tilde y_0^i-\tilde y_1^i)\\ z_m^i=w_{al}^i(\tilde z_0^i-\tilde z_1^i) \end{eqnarray}

\begin{eqnarray} x_p^i=w_{al}^i(x_0^i+x_1^i)\\ y_p^i=w_{al}^i(y_0^i+y_1^i)\\ z_p^i=w_{al}^i(z_0^i+z_1^i) \end{eqnarray}

Then one build the COM-resetted matrix

\begin{displaymath} \mathbf{M}=\left[ \begin{array}{cccc} \sum ( {x}_{m}^{2}+{y}_{m}^{2}+{z}_{m}^{2}) & \sum (y_{p}z_{m} -y_{m}z_{p}) & \sum ( x_{m}z_{p} -x_{p}z_{m}) & \sum (x_{p}y_{m}-x_{m}y_{p} ) \\ \sum ( y_{p}z_{m} -y_{m}z_{p}) & \sum ( {x}_{m}^{2}+{y}_{p}^{2}+{z}_{p}^{2}) & \sum ( x_{m}y_{m} -x_{p}y_{p} ) & \sum (x_{m}z_{m}-x_{p}z_{p} ) \\ \sum (x_m z_p - x_p z_m ) & \sum ( x_m y_m -x_p y_p) & \sum ( {x}_{p}^{2}+{y}_{m}^{2}+{z}_{p}^{2}) & \sum ( y_m z_m -y_p z_p) \\ \sum (x_p y_m -x_m y_p ) & \sum (x_m z_m - x_p z_p ) & \sum (y_m z_m- y_p z_p ) & \sum ( {x}_{p}^{2}+{y}_{p}^{2}+{z}_{m}^{2}) \\ \end{array} \right] \end{displaymath}

by diagonalizing one obtains the mean square deviation by using the lowest eigenvalue \( \lambda_0 \)

\begin{equation} MSD= \frac{\lambda_{0}}{w_{tot}^{al}} \end{equation}

The rotation matrix is obtained from the eigenvector corresponding to \( \lambda_0 \) eigenvalue having components \( q_1, q_2, q_3, q_4 \)

\begin{displaymath} \mathbf{R}=\left[ \begin{array}{ccc} q_1 ^2 + q_2 ^2 - q_3 ^2 - q_4^2 & 2(q_2 q_3 + q_1 q_4) & 2(q_1 q_4 -q_1 q_3 )\\ 2(q_2 q_3 - q_1 q_4) & q_1 ^2 +q_3 ^2 -q_2 ^2 -q_4^2 & 2(q_3 q_4 - q_1 q_2)\\ 2( q_2 q_4 + q_1 q_3) & 2( q_3 q_4 - q_1 q_2) & q_1^2 +q_4 ^2 - q_2^2 - q_3 ^2 \\ \end{array} \right] \end{displaymath}

by using the perturbation theory one can retrieve the various derivatives:

In derivative calculation we exploited the classical Perturbation Theory up to the first order. In extensive manner, we introduce a perturbation over \(\lambda_{0}\) correlated with a pertubation of the states \(\vert q_{0}\rangle \) (in bra-ket notation):

\begin{displaymath} [\mathbf{M}+d\mathbf{M}][\vert q_{0}\rangle + \vert dq_{0}\rangle ]= [\lambda_{0}+d\lambda_{0}][\vert q_{0}\rangle +\vert dq_{0}\rangle ] \end{displaymath}

Grouping the zero order we recollect the unperturbed equation(see before). To the first order:

\begin{displaymath} d\mathbf{M}q_{0}+\mathbf{M}\vert dq_{0}\rangle =d\lambda_{0}\vert q_{0}\rangle +\lambda_{0} \vert dq_{0}\rangle \end{displaymath}

Now we express \(dq_{0}\) as linear combination of the other ortogonal eigenvectors:

\begin{displaymath} \vert dq_{0}\rangle =\sum_{j\neq0}c_{j}\vert q_{j}\rangle \end{displaymath}

thus we have

\begin{displaymath} d\mathbf{M}\vert q_{0}\rangle +\sum_{j\neq0}c_{j}\mathbf{M}\vert q_{j}\rangle= d\lambda_{0}\vert q_{0}\rangle+\lambda_{0}\sum_{j\neq0}c_{j}\vert q_{j}\rangle \end{displaymath}

projecting onto the \(q_{0}\) state and deleting the projection onto \(\vert dq_{0}\rangle\) beacuse of ortogonality:

\begin{displaymath} \langle q_{0}\vert d\mathbf{M}\vert q_{0}\rangle +\sum_{j\neq0}c_{j}\lambda_{j}\langle q_{0} \vert q_{j}\rangle= d\lambda_{0}\langle q_{0}\vert q_{0}\rangle+\lambda_{0}\sum_{j\neq0}c_{j}\langle q_{0}\vert q_{j}\rangle \end{displaymath}

we get

\begin{displaymath} \langle q_{0}\vert d\mathbf{M}\vert q_{0}\rangle=d\lambda_{0} \end{displaymath}

So, using simple chain rules:

\begin{displaymath} \langle q_{0}\vert \frac{d\mathbf{M}}{dr_{k}^{\gamma}}\vert q_{0}\rangle dr_{k}^{\gamma}=d\lambda_{0} \end{displaymath}

where here we used the notation \(r_{k}^{\gamma}\) to denote an arbitrary position which can be \(\tilde x_0 ,\tilde y_0,\tilde z_0\) or \(\tilde x_1 ,\tilde y_1,\tilde z_1\) we get

\begin{displaymath} \langle q_{0}\vert \frac{d\mathbf{M}}{dr_{k}^{\gamma}}\vert q_{0}\rangle =\frac{d\lambda_{0}}{dr_{k}^{\gamma}} \end{displaymath}

The derivatives of the matrix \(\frac{d\mathbf{M}}{dr_{k}^{\gamma}} \) can be readily obtained via the chain rule

\begin{displaymath} \frac{d\mathbf{M}}{dr_{k}^{\gamma}}=\sum_{\l}^{nat}\sum_{\alpha}^{x,y,z} \frac{d\mathbf{M}}{dP_{l}^{\alpha}}\frac{dP_{l}^{\alpha}}{dr_{k}^{\gamma}} +\\ \frac{d\mathbf{M}}{dM_{l}^{\alpha}}\frac{dM_{l}^{\alpha}}{dr_{k}^{\gamma}} \end{displaymath}

where \( M_{l}^{\alpha} \) corresponds to \( x_m^{l},y_m^{l},z_m^{l} \) and \( P_{l}^{\alpha} \) corresponds to \( x_p^{l},y_p^{l},z_p^{l} \) according to the \( \alpha \) component.

Constructor & Destructor Documentation

◆ Kearsley()

PLMD::Kearsley::Kearsley ( const std::vector< Vector > &  p0,
const std::vector< Vector > &  p1,
const std::vector< double > &  align,
Log *&  log 
)

constructor: need the two structure, the alignment vector and the log reference

Member Function Documentation

◆ assignAlign()

void PLMD::Kearsley::assignAlign ( const std::vector< double > &  align)

transfer the alignment vector

◆ assignP0()

void PLMD::Kearsley::assignP0 ( const std::vector< Vector > &  p0)

switch the assignment of the structure p0 (e.g. at each md step)

◆ assignP1()

void PLMD::Kearsley::assignP1 ( const std::vector< Vector > &  p1)

derivatives: derivative of the error respect p1

◆ calculate()

double PLMD::Kearsley::calculate ( bool  rmsd)

◆ finiteDifferenceInterface()

void PLMD::Kearsley::finiteDifferenceInterface ( bool  rmsd)

finite differences of all the relevant quantities: takes a bool which decides if giving back rmsd or not (msd in this case)

Member Data Documentation

◆ align

std::vector<double> PLMD::Kearsley::align
private

alignment weight: the rmsd/msd that it provides is only based on this scalar

◆ com0

Vector PLMD::Kearsley::com0

center of mass of p0

◆ com0_is_removed

bool PLMD::Kearsley::com0_is_removed
private

◆ com1

Vector PLMD::Kearsley::com1

center of mass of p1

◆ com1_is_removed

bool PLMD::Kearsley::com1_is_removed
private

◆ derrdp0

std::vector<Vector> PLMD::Kearsley::derrdp0

derivatives: derivative of the error respect p0

◆ derrdp1

std::vector<Vector> PLMD::Kearsley::derrdp1

derivatives: derivative of the error respect p1

◆ diff0on1

std::vector<Vector> PLMD::Kearsley::diff0on1

displacement: the vector that goes from the p0 onto p1

◆ diff1on0

std::vector<Vector> PLMD::Kearsley::diff1on0

displacement: the vector that goes from the p1 onto p0 (via inverse rotation)

◆ dmatdp0

std::vector<double> PLMD::Kearsley::dmatdp0

derivative of the rotation matrix note the dimension 3x3 x 3 x N

◆ dmatdp1

std::vector<double> PLMD::Kearsley::dmatdp1

◆ err

double PLMD::Kearsley::err

error: the distance between two frames (might be rmsd/msd. See below)

◆ log

Log* PLMD::Kearsley::log
private

general log reference that needs to be initialized when constructed

◆ p0

std::vector<Vector> PLMD::Kearsley::p0
private

position of atoms (first frame. In md is the running frame)

◆ p0reset

std::vector<Vector> PLMD::Kearsley::p0reset

position resetted wrt coms p0

◆ p0rotated

std::vector<Vector> PLMD::Kearsley::p0rotated

position rotated: p0

◆ p1

std::vector<Vector> PLMD::Kearsley::p1
private

position of atoms (second frame. In md is the reference frame)

◆ p1reset

std::vector<Vector> PLMD::Kearsley::p1reset

position resetted wrt coms p1

◆ p1rotated

std::vector<Vector> PLMD::Kearsley::p1rotated

position rotated: p1

◆ rotmat0on1

Tensor PLMD::Kearsley::rotmat0on1

rotation matrices p0 on p1 and reverse (p1 over p0)

◆ rotmat1on0

Tensor PLMD::Kearsley::rotmat1on0

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