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