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#include <array>
8
9#ifndef _noalias
10#if defined(__INTEL_COMPILER) || (defined(__PGI) && !defined(__NVCOMPILER))
11#define _noalias restrict
12#elif defined(__GNUC__) || defined(__INTEL_LLVM_COMPILER) || defined(__NVCOMPILER)
13#define _noalias __restrict
14#else
15#define _noalias
16#endif
17#endif
18
20enum class rotation_derivative_dldq {
22 use_dl = 1 << 0,
24 use_dq = 1 << 1
25};
26
27inline constexpr rotation_derivative_dldq operator|(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs) {
28 return static_cast<rotation_derivative_dldq>(
29 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Lhs) |
30 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Rhs));
31}
32
33inline constexpr bool operator&(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs)
34{
35 return (static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Lhs) &
36 static_cast<std::underlying_type<rotation_derivative_dldq>::type>(Rhs));
37}
38
40// template <typename T1, typename T2, bool soa = false>
41struct rotation_derivative {
43 const cvm::rotation &m_rot;
45 const std::vector<cvm::real> &m_pos1;
47 const std::vector<cvm::real> &m_pos2;
49 size_t m_num_atoms_pos1;
51 size_t m_num_atoms_pos2;
53 cvm::real tmp_Q0Q0[4][4];
54 cvm::real tmp_Q0Q0_L[4][4][4];
65 rotation_derivative(
66 const cvm::rotation &rot,
67 const std::vector<cvm::real> &pos1,
68 const std::vector<cvm::real> &pos2,
69 const size_t num_atoms_pos1,
70 const size_t num_atoms_pos2):
71 m_rot(rot), m_pos1(pos1), m_pos2(pos2),
72 m_num_atoms_pos1(num_atoms_pos1),
73 m_num_atoms_pos2(num_atoms_pos2) {}
80 void prepare_derivative(rotation_derivative_dldq require_dl_dq) {
81 if (require_dl_dq & rotation_derivative_dldq::use_dl) {
82 const auto &Q0 = m_rot.S_eigvec[0];
83 tmp_Q0Q0[0][0] = Q0[0] * Q0[0];
84 tmp_Q0Q0[0][1] = Q0[0] * Q0[1];
85 tmp_Q0Q0[0][2] = Q0[0] * Q0[2];
86 tmp_Q0Q0[0][3] = Q0[0] * Q0[3];
87 tmp_Q0Q0[1][0] = Q0[1] * Q0[0];
88 tmp_Q0Q0[1][1] = Q0[1] * Q0[1];
89 tmp_Q0Q0[1][2] = Q0[1] * Q0[2];
90 tmp_Q0Q0[1][3] = Q0[1] * Q0[3];
91 tmp_Q0Q0[2][0] = Q0[2] * Q0[0];
92 tmp_Q0Q0[2][1] = Q0[2] * Q0[1];
93 tmp_Q0Q0[2][2] = Q0[2] * Q0[2];
94 tmp_Q0Q0[2][3] = Q0[2] * Q0[3];
95 tmp_Q0Q0[3][0] = Q0[3] * Q0[0];
96 tmp_Q0Q0[3][1] = Q0[3] * Q0[1];
97 tmp_Q0Q0[3][2] = Q0[3] * Q0[2];
98 tmp_Q0Q0[3][3] = Q0[3] * Q0[3];
99 }
100 if (require_dl_dq & rotation_derivative_dldq::use_dq) {
101 const auto &Q0 = m_rot.S_eigvec[0];
102 const auto &Q1 = m_rot.S_eigvec[1];
103 const auto &Q2 = m_rot.S_eigvec[2];
104 const auto &Q3 = m_rot.S_eigvec[3];
105 cvm::real const L0 = m_rot.S_eigval[0];
106 cvm::real const L1 = m_rot.S_eigval[1];
107 cvm::real const L2 = m_rot.S_eigval[2];
108 cvm::real const L3 = m_rot.S_eigval[3];
109
110 tmp_Q0Q0_L[0][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[0] +
111 (Q2[0] * Q0[0]) / (L0-L2) * Q2[0] +
112 (Q3[0] * Q0[0]) / (L0-L3) * Q3[0];
113 tmp_Q0Q0_L[1][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[1] +
114 (Q2[0] * Q0[0]) / (L0-L2) * Q2[1] +
115 (Q3[0] * Q0[0]) / (L0-L3) * Q3[1];
116 tmp_Q0Q0_L[2][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[2] +
117 (Q2[0] * Q0[0]) / (L0-L2) * Q2[2] +
118 (Q3[0] * Q0[0]) / (L0-L3) * Q3[2];
119 tmp_Q0Q0_L[3][0][0] = (Q1[0] * Q0[0]) / (L0-L1) * Q1[3] +
120 (Q2[0] * Q0[0]) / (L0-L2) * Q2[3] +
121 (Q3[0] * Q0[0]) / (L0-L3) * Q3[3];
122
123 tmp_Q0Q0_L[0][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[0] +
124 (Q2[0] * Q0[1]) / (L0-L2) * Q2[0] +
125 (Q3[0] * Q0[1]) / (L0-L3) * Q3[0];
126 tmp_Q0Q0_L[1][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[1] +
127 (Q2[0] * Q0[1]) / (L0-L2) * Q2[1] +
128 (Q3[0] * Q0[1]) / (L0-L3) * Q3[1];
129 tmp_Q0Q0_L[2][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[2] +
130 (Q2[0] * Q0[1]) / (L0-L2) * Q2[2] +
131 (Q3[0] * Q0[1]) / (L0-L3) * Q3[2];
132 tmp_Q0Q0_L[3][0][1] = (Q1[0] * Q0[1]) / (L0-L1) * Q1[3] +
133 (Q2[0] * Q0[1]) / (L0-L2) * Q2[3] +
134 (Q3[0] * Q0[1]) / (L0-L3) * Q3[3];
135
136
137 tmp_Q0Q0_L[0][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[0] +
138 (Q2[0] * Q0[2]) / (L0-L2) * Q2[0] +
139 (Q3[0] * Q0[2]) / (L0-L3) * Q3[0];
140 tmp_Q0Q0_L[1][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[1] +
141 (Q2[0] * Q0[2]) / (L0-L2) * Q2[1] +
142 (Q3[0] * Q0[2]) / (L0-L3) * Q3[1];
143 tmp_Q0Q0_L[2][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[2] +
144 (Q2[0] * Q0[2]) / (L0-L2) * Q2[2] +
145 (Q3[0] * Q0[2]) / (L0-L3) * Q3[2];
146 tmp_Q0Q0_L[3][0][2] = (Q1[0] * Q0[2]) / (L0-L1) * Q1[3] +
147 (Q2[0] * Q0[2]) / (L0-L2) * Q2[3] +
148 (Q3[0] * Q0[2]) / (L0-L3) * Q3[3];
149
150 tmp_Q0Q0_L[0][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[0] +
151 (Q2[0] * Q0[3]) / (L0-L2) * Q2[0] +
152 (Q3[0] * Q0[3]) / (L0-L3) * Q3[0];
153 tmp_Q0Q0_L[1][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[1] +
154 (Q2[0] * Q0[3]) / (L0-L2) * Q2[1] +
155 (Q3[0] * Q0[3]) / (L0-L3) * Q3[1];
156 tmp_Q0Q0_L[2][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[2] +
157 (Q2[0] * Q0[3]) / (L0-L2) * Q2[2] +
158 (Q3[0] * Q0[3]) / (L0-L3) * Q3[2];
159 tmp_Q0Q0_L[3][0][3] = (Q1[0] * Q0[3]) / (L0-L1) * Q1[3] +
160 (Q2[0] * Q0[3]) / (L0-L2) * Q2[3] +
161 (Q3[0] * Q0[3]) / (L0-L3) * Q3[3];
162
163 tmp_Q0Q0_L[0][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[0] +
164 (Q2[1] * Q0[0]) / (L0-L2) * Q2[0] +
165 (Q3[1] * Q0[0]) / (L0-L3) * Q3[0];
166 tmp_Q0Q0_L[1][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[1] +
167 (Q2[1] * Q0[0]) / (L0-L2) * Q2[1] +
168 (Q3[1] * Q0[0]) / (L0-L3) * Q3[1];
169 tmp_Q0Q0_L[2][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[2] +
170 (Q2[1] * Q0[0]) / (L0-L2) * Q2[2] +
171 (Q3[1] * Q0[0]) / (L0-L3) * Q3[2];
172 tmp_Q0Q0_L[3][1][0] = (Q1[1] * Q0[0]) / (L0-L1) * Q1[3] +
173 (Q2[1] * Q0[0]) / (L0-L2) * Q2[3] +
174 (Q3[1] * Q0[0]) / (L0-L3) * Q3[3];
175
176 tmp_Q0Q0_L[0][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[0] +
177 (Q2[1] * Q0[1]) / (L0-L2) * Q2[0] +
178 (Q3[1] * Q0[1]) / (L0-L3) * Q3[0];
179 tmp_Q0Q0_L[1][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[1] +
180 (Q2[1] * Q0[1]) / (L0-L2) * Q2[1] +
181 (Q3[1] * Q0[1]) / (L0-L3) * Q3[1];
182 tmp_Q0Q0_L[2][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[2] +
183 (Q2[1] * Q0[1]) / (L0-L2) * Q2[2] +
184 (Q3[1] * Q0[1]) / (L0-L3) * Q3[2];
185 tmp_Q0Q0_L[3][1][1] = (Q1[1] * Q0[1]) / (L0-L1) * Q1[3] +
186 (Q2[1] * Q0[1]) / (L0-L2) * Q2[3] +
187 (Q3[1] * Q0[1]) / (L0-L3) * Q3[3];
188
189 tmp_Q0Q0_L[0][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[0] +
190 (Q2[1] * Q0[2]) / (L0-L2) * Q2[0] +
191 (Q3[1] * Q0[2]) / (L0-L3) * Q3[0];
192 tmp_Q0Q0_L[1][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[1] +
193 (Q2[1] * Q0[2]) / (L0-L2) * Q2[1] +
194 (Q3[1] * Q0[2]) / (L0-L3) * Q3[1];
195 tmp_Q0Q0_L[2][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[2] +
196 (Q2[1] * Q0[2]) / (L0-L2) * Q2[2] +
197 (Q3[1] * Q0[2]) / (L0-L3) * Q3[2];
198 tmp_Q0Q0_L[3][1][2] = (Q1[1] * Q0[2]) / (L0-L1) * Q1[3] +
199 (Q2[1] * Q0[2]) / (L0-L2) * Q2[3] +
200 (Q3[1] * Q0[2]) / (L0-L3) * Q3[3];
201
202 tmp_Q0Q0_L[0][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[0] +
203 (Q2[1] * Q0[3]) / (L0-L2) * Q2[0] +
204 (Q3[1] * Q0[3]) / (L0-L3) * Q3[0];
205 tmp_Q0Q0_L[1][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[1] +
206 (Q2[1] * Q0[3]) / (L0-L2) * Q2[1] +
207 (Q3[1] * Q0[3]) / (L0-L3) * Q3[1];
208 tmp_Q0Q0_L[2][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[2] +
209 (Q2[1] * Q0[3]) / (L0-L2) * Q2[2] +
210 (Q3[1] * Q0[3]) / (L0-L3) * Q3[2];
211 tmp_Q0Q0_L[3][1][3] = (Q1[1] * Q0[3]) / (L0-L1) * Q1[3] +
212 (Q2[1] * Q0[3]) / (L0-L2) * Q2[3] +
213 (Q3[1] * Q0[3]) / (L0-L3) * Q3[3];
214
215
216 tmp_Q0Q0_L[0][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[0] +
217 (Q2[2] * Q0[0]) / (L0-L2) * Q2[0] +
218 (Q3[2] * Q0[0]) / (L0-L3) * Q3[0];
219 tmp_Q0Q0_L[1][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[1] +
220 (Q2[2] * Q0[0]) / (L0-L2) * Q2[1] +
221 (Q3[2] * Q0[0]) / (L0-L3) * Q3[1];
222 tmp_Q0Q0_L[2][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[2] +
223 (Q2[2] * Q0[0]) / (L0-L2) * Q2[2] +
224 (Q3[2] * Q0[0]) / (L0-L3) * Q3[2];
225 tmp_Q0Q0_L[3][2][0] = (Q1[2] * Q0[0]) / (L0-L1) * Q1[3] +
226 (Q2[2] * Q0[0]) / (L0-L2) * Q2[3] +
227 (Q3[2] * Q0[0]) / (L0-L3) * Q3[3];
228
229 tmp_Q0Q0_L[0][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[0] +
230 (Q2[2] * Q0[1]) / (L0-L2) * Q2[0] +
231 (Q3[2] * Q0[1]) / (L0-L3) * Q3[0];
232 tmp_Q0Q0_L[1][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[1] +
233 (Q2[2] * Q0[1]) / (L0-L2) * Q2[1] +
234 (Q3[2] * Q0[1]) / (L0-L3) * Q3[1];
235 tmp_Q0Q0_L[2][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[2] +
236 (Q2[2] * Q0[1]) / (L0-L2) * Q2[2] +
237 (Q3[2] * Q0[1]) / (L0-L3) * Q3[2];
238 tmp_Q0Q0_L[3][2][1] = (Q1[2] * Q0[1]) / (L0-L1) * Q1[3] +
239 (Q2[2] * Q0[1]) / (L0-L2) * Q2[3] +
240 (Q3[2] * Q0[1]) / (L0-L3) * Q3[3];
241
242 tmp_Q0Q0_L[0][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[0] +
243 (Q2[2] * Q0[2]) / (L0-L2) * Q2[0] +
244 (Q3[2] * Q0[2]) / (L0-L3) * Q3[0];
245 tmp_Q0Q0_L[1][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[1] +
246 (Q2[2] * Q0[2]) / (L0-L2) * Q2[1] +
247 (Q3[2] * Q0[2]) / (L0-L3) * Q3[1];
248 tmp_Q0Q0_L[2][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[2] +
249 (Q2[2] * Q0[2]) / (L0-L2) * Q2[2] +
250 (Q3[2] * Q0[2]) / (L0-L3) * Q3[2];
251 tmp_Q0Q0_L[3][2][2] = (Q1[2] * Q0[2]) / (L0-L1) * Q1[3] +
252 (Q2[2] * Q0[2]) / (L0-L2) * Q2[3] +
253 (Q3[2] * Q0[2]) / (L0-L3) * Q3[3];
254
255 tmp_Q0Q0_L[0][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[0] +
256 (Q2[2] * Q0[3]) / (L0-L2) * Q2[0] +
257 (Q3[2] * Q0[3]) / (L0-L3) * Q3[0];
258 tmp_Q0Q0_L[1][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[1] +
259 (Q2[2] * Q0[3]) / (L0-L2) * Q2[1] +
260 (Q3[2] * Q0[3]) / (L0-L3) * Q3[1];
261 tmp_Q0Q0_L[2][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[2] +
262 (Q2[2] * Q0[3]) / (L0-L2) * Q2[2] +
263 (Q3[2] * Q0[3]) / (L0-L3) * Q3[2];
264 tmp_Q0Q0_L[3][2][3] = (Q1[2] * Q0[3]) / (L0-L1) * Q1[3] +
265 (Q2[2] * Q0[3]) / (L0-L2) * Q2[3] +
266 (Q3[2] * Q0[3]) / (L0-L3) * Q3[3];
267
268 tmp_Q0Q0_L[0][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[0] +
269 (Q2[3] * Q0[0]) / (L0-L2) * Q2[0] +
270 (Q3[3] * Q0[0]) / (L0-L3) * Q3[0];
271 tmp_Q0Q0_L[1][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[1] +
272 (Q2[3] * Q0[0]) / (L0-L2) * Q2[1] +
273 (Q3[3] * Q0[0]) / (L0-L3) * Q3[1];
274 tmp_Q0Q0_L[2][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[2] +
275 (Q2[3] * Q0[0]) / (L0-L2) * Q2[2] +
276 (Q3[3] * Q0[0]) / (L0-L3) * Q3[2];
277 tmp_Q0Q0_L[3][3][0] = (Q1[3] * Q0[0]) / (L0-L1) * Q1[3] +
278 (Q2[3] * Q0[0]) / (L0-L2) * Q2[3] +
279 (Q3[3] * Q0[0]) / (L0-L3) * Q3[3];
280
281 tmp_Q0Q0_L[0][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[0] +
282 (Q2[3] * Q0[1]) / (L0-L2) * Q2[0] +
283 (Q3[3] * Q0[1]) / (L0-L3) * Q3[0];
284 tmp_Q0Q0_L[1][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[1] +
285 (Q2[3] * Q0[1]) / (L0-L2) * Q2[1] +
286 (Q3[3] * Q0[1]) / (L0-L3) * Q3[1];
287 tmp_Q0Q0_L[2][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[2] +
288 (Q2[3] * Q0[1]) / (L0-L2) * Q2[2] +
289 (Q3[3] * Q0[1]) / (L0-L3) * Q3[2];
290 tmp_Q0Q0_L[3][3][1] = (Q1[3] * Q0[1]) / (L0-L1) * Q1[3] +
291 (Q2[3] * Q0[1]) / (L0-L2) * Q2[3] +
292 (Q3[3] * Q0[1]) / (L0-L3) * Q3[3];
293
294 tmp_Q0Q0_L[0][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[0] +
295 (Q2[3] * Q0[2]) / (L0-L2) * Q2[0] +
296 (Q3[3] * Q0[2]) / (L0-L3) * Q3[0];
297 tmp_Q0Q0_L[1][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[1] +
298 (Q2[3] * Q0[2]) / (L0-L2) * Q2[1] +
299 (Q3[3] * Q0[2]) / (L0-L3) * Q3[1];
300 tmp_Q0Q0_L[2][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[2] +
301 (Q2[3] * Q0[2]) / (L0-L2) * Q2[2] +
302 (Q3[3] * Q0[2]) / (L0-L3) * Q3[2];
303 tmp_Q0Q0_L[3][3][2] = (Q1[3] * Q0[2]) / (L0-L1) * Q1[3] +
304 (Q2[3] * Q0[2]) / (L0-L2) * Q2[3] +
305 (Q3[3] * Q0[2]) / (L0-L3) * Q3[3];
306
307 tmp_Q0Q0_L[0][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[0] +
308 (Q2[3] * Q0[3]) / (L0-L2) * Q2[0] +
309 (Q3[3] * Q0[3]) / (L0-L3) * Q3[0];
310 tmp_Q0Q0_L[1][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[1] +
311 (Q2[3] * Q0[3]) / (L0-L2) * Q2[1] +
312 (Q3[3] * Q0[3]) / (L0-L3) * Q3[1];
313 tmp_Q0Q0_L[2][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[2] +
314 (Q2[3] * Q0[3]) / (L0-L2) * Q2[2] +
315 (Q3[3] * Q0[3]) / (L0-L3) * Q3[2];
316 tmp_Q0Q0_L[3][3][3] = (Q1[3] * Q0[3]) / (L0-L1) * Q1[3] +
317 (Q2[3] * Q0[3]) / (L0-L2) * Q2[3] +
318 (Q3[3] * Q0[3]) / (L0-L3) * Q3[3];
319 }
320 }
328 template <bool use_dl, bool use_dq, bool use_ds>
329 void calc_derivative_impl(
330 const cvm::rvector (&ds)[4][4],
331 cvm::rvector* _noalias const dl0_out,
332 std::array<cvm::rvector, 4>* _noalias const dq0_out,
333 std::array<std::array<cvm::rvector, 4>, 4>* _noalias const ds_out) const {
334 if (use_ds) {
335 // this code path is for debug_gradients, so not necessary to unroll the loop
336 *ds_out = std::array<std::array<cvm::rvector, 4>, 4>();
337 for (int i = 0; i < 4; ++i) {
338 for (int j = 0; j < 4; ++j) {
339 (*ds_out)[i][j] = ds[i][j];
340 }
341 }
342 }
343 if (use_dl) {
344 /* manually loop unrolling of the following loop:
345 dl0_1.reset();
346 for (size_t i = 0; i < 4; i++) {
347 for (size_t j = 0; j < 4; j++) {
348 dl0_1 += Q0[i] * ds_1[i][j] * Q0[j];
349 }
350 }
351 */
352 *dl0_out = tmp_Q0Q0[0][0] * ds[0][0] +
353 tmp_Q0Q0[0][1] * ds[0][1] +
354 tmp_Q0Q0[0][2] * ds[0][2] +
355 tmp_Q0Q0[0][3] * ds[0][3] +
356 tmp_Q0Q0[1][0] * ds[1][0] +
357 tmp_Q0Q0[1][1] * ds[1][1] +
358 tmp_Q0Q0[1][2] * ds[1][2] +
359 tmp_Q0Q0[1][3] * ds[1][3] +
360 tmp_Q0Q0[2][0] * ds[2][0] +
361 tmp_Q0Q0[2][1] * ds[2][1] +
362 tmp_Q0Q0[2][2] * ds[2][2] +
363 tmp_Q0Q0[2][3] * ds[2][3] +
364 tmp_Q0Q0[3][0] * ds[3][0] +
365 tmp_Q0Q0[3][1] * ds[3][1] +
366 tmp_Q0Q0[3][2] * ds[3][2] +
367 tmp_Q0Q0[3][3] * ds[3][3];
368 }
369 if (use_dq) {
370 // we can skip this check if a fixed-size array is used
371 // if (dq0_out->size() != 4) dq0_out->resize(4);
372 /* manually loop unrolling of the following loop:
373 dq0_1.reset();
374 for (size_t p = 0; p < 4; p++) {
375 for (size_t i = 0; i < 4; i++) {
376 for (size_t j = 0; j < 4; j++) {
377 dq0_1[p] +=
378 (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
379 (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
380 (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p];
381 }
382 }
383 }
384 */
385 (*dq0_out)[0] = tmp_Q0Q0_L[0][0][0] * ds[0][0] +
386 tmp_Q0Q0_L[0][0][1] * ds[0][1] +
387 tmp_Q0Q0_L[0][0][2] * ds[0][2] +
388 tmp_Q0Q0_L[0][0][3] * ds[0][3] +
389 tmp_Q0Q0_L[0][1][0] * ds[1][0] +
390 tmp_Q0Q0_L[0][1][1] * ds[1][1] +
391 tmp_Q0Q0_L[0][1][2] * ds[1][2] +
392 tmp_Q0Q0_L[0][1][3] * ds[1][3] +
393 tmp_Q0Q0_L[0][2][0] * ds[2][0] +
394 tmp_Q0Q0_L[0][2][1] * ds[2][1] +
395 tmp_Q0Q0_L[0][2][2] * ds[2][2] +
396 tmp_Q0Q0_L[0][2][3] * ds[2][3] +
397 tmp_Q0Q0_L[0][3][0] * ds[3][0] +
398 tmp_Q0Q0_L[0][3][1] * ds[3][1] +
399 tmp_Q0Q0_L[0][3][2] * ds[3][2] +
400 tmp_Q0Q0_L[0][3][3] * ds[3][3];
401
402 (*dq0_out)[1] = tmp_Q0Q0_L[1][0][0] * ds[0][0] +
403 tmp_Q0Q0_L[1][0][1] * ds[0][1] +
404 tmp_Q0Q0_L[1][0][2] * ds[0][2] +
405 tmp_Q0Q0_L[1][0][3] * ds[0][3] +
406 tmp_Q0Q0_L[1][1][0] * ds[1][0] +
407 tmp_Q0Q0_L[1][1][1] * ds[1][1] +
408 tmp_Q0Q0_L[1][1][2] * ds[1][2] +
409 tmp_Q0Q0_L[1][1][3] * ds[1][3] +
410 tmp_Q0Q0_L[1][2][0] * ds[2][0] +
411 tmp_Q0Q0_L[1][2][1] * ds[2][1] +
412 tmp_Q0Q0_L[1][2][2] * ds[2][2] +
413 tmp_Q0Q0_L[1][2][3] * ds[2][3] +
414 tmp_Q0Q0_L[1][3][0] * ds[3][0] +
415 tmp_Q0Q0_L[1][3][1] * ds[3][1] +
416 tmp_Q0Q0_L[1][3][2] * ds[3][2] +
417 tmp_Q0Q0_L[1][3][3] * ds[3][3];
418
419 (*dq0_out)[2] = tmp_Q0Q0_L[2][0][0] * ds[0][0] +
420 tmp_Q0Q0_L[2][0][1] * ds[0][1] +
421 tmp_Q0Q0_L[2][0][2] * ds[0][2] +
422 tmp_Q0Q0_L[2][0][3] * ds[0][3] +
423 tmp_Q0Q0_L[2][1][0] * ds[1][0] +
424 tmp_Q0Q0_L[2][1][1] * ds[1][1] +
425 tmp_Q0Q0_L[2][1][2] * ds[1][2] +
426 tmp_Q0Q0_L[2][1][3] * ds[1][3] +
427 tmp_Q0Q0_L[2][2][0] * ds[2][0] +
428 tmp_Q0Q0_L[2][2][1] * ds[2][1] +
429 tmp_Q0Q0_L[2][2][2] * ds[2][2] +
430 tmp_Q0Q0_L[2][2][3] * ds[2][3] +
431 tmp_Q0Q0_L[2][3][0] * ds[3][0] +
432 tmp_Q0Q0_L[2][3][1] * ds[3][1] +
433 tmp_Q0Q0_L[2][3][2] * ds[3][2] +
434 tmp_Q0Q0_L[2][3][3] * ds[3][3];
435
436 (*dq0_out)[3] = tmp_Q0Q0_L[3][0][0] * ds[0][0] +
437 tmp_Q0Q0_L[3][0][1] * ds[0][1] +
438 tmp_Q0Q0_L[3][0][2] * ds[0][2] +
439 tmp_Q0Q0_L[3][0][3] * ds[0][3] +
440 tmp_Q0Q0_L[3][1][0] * ds[1][0] +
441 tmp_Q0Q0_L[3][1][1] * ds[1][1] +
442 tmp_Q0Q0_L[3][1][2] * ds[1][2] +
443 tmp_Q0Q0_L[3][1][3] * ds[1][3] +
444 tmp_Q0Q0_L[3][2][0] * ds[2][0] +
445 tmp_Q0Q0_L[3][2][1] * ds[2][1] +
446 tmp_Q0Q0_L[3][2][2] * ds[2][2] +
447 tmp_Q0Q0_L[3][2][3] * ds[2][3] +
448 tmp_Q0Q0_L[3][3][0] * ds[3][0] +
449 tmp_Q0Q0_L[3][3][1] * ds[3][1] +
450 tmp_Q0Q0_L[3][3][2] * ds[3][2] +
451 tmp_Q0Q0_L[3][3][3] * ds[3][3];
452 }
453 }
464 template <bool use_dl, bool use_dq, bool use_ds>
465 void calc_derivative_wrt_group1(
466 size_t ia, cvm::rvector* _noalias const dl0_1_out = nullptr,
467 std::array<cvm::rvector, 4>* _noalias const dq0_1_out = nullptr,
468 std::array<std::array<cvm::rvector, 4>, 4>* _noalias const ds_1_out = nullptr) const {
469 // if (dl0_1_out == nullptr && dq0_1_out == nullptr) return;
470 const cvm::real a2x = m_pos2[ia];
471 const cvm::real a2y = m_pos2[ia + m_num_atoms_pos2];
472 const cvm::real a2z = m_pos2[ia + 2 * m_num_atoms_pos2];
473 const cvm::rvector ds_1[4][4] = {
474 {{ a2x, a2y, a2z}, { 0.0, a2z, -a2y}, {-a2z, 0.0, a2x}, { a2y, -a2x, 0.0}},
475 {{ 0.0, a2z, -a2y}, { a2x, -a2y, -a2z}, { a2y, a2x, 0.0}, { a2z, 0.0, a2x}},
476 {{-a2z, 0.0, a2x}, { a2y, a2x, 0.0}, {-a2x, a2y, -a2z}, { 0.0, a2z, a2y}},
477 {{ a2y, -a2x, 0.0}, { a2z, 0.0, a2x}, { 0.0, a2z, a2y}, {-a2x, -a2y, a2z}}};
478 calc_derivative_impl<use_dl, use_dq, use_ds>(ds_1, dl0_1_out, dq0_1_out, ds_1_out);
479 }
490 template <bool use_dl, bool use_dq, bool use_ds>
491 void calc_derivative_wrt_group2(
492 size_t ia, cvm::rvector* _noalias const dl0_2_out = nullptr,
493 std::array<cvm::rvector, 4>* _noalias const dq0_2_out = nullptr,
494 std::array<std::array<cvm::rvector, 4>, 4>* _noalias const ds_2_out = nullptr) const {
495 // if (dl0_2_out == nullptr && dq0_2_out == nullptr) return;
496 const cvm::real a1x = m_pos1[ia];
497 const cvm::real a1y = m_pos1[ia + m_num_atoms_pos1];
498 const cvm::real a1z = m_pos1[ia + 2 * m_num_atoms_pos1];
499 const cvm::rvector ds_2[4][4] = {
500 {{ a1x, a1y, a1z}, { 0.0, -a1z, a1y}, { a1z, 0.0, -a1x}, {-a1y, a1x, 0.0}},
501 {{ 0.0, -a1z, a1y}, { a1x, -a1y, -a1z}, { a1y, a1x, 0.0}, { a1z, 0.0, a1x}},
502 {{ a1z, 0.0, -a1x}, { a1y, a1x, 0.0}, {-a1x, a1y, -a1z}, { 0.0, a1z, a1y}},
503 {{-a1y, a1x, 0.0}, { a1z, 0.0, a1x}, { 0.0, a1z, a1y}, {-a1x, -a1y, a1z}}};
504 calc_derivative_impl<use_dl, use_dq, use_ds>(ds_2, dl0_2_out, dq0_2_out, ds_2_out);
505 }
506};
507
508#endif // COLVAR_ROTATION_DERIVATIVE
colvarmodule::rotation
A rotation between two sets of coordinates (for the moment a wrapper for colvarmodule::quaternion)
Definition: colvartypes.h:1313
colvarmodule::rotation::S_eigval
cvm::real S_eigval[4]
Eigenvalues of S.
Definition: colvartypes.h:1322
colvarmodule::rotation::S_eigvec
cvm::real S_eigvec[4][4]
Eigenvectors of S.
Definition: colvartypes.h:1325
colvarmodule::rvector
vector of real numbers with three components
Definition: colvartypes.h:724
colvarmodule::real
double real
Defining an abstract real number allows to switch precision.
Definition: colvarmodule.h:139
rotation_derivative
Helper class for calculating the derivative of rotation.
Definition: colvar_rotation_derivative.h:41
rotation_derivative::rotation_derivative
rotation_derivative(const cvm::rotation &rot, const std::vector< cvm::real > &pos1, const std::vector< cvm::real > &pos2, const size_t num_atoms_pos1, const size_t num_atoms_pos2)
Constructor of the cvm::rotation::derivative class for SOA.
Definition: colvar_rotation_derivative.h:65
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:80
rotation_derivative::calc_derivative_wrt_group1
void calc_derivative_wrt_group1(size_t ia, cvm::rvector *_noalias const dl0_1_out=nullptr, std::array< cvm::rvector, 4 > *_noalias const dq0_1_out=nullptr, std::array< std::array< cvm::rvector, 4 >, 4 > *_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:465
rotation_derivative::m_num_atoms_pos1
size_t m_num_atoms_pos1
Number of atoms in group1 (used in SOA)
Definition: colvar_rotation_derivative.h:49
rotation_derivative::calc_derivative_impl
void calc_derivative_impl(const cvm::rvector(&ds)[4][4], cvm::rvector *_noalias const dl0_out, std::array< cvm::rvector, 4 > *_noalias const dq0_out, std::array< std::array< cvm::rvector, 4 >, 4 > *_noalias const ds_out) const
Actual implementation of the derivative calculation.
Definition: colvar_rotation_derivative.h:329
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:53
rotation_derivative::m_pos1
const std::vector< cvm::real > & m_pos1
Reference to the atom positions of group 1.
Definition: colvar_rotation_derivative.h:45
rotation_derivative::calc_derivative_wrt_group2
void calc_derivative_wrt_group2(size_t ia, cvm::rvector *_noalias const dl0_2_out=nullptr, std::array< cvm::rvector, 4 > *_noalias const dq0_2_out=nullptr, std::array< std::array< cvm::rvector, 4 >, 4 > *_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:491
rotation_derivative::m_rot
const cvm::rotation & m_rot
Reference to the rotation.
Definition: colvar_rotation_derivative.h:43
rotation_derivative::m_num_atoms_pos2
size_t m_num_atoms_pos2
Number of atoms in group1 (used in SOA)
Definition: colvar_rotation_derivative.h:51
rotation_derivative::m_pos2
const std::vector< cvm::real > & m_pos2
Reference to the atom positions of group 2.
Definition: colvar_rotation_derivative.h:47
Generated on Tue Sep 9 2025 17:38:10 for Collective Variables Module - Developer Documentation by doxygen 1.9.5