 |
 |
 |
| |
|
|
#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< Vector > | diff0on1 |
| displacement: the vector that goes from the p0 onto p1 More... | |
| std::vector< Vector > | diff1on0 |
| 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< Vector > | p0reset |
| position resetted wrt coms p0 More... | |
| std::vector< Vector > | p1reset |
| position resetted wrt coms p1 More... | |
| std::vector< Vector > | p0rotated |
| position rotated: p0 More... | |
| std::vector< Vector > | p1rotated |
| position rotated: p1 More... | |
| Tensor | rotmat0on1 |
| rotation matrices p0 on p1 and reverse (p1 over p0) More... | |
| Tensor | rotmat1on0 |
| std::vector< Vector > | derrdp0 |
| derivatives: derivative of the error respect p0 More... | |
| std::vector< Vector > | derrdp1 |
| 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 | |
| Log * | log |
| general log reference that needs to be initialized when constructed More... | |
| std::vector< Vector > | p0 |
| position of atoms (first frame. In md is the running frame) More... | |
| std::vector< Vector > | p1 |
| 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 |
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.
| 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
| void PLMD::Kearsley::assignAlign | ( | const std::vector< double > & | align | ) |
transfer the alignment vector
| void PLMD::Kearsley::assignP0 | ( | const std::vector< Vector > & | p0 | ) |
switch the assignment of the structure p0 (e.g. at each md step)
| void PLMD::Kearsley::assignP1 | ( | const std::vector< Vector > & | p1 | ) |
derivatives: derivative of the error respect p1
| double PLMD::Kearsley::calculate | ( | bool | rmsd | ) |
| 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)
|
private |
alignment weight: the rmsd/msd that it provides is only based on this scalar
| Vector PLMD::Kearsley::com0 |
center of mass of p0
|
private |
| Vector PLMD::Kearsley::com1 |
center of mass of p1
|
private |
| std::vector<Vector> PLMD::Kearsley::derrdp0 |
derivatives: derivative of the error respect p0
| std::vector<Vector> PLMD::Kearsley::derrdp1 |
derivatives: derivative of the error respect p1
| std::vector<Vector> PLMD::Kearsley::diff0on1 |
displacement: the vector that goes from the p0 onto p1
| std::vector<Vector> PLMD::Kearsley::diff1on0 |
displacement: the vector that goes from the p1 onto p0 (via inverse rotation)
| std::vector<double> PLMD::Kearsley::dmatdp0 |
derivative of the rotation matrix note the dimension 3x3 x 3 x N
| std::vector<double> PLMD::Kearsley::dmatdp1 |
| double PLMD::Kearsley::err |
error: the distance between two frames (might be rmsd/msd. See below)
|
private |
general log reference that needs to be initialized when constructed
|
private |
position of atoms (first frame. In md is the running frame)
| std::vector<Vector> PLMD::Kearsley::p0reset |
position resetted wrt coms p0
| std::vector<Vector> PLMD::Kearsley::p0rotated |
position rotated: p0
|
private |
position of atoms (second frame. In md is the reference frame)
| std::vector<Vector> PLMD::Kearsley::p1reset |
position resetted wrt coms p1
| std::vector<Vector> PLMD::Kearsley::p1rotated |
position rotated: p1
| Tensor PLMD::Kearsley::rotmat0on1 |
rotation matrices p0 on p1 and reverse (p1 over p0)
| Tensor PLMD::Kearsley::rotmat1on0 |
Hosted by GitHub
|
1.8.14
|