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
8#ifndef _noalias
9#if defined(__INTEL_COMPILER) || (defined(__PGI) && !defined(__NVCOMPILER))
10#define _noalias restrict
11#elif defined(__GNUC__) || defined(__INTEL_LLVM_COMPILER) || defined(__NVCOMPILER)
12#define _noalias __restrict
13#else
14#define _noalias
15#endif
16#endif
17
19template <typename T, typename std::enable_if<std::is_same<T, cvm::atom_pos>::value, bool>::type = true>
20inline void read_atom_coord(
21 size_t ia, const std::vector<T>& pos,
22 cvm::real* _noalias x, cvm::real* _noalias y, cvm::real* _noalias z) {
23 *x = pos[ia].x;
24 *y = pos[ia].y;
25 *z = pos[ia].z;
26}
27
28template <typename T, typename std::enable_if<std::is_same<T, cvm::atom>::value, bool>::type = true>
29inline void read_atom_coord(
30 size_t ia, const std::vector<T>& pos,
31 cvm::real* _noalias x, cvm::real* _noalias y, cvm::real* _noalias z) {
32 *x = pos[ia].pos.x;
33 *y = pos[ia].pos.y;
34 *z = pos[ia].pos.z;
35}
36
38enum class rotation_derivative_dldq {
40 use_dl = 1 << 0,
42 use_dq = 1 << 1
43};
44
45inline constexpr rotation_derivative_dldq operator|(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs) {
46 return static_cast<rotation_derivative_dldq>(
47 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Lhs) |
48 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Rhs));
49}
50
51inline constexpr bool operator&(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs)
52{
53 return (static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Lhs) &
54 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Rhs));
55}
56
58template <typename T1, typename T2>
59struct rotation_derivative {
60 static_assert(std::is_same<T1, cvm::atom_pos>::value || std::is_same<T1, cvm::atom>::value,
61 "class template rotation_derivative only supports cvm::atom_pos or cvm::atom types.");
62 static_assert(std::is_same<T2, cvm::atom_pos>::value || std::is_same<T2, cvm::atom>::value,
63 "class template rotation_derivative only supports cvm::atom_pos or cvm::atom types.");
65 const cvm::rotation &m_rot;
67 const std::vector<T1> &m_pos1;
69 const std::vector<T2> &m_pos2;
71 cvm::real tmp_Q0Q0[4][4];
72 cvm::real tmp_Q0Q0_L[4][4][4];
81 rotation_derivative(
82 const cvm::rotation &rot,
83 const std::vector<T1> &pos1,
84 const std::vector<T2> &pos2):
85 m_rot(rot), m_pos1(pos1), m_pos2(pos2) {};
92 void prepare_derivative(rotation_derivative_dldq require_dl_dq) {
93 if (require_dl_dq & rotation_derivative_dldq::use_dl) {
94 const auto &Q0 = m_rot.S_eigvec[0];
95 tmp_Q0Q0[0][0] = Q0[0] * Q0[0];
96 tmp_Q0Q0[0][1] = Q0[0] * Q0[1];
97 tmp_Q0Q0[0][2] = Q0[0] * Q0[2];
98 tmp_Q0Q0[0][3] = Q0[0] * Q0[3];
99 tmp_Q0Q0[1][0] = Q0[1] * Q0[0];
100 tmp_Q0Q0[1][1] = Q0[1] * Q0[1];
101 tmp_Q0Q0[1][2] = Q0[1] * Q0[2];
102 tmp_Q0Q0[1][3] = Q0[1] * Q0[3];
103 tmp_Q0Q0[2][0] = Q0[2] * Q0[0];
104 tmp_Q0Q0[2][1] = Q0[2] * Q0[1];
105 tmp_Q0Q0[2][2] = Q0[2] * Q0[2];
106 tmp_Q0Q0[2][3] = Q0[2] * Q0[3];
107 tmp_Q0Q0[3][0] = Q0[3] * Q0[0];
108 tmp_Q0Q0[3][1] = Q0[3] * Q0[1];
109 tmp_Q0Q0[3][2] = Q0[3] * Q0[2];
110 tmp_Q0Q0[3][3] = Q0[3] * Q0[3];
111 }
112 if (require_dl_dq & rotation_derivative_dldq::use_dq) {
113 const auto &Q0 = m_rot.S_eigvec[0];
114 const auto &Q1 = m_rot.S_eigvec[1];
115 const auto &Q2 = m_rot.S_eigvec[2];
116 const auto &Q3 = m_rot.S_eigvec[3];
117 cvm::real const L0 = m_rot.S_eigval[0];
118 cvm::real const L1 = m_rot.S_eigval[1];
119 cvm::real const L2 = m_rot.S_eigval[2];
120 cvm::real const L3 = m_rot.S_eigval[3];
121
122 tmp_Q0Q0_L[0][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[0] +
123 (Q2[0] * Q0[0]) / (L0-L2) * Q2[0] +
124 (Q3[0] * Q0[0]) / (L0-L3) * Q3[0];
125 tmp_Q0Q0_L[1][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[1] +
126 (Q2[0] * Q0[0]) / (L0-L2) * Q2[1] +
127 (Q3[0] * Q0[0]) / (L0-L3) * Q3[1];
128 tmp_Q0Q0_L[2][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[2] +
129 (Q2[0] * Q0[0]) / (L0-L2) * Q2[2] +
130 (Q3[0] * Q0[0]) / (L0-L3) * Q3[2];
131 tmp_Q0Q0_L[3][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[3] +
132 (Q2[0] * Q0[0]) / (L0-L2) * Q2[3] +
133 (Q3[0] * Q0[0]) / (L0-L3) * Q3[3];
134
135 tmp_Q0Q0_L[0][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[0] +
136 (Q2[0] * Q0[1]) / (L0-L2) * Q2[0] +
137 (Q3[0] * Q0[1]) / (L0-L3) * Q3[0];
138 tmp_Q0Q0_L[1][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[1] +
139 (Q2[0] * Q0[1]) / (L0-L2) * Q2[1] +
140 (Q3[0] * Q0[1]) / (L0-L3) * Q3[1];
141 tmp_Q0Q0_L[2][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[2] +
142 (Q2[0] * Q0[1]) / (L0-L2) * Q2[2] +
143 (Q3[0] * Q0[1]) / (L0-L3) * Q3[2];
144 tmp_Q0Q0_L[3][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[3] +
145 (Q2[0] * Q0[1]) / (L0-L2) * Q2[3] +
146 (Q3[0] * Q0[1]) / (L0-L3) * Q3[3];
147
148
149 tmp_Q0Q0_L[0][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[0] +
150 (Q2[0] * Q0[2]) / (L0-L2) * Q2[0] +
151 (Q3[0] * Q0[2]) / (L0-L3) * Q3[0];
152 tmp_Q0Q0_L[1][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[1] +
153 (Q2[0] * Q0[2]) / (L0-L2) * Q2[1] +
154 (Q3[0] * Q0[2]) / (L0-L3) * Q3[1];
155 tmp_Q0Q0_L[2][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[2] +
156 (Q2[0] * Q0[2]) / (L0-L2) * Q2[2] +
157 (Q3[0] * Q0[2]) / (L0-L3) * Q3[2];
158 tmp_Q0Q0_L[3][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[3] +
159 (Q2[0] * Q0[2]) / (L0-L2) * Q2[3] +
160 (Q3[0] * Q0[2]) / (L0-L3) * Q3[3];
161
162 tmp_Q0Q0_L[0][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[0] +
163 (Q2[0] * Q0[3]) / (L0-L2) * Q2[0] +
164 (Q3[0] * Q0[3]) / (L0-L3) * Q3[0];
165 tmp_Q0Q0_L[1][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[1] +
166 (Q2[0] * Q0[3]) / (L0-L2) * Q2[1] +
167 (Q3[0] * Q0[3]) / (L0-L3) * Q3[1];
168 tmp_Q0Q0_L[2][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[2] +
169 (Q2[0] * Q0[3]) / (L0-L2) * Q2[2] +
170 (Q3[0] * Q0[3]) / (L0-L3) * Q3[2];
171 tmp_Q0Q0_L[3][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[3] +
172 (Q2[0] * Q0[3]) / (L0-L2) * Q2[3] +
173 (Q3[0] * Q0[3]) / (L0-L3) * Q3[3];
174
175 tmp_Q0Q0_L[0][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[0] +
176 (Q2[1] * Q0[0]) / (L0-L2) * Q2[0] +
177 (Q3[1] * Q0[0]) / (L0-L3) * Q3[0];
178 tmp_Q0Q0_L[1][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[1] +
179 (Q2[1] * Q0[0]) / (L0-L2) * Q2[1] +
180 (Q3[1] * Q0[0]) / (L0-L3) * Q3[1];
181 tmp_Q0Q0_L[2][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[2] +
182 (Q2[1] * Q0[0]) / (L0-L2) * Q2[2] +
183 (Q3[1] * Q0[0]) / (L0-L3) * Q3[2];
184 tmp_Q0Q0_L[3][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[3] +
185 (Q2[1] * Q0[0]) / (L0-L2) * Q2[3] +
186 (Q3[1] * Q0[0]) / (L0-L3) * Q3[3];
187
188 tmp_Q0Q0_L[0][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[0] +
189 (Q2[1] * Q0[1]) / (L0-L2) * Q2[0] +
190 (Q3[1] * Q0[1]) / (L0-L3) * Q3[0];
191 tmp_Q0Q0_L[1][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[1] +
192 (Q2[1] * Q0[1]) / (L0-L2) * Q2[1] +
193 (Q3[1] * Q0[1]) / (L0-L3) * Q3[1];
194 tmp_Q0Q0_L[2][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[2] +
195 (Q2[1] * Q0[1]) / (L0-L2) * Q2[2] +
196 (Q3[1] * Q0[1]) / (L0-L3) * Q3[2];
197 tmp_Q0Q0_L[3][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[3] +
198 (Q2[1] * Q0[1]) / (L0-L2) * Q2[3] +
199 (Q3[1] * Q0[1]) / (L0-L3) * Q3[3];
200
201 tmp_Q0Q0_L[0][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[0] +
202 (Q2[1] * Q0[2]) / (L0-L2) * Q2[0] +
203 (Q3[1] * Q0[2]) / (L0-L3) * Q3[0];
204 tmp_Q0Q0_L[1][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[1] +
205 (Q2[1] * Q0[2]) / (L0-L2) * Q2[1] +
206 (Q3[1] * Q0[2]) / (L0-L3) * Q3[1];
207 tmp_Q0Q0_L[2][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[2] +
208 (Q2[1] * Q0[2]) / (L0-L2) * Q2[2] +
209 (Q3[1] * Q0[2]) / (L0-L3) * Q3[2];
210 tmp_Q0Q0_L[3][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[3] +
211 (Q2[1] * Q0[2]) / (L0-L2) * Q2[3] +
212 (Q3[1] * Q0[2]) / (L0-L3) * Q3[3];
213
214 tmp_Q0Q0_L[0][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[0] +
215 (Q2[1] * Q0[3]) / (L0-L2) * Q2[0] +
216 (Q3[1] * Q0[3]) / (L0-L3) * Q3[0];
217 tmp_Q0Q0_L[1][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[1] +
218 (Q2[1] * Q0[3]) / (L0-L2) * Q2[1] +
219 (Q3[1] * Q0[3]) / (L0-L3) * Q3[1];
220 tmp_Q0Q0_L[2][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[2] +
221 (Q2[1] * Q0[3]) / (L0-L2) * Q2[2] +
222 (Q3[1] * Q0[3]) / (L0-L3) * Q3[2];
223 tmp_Q0Q0_L[3][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[3] +
224 (Q2[1] * Q0[3]) / (L0-L2) * Q2[3] +
225 (Q3[1] * Q0[3]) / (L0-L3) * Q3[3];
226
227
228 tmp_Q0Q0_L[0][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[0] +
229 (Q2[2] * Q0[0]) / (L0-L2) * Q2[0] +
230 (Q3[2] * Q0[0]) / (L0-L3) * Q3[0];
231 tmp_Q0Q0_L[1][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[1] +
232 (Q2[2] * Q0[0]) / (L0-L2) * Q2[1] +
233 (Q3[2] * Q0[0]) / (L0-L3) * Q3[1];
234 tmp_Q0Q0_L[2][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[2] +
235 (Q2[2] * Q0[0]) / (L0-L2) * Q2[2] +
236 (Q3[2] * Q0[0]) / (L0-L3) * Q3[2];
237 tmp_Q0Q0_L[3][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[3] +
238 (Q2[2] * Q0[0]) / (L0-L2) * Q2[3] +
239 (Q3[2] * Q0[0]) / (L0-L3) * Q3[3];
240
241 tmp_Q0Q0_L[0][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[0] +
242 (Q2[2] * Q0[1]) / (L0-L2) * Q2[0] +
243 (Q3[2] * Q0[1]) / (L0-L3) * Q3[0];
244 tmp_Q0Q0_L[1][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[1] +
245 (Q2[2] * Q0[1]) / (L0-L2) * Q2[1] +
246 (Q3[2] * Q0[1]) / (L0-L3) * Q3[1];
247 tmp_Q0Q0_L[2][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[2] +
248 (Q2[2] * Q0[1]) / (L0-L2) * Q2[2] +
249 (Q3[2] * Q0[1]) / (L0-L3) * Q3[2];
250 tmp_Q0Q0_L[3][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[3] +
251 (Q2[2] * Q0[1]) / (L0-L2) * Q2[3] +
252 (Q3[2] * Q0[1]) / (L0-L3) * Q3[3];
253
254 tmp_Q0Q0_L[0][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[0] +
255 (Q2[2] * Q0[2]) / (L0-L2) * Q2[0] +
256 (Q3[2] * Q0[2]) / (L0-L3) * Q3[0];
257 tmp_Q0Q0_L[1][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[1] +
258 (Q2[2] * Q0[2]) / (L0-L2) * Q2[1] +
259 (Q3[2] * Q0[2]) / (L0-L3) * Q3[1];
260 tmp_Q0Q0_L[2][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[2] +
261 (Q2[2] * Q0[2]) / (L0-L2) * Q2[2] +
262 (Q3[2] * Q0[2]) / (L0-L3) * Q3[2];
263 tmp_Q0Q0_L[3][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[3] +
264 (Q2[2] * Q0[2]) / (L0-L2) * Q2[3] +
265 (Q3[2] * Q0[2]) / (L0-L3) * Q3[3];
266
267 tmp_Q0Q0_L[0][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[0] +
268 (Q2[2] * Q0[3]) / (L0-L2) * Q2[0] +
269 (Q3[2] * Q0[3]) / (L0-L3) * Q3[0];
270 tmp_Q0Q0_L[1][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[1] +
271 (Q2[2] * Q0[3]) / (L0-L2) * Q2[1] +
272 (Q3[2] * Q0[3]) / (L0-L3) * Q3[1];
273 tmp_Q0Q0_L[2][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[2] +
274 (Q2[2] * Q0[3]) / (L0-L2) * Q2[2] +
275 (Q3[2] * Q0[3]) / (L0-L3) * Q3[2];
276 tmp_Q0Q0_L[3][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[3] +
277 (Q2[2] * Q0[3]) / (L0-L2) * Q2[3] +
278 (Q3[2] * Q0[3]) / (L0-L3) * Q3[3];
279
280 tmp_Q0Q0_L[0][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[0] +
281 (Q2[3] * Q0[0]) / (L0-L2) * Q2[0] +
282 (Q3[3] * Q0[0]) / (L0-L3) * Q3[0];
283 tmp_Q0Q0_L[1][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[1] +
284 (Q2[3] * Q0[0]) / (L0-L2) * Q2[1] +
285 (Q3[3] * Q0[0]) / (L0-L3) * Q3[1];
286 tmp_Q0Q0_L[2][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[2] +
287 (Q2[3] * Q0[0]) / (L0-L2) * Q2[2] +
288 (Q3[3] * Q0[0]) / (L0-L3) * Q3[2];
289 tmp_Q0Q0_L[3][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[3] +
290 (Q2[3] * Q0[0]) / (L0-L2) * Q2[3] +
291 (Q3[3] * Q0[0]) / (L0-L3) * Q3[3];
292
293 tmp_Q0Q0_L[0][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[0] +
294 (Q2[3] * Q0[1]) / (L0-L2) * Q2[0] +
295 (Q3[3] * Q0[1]) / (L0-L3) * Q3[0];
296 tmp_Q0Q0_L[1][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[1] +
297 (Q2[3] * Q0[1]) / (L0-L2) * Q2[1] +
298 (Q3[3] * Q0[1]) / (L0-L3) * Q3[1];
299 tmp_Q0Q0_L[2][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[2] +
300 (Q2[3] * Q0[1]) / (L0-L2) * Q2[2] +
301 (Q3[3] * Q0[1]) / (L0-L3) * Q3[2];
302 tmp_Q0Q0_L[3][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[3] +
303 (Q2[3] * Q0[1]) / (L0-L2) * Q2[3] +
304 (Q3[3] * Q0[1]) / (L0-L3) * Q3[3];
305
306 tmp_Q0Q0_L[0][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[0] +
307 (Q2[3] * Q0[2]) / (L0-L2) * Q2[0] +
308 (Q3[3] * Q0[2]) / (L0-L3) * Q3[0];
309 tmp_Q0Q0_L[1][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[1] +
310 (Q2[3] * Q0[2]) / (L0-L2) * Q2[1] +
311 (Q3[3] * Q0[2]) / (L0-L3) * Q3[1];
312 tmp_Q0Q0_L[2][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[2] +
313 (Q2[3] * Q0[2]) / (L0-L2) * Q2[2] +
314 (Q3[3] * Q0[2]) / (L0-L3) * Q3[2];
315 tmp_Q0Q0_L[3][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[3] +
316 (Q2[3] * Q0[2]) / (L0-L2) * Q2[3] +
317 (Q3[3] * Q0[2]) / (L0-L3) * Q3[3];
318
319 tmp_Q0Q0_L[0][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[0] +
320 (Q2[3] * Q0[3]) / (L0-L2) * Q2[0] +
321 (Q3[3] * Q0[3]) / (L0-L3) * Q3[0];
322 tmp_Q0Q0_L[1][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[1] +
323 (Q2[3] * Q0[3]) / (L0-L2) * Q2[1] +
324 (Q3[3] * Q0[3]) / (L0-L3) * Q3[1];
325 tmp_Q0Q0_L[2][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[2] +
326 (Q2[3] * Q0[3]) / (L0-L2) * Q2[2] +
327 (Q3[3] * Q0[3]) / (L0-L3) * Q3[2];
328 tmp_Q0Q0_L[3][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[3] +
329 (Q2[3] * Q0[3]) / (L0-L2) * Q2[3] +
330 (Q3[3] * Q0[3]) / (L0-L3) * Q3[3];
331 }
332 }
340 template <bool use_dl, bool use_dq, bool use_ds>
341 void calc_derivative_impl(
342 const cvm::rvector (&ds)[4][4],
343 cvm::rvector* _noalias const dl0_out,
344 cvm::vector1d<cvm::rvector>* _noalias const dq0_out,
345 cvm::matrix2d<cvm::rvector>* _noalias const ds_out) const {
346 if (use_ds) {
347 // this code path is for debug_gradients, so not necessary to unroll the loop
348 *ds_out = cvm::matrix2d<cvm::rvector>(4, 4);
349 for (int i = 0; i < 4; ++i) {
350 for (int j = 0; j < 4; ++j) {
351 (*ds_out)[i][j] = ds[i][j];
352 }
353 }
354 }
355 if (use_dl) {
356 /* manually loop unrolling of the following loop:
357 dl0_1.reset();
358 for (size_t i = 0; i < 4; i++) {
359 for (size_t j = 0; j < 4; j++) {
360 dl0_1 += Q0[i] * ds_1[i][j] * Q0[j];
361 }
362 }
363 */
364 *dl0_out = tmp_Q0Q0[0][0] * ds[0][0] +
365 tmp_Q0Q0[0][1] * ds[0][1] +
366 tmp_Q0Q0[0][2] * ds[0][2] +
367 tmp_Q0Q0[0][3] * ds[0][3] +
368 tmp_Q0Q0[1][0] * ds[1][0] +
369 tmp_Q0Q0[1][1] * ds[1][1] +
370 tmp_Q0Q0[1][2] * ds[1][2] +
371 tmp_Q0Q0[1][3] * ds[1][3] +
372 tmp_Q0Q0[2][0] * ds[2][0] +
373 tmp_Q0Q0[2][1] * ds[2][1] +
374 tmp_Q0Q0[2][2] * ds[2][2] +
375 tmp_Q0Q0[2][3] * ds[2][3] +
376 tmp_Q0Q0[3][0] * ds[3][0] +
377 tmp_Q0Q0[3][1] * ds[3][1] +
378 tmp_Q0Q0[3][2] * ds[3][2] +
379 tmp_Q0Q0[3][3] * ds[3][3];
380 }
381 if (use_dq) {
382 // we can skip this check if a fixed-size array is used
383 if (dq0_out->size() != 4) dq0_out->resize(4);
384 /* manually loop unrolling of the following loop:
385 dq0_1.reset();
386 for (size_t p = 0; p < 4; p++) {
387 for (size_t i = 0; i < 4; i++) {
388 for (size_t j = 0; j < 4; j++) {
389 dq0_1[p] +=
390 (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
391 (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
392 (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p];
393 }
394 }
395 }
396 */
397 (*dq0_out)[0] = tmp_Q0Q0_L[0][0][0] * ds[0][0] +
398 tmp_Q0Q0_L[0][0][1] * ds[0][1] +
399 tmp_Q0Q0_L[0][0][2] * ds[0][2] +
400 tmp_Q0Q0_L[0][0][3] * ds[0][3] +
401 tmp_Q0Q0_L[0][1][0] * ds[1][0] +
402 tmp_Q0Q0_L[0][1][1] * ds[1][1] +
403 tmp_Q0Q0_L[0][1][2] * ds[1][2] +
404 tmp_Q0Q0_L[0][1][3] * ds[1][3] +
405 tmp_Q0Q0_L[0][2][0] * ds[2][0] +
406 tmp_Q0Q0_L[0][2][1] * ds[2][1] +
407 tmp_Q0Q0_L[0][2][2] * ds[2][2] +
408 tmp_Q0Q0_L[0][2][3] * ds[2][3] +
409 tmp_Q0Q0_L[0][3][0] * ds[3][0] +
410 tmp_Q0Q0_L[0][3][1] * ds[3][1] +
411 tmp_Q0Q0_L[0][3][2] * ds[3][2] +
412 tmp_Q0Q0_L[0][3][3] * ds[3][3];
413
414 (*dq0_out)[1] = tmp_Q0Q0_L[1][0][0] * ds[0][0] +
415 tmp_Q0Q0_L[1][0][1] * ds[0][1] +
416 tmp_Q0Q0_L[1][0][2] * ds[0][2] +
417 tmp_Q0Q0_L[1][0][3] * ds[0][3] +
418 tmp_Q0Q0_L[1][1][0] * ds[1][0] +
419 tmp_Q0Q0_L[1][1][1] * ds[1][1] +
420 tmp_Q0Q0_L[1][1][2] * ds[1][2] +
421 tmp_Q0Q0_L[1][1][3] * ds[1][3] +
422 tmp_Q0Q0_L[1][2][0] * ds[2][0] +
423 tmp_Q0Q0_L[1][2][1] * ds[2][1] +
424 tmp_Q0Q0_L[1][2][2] * ds[2][2] +
425 tmp_Q0Q0_L[1][2][3] * ds[2][3] +
426 tmp_Q0Q0_L[1][3][0] * ds[3][0] +
427 tmp_Q0Q0_L[1][3][1] * ds[3][1] +
428 tmp_Q0Q0_L[1][3][2] * ds[3][2] +
429 tmp_Q0Q0_L[1][3][3] * ds[3][3];
430
431 (*dq0_out)[2] = tmp_Q0Q0_L[2][0][0] * ds[0][0] +
432 tmp_Q0Q0_L[2][0][1] * ds[0][1] +
433 tmp_Q0Q0_L[2][0][2] * ds[0][2] +
434 tmp_Q0Q0_L[2][0][3] * ds[0][3] +
435 tmp_Q0Q0_L[2][1][0] * ds[1][0] +
436 tmp_Q0Q0_L[2][1][1] * ds[1][1] +
437 tmp_Q0Q0_L[2][1][2] * ds[1][2] +
438 tmp_Q0Q0_L[2][1][3] * ds[1][3] +
439 tmp_Q0Q0_L[2][2][0] * ds[2][0] +
440 tmp_Q0Q0_L[2][2][1] * ds[2][1] +
441 tmp_Q0Q0_L[2][2][2] * ds[2][2] +
442 tmp_Q0Q0_L[2][2][3] * ds[2][3] +
443 tmp_Q0Q0_L[2][3][0] * ds[3][0] +
444 tmp_Q0Q0_L[2][3][1] * ds[3][1] +
445 tmp_Q0Q0_L[2][3][2] * ds[3][2] +
446 tmp_Q0Q0_L[2][3][3] * ds[3][3];
447
448 (*dq0_out)[3] = tmp_Q0Q0_L[3][0][0] * ds[0][0] +
449 tmp_Q0Q0_L[3][0][1] * ds[0][1] +
450 tmp_Q0Q0_L[3][0][2] * ds[0][2] +
451 tmp_Q0Q0_L[3][0][3] * ds[0][3] +
452 tmp_Q0Q0_L[3][1][0] * ds[1][0] +
453 tmp_Q0Q0_L[3][1][1] * ds[1][1] +
454 tmp_Q0Q0_L[3][1][2] * ds[1][2] +
455 tmp_Q0Q0_L[3][1][3] * ds[1][3] +
456 tmp_Q0Q0_L[3][2][0] * ds[2][0] +
457 tmp_Q0Q0_L[3][2][1] * ds[2][1] +
458 tmp_Q0Q0_L[3][2][2] * ds[2][2] +
459 tmp_Q0Q0_L[3][2][3] * ds[2][3] +
460 tmp_Q0Q0_L[3][3][0] * ds[3][0] +
461 tmp_Q0Q0_L[3][3][1] * ds[3][1] +
462 tmp_Q0Q0_L[3][3][2] * ds[3][2] +
463 tmp_Q0Q0_L[3][3][3] * ds[3][3];
464 }
465 }
476 template <bool use_dl, bool use_dq, bool use_ds>
477 void calc_derivative_wrt_group1(
478 size_t ia, cvm::rvector* _noalias const dl0_1_out = nullptr,
479 cvm::vector1d<cvm::rvector>* _noalias const dq0_1_out = nullptr,
480 cvm::matrix2d<cvm::rvector>* _noalias const ds_1_out = nullptr) const {
481 // if (dl0_1_out == nullptr && dq0_1_out == nullptr) return;
482 cvm::real a2x, a2y, a2z;
483 // we can get rid of the helper function read_atom_coord if C++17 (constexpr) is available
484 read_atom_coord(ia, m_pos2, &a2x, &a2y, &a2z);
485 const cvm::rvector ds_1[4][4] = {
486 {{ a2x, a2y, a2z}, { 0.0, a2z, -a2y}, {-a2z, 0.0, a2x}, { a2y, -a2x, 0.0}},
487 {{ 0.0, a2z, -a2y}, { a2x, -a2y, -a2z}, { a2y, a2x, 0.0}, { a2z, 0.0, a2x}},
488 {{-a2z, 0.0, a2x}, { a2y, a2x, 0.0}, {-a2x, a2y, -a2z}, { 0.0, a2z, a2y}},
489 {{ a2y, -a2x, 0.0}, { a2z, 0.0, a2x}, { 0.0, a2z, a2y}, {-a2x, -a2y, a2z}}};
490 calc_derivative_impl<use_dl, use_dq, use_ds>(ds_1, dl0_1_out, dq0_1_out, ds_1_out);
491 }
502 template <bool use_dl, bool use_dq, bool use_ds>
503 void calc_derivative_wrt_group2(
504 size_t ia, cvm::rvector* _noalias const dl0_2_out = nullptr,
505 cvm::vector1d<cvm::rvector>* _noalias const dq0_2_out = nullptr,
506 cvm::matrix2d<cvm::rvector>* _noalias const ds_2_out = nullptr) const {
507 // if (dl0_2_out == nullptr && dq0_2_out == nullptr) return;
508 cvm::real a1x, a1y, a1z;
509 // we can get rid of the helper function read_atom_coord if C++17 (constexpr) is available
510 read_atom_coord(ia, m_pos1, &a1x, &a1y, &a1z);
511 const cvm::rvector ds_2[4][4] = {
512 {{ a1x, a1y, a1z}, { 0.0, -a1z, a1y}, { a1z, 0.0, -a1x}, {-a1y, a1x, 0.0}},
513 {{ 0.0, -a1z, a1y}, { a1x, -a1y, -a1z}, { a1y, a1x, 0.0}, { a1z, 0.0, a1x}},
514 {{ a1z, 0.0, -a1x}, { a1y, a1x, 0.0}, {-a1x, a1y, -a1z}, { 0.0, a1z, a1y}},
515 {{-a1y, a1x, 0.0}, { a1z, 0.0, a1x}, { 0.0, a1z, a1y}, {-a1x, -a1y, a1z}}};
516 calc_derivative_impl<use_dl, use_dq, use_ds>(ds_2, dl0_2_out, dq0_2_out, ds_2_out);
517 }
518};
519
526template<typename T1, typename T2>
527void debug_gradients(
528 cvm::rotation &rot,
529 const std::vector<T1> &pos1,
530 const std::vector<T2> &pos2) {
531 static_assert(std::is_same<T1, cvm::atom_pos>::value || std::is_same<T1, cvm::atom>::value, "");
532 static_assert(std::is_same<T2, cvm::atom_pos>::value || std::is_same<T2, cvm::atom>::value, "");
533 // eigenvalues and eigenvectors
534 cvm::real const L0 = rot.S_eigval[0];
535 cvm::real const L1 = rot.S_eigval[1];
536 cvm::real const L2 = rot.S_eigval[2];
537 cvm::real const L3 = rot.S_eigval[3];
538 cvm::quaternion const Q0(rot.S_eigvec[0]);
539 cvm::quaternion const Q1(rot.S_eigvec[1]);
540 cvm::quaternion const Q2(rot.S_eigvec[2]);
541 cvm::quaternion const Q3(rot.S_eigvec[3]);
542
543 cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
544 ", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
545 ", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
546 "\n");
547 cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
548 ", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
549 ", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
550 "\n");
551 cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
552 ", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
553 ", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
554 "\n");
555 cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
556 ", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
557 ", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
558 "\n");
559
560 rotation_derivative<T1, T2> deriv(rot, pos1, pos2);
561 cvm::rvector dl0_2;
562 cvm::vector1d<cvm::rvector> dq0_2(4);
563 cvm::matrix2d<cvm::rvector> ds_2;
564#ifdef COLVARS_LAMMPS
565 MathEigen::Jacobi<cvm::real,
566 cvm::real[4],
567 cvm::real[4][4]> *ecalc =
568 reinterpret_cast<MathEigen::Jacobi<cvm::real,
569 cvm::real[4],
570 cvm::real[4][4]> *>(rot.jacobi);
571#endif
572 deriv.prepare_derivative(rotation_derivative_dldq::use_dl | rotation_derivative_dldq::use_dq);
573 cvm::real S_new[4][4];
574 cvm::real S_new_eigval[4];
575 cvm::real S_new_eigvec[4][4];
576 for (size_t ia = 0; ia < pos2.size(); ++ia) {
577 deriv.template calc_derivative_wrt_group2<true, true, true>(ia, &dl0_2, &dq0_2, &ds_2);
578 // make an infitesimal move along each cartesian coordinate of
579 // this atom, and solve again the eigenvector problem
580 for (size_t comp = 0; comp < 3; comp++) {
581 std::memcpy(S_new, rot.S_backup, sizeof(cvm::real) * 4 * 4);
582 std::memset(S_new_eigval, 0, sizeof(cvm::real) * 4);
583 std::memset(S_new_eigvec, 0, sizeof(cvm::real) * 4 * 4);
584 for (size_t i = 0; i < 4; i++) {
585 for (size_t j = 0; j < 4; j++) {
586 S_new[i][j] +=
587 colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
588 }
589 }
590#ifdef COLVARS_LAMMPS
591 ecalc->Diagonalize(S_new, S_new_eigval, S_new_eigvec);
592#else
593 NR::diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
594#endif
595 cvm::real const &L0_new = S_new_eigval[0];
596 cvm::quaternion const Q0_new(S_new_eigvec[0]);
597
598 cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
599 cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
600 dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
601 dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
602 dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
603
604 cvm::log( "|(l_0+dl_0) - l_0^new|/l_0 = "+
605 cvm::to_str(cvm::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
606 ", |(q_0+dq_0) - q_0^new| = "+
607 cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
608 "\n");
609 }
610 }
611}
612
613#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:1341
colvarmodule::rotation
A rotation between two sets of coordinates (for the moment a wrapper for colvarmodule::quaternion)
Definition: colvartypes.h:1363
colvarmodule::rotation::S_eigval
cvm::real S_eigval[4]
Eigenvalues of S.
Definition: colvartypes.h:1372
colvarmodule::rotation::S_backup
cvm::real S_backup[4][4]
Used for debugging gradients.
Definition: colvartypes.h:1378
colvarmodule::rotation::S_eigvec
cvm::real S_eigvec[4][4]
Eigenvectors of S.
Definition: colvartypes.h:1375
colvarmodule::rotation::jacobi
void * jacobi
Pointer to instance of Jacobi solver.
Definition: colvartypes.h:1554
colvarmodule::rvector
vector of real numbers with three components
Definition: colvartypes.h:727
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:1969
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:2392
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:59
rotation_derivative::calc_derivative_wrt_group1
void calc_derivative_wrt_group1(size_t ia, cvm::rvector *_noalias const dl0_1_out=nullptr, cvm::vector1d< cvm::rvector > *_noalias const dq0_1_out=nullptr, cvm::matrix2d< cvm::rvector > *_noalias 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:477
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:92
rotation_derivative::m_pos1
const std::vector< T1 > & m_pos1
Reference to the atom positions of group 1.
Definition: colvar_rotation_derivative.h:67
rotation_derivative::m_rot
const cvm::rotation & m_rot
Reference to the rotation.
Definition: colvar_rotation_derivative.h:65
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:71
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:81
rotation_derivative::calc_derivative_wrt_group2
void calc_derivative_wrt_group2(size_t ia, cvm::rvector *_noalias const dl0_2_out=nullptr, cvm::vector1d< cvm::rvector > *_noalias const dq0_2_out=nullptr, cvm::matrix2d< cvm::rvector > *_noalias 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:503
rotation_derivative::m_pos2
const std::vector< T2 > & m_pos2
Reference to the atom positions of group 2.
Definition: colvar_rotation_derivative.h:69
rotation_derivative::calc_derivative_impl
void calc_derivative_impl(const cvm::rvector(&ds)[4][4], cvm::rvector *_noalias const dl0_out, cvm::vector1d< cvm::rvector > *_noalias const dq0_out, cvm::matrix2d< cvm::rvector > *_noalias const ds_out) const
Actual implementation of the derivative calculation.
Definition: colvar_rotation_derivative.h:341
Generated on Mon Dec 23 2024 17:18:10 for Collective Variables Module - Developer Documentation by doxygen 1.9.5