Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2025 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: BSD-3-Clause */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #include <assert.h>
9 : #include <limits.h>
10 : #include <math.h>
11 : #include <stdio.h>
12 : #include <stdlib.h>
13 : #include <string.h>
14 :
15 : #define GRID_DO_COLLOCATE 0
16 : #include "../common/grid_common.h"
17 : #include "grid_cpu_collint.h"
18 : #include "grid_cpu_integrate.h"
19 :
20 : /*******************************************************************************
21 : * \brief Cab matrix container to be passed through get_force/virial to cab_get.
22 : * \author Ole Schuett
23 : ******************************************************************************/
24 : typedef struct {
25 : const double *data;
26 : const int m1;
27 : } cab_store;
28 :
29 : /*******************************************************************************
30 : * \brief Returns matrix element cab[idx(b)][idx(a)].
31 : * \author Ole Schuett
32 : ******************************************************************************/
33 5362459222 : static inline double cab_get(const cab_store *cab, const orbital a,
34 : const orbital b) {
35 5362459222 : return cab->data[idx(b) * cab->m1 + idx(a)];
36 : }
37 :
38 : #include "../common/grid_process_vab.h"
39 :
40 : /*******************************************************************************
41 : * \brief Integrates a single task. See grid_cpu_integrate.h for details.
42 : * \author Ole Schuett
43 : ******************************************************************************/
44 90803791 : void grid_cpu_integrate_pgf_product(
45 : const bool orthorhombic, const bool compute_tau, const int border_mask,
46 : const int la_max, const int la_min, const int lb_max, const int lb_min,
47 : const double zeta, const double zetb, const double dh[3][3],
48 : const double dh_inv[3][3], const double ra[3], const double rab[3],
49 : const int npts_global[3], const int npts_local[3], const int shift_local[3],
50 : const int border_width[3], const double radius, const int o1, const int o2,
51 90803791 : const int n1, const int n2, const double *grid, double hab[n2][n1],
52 : const double pab[n2][n1], double forces[2][3], double virials[2][3][3],
53 : double hdab[n2][n1][3], double hadb[n2][n1][3],
54 90803791 : double a_hdab[n2][n1][3][3]) {
55 :
56 90803791 : const bool calculate_forces = (forces != NULL || hdab != NULL);
57 90803791 : const bool calculate_virial = (virials != NULL || a_hdab != NULL);
58 90803791 : assert(!calculate_virial || calculate_forces);
59 90803791 : const process_ldiffs ldiffs =
60 90803791 : process_get_ldiffs(calculate_forces, calculate_virial, compute_tau);
61 90803791 : int la_max_local = la_max + ldiffs.la_max_diff;
62 90803791 : int lb_max_local = lb_max + ldiffs.lb_max_diff;
63 90803791 : int la_min_local = imax(0, la_min + ldiffs.la_min_diff);
64 90803791 : int lb_min_local = imax(0, lb_min + ldiffs.lb_min_diff);
65 :
66 90803791 : const int m1 = ncoset(la_max_local);
67 90803791 : const int m2 = ncoset(lb_max_local);
68 90803791 : double cab[m2 * m1];
69 90803791 : memset(cab, 0, m2 * m1 * sizeof(double));
70 :
71 90803791 : const cab_store cab_obj = {.data = cab, .m1 = m1};
72 :
73 90803791 : const double rscale = 1.0; // TODO: remove rscale from cab_to_grid
74 90803791 : cab_to_grid(orthorhombic, border_mask, la_max_local, la_min_local,
75 : lb_max_local, lb_min_local, zeta, zetb, rscale, dh, dh_inv, ra,
76 : rab, npts_global, npts_local, shift_local, border_width, radius,
77 : cab, grid);
78 :
79 : // cab contains all the information needed to find the elements of hab
80 : // and optionally of derivatives of these elements
81 219450132 : for (int la = la_min; la <= la_max; la++) {
82 318438829 : for (int ax = 0; ax <= la; ax++) {
83 450748629 : for (int ay = 0; ay <= la - ax; ay++) {
84 260956141 : const int az = la - ax - ay;
85 260956141 : const orbital a = {{ax, ay, az}};
86 701936681 : for (int lb = lb_min; lb <= lb_max; lb++) {
87 1155605435 : for (int bx = 0; bx <= lb; bx++) {
88 1767598383 : for (int by = 0; by <= lb - bx; by++) {
89 1052973488 : const int bz = lb - bx - by;
90 1052973488 : const orbital b = {{bx, by, bz}};
91 :
92 : // Update hab block.
93 1052973488 : hab[o2 + idx(b)][o1 + idx(a)] +=
94 1052973488 : get_hab(a, b, zeta, zetb, &cab_obj, compute_tau);
95 :
96 : // Update forces.
97 1052973488 : if (forces != NULL) {
98 143708462 : const double pabval = pab[o2 + idx(b)][o1 + idx(a)];
99 574833848 : for (int i = 0; i < 3; i++) {
100 431125386 : forces[0][i] += pabval * get_force_a(a, b, i, zeta, zetb,
101 : &cab_obj, compute_tau);
102 431125386 : forces[1][i] += pabval * get_force_b(a, b, i, zeta, zetb, rab,
103 : &cab_obj, compute_tau);
104 : }
105 : }
106 :
107 : // Update virials.
108 1052973488 : if (virials != NULL) {
109 22168613 : const double pabval = pab[o2 + idx(b)][o1 + idx(a)];
110 88674452 : for (int i = 0; i < 3; i++) {
111 266023356 : for (int j = 0; j < 3; j++) {
112 199517517 : virials[0][i][j] +=
113 199517517 : pabval * get_virial_a(a, b, i, j, zeta, zetb, &cab_obj,
114 : compute_tau);
115 199517517 : virials[1][i][j] +=
116 199517517 : pabval * get_virial_b(a, b, i, j, zeta, zetb, rab,
117 : &cab_obj, compute_tau);
118 : }
119 : }
120 : }
121 :
122 : // Update hdab, hadb, and a_hdab (not used in batch mode).
123 1052973488 : if (hdab != NULL) {
124 17354 : assert(!compute_tau);
125 69416 : for (int i = 0; i < 3; i++) {
126 52062 : hdab[o2 + idx(b)][o1 + idx(a)][i] +=
127 52062 : get_force_a(a, b, i, zeta, zetb, &cab_obj, false);
128 : }
129 : }
130 1052973488 : if (hadb != NULL) {
131 7902 : assert(!compute_tau);
132 31608 : for (int i = 0; i < 3; i++) {
133 23706 : hadb[o2 + idx(b)][o1 + idx(a)][i] +=
134 23706 : get_force_b(a, b, i, zeta, zetb, rab, &cab_obj, false);
135 : }
136 : }
137 1052973488 : if (a_hdab != NULL) {
138 360 : assert(!compute_tau);
139 1440 : for (int i = 0; i < 3; i++) {
140 4320 : for (int j = 0; j < 3; j++) {
141 3240 : a_hdab[o2 + idx(b)][o1 + idx(a)][i][j] +=
142 3240 : get_virial_a(a, b, i, j, zeta, zetb, &cab_obj, false);
143 : }
144 : }
145 : }
146 : }
147 : }
148 : }
149 : }
150 : }
151 : }
152 90803791 : }
153 :
154 : // EOF
|
