Collective Variables Module - Developer Documentation
Loading...
Searching...
No Matches
colvar_rotation_derivative.h
1#ifndef COLVAR_ROTATION_DERIVATIVE
2#define COLVAR_ROTATION_DERIVATIVE
3
4#include "colvartypes.h"
5#include <type_traits>
6#include <cstring>
7
9template <typename T, typename std::enable_if<std::is_same<T, cvm::atom_pos>::value, bool>::type = true>
10inline void read_atom_coord(
11 size_t ia, const std::vector<T>& pos,
12 cvm::real* x, cvm::real* y, cvm::real* z) {
13 *x = pos[ia].x;
14 *y = pos[ia].y;
15 *z = pos[ia].z;
16}
17
18template <typename T, typename std::enable_if<std::is_same<T, cvm::atom>::value, bool>::type = true>
19inline void read_atom_coord(
20 size_t ia, const std::vector<T>& pos,
21 cvm::real* x, cvm::real* y, cvm::real* z) {
22 *x = pos[ia].pos.x;
23 *y = pos[ia].pos.y;
24 *z = pos[ia].pos.z;
25}
26
28enum class rotation_derivative_dldq {
30 use_dl = 1 << 0,
32 use_dq = 1 << 1
33};
34
35inline constexpr rotation_derivative_dldq operator|(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs) {
36 return static_cast<rotation_derivative_dldq>(
37 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Lhs) |
38 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Rhs));
39}
40
41inline constexpr bool operator&(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs)
42{
43 return (static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Lhs) &
44 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Rhs));
45}
46
48template <typename T1, typename T2>
49struct rotation_derivative {
50 static_assert(std::is_same<T1, cvm::atom_pos>::value || std::is_same<T1, cvm::atom>::value,
51 "class template rotation_derivative only supports cvm::atom_pos or cvm::atom types.");
52 static_assert(std::is_same<T2, cvm::atom_pos>::value || std::is_same<T2, cvm::atom>::value,
53 "class template rotation_derivative only supports cvm::atom_pos or cvm::atom types.");
55 const cvm::rotation &m_rot;
57 const std::vector<T1> &m_pos1;
59 const std::vector<T2> &m_pos2;
61 cvm::real tmp_Q0Q0[4][4];
62 cvm::real tmp_Q0Q0_L[4][4][4];
71 rotation_derivative(
72 const cvm::rotation &rot,
73 const std::vector<T1> &pos1,
74 const std::vector<T2> &pos2):
75 m_rot(rot), m_pos1(pos1), m_pos2(pos2) {};
82 void prepare_derivative(rotation_derivative_dldq require_dl_dq) {
83 if (require_dl_dq & rotation_derivative_dldq::use_dl) {
84 const auto &Q0 = m_rot.S_eigvec[0];
85 tmp_Q0Q0[0][0] = Q0[0] * Q0[0];
86 tmp_Q0Q0[0][1] = Q0[0] * Q0[1];
87 tmp_Q0Q0[0][2] = Q0[0] * Q0[2];
88 tmp_Q0Q0[0][3] = Q0[0] * Q0[3];
89 tmp_Q0Q0[1][0] = Q0[1] * Q0[0];
90 tmp_Q0Q0[1][1] = Q0[1] * Q0[1];
91 tmp_Q0Q0[1][2] = Q0[1] * Q0[2];
92 tmp_Q0Q0[1][3] = Q0[1] * Q0[3];
93 tmp_Q0Q0[2][0] = Q0[2] * Q0[0];
94 tmp_Q0Q0[2][1] = Q0[2] * Q0[1];
95 tmp_Q0Q0[2][2] = Q0[2] * Q0[2];
96 tmp_Q0Q0[2][3] = Q0[2] * Q0[3];
97 tmp_Q0Q0[3][0] = Q0[3] * Q0[0];
98 tmp_Q0Q0[3][1] = Q0[3] * Q0[1];
99 tmp_Q0Q0[3][2] = Q0[3] * Q0[2];
100 tmp_Q0Q0[3][3] = Q0[3] * Q0[3];
101 }
102 if (require_dl_dq & rotation_derivative_dldq::use_dq) {
103 const auto &Q0 = m_rot.S_eigvec[0];
104 const auto &Q1 = m_rot.S_eigvec[1];
105 const auto &Q2 = m_rot.S_eigvec[2];
106 const auto &Q3 = m_rot.S_eigvec[3];
107 cvm::real const L0 = m_rot.S_eigval[0];
108 cvm::real const L1 = m_rot.S_eigval[1];
109 cvm::real const L2 = m_rot.S_eigval[2];
110 cvm::real const L3 = m_rot.S_eigval[3];
111
112 tmp_Q0Q0_L[0][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[0] +
113 (Q2[0] * Q0[0]) / (L0-L2) * Q2[0] +
114 (Q3[0] * Q0[0]) / (L0-L3) * Q3[0];
115 tmp_Q0Q0_L[1][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[1] +
116 (Q2[0] * Q0[0]) / (L0-L2) * Q2[1] +
117 (Q3[0] * Q0[0]) / (L0-L3) * Q3[1];
118 tmp_Q0Q0_L[2][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[2] +
119 (Q2[0] * Q0[0]) / (L0-L2) * Q2[2] +
120 (Q3[0] * Q0[0]) / (L0-L3) * Q3[2];
121 tmp_Q0Q0_L[3][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[3] +
122 (Q2[0] * Q0[0]) / (L0-L2) * Q2[3] +
123 (Q3[0] * Q0[0]) / (L0-L3) * Q3[3];
124
125 tmp_Q0Q0_L[0][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[0] +
126 (Q2[0] * Q0[1]) / (L0-L2) * Q2[0] +
127 (Q3[0] * Q0[1]) / (L0-L3) * Q3[0];
128 tmp_Q0Q0_L[1][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[1] +
129 (Q2[0] * Q0[1]) / (L0-L2) * Q2[1] +
130 (Q3[0] * Q0[1]) / (L0-L3) * Q3[1];
131 tmp_Q0Q0_L[2][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[2] +
132 (Q2[0] * Q0[1]) / (L0-L2) * Q2[2] +
133 (Q3[0] * Q0[1]) / (L0-L3) * Q3[2];
134 tmp_Q0Q0_L[3][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[3] +
135 (Q2[0] * Q0[1]) / (L0-L2) * Q2[3] +
136 (Q3[0] * Q0[1]) / (L0-L3) * Q3[3];
137
138
139 tmp_Q0Q0_L[0][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[0] +
140 (Q2[0] * Q0[2]) / (L0-L2) * Q2[0] +
141 (Q3[0] * Q0[2]) / (L0-L3) * Q3[0];
142 tmp_Q0Q0_L[1][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[1] +
143 (Q2[0] * Q0[2]) / (L0-L2) * Q2[1] +
144 (Q3[0] * Q0[2]) / (L0-L3) * Q3[1];
145 tmp_Q0Q0_L[2][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[2] +
146 (Q2[0] * Q0[2]) / (L0-L2) * Q2[2] +
147 (Q3[0] * Q0[2]) / (L0-L3) * Q3[2];
148 tmp_Q0Q0_L[3][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[3] +
149 (Q2[0] * Q0[2]) / (L0-L2) * Q2[3] +
150 (Q3[0] * Q0[2]) / (L0-L3) * Q3[3];
151
152 tmp_Q0Q0_L[0][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[0] +
153 (Q2[0] * Q0[3]) / (L0-L2) * Q2[0] +
154 (Q3[0] * Q0[3]) / (L0-L3) * Q3[0];
155 tmp_Q0Q0_L[1][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[1] +
156 (Q2[0] * Q0[3]) / (L0-L2) * Q2[1] +
157 (Q3[0] * Q0[3]) / (L0-L3) * Q3[1];
158 tmp_Q0Q0_L[2][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[2] +
159 (Q2[0] * Q0[3]) / (L0-L2) * Q2[2] +
160 (Q3[0] * Q0[3]) / (L0-L3) * Q3[2];
161 tmp_Q0Q0_L[3][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[3] +
162 (Q2[0] * Q0[3]) / (L0-L2) * Q2[3] +
163 (Q3[0] * Q0[3]) / (L0-L3) * Q3[3];
164
165 tmp_Q0Q0_L[0][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[0] +
166 (Q2[1] * Q0[0]) / (L0-L2) * Q2[0] +
167 (Q3[1] * Q0[0]) / (L0-L3) * Q3[0];
168 tmp_Q0Q0_L[1][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[1] +
169 (Q2[1] * Q0[0]) / (L0-L2) * Q2[1] +
170 (Q3[1] * Q0[0]) / (L0-L3) * Q3[1];
171 tmp_Q0Q0_L[2][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[2] +
172 (Q2[1] * Q0[0]) / (L0-L2) * Q2[2] +
173 (Q3[1] * Q0[0]) / (L0-L3) * Q3[2];
174 tmp_Q0Q0_L[3][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[3] +
175 (Q2[1] * Q0[0]) / (L0-L2) * Q2[3] +
176 (Q3[1] * Q0[0]) / (L0-L3) * Q3[3];
177
178 tmp_Q0Q0_L[0][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[0] +
179 (Q2[1] * Q0[1]) / (L0-L2) * Q2[0] +
180 (Q3[1] * Q0[1]) / (L0-L3) * Q3[0];
181 tmp_Q0Q0_L[1][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[1] +
182 (Q2[1] * Q0[1]) / (L0-L2) * Q2[1] +
183 (Q3[1] * Q0[1]) / (L0-L3) * Q3[1];
184 tmp_Q0Q0_L[2][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[2] +
185 (Q2[1] * Q0[1]) / (L0-L2) * Q2[2] +
186 (Q3[1] * Q0[1]) / (L0-L3) * Q3[2];
187 tmp_Q0Q0_L[3][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[3] +
188 (Q2[1] * Q0[1]) / (L0-L2) * Q2[3] +
189 (Q3[1] * Q0[1]) / (L0-L3) * Q3[3];
190
191 tmp_Q0Q0_L[0][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[0] +
192 (Q2[1] * Q0[2]) / (L0-L2) * Q2[0] +
193 (Q3[1] * Q0[2]) / (L0-L3) * Q3[0];
194 tmp_Q0Q0_L[1][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[1] +
195 (Q2[1] * Q0[2]) / (L0-L2) * Q2[1] +
196 (Q3[1] * Q0[2]) / (L0-L3) * Q3[1];
197 tmp_Q0Q0_L[2][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[2] +
198 (Q2[1] * Q0[2]) / (L0-L2) * Q2[2] +
199 (Q3[1] * Q0[2]) / (L0-L3) * Q3[2];
200 tmp_Q0Q0_L[3][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[3] +
201 (Q2[1] * Q0[2]) / (L0-L2) * Q2[3] +
202 (Q3[1] * Q0[2]) / (L0-L3) * Q3[3];
203
204 tmp_Q0Q0_L[0][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[0] +
205 (Q2[1] * Q0[3]) / (L0-L2) * Q2[0] +
206 (Q3[1] * Q0[3]) / (L0-L3) * Q3[0];
207 tmp_Q0Q0_L[1][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[1] +
208 (Q2[1] * Q0[3]) / (L0-L2) * Q2[1] +
209 (Q3[1] * Q0[3]) / (L0-L3) * Q3[1];
210 tmp_Q0Q0_L[2][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[2] +
211 (Q2[1] * Q0[3]) / (L0-L2) * Q2[2] +
212 (Q3[1] * Q0[3]) / (L0-L3) * Q3[2];
213 tmp_Q0Q0_L[3][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[3] +
214 (Q2[1] * Q0[3]) / (L0-L2) * Q2[3] +
215 (Q3[1] * Q0[3]) / (L0-L3) * Q3[3];
216
217
218 tmp_Q0Q0_L[0][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[0] +
219 (Q2[2] * Q0[0]) / (L0-L2) * Q2[0] +
220 (Q3[2] * Q0[0]) / (L0-L3) * Q3[0];
221 tmp_Q0Q0_L[1][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[1] +
222 (Q2[2] * Q0[0]) / (L0-L2) * Q2[1] +
223 (Q3[2] * Q0[0]) / (L0-L3) * Q3[1];
224 tmp_Q0Q0_L[2][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[2] +
225 (Q2[2] * Q0[0]) / (L0-L2) * Q2[2] +
226 (Q3[2] * Q0[0]) / (L0-L3) * Q3[2];
227 tmp_Q0Q0_L[3][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[3] +
228 (Q2[2] * Q0[0]) / (L0-L2) * Q2[3] +
229 (Q3[2] * Q0[0]) / (L0-L3) * Q3[3];
230
231 tmp_Q0Q0_L[0][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[0] +
232 (Q2[2] * Q0[1]) / (L0-L2) * Q2[0] +
233 (Q3[2] * Q0[1]) / (L0-L3) * Q3[0];
234 tmp_Q0Q0_L[1][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[1] +
235 (Q2[2] * Q0[1]) / (L0-L2) * Q2[1] +
236 (Q3[2] * Q0[1]) / (L0-L3) * Q3[1];
237 tmp_Q0Q0_L[2][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[2] +
238 (Q2[2] * Q0[1]) / (L0-L2) * Q2[2] +
239 (Q3[2] * Q0[1]) / (L0-L3) * Q3[2];
240 tmp_Q0Q0_L[3][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[3] +
241 (Q2[2] * Q0[1]) / (L0-L2) * Q2[3] +
242 (Q3[2] * Q0[1]) / (L0-L3) * Q3[3];
243
244 tmp_Q0Q0_L[0][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[0] +
245 (Q2[2] * Q0[2]) / (L0-L2) * Q2[0] +
246 (Q3[2] * Q0[2]) / (L0-L3) * Q3[0];
247 tmp_Q0Q0_L[1][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[1] +
248 (Q2[2] * Q0[2]) / (L0-L2) * Q2[1] +
249 (Q3[2] * Q0[2]) / (L0-L3) * Q3[1];
250 tmp_Q0Q0_L[2][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[2] +
251 (Q2[2] * Q0[2]) / (L0-L2) * Q2[2] +
252 (Q3[2] * Q0[2]) / (L0-L3) * Q3[2];
253 tmp_Q0Q0_L[3][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[3] +
254 (Q2[2] * Q0[2]) / (L0-L2) * Q2[3] +
255 (Q3[2] * Q0[2]) / (L0-L3) * Q3[3];
256
257 tmp_Q0Q0_L[0][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[0] +
258 (Q2[2] * Q0[3]) / (L0-L2) * Q2[0] +
259 (Q3[2] * Q0[3]) / (L0-L3) * Q3[0];
260 tmp_Q0Q0_L[1][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[1] +
261 (Q2[2] * Q0[3]) / (L0-L2) * Q2[1] +
262 (Q3[2] * Q0[3]) / (L0-L3) * Q3[1];
263 tmp_Q0Q0_L[2][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[2] +
264 (Q2[2] * Q0[3]) / (L0-L2) * Q2[2] +
265 (Q3[2] * Q0[3]) / (L0-L3) * Q3[2];
266 tmp_Q0Q0_L[3][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[3] +
267 (Q2[2] * Q0[3]) / (L0-L2) * Q2[3] +
268 (Q3[2] * Q0[3]) / (L0-L3) * Q3[3];
269
270 tmp_Q0Q0_L[0][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[0] +
271 (Q2[3] * Q0[0]) / (L0-L2) * Q2[0] +
272 (Q3[3] * Q0[0]) / (L0-L3) * Q3[0];
273 tmp_Q0Q0_L[1][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[1] +
274 (Q2[3] * Q0[0]) / (L0-L2) * Q2[1] +
275 (Q3[3] * Q0[0]) / (L0-L3) * Q3[1];
276 tmp_Q0Q0_L[2][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[2] +
277 (Q2[3] * Q0[0]) / (L0-L2) * Q2[2] +
278 (Q3[3] * Q0[0]) / (L0-L3) * Q3[2];
279 tmp_Q0Q0_L[3][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[3] +
280 (Q2[3] * Q0[0]) / (L0-L2) * Q2[3] +
281 (Q3[3] * Q0[0]) / (L0-L3) * Q3[3];
282
283 tmp_Q0Q0_L[0][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[0] +
284 (Q2[3] * Q0[1]) / (L0-L2) * Q2[0] +
285 (Q3[3] * Q0[1]) / (L0-L3) * Q3[0];
286 tmp_Q0Q0_L[1][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[1] +
287 (Q2[3] * Q0[1]) / (L0-L2) * Q2[1] +
288 (Q3[3] * Q0[1]) / (L0-L3) * Q3[1];
289 tmp_Q0Q0_L[2][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[2] +
290 (Q2[3] * Q0[1]) / (L0-L2) * Q2[2] +
291 (Q3[3] * Q0[1]) / (L0-L3) * Q3[2];
292 tmp_Q0Q0_L[3][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[3] +
293 (Q2[3] * Q0[1]) / (L0-L2) * Q2[3] +
294 (Q3[3] * Q0[1]) / (L0-L3) * Q3[3];
295
296 tmp_Q0Q0_L[0][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[0] +
297 (Q2[3] * Q0[2]) / (L0-L2) * Q2[0] +
298 (Q3[3] * Q0[2]) / (L0-L3) * Q3[0];
299 tmp_Q0Q0_L[1][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[1] +
300 (Q2[3] * Q0[2]) / (L0-L2) * Q2[1] +
301 (Q3[3] * Q0[2]) / (L0-L3) * Q3[1];
302 tmp_Q0Q0_L[2][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[2] +
303 (Q2[3] * Q0[2]) / (L0-L2) * Q2[2] +
304 (Q3[3] * Q0[2]) / (L0-L3) * Q3[2];
305 tmp_Q0Q0_L[3][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[3] +
306 (Q2[3] * Q0[2]) / (L0-L2) * Q2[3] +
307 (Q3[3] * Q0[2]) / (L0-L3) * Q3[3];
308
309 tmp_Q0Q0_L[0][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[0] +
310 (Q2[3] * Q0[3]) / (L0-L2) * Q2[0] +
311 (Q3[3] * Q0[3]) / (L0-L3) * Q3[0];
312 tmp_Q0Q0_L[1][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[1] +
313 (Q2[3] * Q0[3]) / (L0-L2) * Q2[1] +
314 (Q3[3] * Q0[3]) / (L0-L3) * Q3[1];
315 tmp_Q0Q0_L[2][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[2] +
316 (Q2[3] * Q0[3]) / (L0-L2) * Q2[2] +
317 (Q3[3] * Q0[3]) / (L0-L3) * Q3[2];
318 tmp_Q0Q0_L[3][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[3] +
319 (Q2[3] * Q0[3]) / (L0-L2) * Q2[3] +
320 (Q3[3] * Q0[3]) / (L0-L3) * Q3[3];
321 }
322 }
330 void calc_derivative_impl(
331 const cvm::rvector (&ds)[4][4],
332 cvm::rvector* const dl0_out,
333 cvm::vector1d<cvm::rvector>* const dq0_out,
334 cvm::matrix2d<cvm::rvector>* const ds_out) const {
335 if (ds_out != nullptr) {
336 // this code path is for debug_gradients, so not necessary to unroll the loop
337 *ds_out = cvm::matrix2d<cvm::rvector>(4, 4);
338 for (int i = 0; i < 4; ++i) {
339 for (int j = 0; j < 4; ++j) {
340 (*ds_out)[i][j] = ds[i][j];
341 }
342 }
343 }
344 if (dl0_out != nullptr) {
345 /* manually loop unrolling of the following loop:
346 dl0_1.reset();
347 for (size_t i = 0; i < 4; i++) {
348 for (size_t j = 0; j < 4; j++) {
349 dl0_1 += Q0[i] * ds_1[i][j] * Q0[j];
350 }
351 }
352 */
353 *dl0_out = tmp_Q0Q0[0][0] * ds[0][0] +
354 tmp_Q0Q0[0][1] * ds[0][1] +
355 tmp_Q0Q0[0][2] * ds[0][2] +
356 tmp_Q0Q0[0][3] * ds[0][3] +
357 tmp_Q0Q0[1][0] * ds[1][0] +
358 tmp_Q0Q0[1][1] * ds[1][1] +
359 tmp_Q0Q0[1][2] * ds[1][2] +
360 tmp_Q0Q0[1][3] * ds[1][3] +
361 tmp_Q0Q0[2][0] * ds[2][0] +
362 tmp_Q0Q0[2][1] * ds[2][1] +
363 tmp_Q0Q0[2][2] * ds[2][2] +
364 tmp_Q0Q0[2][3] * ds[2][3] +
365 tmp_Q0Q0[3][0] * ds[3][0] +
366 tmp_Q0Q0[3][1] * ds[3][1] +
367 tmp_Q0Q0[3][2] * ds[3][2] +
368 tmp_Q0Q0[3][3] * ds[3][3];
369 }
370 if (dq0_out != nullptr) {
371 // we can skip this check if a fixed-size array is used
372 if (dq0_out->size() != 4) dq0_out->resize(4);
373 /* manually loop unrolling of the following loop:
374 dq0_1.reset();
375 for (size_t p = 0; p < 4; p++) {
376 for (size_t i = 0; i < 4; i++) {
377 for (size_t j = 0; j < 4; j++) {
378 dq0_1[p] +=
379 (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
380 (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
381 (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p];
382 }
383 }
384 }
385 */
386 (*dq0_out)[0] = tmp_Q0Q0_L[0][0][0] * ds[0][0] +
387 tmp_Q0Q0_L[0][0][1] * ds[0][1] +
388 tmp_Q0Q0_L[0][0][2] * ds[0][2] +
389 tmp_Q0Q0_L[0][0][3] * ds[0][3] +
390 tmp_Q0Q0_L[0][1][0] * ds[1][0] +
391 tmp_Q0Q0_L[0][1][1] * ds[1][1] +
392 tmp_Q0Q0_L[0][1][2] * ds[1][2] +
393 tmp_Q0Q0_L[0][1][3] * ds[1][3] +
394 tmp_Q0Q0_L[0][2][0] * ds[2][0] +
395 tmp_Q0Q0_L[0][2][1] * ds[2][1] +
396 tmp_Q0Q0_L[0][2][2] * ds[2][2] +
397 tmp_Q0Q0_L[0][2][3] * ds[2][3] +
398 tmp_Q0Q0_L[0][3][0] * ds[3][0] +
399 tmp_Q0Q0_L[0][3][1] * ds[3][1] +
400 tmp_Q0Q0_L[0][3][2] * ds[3][2] +
401 tmp_Q0Q0_L[0][3][3] * ds[3][3];
402
403 (*dq0_out)[1] = tmp_Q0Q0_L[1][0][0] * ds[0][0] +
404 tmp_Q0Q0_L[1][0][1] * ds[0][1] +
405 tmp_Q0Q0_L[1][0][2] * ds[0][2] +
406 tmp_Q0Q0_L[1][0][3] * ds[0][3] +
407 tmp_Q0Q0_L[1][1][0] * ds[1][0] +
408 tmp_Q0Q0_L[1][1][1] * ds[1][1] +
409 tmp_Q0Q0_L[1][1][2] * ds[1][2] +
410 tmp_Q0Q0_L[1][1][3] * ds[1][3] +
411 tmp_Q0Q0_L[1][2][0] * ds[2][0] +
412 tmp_Q0Q0_L[1][2][1] * ds[2][1] +
413 tmp_Q0Q0_L[1][2][2] * ds[2][2] +
414 tmp_Q0Q0_L[1][2][3] * ds[2][3] +
415 tmp_Q0Q0_L[1][3][0] * ds[3][0] +
416 tmp_Q0Q0_L[1][3][1] * ds[3][1] +
417 tmp_Q0Q0_L[1][3][2] * ds[3][2] +
418 tmp_Q0Q0_L[1][3][3] * ds[3][3];
419
420 (*dq0_out)[2] = tmp_Q0Q0_L[2][0][0] * ds[0][0] +
421 tmp_Q0Q0_L[2][0][1] * ds[0][1] +
422 tmp_Q0Q0_L[2][0][2] * ds[0][2] +
423 tmp_Q0Q0_L[2][0][3] * ds[0][3] +
424 tmp_Q0Q0_L[2][1][0] * ds[1][0] +
425 tmp_Q0Q0_L[2][1][1] * ds[1][1] +
426 tmp_Q0Q0_L[2][1][2] * ds[1][2] +
427 tmp_Q0Q0_L[2][1][3] * ds[1][3] +
428 tmp_Q0Q0_L[2][2][0] * ds[2][0] +
429 tmp_Q0Q0_L[2][2][1] * ds[2][1] +
430 tmp_Q0Q0_L[2][2][2] * ds[2][2] +
431 tmp_Q0Q0_L[2][2][3] * ds[2][3] +
432 tmp_Q0Q0_L[2][3][0] * ds[3][0] +
433 tmp_Q0Q0_L[2][3][1] * ds[3][1] +
434 tmp_Q0Q0_L[2][3][2] * ds[3][2] +
435 tmp_Q0Q0_L[2][3][3] * ds[3][3];
436
437 (*dq0_out)[3] = tmp_Q0Q0_L[3][0][0] * ds[0][0] +
438 tmp_Q0Q0_L[3][0][1] * ds[0][1] +
439 tmp_Q0Q0_L[3][0][2] * ds[0][2] +
440 tmp_Q0Q0_L[3][0][3] * ds[0][3] +
441 tmp_Q0Q0_L[3][1][0] * ds[1][0] +
442 tmp_Q0Q0_L[3][1][1] * ds[1][1] +
443 tmp_Q0Q0_L[3][1][2] * ds[1][2] +
444 tmp_Q0Q0_L[3][1][3] * ds[1][3] +
445 tmp_Q0Q0_L[3][2][0] * ds[2][0] +
446 tmp_Q0Q0_L[3][2][1] * ds[2][1] +
447 tmp_Q0Q0_L[3][2][2] * ds[2][2] +
448 tmp_Q0Q0_L[3][2][3] * ds[2][3] +
449 tmp_Q0Q0_L[3][3][0] * ds[3][0] +
450 tmp_Q0Q0_L[3][3][1] * ds[3][1] +
451 tmp_Q0Q0_L[3][3][2] * ds[3][2] +
452 tmp_Q0Q0_L[3][3][3] * ds[3][3];
453 }
454 }
465 void calc_derivative_wrt_group1(
466 size_t ia, cvm::rvector* const dl0_1_out = nullptr,
467 cvm::vector1d<cvm::rvector>* const dq0_1_out = nullptr,
468 cvm::matrix2d<cvm::rvector>* const ds_1_out = nullptr) const {
469 if (dl0_1_out == nullptr && dq0_1_out == nullptr) return;
470 cvm::real a2x, a2y, a2z;
471 // we can get rid of the helper function read_atom_coord if C++17 (constexpr) is available
472 read_atom_coord(ia, m_pos2, &a2x, &a2y, &a2z);
473 cvm::rvector ds_1[4][4];
474 ds_1[0][0].set( a2x, a2y, a2z);
475 ds_1[1][0].set( 0.0, a2z, -a2y);
476 ds_1[0][1] = ds_1[1][0];
477 ds_1[2][0].set(-a2z, 0.0, a2x);
478 ds_1[0][2] = ds_1[2][0];
479 ds_1[3][0].set( a2y, -a2x, 0.0);
480 ds_1[0][3] = ds_1[3][0];
481 ds_1[1][1].set( a2x, -a2y, -a2z);
482 ds_1[2][1].set( a2y, a2x, 0.0);
483 ds_1[1][2] = ds_1[2][1];
484 ds_1[3][1].set( a2z, 0.0, a2x);
485 ds_1[1][3] = ds_1[3][1];
486 ds_1[2][2].set(-a2x, a2y, -a2z);
487 ds_1[3][2].set( 0.0, a2z, a2y);
488 ds_1[2][3] = ds_1[3][2];
489 ds_1[3][3].set(-a2x, -a2y, a2z);
490 calc_derivative_impl(ds_1, dl0_1_out, dq0_1_out, ds_1_out);
491 }
502 void calc_derivative_wrt_group2(
503 size_t ia, cvm::rvector* const dl0_2_out = nullptr,
504 cvm::vector1d<cvm::rvector>* const dq0_2_out = nullptr,
505 cvm::matrix2d<cvm::rvector>* const ds_2_out = nullptr) const {
506 if (dl0_2_out == nullptr && dq0_2_out == nullptr) return;
507 cvm::real a1x, a1y, a1z;
508 // we can get rid of the helper function read_atom_coord if C++17 (constexpr) is available
509 read_atom_coord(ia, m_pos1, &a1x, &a1y, &a1z);
510 cvm::rvector ds_2[4][4];
511 ds_2[0][0].set( a1x, a1y, a1z);
512 ds_2[1][0].set( 0.0, -a1z, a1y);
513 ds_2[0][1] = ds_2[1][0];
514 ds_2[2][0].set( a1z, 0.0, -a1x);
515 ds_2[0][2] = ds_2[2][0];
516 ds_2[3][0].set(-a1y, a1x, 0.0);
517 ds_2[0][3] = ds_2[3][0];
518 ds_2[1][1].set( a1x, -a1y, -a1z);
519 ds_2[2][1].set( a1y, a1x, 0.0);
520 ds_2[1][2] = ds_2[2][1];
521 ds_2[3][1].set( a1z, 0.0, a1x);
522 ds_2[1][3] = ds_2[3][1];
523 ds_2[2][2].set(-a1x, a1y, -a1z);
524 ds_2[3][2].set( 0.0, a1z, a1y);
525 ds_2[2][3] = ds_2[3][2];
526 ds_2[3][3].set(-a1x, -a1y, a1z);
527 calc_derivative_impl(ds_2, dl0_2_out, dq0_2_out, ds_2_out);
528 }
529};
530
537template<typename T1, typename T2>
538void debug_gradients(
539 cvm::rotation &rot,
540 const std::vector<T1> &pos1,
541 const std::vector<T2> &pos2) {
542 static_assert(std::is_same<T1, cvm::atom_pos>::value || std::is_same<T1, cvm::atom>::value, "");
543 static_assert(std::is_same<T2, cvm::atom_pos>::value || std::is_same<T2, cvm::atom>::value, "");
544 // eigenvalues and eigenvectors
545 cvm::real const L0 = rot.S_eigval[0];
546 cvm::real const L1 = rot.S_eigval[1];
547 cvm::real const L2 = rot.S_eigval[2];
548 cvm::real const L3 = rot.S_eigval[3];
549 cvm::quaternion const Q0(rot.S_eigvec[0]);
550 cvm::quaternion const Q1(rot.S_eigvec[1]);
551 cvm::quaternion const Q2(rot.S_eigvec[2]);
552 cvm::quaternion const Q3(rot.S_eigvec[3]);
553
554 cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
555 ", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
556 ", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
557 "\n");
558 cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
559 ", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
560 ", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
561 "\n");
562 cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
563 ", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
564 ", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
565 "\n");
566 cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
567 ", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
568 ", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
569 "\n");
570
571 rotation_derivative<T1, T2> deriv(rot, pos1, pos2);
572 cvm::rvector dl0_2;
573 cvm::vector1d<cvm::rvector> dq0_2(4);
574 cvm::matrix2d<cvm::rvector> ds_2;
575#ifdef COLVARS_LAMMPS
576 MathEigen::Jacobi<cvm::real,
577 cvm::real[4],
578 cvm::real[4][4]> *ecalc =
579 reinterpret_cast<MathEigen::Jacobi<cvm::real,
580 cvm::real[4],
581 cvm::real[4][4]> *>(rot.jacobi);
582#endif
583 deriv.prepare_derivative(rotation_derivative_dldq::use_dl | rotation_derivative_dldq::use_dq);
584 cvm::real S_new[4][4];
585 cvm::real S_new_eigval[4];
586 cvm::real S_new_eigvec[4][4];
587 for (size_t ia = 0; ia < pos2.size(); ++ia) {
588 // cvm::real const &a1x = pos1[ia].x;
589 // cvm::real const &a1y = pos1[ia].y;
590 // cvm::real const &a1z = pos1[ia].z;
591 deriv.calc_derivative_wrt_group2(ia, &dl0_2, &dq0_2, &ds_2);
592 // make an infitesimal move along each cartesian coordinate of
593 // this atom, and solve again the eigenvector problem
594 for (size_t comp = 0; comp < 3; comp++) {
595 std::memcpy(S_new, rot.S_backup, sizeof(cvm::real) * 4 * 4);
596 std::memset(S_new_eigval, 0, sizeof(cvm::real) * 4);
597 std::memset(S_new_eigvec, 0, sizeof(cvm::real) * 4 * 4);
598 for (size_t i = 0; i < 4; i++) {
599 for (size_t j = 0; j < 4; j++) {
600 S_new[i][j] +=
601 colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
602 }
603 }
604#ifdef COLVARS_LAMMPS
605 ecalc->Diagonalize(S_new, S_new_eigval, S_new_eigvec);
606#else
607 NR::diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
608#endif
609 cvm::real const &L0_new = S_new_eigval[0];
610 cvm::quaternion const Q0_new(S_new_eigvec[0]);
611
612 cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
613 cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
614 dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
615 dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
616 dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
617
618 cvm::log( "|(l_0+dl_0) - l_0^new|/l_0 = "+
619 cvm::to_str(cvm::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
620 ", |(q_0+dq_0) - q_0^new| = "+
621 cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
622 "\n");
623 }
624 }
625}
626
627#endif // COLVAR_ROTATION_DERIVATIVE
colvarmodule::matrix2d
Arbitrary size array (two dimensions) suitable for linear algebra operations (i.e....
Definition: colvartypes.h:376
colvarmodule::quaternion
1-dimensional vector of real numbers with four components and a quaternion algebra
Definition: colvartypes.h:958
colvarmodule::quaternion::inner
cvm::real inner(cvm::quaternion const &Q2) const
Inner product (as a 4-d vector) with Q2; requires match() if the largest overlap is looked for.
Definition: colvartypes.h:1292
colvarmodule::rotation
A rotation between two sets of coordinates (for the moment a wrapper for colvarmodule::quaternion)
Definition: colvartypes.h:1314
colvarmodule::rotation::S_eigval
cvm::real S_eigval[4]
Eigenvalues of S.
Definition: colvartypes.h:1323
colvarmodule::rotation::S_backup
cvm::real S_backup[4][4]
Used for debugging gradients.
Definition: colvartypes.h:1329
colvarmodule::rotation::S_eigvec
cvm::real S_eigvec[4][4]
Eigenvectors of S.
Definition: colvartypes.h:1326
colvarmodule::rotation::jacobi
void * jacobi
Pointer to instance of Jacobi solver.
Definition: colvartypes.h:1505
colvarmodule::rvector
vector of real numbers with three components
Definition: colvartypes.h:727
colvarmodule::rvector::set
void set(cvm::real value)
Set all components to a scalar.
Definition: colvartypes.h:760
colvarmodule::vector1d
Arbitrary size array (one dimensions) suitable for linear algebra operations (i.e....
Definition: colvartypes.h:37
colvarmodule::real
double real
Defining an abstract real number allows to switch precision.
Definition: colvarmodule.h:95
colvarmodule::log
static void log(std::string const &message, int min_log_level=10)
Definition: colvarmodule.cpp:1950
colvarmodule::cv_prec
static size_t const cv_prec
Number of digits to represent a collective variables value(s)
Definition: colvarmodule.h:664
colvarmodule::cv_width
static size_t const cv_width
Number of characters to represent a collective variables value(s)
Definition: colvarmodule.h:666
colvarmodule::fabs
static real fabs(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:127
colvarmodule::to_str
static std::string to_str(char const *s)
Convert to string for output purposes.
Definition: colvarmodule.cpp:2373
colvarmodule::debug_gradients_step_size
static real debug_gradients_step_size
Finite difference step size (if there is no dynamics, or if gradients need to be tested independently...
Definition: colvarmodule.h:250
rotation_derivative
Helper class for calculating the derivative of rotation.
Definition: colvar_rotation_derivative.h:49
rotation_derivative::calc_derivative_impl
void calc_derivative_impl(const cvm::rvector(&ds)[4][4], cvm::rvector *const dl0_out, cvm::vector1d< cvm::rvector > *const dq0_out, cvm::matrix2d< cvm::rvector > *const ds_out) const
Actual implementation of the derivative calculation.
Definition: colvar_rotation_derivative.h:330
rotation_derivative::prepare_derivative
void prepare_derivative(rotation_derivative_dldq require_dl_dq)
This function must be called before calc_derivative_wrt_group1 and calc_derivative_wrt_group2 in orde...
Definition: colvar_rotation_derivative.h:82
rotation_derivative::m_pos1
const std::vector< T1 > & m_pos1
Reference to the atom positions of group 1.
Definition: colvar_rotation_derivative.h:57
rotation_derivative::calc_derivative_wrt_group2
void calc_derivative_wrt_group2(size_t ia, cvm::rvector *const dl0_2_out=nullptr, cvm::vector1d< cvm::rvector > *const dq0_2_out=nullptr, cvm::matrix2d< cvm::rvector > *const ds_2_out=nullptr) const
Calculate the derivatives of S, the leading eigenvalue L and the leading eigenvector Q with respect t...
Definition: colvar_rotation_derivative.h:502
rotation_derivative::m_rot
const cvm::rotation & m_rot
Reference to the rotation.
Definition: colvar_rotation_derivative.h:55
rotation_derivative::tmp_Q0Q0
cvm::real tmp_Q0Q0[4][4]
Temporary variable that will be updated if prepare_derivative called.
Definition: colvar_rotation_derivative.h:61
rotation_derivative::rotation_derivative
rotation_derivative(const cvm::rotation &rot, const std::vector< T1 > &pos1, const std::vector< T2 > &pos2)
Constructor of the cvm::rotation::derivative class.
Definition: colvar_rotation_derivative.h:71
rotation_derivative::calc_derivative_wrt_group1
void calc_derivative_wrt_group1(size_t ia, cvm::rvector *const dl0_1_out=nullptr, cvm::vector1d< cvm::rvector > *const dq0_1_out=nullptr, cvm::matrix2d< cvm::rvector > *const ds_1_out=nullptr) const
Calculate the derivatives of S, the leading eigenvalue L and the leading eigenvector Q with respect t...
Definition: colvar_rotation_derivative.h:465
rotation_derivative::m_pos2
const std::vector< T2 > & m_pos2
Reference to the atom positions of group 2.
Definition: colvar_rotation_derivative.h:59
Generated on Tue Apr 16 2024 15:29:55 for Collective Variables Module - Developer Documentation by doxygen 1.9.5