Line data    Source code

       1             : /*----------------------------------------------------------------------------*/
       2             : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3             : /*  Copyright 2000-2024 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     1573848 : 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     1573848 :   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     1573848 :   const int kgmin = ceil(-1e-8 - disr_radius * dh_inv[2][2]);
      31     1573848 :   if (bounds != NULL) {
      32      786924 :     bounds[ibound] = kgmin;
      33             :   }
      34             :   ibound++;
      35    15481452 :   for (int kg = kgmin; kg <= 0; kg++) {
      36    13907604 :     const int kd = (2 * kg - 1) / 2; // distance from center in grid points
      37    13907604 :     const double kr = kd * dh[2][2]; // distance from center in a.u.
      38    13907604 :     const double kremain = disr_radius * disr_radius - kr * kr;
      39    13907604 :     const int jgmin = ceil(-1e-8 - sqrt(fmax(0.0, kremain)) * dh_inv[1][1]);
      40    13907604 :     if (bounds != NULL) {
      41     6953802 :       bounds[ibound] = jgmin;
      42             :     }
      43    13907604 :     ibound++;
      44   144891512 :     for (int jg = jgmin; jg <= 0; jg++) {
      45   130983908 :       const int jd = (2 * jg - 1) / 2; // distance from center in grid points
      46   130983908 :       const double jr = jd * dh[1][1]; // distance from center in a.u.
      47   130983908 :       const double jremain = kremain - jr * jr;
      48   130983908 :       const int igmin = ceil(-1e-8 - sqrt(fmax(0.0, jremain)) * dh_inv[0][0]);
      49   130983908 :       if (bounds != NULL) {
      50    65491954 :         bounds[ibound] = igmin;
      51             :       }
      52   130983908 :       ibound++;
      53             :     }
      54             :   }
      55     1573848 :   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       62709 : 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       62709 :   if (entry->max_imr > 0) {
      67       49197 :     free(entry->offsets);
      68       49197 :     free(entry->storage);
      69             :   }
      70       62709 :   entry->max_imr = max_imr;
      71             : 
      72             :   // Compute required storage size.
      73       62709 :   entry->offsets = malloc(max_imr * sizeof(int));
      74       62709 :   int nbounds_total = 0;
      75      849633 :   for (int imr = 1; imr <= max_imr; imr++) {
      76      786924 :     const double radius = imr * drmin;
      77      786924 :     const int nbounds = single_sphere_bounds(radius, dh, dh_inv, NULL);
      78      786924 :     entry->offsets[imr - 1] = nbounds_total;
      79      786924 :     nbounds_total += nbounds;
      80             :   }
      81             : 
      82             :   // Allocate and fill storage.
      83       62709 :   entry->storage = malloc(nbounds_total * sizeof(int));
      84      849633 :   for (int imr = 1; imr <= max_imr; imr++) {
      85      786924 :     const double radius = imr * drmin;
      86      786924 :     const int offset = entry->offsets[imr - 1];
      87      786924 :     single_sphere_bounds(radius, dh, dh_inv, &entry->storage[offset]);
      88             :   }
      89       62709 : }
      90             : 
      91             : /*******************************************************************************
      92             :  * \brief Lookup the sphere bound from cache and compute them as needed.
      93             :  *        See grid_sphere_cache.h for details.
      94             :  * \author Ole Schuett
      95             :  ******************************************************************************/
      96   119703585 : void grid_sphere_cache_lookup(const double radius, const double dh[3][3],
      97             :                               const double dh_inv[3][3], int **sphere_bounds,
      98             :                               double *discr_radius) {
      99             : 
     100             :   // Prepare the cache.
     101   119703585 :   grid_sphere_cache *cache = grid_library_get_sphere_cache();
     102             : 
     103             :   // Find or create cache entry for given grid.
     104   119703585 :   const double dr0 = dh[0][0], dr1 = dh[1][1], dr2 = dh[2][2];
     105   119703585 :   grid_sphere_cache_entry *entry;
     106   119703585 :   bool found = false;
     107             : 
     108             :   // Fast path: check prev match.
     109   119703585 :   if (cache->prev_match < cache->size) {
     110   119699138 :     entry = &cache->entries[cache->prev_match];
     111   119699138 :     if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) {
     112   118156832 :       found = true;
     113             :     }
     114             :   }
     115             : 
     116             :   // Full search.
     117   119703585 :   if (!found) {
     118     4052030 :     for (int i = 0; i < cache->size; i++) {
     119     4038518 :       entry = &cache->entries[i];
     120     4038518 :       if (entry->dr[0] == dr0 && entry->dr[1] == dr1 && entry->dr[2] == dr2) {
     121     1533241 :         cache->prev_match = i;
     122     1533241 :         found = true;
     123     1533241 :         break;
     124             :       }
     125             :     }
     126             :   }
     127             : 
     128             :   // If no existing cache entry was found then create a new one.
     129   119703585 :   if (!found) {
     130       13512 :     cache->size++;
     131       13512 :     grid_sphere_cache_entry *old_entries = cache->entries;
     132       13512 :     const size_t entry_size = sizeof(grid_sphere_cache_entry);
     133       13512 :     cache->entries = malloc(cache->size * entry_size);
     134       13512 :     memcpy(cache->entries, old_entries, (cache->size - 1) * entry_size);
     135       13512 :     free(old_entries);
     136       13512 :     cache->prev_match = cache->size - 1;
     137       13512 :     entry = &cache->entries[cache->size - 1];
     138             :     // Initialize new cache entry
     139       13512 :     entry->max_imr = 0;
     140       13512 :     entry->dr[0] = dr0;
     141       13512 :     entry->dr[1] = dr1;
     142       13512 :     entry->dr[2] = dr2;
     143       13512 :     entry->drmin = fmin(dr0, fmin(dr1, dr2));
     144       13512 :     entry->drmin_inv = 1.0 / entry->drmin;
     145             :   }
     146             : 
     147             :   // Discretize the radius.
     148   119703585 :   const int imr = imax(1, (int)ceil(radius * entry->drmin_inv));
     149   119703585 :   *discr_radius = entry->drmin * imr;
     150             : 
     151             :   // Rebuild cache entry if requested radius is too large.
     152   119703585 :   if (entry->max_imr < imr) {
     153       62709 :     rebuild_cache_entry(imr, entry->drmin, dh, dh_inv, entry);
     154             :   }
     155   119703585 :   const int offset = entry->offsets[imr - 1];
     156   119703585 :   *sphere_bounds = &entry->storage[offset];
     157   119703585 : }
     158             : 
     159             : /*******************************************************************************
     160             :  * \brief Free the memory of the sphere cache.
     161             :  * \author Ole Schuett
     162             :  ******************************************************************************/
     163        8394 : void grid_sphere_cache_free(grid_sphere_cache *cache) {
     164       21906 :   for (int i = 0; i < cache->size; i++) {
     165       13512 :     if (cache->entries[i].max_imr > 0) {
     166       13512 :       free(cache->entries[i].offsets);
     167       13512 :       free(cache->entries[i].storage);
     168             :     }
     169             :   }
     170        8394 :   free(cache->entries);
     171        8394 :   cache->size = 0;
     172        8394 : }
     173             : 
     174             : // EOF