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 <float.h>
10 : #include <limits.h>
11 : #include <math.h>
12 : #include <stdio.h>
13 : #include <stdlib.h>
14 : #include <string.h>
15 :
16 : #define GRID_DO_COLLOCATE 1
17 : #include "../common/grid_common.h"
18 : #include "grid_cpu_collint.h"
19 : #include "grid_cpu_collocate.h"
20 : #include "grid_cpu_integrate.h"
21 : #include "grid_cpu_prepare_pab.h"
22 :
23 : /*******************************************************************************
24 : * \brief Writes the given arguments into a .task file.
25 : * See grid_replay.h for details.
26 : * \author Ole Schuett
27 : ******************************************************************************/
28 0 : static void write_task_file(
29 : const bool orthorhombic, const int border_mask, const enum grid_func func,
30 : const int la_max, const int la_min, const int lb_max, const int lb_min,
31 : const double zeta, const double zetb, const double rscale,
32 : const double dh[3][3], const double dh_inv[3][3], const double ra[3],
33 : const double rab[3], const int npts_global[3], const int npts_local[3],
34 : const int shift_local[3], const int border_width[3], const double radius,
35 : const int o1, const int o2, const int n1, const int n2,
36 0 : const double pab[n2][n1], const double *grid) {
37 :
38 0 : static int counter = 1;
39 0 : int my_number;
40 :
41 0 : #pragma omp critical
42 0 : my_number = counter++;
43 :
44 0 : char filename[100];
45 0 : snprintf(filename, sizeof(filename), "grid_collocate_%05i.task", my_number);
46 :
47 0 : const int D = DBL_DECIMAL_DIG;
48 0 : FILE *fp = fopen(filename, "w+");
49 0 : fprintf(fp, "#Grid task v10\n");
50 0 : fprintf(fp, "orthorhombic %i\n", orthorhombic);
51 0 : fprintf(fp, "border_mask %i\n", border_mask);
52 0 : fprintf(fp, "func %i\n", func);
53 0 : fprintf(fp, "la_max %i\n", la_max);
54 0 : fprintf(fp, "la_min %i\n", la_min);
55 0 : fprintf(fp, "lb_max %i\n", lb_max);
56 0 : fprintf(fp, "lb_min %i\n", lb_min);
57 0 : fprintf(fp, "zeta %.*e\n", D, zeta);
58 0 : fprintf(fp, "zetb %.*e\n", D, zetb);
59 0 : fprintf(fp, "rscale %.*e\n", D, rscale);
60 0 : for (int i = 0; i < 3; i++)
61 0 : fprintf(fp, "dh %i %.*e %.*e %.*e\n", i, D, dh[i][0], D, dh[i][1], D,
62 0 : dh[i][2]);
63 0 : for (int i = 0; i < 3; i++)
64 0 : fprintf(fp, "dh_inv %i %.*e %.*e %.*e\n", i, D, dh_inv[i][0], D,
65 0 : dh_inv[i][1], D, dh_inv[i][2]);
66 0 : fprintf(fp, "ra %.*e %.*e %.*e\n", D, ra[0], D, ra[1], D, ra[2]);
67 0 : fprintf(fp, "rab %.*e %.*e %.*e\n", D, rab[0], D, rab[1], D, rab[2]);
68 0 : fprintf(fp, "npts_global %i %i %i\n", npts_global[0], npts_global[1],
69 : npts_global[2]);
70 0 : fprintf(fp, "npts_local %i %i %i\n", npts_local[0], npts_local[1],
71 : npts_local[2]);
72 0 : fprintf(fp, "shift_local %i %i %i\n", shift_local[0], shift_local[1],
73 : shift_local[2]);
74 0 : fprintf(fp, "border_width %i %i %i\n", border_width[0], border_width[1],
75 : border_width[2]);
76 0 : fprintf(fp, "radius %.*e\n", D, radius);
77 0 : fprintf(fp, "o1 %i\n", o1);
78 0 : fprintf(fp, "o2 %i\n", o2);
79 0 : fprintf(fp, "n1 %i\n", n1);
80 0 : fprintf(fp, "n2 %i\n", n2);
81 :
82 0 : for (int i = 0; i < n2; i++) {
83 0 : for (int j = 0; j < n1; j++) {
84 0 : fprintf(fp, "pab %i %i %.*e\n", i, j, D, pab[i][j]);
85 : }
86 : }
87 :
88 0 : const int npts_local_total = npts_local[0] * npts_local[1] * npts_local[2];
89 :
90 0 : int ngrid_nonzero = 0;
91 0 : for (int i = 0; i < npts_local_total; i++) {
92 0 : if (grid[i] != 0.0) {
93 0 : ngrid_nonzero++;
94 : }
95 : }
96 0 : fprintf(fp, "ngrid_nonzero %i\n", ngrid_nonzero);
97 :
98 0 : for (int k = 0; k < npts_local[2]; k++) {
99 0 : for (int j = 0; j < npts_local[1]; j++) {
100 0 : for (int i = 0; i < npts_local[0]; i++) {
101 0 : const double val =
102 0 : grid[k * npts_local[1] * npts_local[0] + j * npts_local[0] + i];
103 0 : if (val != 0.0) {
104 0 : fprintf(fp, "grid %i %i %i %.*e\n", i, j, k, D, val);
105 : }
106 : }
107 : }
108 : }
109 :
110 0 : double hab[n2][n1];
111 0 : double forces[2][3] = {0};
112 0 : double virials[2][3][3] = {0};
113 0 : memset(hab, 0, n2 * n1 * sizeof(double));
114 :
115 0 : const bool compute_tau = (func == GRID_FUNC_DADB);
116 :
117 0 : grid_cpu_integrate_pgf_product(orthorhombic, compute_tau, border_mask, la_max,
118 : la_min, lb_max, lb_min, zeta, zetb, dh, dh_inv,
119 : ra, rab, npts_global, npts_local, shift_local,
120 : border_width, radius, o1, o2, n1, n2, grid,
121 : hab, pab, forces, virials, NULL, NULL, NULL);
122 :
123 0 : for (int i = o2; i < ncoset(lb_max) + o2; i++) {
124 0 : for (int j = o1; j < ncoset(la_max) + o1; j++) {
125 0 : fprintf(fp, "hab %i %i %.*e\n", i, j, D, hab[i][j]);
126 : }
127 : }
128 0 : fprintf(fp, "force_a %.*e %.*e %.*e\n", D, forces[0][0], D, forces[0][1], D,
129 : forces[0][2]);
130 0 : fprintf(fp, "force_b %.*e %.*e %.*e\n", D, forces[1][0], D, forces[1][1], D,
131 : forces[1][2]);
132 0 : for (int i = 0; i < 3; i++)
133 0 : fprintf(fp, "virial %i %.*e %.*e %.*e\n", i, D,
134 0 : virials[0][i][0] + virials[1][i][0], D,
135 0 : virials[0][i][1] + virials[1][i][1], D,
136 0 : virials[0][i][2] + virials[1][i][2]);
137 :
138 0 : fprintf(fp, "#THE_END\n");
139 0 : fclose(fp);
140 0 : printf("Wrote %s\n", filename);
141 0 : }
142 :
143 : /*******************************************************************************
144 : * \brief Collocates a single product of primitiv Gaussians.
145 : * See grid_cpu_collocate.h for details.
146 : * \author Ole Schuett
147 : ******************************************************************************/
148 90379116 : static void collocate_internal(
149 : const bool orthorhombic, const int border_mask, const enum grid_func func,
150 : const int la_max, const int la_min, const int lb_max, const int lb_min,
151 : const double zeta, const double zetb, const double rscale,
152 : const double dh[3][3], const double dh_inv[3][3], const double ra[3],
153 : const double rab[3], const int npts_global[3], const int npts_local[3],
154 : const int shift_local[3], const int border_width[3], const double radius,
155 : const int o1, const int o2, const int n1, const int n2,
156 90379116 : const double pab[n2][n1], double *grid) {
157 :
158 90379116 : int la_min_diff, la_max_diff, lb_min_diff, lb_max_diff;
159 90379116 : grid_cpu_prepare_get_ldiffs(func, &la_min_diff, &la_max_diff, &lb_min_diff,
160 : &lb_max_diff);
161 :
162 90379116 : const int la_min_cab = imax(la_min + la_min_diff, 0);
163 90379116 : const int lb_min_cab = imax(lb_min + lb_min_diff, 0);
164 90379116 : const int la_max_cab = la_max + la_max_diff;
165 90379116 : const int lb_max_cab = lb_max + lb_max_diff;
166 90379116 : const int n1_cab = ncoset(la_max_cab);
167 90379116 : const int n2_cab = ncoset(lb_max_cab);
168 :
169 90379116 : const size_t cab_size = n2_cab * n1_cab;
170 90379116 : double cab[cab_size];
171 90379116 : memset(cab, 0, cab_size * sizeof(double));
172 :
173 90379116 : grid_cpu_prepare_pab(func, o1, o2, la_max, la_min, lb_max, lb_min, zeta, zetb,
174 90379116 : n1, n2, pab, n1_cab, n2_cab, (double(*)[n1_cab])cab);
175 90379116 : cab_to_grid(orthorhombic, border_mask, la_max_cab, la_min_cab, lb_max_cab,
176 : lb_min_cab, zeta, zetb, rscale, dh, dh_inv, ra, rab, npts_global,
177 : npts_local, shift_local, border_width, radius, cab, grid);
178 90379116 : }
179 :
180 : /*******************************************************************************
181 : * \brief Public entry point. A thin wrapper with the only purpose of calling
182 : * write_task_file when DUMP_TASKS = true.
183 : * \author Ole Schuett
184 : ******************************************************************************/
185 90379116 : void grid_cpu_collocate_pgf_product(
186 : const bool orthorhombic, const int border_mask, const enum grid_func func,
187 : const int la_max, const int la_min, const int lb_max, const int lb_min,
188 : const double zeta, const double zetb, const double rscale,
189 : const double dh[3][3], const double dh_inv[3][3], const double ra[3],
190 : const double rab[3], const int npts_global[3], const int npts_local[3],
191 : const int shift_local[3], const int border_width[3], const double radius,
192 : const int o1, const int o2, const int n1, const int n2,
193 90379116 : const double pab[n2][n1], double *grid) {
194 :
195 : // Set this to true to write each task to a file.
196 90379116 : const bool DUMP_TASKS = false;
197 :
198 90379116 : double *grid_before = NULL;
199 90379116 : const size_t npts_local_total = npts_local[0] * npts_local[1] * npts_local[2];
200 :
201 90379116 : if (DUMP_TASKS) {
202 : const size_t sizeof_grid = sizeof(double) * npts_local_total;
203 : grid_before = malloc(sizeof_grid);
204 : assert(grid_before != NULL);
205 : memcpy(grid_before, grid, sizeof_grid);
206 : memset(grid, 0, sizeof_grid);
207 : }
208 :
209 90379116 : collocate_internal(orthorhombic, border_mask, func, la_max, la_min, lb_max,
210 : lb_min, zeta, zetb, rscale, dh, dh_inv, ra, rab,
211 : npts_global, npts_local, shift_local, border_width, radius,
212 : o1, o2, n1, n2, pab, grid);
213 :
214 90379116 : if (DUMP_TASKS) {
215 : write_task_file(orthorhombic, border_mask, func, la_max, la_min, lb_max,
216 : lb_min, zeta, zetb, rscale, dh, dh_inv, ra, rab,
217 : npts_global, npts_local, shift_local, border_width, radius,
218 : o1, o2, n1, n2, pab, grid);
219 :
220 : for (size_t i = 0; i < npts_local_total; i++) {
221 : grid[i] += grid_before[i];
222 : }
223 : free(grid_before);
224 : }
225 90379116 : }
226 :
227 : // EOF
|
