Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2026 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 1690640 : 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 1690640 : 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 1690640 : const int kgmin = ceil(-1e-8 - disr_radius * dh_inv[2][2]);
31 1690640 : if (bounds != NULL) {
32 845320 : bounds[ibound] = kgmin;
33 : }
34 : ibound++;
35 16709630 : for (int kg = kgmin; kg <= 0; kg++) {
36 15018990 : const int kd = (2 * kg - 1) / 2; // distance from center in grid points
37 15018990 : const double kr = kd * dh[2][2]; // distance from center in a.u.
38 15018990 : const double kremain = disr_radius * disr_radius - kr * kr;
39 15018990 : const int jgmin = ceil(-1e-8 - sqrt(fmax(0.0, kremain)) * dh_inv[1][1]);
40 15018990 : if (bounds != NULL) {
41 7509495 : bounds[ibound] = jgmin;
42 : }
43 15018990 : ibound++;
44 157013154 : for (int jg = jgmin; jg <= 0; jg++) {
45 141994164 : const int jd = (2 * jg - 1) / 2; // distance from center in grid points
46 141994164 : const double jr = jd * dh[1][1]; // distance from center in a.u.
47 141994164 : const double jremain = kremain - jr * jr;
48 141994164 : const int igmin = ceil(-1e-8 - sqrt(fmax(0.0, jremain)) * dh_inv[0][0]);
49 141994164 : if (bounds != NULL) {
50 70997082 : bounds[ibound] = igmin;
51 : }
52 141994164 : ibound++;
53 : }
54 : }
55 1690640 : 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 66953 : 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 66953 : if (entry->max_imr > 0) {
67 52347 : free(entry->offsets);
68 52347 : free(entry->storage);
69 : }
70 66953 : entry->max_imr = max_imr;
71 :
72 : // Compute required storage size.
73 66953 : entry->offsets = malloc(max_imr * sizeof(int));
74 66953 : assert(entry->offsets != NULL || max_imr == 0);
75 : int nbounds_total = 0;
76 912273 : for (int imr = 1; imr <= max_imr; imr++) {
77 845320 : const double radius = imr * drmin;
78 845320 : const int nbounds = single_sphere_bounds(radius, dh, dh_inv, NULL);
79 845320 : entry->offsets[imr - 1] = nbounds_total;
80 845320 : nbounds_total += nbounds;
81 : }
82 :
83 : // Allocate and fill storage.
84 66953 : entry->storage = malloc(nbounds_total * sizeof(int));
85 66953 : assert(entry->storage != NULL || nbounds_total == 0);
86 912273 : for (int imr = 1; imr <= max_imr; imr++) {
87 845320 : const double radius = imr * drmin;
88 845320 : const int offset = entry->offsets[imr - 1];
89 845320 : single_sphere_bounds(radius, dh, dh_inv, &entry->storage[offset]);
90 : }
91 66953 : }
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 130927310 : 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 130927310 : grid_sphere_cache *cache = grid_library_get_sphere_cache();
104 :
105 : // Find or create cache entry for given grid.
106 130927310 : const double dr0 = dh[0][0], dr1 = dh[1][1], dr2 = dh[2][2];
107 130927310 : grid_sphere_cache_entry *entry = NULL;
108 130927310 : bool found = false;
109 :
110 : // Fast path: check prev match.
111 130927310 : if (cache->prev_match < cache->size) {
112 130922492 : entry = &cache->entries[cache->prev_match];
113 130922492 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) {
114 129201177 : found = true;
115 : }
116 : }
117 :
118 : // Full search.
119 130927310 : if (!found) {
120 4454053 : for (int i = 0; i < cache->size; i++) {
121 4439447 : entry = &cache->entries[i];
122 4439447 : if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) {
123 1711527 : cache->prev_match = i;
124 1711527 : found = true;
125 1711527 : break;
126 : }
127 : }
128 : }
129 :
130 : // If no existing cache entry was found then create a new one.
131 130927310 : if (!found) {
132 14606 : cache->size++;
133 14606 : grid_sphere_cache_entry *old_entries = cache->entries;
134 14606 : const size_t entry_size = sizeof(grid_sphere_cache_entry);
135 14606 : cache->entries = malloc(cache->size * entry_size);
136 14606 : assert(cache->entries != NULL);
137 14606 : if (old_entries != NULL) {
138 9788 : memcpy(cache->entries, old_entries, (cache->size - 1) * entry_size);
139 9788 : free(old_entries);
140 : }
141 14606 : cache->prev_match = cache->size - 1;
142 14606 : entry = &cache->entries[cache->size - 1];
143 : // Initialize new cache entry
144 14606 : entry->max_imr = 0;
145 14606 : entry->dr[0] = dr0;
146 14606 : entry->dr[1] = dr1;
147 14606 : entry->dr[2] = dr2;
148 14606 : entry->drmin = fmin(dr0, fmin(dr1, dr2));
149 14606 : entry->drmin_inv = 1.0 / entry->drmin;
150 : }
151 :
152 : // Discretize the radius.
153 130927310 : const int imr = imax(1, (int)ceil(radius * entry->drmin_inv));
154 130927310 : *discr_radius = entry->drmin * imr;
155 :
156 : // Rebuild cache entry if requested radius is too large.
157 130927310 : if (entry->max_imr < imr) {
158 66953 : rebuild_cache_entry(imr, entry->drmin, dh, dh_inv, entry);
159 : }
160 130927310 : const int offset = entry->offsets[imr - 1];
161 130927310 : *sphere_bounds = &entry->storage[offset];
162 130927310 : }
163 :
164 : /*******************************************************************************
165 : * \brief Free the memory of the sphere cache.
166 : * \author Ole Schuett
167 : ******************************************************************************/
168 9522 : void grid_sphere_cache_free(grid_sphere_cache *cache) {
169 24128 : for (int i = 0; i < cache->size; i++) {
170 14606 : if (cache->entries[i].max_imr > 0) {
171 14606 : free(cache->entries[i].offsets);
172 14606 : free(cache->entries[i].storage);
173 : }
174 : }
175 9522 : free(cache->entries);
176 9522 : cache->size = 0;
177 9522 : }
178 :
179 : // EOF
|
