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 "grid_sphere_cache.h"
9 : #include "grid_common.h"
10 : #include "grid_library.h"
11 : #include <assert.h>
12 : #include <math.h>
13 : #include <stdio.h>
14 : #include <stdlib.h>
15 : #include <string.h>
16 :
17 : /*******************************************************************************
18 : * \brief Compute the sphere bounds for a given single radius.
19 : * \author Ole Schuett
20 : ******************************************************************************/
21 1641758 : static int single_sphere_bounds(const double disr_radius, const double dh[3][3],
22 : const double dh_inv[3][3], int *bounds) {
23 :
24 1641758 : int ibound = 0;
25 :
26 : // The cube contains an even number of grid points in each direction and
27 : // collocation is always performed on a pair of two opposing grid points.
28 : // Hence, the points with index 0 and 1 are both assigned distance zero via
29 : // the formular distance=(2*index-1)/2.
30 1641758 : const int kgmin = ceil(-1e-8 - disr_radius * dh_inv[2][2]);
31 1641758 : if (bounds != NULL) {
32 820879 : bounds[ibound] = kgmin;
33 : }
34 : ibound++;
35 16210420 : for (int kg = kgmin; kg <= 0; kg++) {
36 14568662 : const int kd = (2 * kg - 1) / 2; // distance from center in grid points
37 14568662 : const double kr = kd * dh[2][2]; // distance from center in a.u.
38 14568662 : const double kremain = disr_radius * disr_radius - kr * kr;
39 14568662 : const int jgmin = ceil(-1e-8 - sqrt(fmax(0.0, kremain)) * dh_inv[1][1]);
40 14568662 : if (bounds != NULL) {
41 7284331 : bounds[ibound] = jgmin;
42 : }
43 14568662 : ibound++;
44 152528978 : for (int jg = jgmin; jg <= 0; jg++) {
45 137960316 : const int jd = (2 * jg - 1) / 2; // distance from center in grid points
46 137960316 : const double jr = jd * dh[1][1]; // distance from center in a.u.
47 137960316 : const double jremain = kremain - jr * jr;
48 137960316 : const int igmin = ceil(-1e-8 - sqrt(fmax(0.0, jremain)) * dh_inv[0][0]);
49 137960316 : if (bounds != NULL) {
50 68980158 : bounds[ibound] = igmin;
51 : }
52 137960316 : ibound++;
53 : }
54 : }
55 1641758 : return ibound; // Number of bounds - needed to allocate array.
56 : }
57 :
58 : /*******************************************************************************
59 : * \brief Rebuild a cache entry for a given cell and max radius.
60 : * \author Ole Schuett
61 : ******************************************************************************/
62 65162 : static void rebuild_cache_entry(const int max_imr, const double drmin,
63 : const double dh[3][3],
64 : const double dh_inv[3][3],
65 : grid_sphere_cache_entry *entry) {
66 65162 : if (entry->max_imr > 0) {
67 51085 : free(entry->offsets);
68 51085 : free(entry->storage);
69 : }
70 65162 : entry->max_imr = max_imr;
71 :
72 : // Compute required storage size.
73 65162 : entry->offsets = malloc(max_imr * sizeof(int));
74 65162 : assert(entry->offsets != NULL || max_imr == 0);
75 : int nbounds_total = 0;
76 886041 : for (int imr = 1; imr <= max_imr; imr++) {
77 820879 : const double radius = imr * drmin;
78 820879 : const int nbounds = single_sphere_bounds(radius, dh, dh_inv, NULL);
79 820879 : entry->offsets[imr - 1] = nbounds_total;
80 820879 : nbounds_total += nbounds;
81 : }
82 :
83 : // Allocate and fill storage.
84 65162 : entry->storage = malloc(nbounds_total * sizeof(int));
85 65162 : assert(entry->storage != NULL || nbounds_total == 0);
86 886041 : for (int imr = 1; imr <= max_imr; imr++) {
87 820879 : const double radius = imr * drmin;
88 820879 : const int offset = entry->offsets[imr - 1];
89 820879 : single_sphere_bounds(radius, dh, dh_inv, &entry->storage[offset]);
90 : }
91 65162 : }
92 :
93 : /*******************************************************************************
94 : * \brief Lookup the sphere bound from cache and compute them as needed.
95 : * See grid_sphere_cache.h for details.
96 : * \author Ole Schuett
97 : ******************************************************************************/
98 128252625 : void grid_sphere_cache_lookup(const double radius, const double dh[3][3],
99 : const double dh_inv[3][3], int **sphere_bounds,
100 : double *discr_radius) {
101 :
102 : // Prepare the cache.
103 128252625 : grid_sphere_cache *cache = grid_library_get_sphere_cache();
104 :
105 : // Find or create cache entry for given grid.
106 128252625 : const double dr0 = dh[0][0], dr1 = dh[1][1], dr2 = dh[2][2];
107 128252625 : grid_sphere_cache_entry *entry = NULL;
108 128252625 : bool found = false;
109 :
110 : // Fast path: check prev match.
111 128252625 : if (cache->prev_match < cache->size) {
112 128248001 : entry = &cache->entries[cache->prev_match];
113 128248001 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) {
114 126586333 : found = true;
115 : }
116 : }
117 :
118 : // Full search.
119 128252625 : if (!found) {
120 4335993 : for (int i = 0; i < cache->size; i++) {
121 4321916 : entry = &cache->entries[i];
122 4321916 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) {
123 1652215 : cache->prev_match = i;
124 1652215 : found = true;
125 1652215 : break;
126 : }
127 : }
128 : }
129 :
130 : // If no existing cache entry was found then create a new one.
131 128252625 : if (!found) {
132 14077 : cache->size++;
133 14077 : grid_sphere_cache_entry *old_entries = cache->entries;
134 14077 : const size_t entry_size = sizeof(grid_sphere_cache_entry);
135 14077 : cache->entries = malloc(cache->size * entry_size);
136 14077 : assert(cache->entries != NULL);
137 14077 : if (old_entries != NULL) {
138 9453 : memcpy(cache->entries, old_entries, (cache->size - 1) * entry_size);
139 9453 : free(old_entries);
140 : }
141 14077 : cache->prev_match = cache->size - 1;
142 14077 : entry = &cache->entries[cache->size - 1];
143 : // Initialize new cache entry
144 14077 : entry->max_imr = 0;
145 14077 : entry->dr[0] = dr0;
146 14077 : entry->dr[1] = dr1;
147 14077 : entry->dr[2] = dr2;
148 14077 : entry->drmin = fmin(dr0, fmin(dr1, dr2));
149 14077 : entry->drmin_inv = 1.0 / entry->drmin;
150 : }
151 :
152 : // Discretize the radius.
153 128252625 : const int imr = imax(1, (int)ceil(radius * entry->drmin_inv));
154 128252625 : *discr_radius = entry->drmin * imr;
155 :
156 : // Rebuild cache entry if requested radius is too large.
157 128252625 : if (entry->max_imr < imr) {
158 65162 : rebuild_cache_entry(imr, entry->drmin, dh, dh_inv, entry);
159 : }
160 128252625 : const int offset = entry->offsets[imr - 1];
161 128252625 : *sphere_bounds = &entry->storage[offset];
162 128252625 : }
163 :
164 : /*******************************************************************************
165 : * \brief Free the memory of the sphere cache.
166 : * \author Ole Schuett
167 : ******************************************************************************/
168 9286 : void grid_sphere_cache_free(grid_sphere_cache *cache) {
169 23363 : for (int i = 0; i < cache->size; i++) {
170 14077 : if (cache->entries[i].max_imr > 0) {
171 14077 : free(cache->entries[i].offsets);
172 14077 : free(cache->entries[i].storage);
173 : }
174 : }
175 9286 : free(cache->entries);
176 9286 : cache->size = 0;
177 9286 : }
178 :
179 : // EOF
|
