(git:58e3e09)
grid_ref_integrate.c
Go to the documentation of this file.
1 /*----------------------------------------------------------------------------*/
2 /* CP2K: A general program to perform molecular dynamics simulations */
3 /* Copyright 2000-2024 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_ref_collint.h"
18 #include "grid_ref_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 static inline double cab_get(const cab_store *cab, const orbital a,
34  const orbital b) {
35  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_ref_integrate.h for details.
42  * \author Ole Schuett
43  ******************************************************************************/
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  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  double a_hdab[n2][n1][3][3]) {
55 
56  const bool calculate_forces = (forces != NULL || hdab != NULL);
57  const bool calculate_virial = (virials != NULL || a_hdab != NULL);
58  assert(!calculate_virial || calculate_forces);
59  const process_ldiffs ldiffs =
60  process_get_ldiffs(calculate_forces, calculate_virial, compute_tau);
61  int la_max_local = la_max + ldiffs.la_max_diff;
62  int lb_max_local = lb_max + ldiffs.lb_max_diff;
63  int la_min_local = imax(0, la_min + ldiffs.la_min_diff);
64  int lb_min_local = imax(0, lb_min + ldiffs.lb_min_diff);
65 
66  const int m1 = ncoset(la_max_local);
67  const int m2 = ncoset(lb_max_local);
68  double cab[m2 * m1];
69  memset(cab, 0, m2 * m1 * sizeof(double));
70 
71  const cab_store cab_obj = {.data = cab, .m1 = m1};
72 
73  const double rscale = 1.0; // TODO: remove rscale from cab_to_grid
74  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  for (int la = la_min; la <= la_max; la++) {
82  for (int ax = 0; ax <= la; ax++) {
83  for (int ay = 0; ay <= la - ax; ay++) {
84  const int az = la - ax - ay;
85  const orbital a = {{ax, ay, az}};
86  for (int lb = lb_min; lb <= lb_max; lb++) {
87  for (int bx = 0; bx <= lb; bx++) {
88  for (int by = 0; by <= lb - bx; by++) {
89  const int bz = lb - bx - by;
90  const orbital b = {{bx, by, bz}};
91 
92  // Update hab block.
93  hab[o2 + idx(b)][o1 + idx(a)] +=
94  get_hab(a, b, zeta, zetb, &cab_obj, compute_tau);
95 
96  // Update forces.
97  if (forces != NULL) {
98  const double pabval = pab[o2 + idx(b)][o1 + idx(a)];
99  for (int i = 0; i < 3; i++) {
100  forces[0][i] += pabval * get_force_a(a, b, i, zeta, zetb,
101  &cab_obj, compute_tau);
102  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  if (virials != NULL) {
109  const double pabval = pab[o2 + idx(b)][o1 + idx(a)];
110  for (int i = 0; i < 3; i++) {
111  for (int j = 0; j < 3; j++) {
112  virials[0][i][j] +=
113  pabval * get_virial_a(a, b, i, j, zeta, zetb, &cab_obj,
114  compute_tau);
115  virials[1][i][j] +=
116  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  if (hdab != NULL) {
124  assert(!compute_tau);
125  for (int i = 0; i < 3; i++) {
126  hdab[o2 + idx(b)][o1 + idx(a)][i] +=
127  get_force_a(a, b, i, zeta, zetb, &cab_obj, false);
128  }
129  }
130  if (hadb != NULL) {
131  assert(!compute_tau);
132  for (int i = 0; i < 3; i++) {
133  hadb[o2 + idx(b)][o1 + idx(a)][i] +=
134  get_force_b(a, b, i, zeta, zetb, rab, &cab_obj, false);
135  }
136  }
137  if (a_hdab != NULL) {
138  assert(!compute_tau);
139  for (int i = 0; i < 3; i++) {
140  for (int j = 0; j < 3; j++) {
141  a_hdab[o2 + idx(b)][o1 + idx(a)][i][j] +=
142  get_virial_a(a, b, i, j, zeta, zetb, &cab_obj, false);
143  }
144  }
145  }
146  }
147  }
148  }
149  }
150  }
151  }
152 }
153 
154 // EOF
static int imax(int x, int y)
Returns the larger of two given integer (missing from the C standard)
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition: grid_common.h:73
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition: grid_common.h:153
static void const int const int const int const int const int const double const int const int const int int GRID_CONST_WHEN_COLLOCATE double GRID_CONST_WHEN_INTEGRATE double * grid
static void const int const int i
static void const int const int const int const int const int const double const int const int const int npts_local[3]
static GRID_DEVICE double get_hab(const orbital a, const orbital b, const double zeta, const double zetb, const cab_store *cab, const bool compute_tau)
Returns element i,j of hab matrix.
static process_ldiffs process_get_ldiffs(bool calculate_forces, bool calculate_virial, bool compute_tau)
Returns difference in angular momentum range for given flags.
static GRID_DEVICE double get_force_b(const orbital a, const orbital b, const int i, const double zeta, const double zetb, const double rab[3], const cab_store *cab, const bool compute_tau)
Returns i'th component of force on atom b.
static GRID_DEVICE double get_virial_b(const orbital a, const orbital b, const int i, const int j, const double zeta, const double zetb, const double rab[3], const cab_store *cab, const bool compute_tau)
Returns element i,j of virial on atom b.
static GRID_DEVICE double get_virial_a(const orbital a, const orbital b, const int i, const int j, const double zeta, const double zetb, const cab_store *cab, const bool compute_tau)
Returns element i,j of virial on atom a.
static GRID_DEVICE double get_force_a(const orbital a, const orbital b, const int i, const double zeta, const double zetb, const cab_store *cab, const bool compute_tau)
Returns i'th component of force on atom a.
static void cab_to_grid(const bool orthorhombic, const int border_mask, const int la_max, const int la_min, const int lb_max, const int lb_min, const double zeta, const double zetb, const double rscale, const double dh[3][3], const double dh_inv[3][3], const double ra[3], const double rab[3], const int npts_global[3], const int npts_local[3], const int shift_local[3], const int border_width[3], const double radius, GRID_CONST_WHEN_COLLOCATE double *cab, GRID_CONST_WHEN_INTEGRATE double *grid)
Collocates coefficients C_ab onto the grid.
void grid_ref_integrate_pgf_product(const bool orthorhombic, const bool compute_tau, const int border_mask, const int la_max, const int la_min, const int lb_max, const int lb_min, const double zeta, const double zetb, const double dh[3][3], const double dh_inv[3][3], const double ra[3], const double rab[3], const int npts_global[3], const int npts_local[3], const int shift_local[3], const int border_width[3], const double radius, const int o1, const int o2, const int n1, const int n2, const double *grid, double hab[n2][n1], const double pab[n2][n1], double forces[2][3], double virials[2][3][3], double hdab[n2][n1][3], double hadb[n2][n1][3], double a_hdab[n2][n1][3][3])
Integrates a single task. See grid_ref_integrate.h for details.
static double cab_get(const cab_store *cab, const orbital a, const orbital b)
Returns matrix element cab[idx(b)][idx(a)].
real(dp), dimension(3) a
Definition: ai_eri_debug.F:31
real(dp), dimension(3) b
Definition: ai_eri_debug.F:31
Cab matrix container to be passed through get_force/virial to cab_get.
const int m1
const double * data
Orbital angular momentum.
Definition: grid_common.h:125
Differences in angular momentum.