Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #ifndef __PLUMED_tools_Kearsley_h
23 : #define __PLUMED_tools_Kearsley_h
24 :
25 : #include "Vector.h"
26 : #include "Tensor.h"
27 : #include <vector>
28 :
29 : namespace PLMD {
30 :
31 : class Log;
32 :
33 : /*! A class that implements Kearsley's calculation
34 : which is optimal alignment via quaternion and
35 : analytical derivatives via perturbation theory
36 :
37 : In Kearsley algorithm (see, S. K. Kearsley, Acta Crystallogr., Sect. A: Found. Crystallogr. 45, 208 1989 ),
38 : here adapted to included the COM reset,
39 : one first calculates the COM reset position respect to frame \f$ \{x_0,y_0,z_0\} \f$ (running frame),
40 : being \f$ w_{tot}^{al}=\sum_i w_{al}^i \f$
41 :
42 : \f{eqnarray}
43 : \tilde x_0^i=x_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} x_0^i\\
44 : \tilde y_0^i=y_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} y_0^i\\
45 : \tilde z_0^i=z_0^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} z_0^i
46 : \f}
47 :
48 : and the same is done with the reference one \f$ \{x_1,y_1,z_1\} \f$
49 : \f{eqnarray}
50 : \tilde x_1^i=x_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} x_1^i\\
51 : \tilde y_1^i=y_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} y_1^i\\
52 : \tilde z_1^i=z_1^i-\sum_i \frac{w_{al}^i}{w_{tot}^{al}} z_1^i
53 : \f}
54 :
55 : Then one can build the \f$ \{x_p,y_p,z_p\} \f$ and \f$ \{x_m,y_m,z_m\} \f$
56 : vectors of weighted summation and differences:
57 : \f{eqnarray}
58 : x_m^i=w_{al}^i(\tilde x_0^i-\tilde x_1^i)\\
59 : y_m^i=w_{al}^i(\tilde y_0^i-\tilde y_1^i)\\
60 : z_m^i=w_{al}^i(\tilde z_0^i-\tilde z_1^i)
61 : \f}
62 :
63 : \f{eqnarray}
64 : x_p^i=w_{al}^i(x_0^i+x_1^i)\\
65 : y_p^i=w_{al}^i(y_0^i+y_1^i)\\
66 : z_p^i=w_{al}^i(z_0^i+z_1^i)
67 : \f}
68 :
69 :
70 : Then one build the COM-resetted matrix
71 :
72 : \f{displaymath}
73 : \mathbf{M}=\left[
74 : \begin{array}{cccc}
75 : \sum ( {x}_{m}^{2}+{y}_{m}^{2}+{z}_{m}^{2}) &
76 : \sum (y_{p}z_{m} -y_{m}z_{p}) &
77 : \sum ( x_{m}z_{p} -x_{p}z_{m}) &
78 : \sum (x_{p}y_{m}-x_{m}y_{p} ) \\
79 : \sum ( y_{p}z_{m} -y_{m}z_{p}) &
80 : \sum ( {x}_{m}^{2}+{y}_{p}^{2}+{z}_{p}^{2}) &
81 : \sum ( x_{m}y_{m} -x_{p}y_{p} ) &
82 : \sum (x_{m}z_{m}-x_{p}z_{p} ) \\
83 : \sum (x_m z_p - x_p z_m ) &
84 : \sum ( x_m y_m -x_p y_p) &
85 : \sum ( {x}_{p}^{2}+{y}_{m}^{2}+{z}_{p}^{2}) &
86 : \sum ( y_m z_m -y_p z_p) \\
87 : \sum (x_p y_m -x_m y_p ) &
88 : \sum (x_m z_m - x_p z_p ) &
89 : \sum (y_m z_m- y_p z_p ) &
90 : \sum ( {x}_{p}^{2}+{y}_{p}^{2}+{z}_{m}^{2}) \\
91 : \end{array}
92 : \right]
93 : \f}
94 :
95 : by diagonalizing one obtains the mean square deviation by using the lowest eigenvalue \f$ \lambda_0 \f$
96 : \f{equation}
97 : MSD= \frac{\lambda_{0}}{w_{tot}^{al}}
98 : \f}
99 :
100 : The rotation matrix is obtained from the eigenvector corresponding to \f$ \lambda_0 \f$ eigenvalue
101 : having components \f$ q_1, q_2, q_3, q_4 \f$
102 :
103 : \f{displaymath}
104 : \mathbf{R}=\left[
105 : \begin{array}{ccc}
106 : q_1 ^2 + q_2 ^2 - q_3 ^2 - q_4^2 &
107 : 2(q_2 q_3 + q_1 q_4) &
108 : 2(q_1 q_4 -q_1 q_3 )\\
109 : 2(q_2 q_3 - q_1 q_4) &
110 : q_1 ^2 +q_3 ^2 -q_2 ^2 -q_4^2 &
111 : 2(q_3 q_4 - q_1 q_2)\\
112 : 2( q_2 q_4 + q_1 q_3) &
113 : 2( q_3 q_4 - q_1 q_2) &
114 : q_1^2 +q_4 ^2 - q_2^2 - q_3 ^2 \\
115 : \end{array}
116 : \right]
117 : \f}
118 :
119 : by using the perturbation theory one can retrieve the various derivatives:
120 :
121 : In derivative calculation we exploited the classical Perturbation Theory up to the first order.
122 : In extensive manner, we introduce a perturbation over \f$\lambda_{0}\f$ correlated with
123 : a pertubation of the states \f$\vert q_{0}\rangle \f$ (in bra-ket notation):
124 : \f{displaymath}
125 : [\mathbf{M}+d\mathbf{M}][\vert q_{0}\rangle + \vert dq_{0}\rangle ]=
126 : [\lambda_{0}+d\lambda_{0}][\vert q_{0}\rangle +\vert dq_{0}\rangle ]
127 : \f}
128 : Grouping the zero order we recollect the unperturbed equation(see before).
129 : To the first order:
130 : \f{displaymath}
131 : d\mathbf{M}q_{0}+\mathbf{M}\vert dq_{0}\rangle =d\lambda_{0}\vert q_{0}\rangle +\lambda_{0} \vert dq_{0}\rangle
132 : \f}
133 : Now we express \f$dq_{0}\f$ as linear combination of the other ortogonal eigenvectors:
134 : \f{displaymath}
135 : \vert dq_{0}\rangle =\sum_{j\neq0}c_{j}\vert q_{j}\rangle
136 : \f}
137 : thus we have
138 : \f{displaymath}
139 : d\mathbf{M}\vert q_{0}\rangle +\sum_{j\neq0}c_{j}\mathbf{M}\vert q_{j}\rangle=
140 : d\lambda_{0}\vert q_{0}\rangle+\lambda_{0}\sum_{j\neq0}c_{j}\vert q_{j}\rangle
141 : \f}
142 : projecting onto the \f$q_{0}\f$ state and deleting the projection onto \f$\vert dq_{0}\rangle\f$ beacuse
143 : of ortogonality:
144 : \f{displaymath}
145 : \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=
146 : d\lambda_{0}\langle q_{0}\vert q_{0}\rangle+\lambda_{0}\sum_{j\neq0}c_{j}\langle q_{0}\vert q_{j}\rangle
147 : \f}
148 : we get
149 : \f{displaymath}
150 : \langle q_{0}\vert d\mathbf{M}\vert q_{0}\rangle=d\lambda_{0}
151 : \f}
152 : So, using simple chain rules:
153 : \f{displaymath}
154 : \langle q_{0}\vert \frac{d\mathbf{M}}{dr_{k}^{\gamma}}\vert q_{0}\rangle
155 : dr_{k}^{\gamma}=d\lambda_{0}
156 : \f}
157 : where here we used the notation \f$r_{k}^{\gamma}\f$ to denote an arbitrary position which can be
158 : \f$\tilde x_0 ,\tilde y_0,\tilde z_0\f$ or \f$\tilde x_1 ,\tilde y_1,\tilde z_1\f$
159 : we get
160 : \f{displaymath}
161 : \langle q_{0}\vert \frac{d\mathbf{M}}{dr_{k}^{\gamma}}\vert q_{0}\rangle
162 : =\frac{d\lambda_{0}}{dr_{k}^{\gamma}}
163 : \f}
164 :
165 : The derivatives of the matrix \f$\frac{d\mathbf{M}}{dr_{k}^{\gamma}} \f$ can be readily obtained via the
166 : chain rule
167 : \f{displaymath}
168 : \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}} +\\
169 : \frac{d\mathbf{M}}{dM_{l}^{\alpha}}\frac{dM_{l}^{\alpha}}{dr_{k}^{\gamma}}
170 : \f}
171 :
172 : where \f$ M_{l}^{\alpha} \f$ corresponds to \f$ x_m^{l},y_m^{l},z_m^{l} \f$ and
173 : \f$ P_{l}^{\alpha} \f$ corresponds to \f$ x_p^{l},y_p^{l},z_p^{l} \f$ according to the \f$ \alpha \f$ component.
174 : */
175 :
176 :
177 0 : class Kearsley
178 : {
179 : /// general log reference that needs to be initialized when constructed
180 : Log* log;
181 : /// position of atoms (first frame. In md is the running frame)
182 : std::vector<Vector> p0;
183 : /// position of atoms (second frame. In md is the reference frame)
184 : std::vector<Vector> p1;
185 : /// alignment weight: the rmsd/msd that it provides is only based on this scalar
186 : std::vector<double> align;
187 :
188 : bool com0_is_removed;
189 : bool com1_is_removed;
190 :
191 : public:
192 : /// error: the distance between two frames (might be rmsd/msd. See below)
193 : double err;
194 : /// displacement: the vector that goes from the p0 onto p1
195 : std::vector<Vector> diff0on1;
196 : /// displacement: the vector that goes from the p1 onto p0 (via inverse rotation)
197 : std::vector<Vector> diff1on0;
198 :
199 : /// center of mass of p0
200 : Vector com0;
201 : /// center of mass of p1
202 : Vector com1;
203 : /// position resetted wrt coms p0
204 : std::vector<Vector> p0reset;
205 : /// position resetted wrt coms p1
206 : std::vector<Vector> p1reset;
207 : /// position rotated: p0
208 : std::vector<Vector> p0rotated;
209 : /// position rotated: p1
210 : std::vector<Vector> p1rotated;
211 : /// rotation matrices p0 on p1 and reverse (p1 over p0)
212 : Tensor rotmat0on1,rotmat1on0;
213 : /// derivatives: derivative of the error respect p0
214 : std::vector<Vector> derrdp0;
215 : /// derivatives: derivative of the error respect p1
216 : std::vector<Vector> derrdp1;
217 : /// derivative of the rotation matrix
218 : /// note the dimension 3x3 x 3 x N
219 : std::vector<double> dmatdp0;
220 : std::vector<double> dmatdp1;
221 :
222 : /// constructor: need the two structure, the alignment vector and the log reference
223 : Kearsley( const std::vector<Vector> &p0, const std::vector<Vector> &p1, const std::vector<double> &align, Log* &log);
224 : /// switch the assignment of the structure p0 (e.g. at each md step)
225 : void assignP0(const std::vector<Vector> & p0);
226 : /// derivatives: derivative of the error respect p1
227 : void assignP1(const std::vector<Vector> & p1);
228 : /// transfer the alignment vector
229 : void assignAlign(const std::vector<double> & align);
230 : /// finite differences of all the relevant quantities: takes a bool which decides if giving back rmsd or not (msd in this case)
231 : void finiteDifferenceInterface(bool rmsd);
232 : // this makes the real calculation: the rmsd bool decides wether doing rmsd or msd
233 : double calculate( bool rmsd );
234 : };
235 :
236 : }
237 :
238 : #endif
239 :
|
