=========================================================================================== * Used to calculate the partial power spectrum derivatives. */
| 2355 | * Used to calculate the partial power spectrum derivatives. |
| 2356 | */ |
| 2357 | void getPDev( |
| 2358 | py::detail::unchecked_mutable_reference<double, 4> &derivatives_mu, |
| 2359 | py::detail::unchecked_reference<double, 2> &positions_u, |
| 2360 | py::detail::unchecked_reference<int, 1> &indices_u, |
| 2361 | CellList &cell_list, |
| 2362 | py::detail::unchecked_reference<double, 5> &CdevX_u, |
| 2363 | py::detail::unchecked_reference<double, 5> &CdevY_u, |
| 2364 | py::detail::unchecked_reference<double, 5> &CdevZ_u, |
| 2365 | py::detail::unchecked_reference<double, 4> &Cnnd_u, |
| 2366 | int Ns, |
| 2367 | int Ts, |
| 2368 | int nCenters, |
| 2369 | int lMax, |
| 2370 | bool crossover |
| 2371 | ) { |
| 2372 | |
| 2373 | // Loop over all given atomic indices for which the derivatives should be |
| 2374 | // calculated for. |
| 2375 | for (int i_idx = 0; i_idx < indices_u.size(); ++i_idx) { |
| 2376 | int i_atom = indices_u(i_idx); |
| 2377 | |
| 2378 | // Get all neighbouring centers for the current atom |
| 2379 | double ix = positions_u(i_atom, 0); |
| 2380 | double iy = positions_u(i_atom, 1); |
| 2381 | double iz = positions_u(i_atom, 2); |
| 2382 | CellListResult result = cell_list.getNeighboursForPosition(ix, iy, iz); |
| 2383 | vector<int> indices = result.indices; |
| 2384 | |
| 2385 | // Loop through all neighbouring centers |
| 2386 | for (size_t j_idx = 0; j_idx < indices.size(); ++j_idx) { |
| 2387 | int i_center = indices[j_idx]; |
| 2388 | int shiftAll = 0; |
| 2389 | for(int j = 0; j < Ts; j++) { |
| 2390 | int jdLimit = crossover ? Ts : j+1; |
| 2391 | for(int jd = j; jd < jdLimit; jd++) { |
| 2392 | for(int m=0; m <= lMax; m++) { |
| 2393 | double prel = m > 1 ? PI*sqrt(8.0/(2.0*m+1.0))*PI3 : PI*sqrt(8.0/(2.0*m+1.0)); |
| 2394 | if (j == jd) { |
| 2395 | for(int k = 0; k < Ns; k++){ |
| 2396 | for(int kd = k; kd < Ns; kd++){ |
| 2397 | for(int buffShift = m*m; buffShift < (m +1)*(m +1); buffShift++){ |
| 2398 | if( abs(Cnnd_u(i_center,j,k,buffShift)) > 1e-8 || abs(Cnnd_u(i_center,j,kd,buffShift)) > 1e-8 ){ |
| 2399 | derivatives_mu(i_center, i_idx, 0, shiftAll) += prel*(Cnnd_u(i_center,j,k,buffShift)*CdevX_u(i_atom, i_center, jd, kd, buffShift) |
| 2400 | +Cnnd_u(i_center,j,kd,buffShift)*CdevX_u(i_atom, i_center, jd, k, buffShift)); |
| 2401 | derivatives_mu(i_center, i_idx, 1, shiftAll) += prel*(Cnnd_u(i_center,j,k,buffShift)*CdevY_u(i_atom, i_center, jd, kd, buffShift) |
| 2402 | +Cnnd_u(i_center,j,kd,buffShift)*CdevY_u(i_atom, i_center, jd, k, buffShift)); |
| 2403 | derivatives_mu(i_center, i_idx, 2, shiftAll) += prel*(Cnnd_u(i_center,j,k,buffShift)*CdevZ_u(i_atom, i_center, jd, kd, buffShift) |
| 2404 | +Cnnd_u(i_center,j,kd,buffShift)*CdevZ_u(i_atom, i_center, jd, k, buffShift)); |
| 2405 | |
| 2406 | }} |
| 2407 | shiftAll++; |
| 2408 | } |
| 2409 | } |
| 2410 | } else { |
| 2411 | for(int k = 0; k < Ns; k++){ |
| 2412 | for(int kd = 0; kd < Ns; kd++){ |
| 2413 | for(int buffShift = m*m; buffShift < (m +1)*(m +1); buffShift++){ |
| 2414 | if( abs(Cnnd_u(i_center,j,k,buffShift)) > 1e-8 || abs(Cnnd_u(i_center,jd,kd,buffShift)) > 1e-8) { |
no test coverage detected