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 <assert.h>
       9             : #include <stdio.h>
      10             : #include <stdlib.h>
      11             : 
      12             : #include "../common/grid_common.h"
      13             : #include "grid_ref_prepare_pab.h"
      14             : 
      15             : /*******************************************************************************
      16             :  * \brief Cab matrix container to be passed through prepare_pab to cab_add.
      17             :  * \author Ole Schuett
      18             :  ******************************************************************************/
      19             : typedef struct {
      20             :   double *data;
      21             :   const int n1;
      22             : } cab_store;
      23             : 
      24             : /*******************************************************************************
      25             :  * \brief Adds given value to matrix element cab[idx(b)][idx(a)].
      26             :  * \author Ole Schuett
      27             :  ******************************************************************************/
      28      129032 : static inline void cab_add(cab_store *cab, const orbital a, const orbital b,
      29             :                            const double value) {
      30      129032 :   cab->data[idx(b) * cab->n1 + idx(a)] += value;
      31      129032 : }
      32             : 
      33             : #include "../common/grid_prepare_pab.h"
      34             : 
      35             : /*******************************************************************************
      36             :  * \brief Returns block size changes due to transformation grid_prepare_pab.
      37             :  * \author Ole Schuett
      38             :  ******************************************************************************/
      39        7862 : void grid_ref_prepare_get_ldiffs(const enum grid_func func, int *la_min_diff,
      40             :                                  int *la_max_diff, int *lb_min_diff,
      41             :                                  int *lb_max_diff) {
      42        7862 :   const prepare_ldiffs ldiffs = prepare_get_ldiffs(func);
      43        7862 :   *la_min_diff = ldiffs.la_min_diff;
      44        7862 :   *la_max_diff = ldiffs.la_max_diff;
      45        7862 :   *lb_min_diff = ldiffs.lb_min_diff;
      46        7862 :   *lb_max_diff = ldiffs.lb_max_diff;
      47        7862 : }
      48             : 
      49             : /*******************************************************************************
      50             :  * \brief Selects and transforms a sub-block of the given density matrix block.
      51             :  *        See grid_ref_prepare_pab.h for details.
      52             :  * \author Ole Schuett
      53             :  ******************************************************************************/
      54        7862 : void grid_ref_prepare_pab(const enum grid_func func, const int o1, const int o2,
      55             :                           const int la_max, const int la_min, const int lb_max,
      56             :                           const int lb_min, const double zeta,
      57             :                           const double zetb, const int n1, const int n2,
      58        7862 :                           const double pab[n2][n1], const int n1_prep,
      59             :                           const int n2_prep,
      60        7862 :                           double pab_prep[n2_prep][n1_prep]) {
      61             : 
      62        7862 :   cab_store cab = {.data = (double *)pab_prep, .n1 = n1_prep};
      63             : 
      64       23784 :   for (int lxa = 0; lxa <= la_max; lxa++) {
      65       48180 :     for (int lxb = 0; lxb <= lb_max; lxb++) {
      66       81266 :       for (int lya = 0; lya <= la_max - lxa; lya++) {
      67      123492 :         for (int lyb = 0; lyb <= lb_max - lxb; lyb++) {
      68       74484 :           const int lza_start = imax(la_min - lxa - lya, 0);
      69      172500 :           for (int lza = lza_start; lza <= la_max - lxa - lya; lza++) {
      70       98016 :             const int lzb_start = imax(lb_min - lxb - lyb, 0);
      71      227048 :             for (int lzb = lzb_start; lzb <= lb_max - lxb - lyb; lzb++) {
      72      129032 :               const orbital a = {{lxa, lya, lza}};
      73      129032 :               const orbital b = {{lxb, lyb, lzb}};
      74      129032 :               const double pab_val = pab[o2 + idx(b)][o1 + idx(a)];
      75      129032 :               prepare_pab(func, a, b, zeta, zetb, pab_val, &cab);
      76             :             }
      77             :           }
      78             :         }
      79             :       }
      80             :     }
      81             :   }
      82        7862 : }
      83             : 
      84             : // EOF