1#ifndef COLVAR_ROTATION_DERIVATIVE
2#define COLVAR_ROTATION_DERIVATIVE
4#include "colvartypes.h"
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
20enum class rotation_derivative_dldq {
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));
33inline constexpr bool operator&(rotation_derivative_dldq Lhs, rotation_derivative_dldq Rhs)
35 return (
static_cast<std::underlying_type<rotation_derivative_dldq>::type
>(Lhs) &
36 static_cast<std::underlying_type<rotation_derivative_dldq>::type
>(Rhs));
46 std::vector<cvm::real>::const_iterator
pos1x;
47 std::vector<cvm::real>::const_iterator pos1y;
48 std::vector<cvm::real>::const_iterator pos1z;
51 std::vector<cvm::real>::const_iterator
pos2x;
52 std::vector<cvm::real>::const_iterator pos2y;
53 std::vector<cvm::real>::const_iterator pos2z;
73 const std::vector<cvm::real> &pos1,
74 const std::vector<cvm::real> &pos2,
75 const size_t num_atoms_pos1,
76 const size_t num_atoms_pos2):
79 pos1y(
pos1x + num_atoms_pos1),
80 pos1z(pos1y + num_atoms_pos1),
82 pos2y(
pos2x + num_atoms_pos2),
83 pos2z(pos2y + num_atoms_pos2),
93 if (require_dl_dq & rotation_derivative_dldq::use_dl) {
112 if (require_dl_dq & rotation_derivative_dldq::use_dq) {
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
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];
340 template <
bool use_dl,
bool use_dq,
bool use_ds>
344 std::array<cvm::rvector, 4>* _noalias
const dq0_out,
345 std::array<std::array<cvm::rvector, 4>, 4>* _noalias
const ds_out)
const {
348 *ds_out = std::array<std::array<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];
364 *dl0_out =
tmp_Q0Q0[0][0] * ds[0][0] +
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];
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];
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];
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];
476 template <
bool use_dl,
bool use_dq,
bool use_ds>
478 size_t ia,
cvm::rvector* _noalias
const dl0_1_out =
nullptr,
479 std::array<cvm::rvector, 4>* _noalias
const dq0_1_out =
nullptr,
480 std::array<std::array<cvm::rvector, 4>, 4>* _noalias
const ds_1_out =
nullptr)
const {
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);
502 template <
bool use_dl,
bool use_dq,
bool use_ds>
504 size_t ia,
cvm::rvector* _noalias
const dl0_2_out =
nullptr,
505 std::array<cvm::rvector, 4>* _noalias
const dq0_2_out =
nullptr,
506 std::array<std::array<cvm::rvector, 4>, 4>* _noalias
const ds_2_out =
nullptr)
const {
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);
537 static_assert((i < 4) && (i >= 0),
"i must be in [0, 3] in project_force_to_C_from_dxdqi.");
539 result.xx = f_on_q * ( tmp_Q0Q0_L[i][0][0] + tmp_Q0Q0_L[i][1][1] - tmp_Q0Q0_L[i][2][2] - tmp_Q0Q0_L[i][3][3] );
540 result.xy = f_on_q * ( tmp_Q0Q0_L[i][0][3] + tmp_Q0Q0_L[i][1][2] + tmp_Q0Q0_L[i][2][1] + tmp_Q0Q0_L[i][3][0] );
541 result.xz = f_on_q * (-tmp_Q0Q0_L[i][0][2] + tmp_Q0Q0_L[i][1][3] - tmp_Q0Q0_L[i][2][0] + tmp_Q0Q0_L[i][3][1] );
542 result.yx = f_on_q * (-tmp_Q0Q0_L[i][0][3] + tmp_Q0Q0_L[i][1][2] + tmp_Q0Q0_L[i][2][1] - tmp_Q0Q0_L[i][3][0] );
543 result.yy = f_on_q * ( tmp_Q0Q0_L[i][0][0] - tmp_Q0Q0_L[i][1][1] + tmp_Q0Q0_L[i][2][2] - tmp_Q0Q0_L[i][3][3] );
544 result.yz = f_on_q * ( tmp_Q0Q0_L[i][0][1] + tmp_Q0Q0_L[i][1][0] + tmp_Q0Q0_L[i][2][3] + tmp_Q0Q0_L[i][3][2] );
545 result.zx = f_on_q * ( tmp_Q0Q0_L[i][0][2] + tmp_Q0Q0_L[i][1][3] + tmp_Q0Q0_L[i][2][0] + tmp_Q0Q0_L[i][3][1] );
546 result.zy = f_on_q * (-tmp_Q0Q0_L[i][0][1] - tmp_Q0Q0_L[i][1][0] + tmp_Q0Q0_L[i][2][3] + tmp_Q0Q0_L[i][3][2] );
547 result.zz = f_on_q * ( tmp_Q0Q0_L[i][0][0] - tmp_Q0Q0_L[i][1][1] - tmp_Q0Q0_L[i][2][2] + tmp_Q0Q0_L[i][3][3] );
568 template <
typename dim4_array_t>
571 result += project_force_to_C_from_dxdqi<0>(sum_dxdq[0]);
572 result += project_force_to_C_from_dxdqi<1>(sum_dxdq[1]);
573 result += project_force_to_C_from_dxdqi<2>(sum_dxdq[2]);
574 result += project_force_to_C_from_dxdqi<3>(sum_dxdq[3]);
593 dxdC.xx * a2x + dxdC.xy * a2y + dxdC.xz * a2z,
594 dxdC.yx * a2x + dxdC.yy * a2y + dxdC.yz * a2z,
595 dxdC.zx * a2x + dxdC.zy * a2y + dxdC.zz * a2z};
614 dxdC.xx * a1x + dxdC.yx * a1y + dxdC.zx * a1z,
615 dxdC.xy * a1x + dxdC.yy * a1y + dxdC.zy * a1z,
616 dxdC.xz * a1x + dxdC.yz * a1y + dxdC.zz * a1z};
2-dimensional array of real numbers with three components along each dimension (works with colvarmodu...
Definition: colvartypes.h:892
A rotation between two sets of coordinates (for the moment a wrapper for colvarmodule::quaternion)
Definition: colvartypes.h:1326
cvm::real S_eigval[4]
Eigenvalues of S.
Definition: colvartypes.h:1335
cvm::real S_eigvec[4][4]
Eigenvectors of S.
Definition: colvartypes.h:1338
vector of real numbers with three components
Definition: colvartypes.h:724
double real
Defining an abstract real number allows to switch precision.
Definition: colvarmodule.h:139
Helper class for calculating the derivative of rotation.
Definition: colvar_rotation_derivative.h:41
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:71
cvm::rmatrix project_force_to_C_from_dxdqi(cvm::real f_on_q) const
Project the force on (or the gradient of ) to the force on the correlation matrix ,...
Definition: colvar_rotation_derivative.h:536
std::vector< cvm::real >::const_iterator pos1x
Reference to the atom positions of group 1.
Definition: colvar_rotation_derivative.h:46
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
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:477
size_t m_num_atoms_pos1
Number of atoms in group1 (used in SOA)
Definition: colvar_rotation_derivative.h:55
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:341
cvm::real tmp_Q0Q0[4][4]
Temporary variable that will be updated if prepare_derivative called.
Definition: colvar_rotation_derivative.h:59
std::vector< cvm::real >::const_iterator pos2x
Reference to the atom positions of group 2.
Definition: colvar_rotation_derivative.h:51
cvm::rvector project_force_to_group2(size_t ia, const cvm::rmatrix &dxdC) const
Project the force on the correlation matrix to .
Definition: colvar_rotation_derivative.h:609
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:503
const cvm::rotation & m_rot
Reference to the rotation.
Definition: colvar_rotation_derivative.h:43
cvm::rmatrix project_force_to_C_from_dxdq(const dim4_array_t &sum_dxdq) const
Project the force on (or the gradient of ) to the force on the correlation matrix ,...
Definition: colvar_rotation_derivative.h:569
size_t m_num_atoms_pos2
Number of atoms in group1 (used in SOA)
Definition: colvar_rotation_derivative.h:57
cvm::rvector project_force_to_group1(size_t ia, const cvm::rmatrix &dxdC) const
Project the force on the correlation matrix to .
Definition: colvar_rotation_derivative.h:588