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