Line data    Source code

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