LCOV - code coverage report
Current view: top level - src/eri_mme - eri_mme_gaussian.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:d0bd076) Lines: 68 70 97.1 %
Date: 2021-09-15 13:52:28 Functions: 4 4 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2021 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Methods related to properties of Hermite and Cartesian Gaussian functions.
      10             : !> \par History
      11             : !>       2015 09 created
      12             : !> \author Patrick Seewald
      13             : ! **************************************************************************************************
      14             : 
      15             : MODULE eri_mme_gaussian
      16             :    USE kinds,                           ONLY: dp
      17             :    USE mathconstants,                   ONLY: gamma1
      18             :    USE minimax_exp,                     ONLY: get_exp_minimax_coeff
      19             : #include "../base/base_uses.f90"
      20             : 
      21             :    IMPLICIT NONE
      22             : 
      23             :    PRIVATE
      24             : 
      25             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
      26             : 
      27             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_gaussian'
      28             : 
      29             :    INTEGER, PARAMETER, PUBLIC :: eri_mme_coulomb = 1, &
      30             :                                  eri_mme_yukawa = 2, &
      31             :                                  eri_mme_longrange = 3
      32             : 
      33             :    PUBLIC :: &
      34             :       create_gaussian_overlap_dist_to_hermite, &
      35             :       create_hermite_to_cartesian, &
      36             :       get_minimax_coeff_v_gspace, &
      37             :       hermite_gauss_norm
      38             : 
      39             : CONTAINS
      40             : 
      41             : ! **************************************************************************************************
      42             : !> \brief Create matrix to transform between cartesian and hermite gaussian
      43             : !>        basis functions.
      44             : !> \param zet    exponent
      45             : !> \param l_max ...
      46             : !> \param h_to_c transformation matrix with dimensions (0:l_max, 0:l_max)
      47             : !> \note  is idempotent, so transformation is the same
      48             : !>        in both directions.
      49             : ! **************************************************************************************************
      50     1530431 :    PURE SUBROUTINE create_hermite_to_cartesian(zet, l_max, h_to_c)
      51             :       REAL(KIND=dp), INTENT(IN)                          :: zet
      52             :       INTEGER, INTENT(IN)                                :: l_max
      53             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
      54             :          INTENT(OUT)                                     :: h_to_c
      55             : 
      56             :       INTEGER                                            :: k, l
      57             : 
      58     6121724 :       ALLOCATE (h_to_c(-1:l_max + 1, 0:l_max))
      59    37045095 :       h_to_c(:, :) = 0.0_dp
      60     1530431 :       h_to_c(0, 0) = 1.0_dp
      61     5023808 :       DO l = 0, l_max - 1
      62    16226901 :          DO k = 0, l + 1
      63    14696470 :             h_to_c(k, l + 1) = -(k + 1)*h_to_c(k + 1, l) + 2.0_dp*zet*h_to_c(k - 1, l)
      64             :          ENDDO
      65             :       ENDDO
      66             : 
      67     1530431 :    END SUBROUTINE create_hermite_to_cartesian
      68             : 
      69             : ! **************************************************************************************************
      70             : !> \brief Norm of 1d Hermite-Gauss functions
      71             : !> \param zet ...
      72             : !> \param l ...
      73             : !> \return ...
      74             : ! **************************************************************************************************
      75     5157124 :    PURE FUNCTION hermite_gauss_norm(zet, l) RESULT(norm)
      76             :       REAL(KIND=dp), INTENT(IN)                          :: zet
      77             :       INTEGER, DIMENSION(3), INTENT(IN)                  :: l
      78             :       REAL(KIND=dp)                                      :: norm
      79             : 
      80    20628496 :       norm = 1.0_dp/SQRT((2.0_dp*zet)**(SUM(l) - 1.5_dp)*(gamma1(l(1))*gamma1(l(2))*gamma1(l(3))))
      81             : 
      82     5157124 :    END FUNCTION hermite_gauss_norm
      83             : 
      84             : ! **************************************************************************************************
      85             : !> \brief Get minimax coefficient a_i and w_i for approximating
      86             : !>        1/G^2 by sum_i w_i exp(-a_i G^2)
      87             : !> \param n_minimax   Number of minimax terms
      88             : !> \param cutoff      Plane Wave cutoff
      89             : !> \param G_min       Minimum absolute value of G
      90             : !> \param minimax_aw  Minimax coefficients a_i, w_i
      91             : !> \param potential   potential to use. Accepts the following values:
      92             : !>                    1: coulomb potential V(r)=1/r
      93             : !>                    2: yukawa potential V(r)=e(-a*r)/r
      94             : !>                    3: long-range coulomb erf(a*r)/r
      95             : !> \param pot_par     potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r
      96             : !> \param err_minimax Maximum error MAX (|1/G^2-\sum_i w_i exp(-a_i G^2)|)
      97             : ! **************************************************************************************************
      98      112692 :    SUBROUTINE get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw, potential, pot_par, err_minimax)
      99             :       INTEGER, INTENT(IN)                                :: n_minimax
     100             :       REAL(KIND=dp), INTENT(IN)                          :: cutoff, G_min
     101             :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: minimax_aw
     102             :       INTEGER, INTENT(IN), OPTIONAL                      :: potential
     103             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pot_par
     104             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: err_minimax
     105             : 
     106             :       INTEGER                                            :: potential_prv
     107             :       REAL(KIND=dp)                                      :: dG, G_max, minimax_Rc
     108      112692 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: a, w
     109             : 
     110      112692 :       IF (PRESENT(potential)) THEN
     111      111236 :          potential_prv = potential
     112             :       ELSE
     113             :          potential_prv = eri_mme_coulomb
     114             :       ENDIF
     115             : 
     116      111236 :       IF (potential_prv > 3) THEN
     117           0 :          CPABORT("unknown potential")
     118             :       ENDIF
     119             : 
     120      112692 :       IF ((potential_prv >= 2) .AND. .NOT. PRESENT(pot_par)) THEN
     121           0 :          CPABORT("potential parameter pot_par required for yukawa or long-range Coulomb")
     122             :       ENDIF
     123             : 
     124      112692 :       dG = 1.0E-3 ! Resolution in G to determine error of minimax approximation
     125             : 
     126             :       ! Note: G_c = SQRT(2*cutoff) cutoff in 1 cartesian direction
     127             :       ! G_max = SQRT(3*G_c**2) maximum absolute value of G vector
     128             :       ! Minimax approx. needs to be valid in range [G_min, G_max]
     129             : 
     130             :       ! 1) compute minimax coefficients
     131             : 
     132      112692 :       G_max = SQRT(3.0_dp*2.0_dp*cutoff)
     133      112692 :       CPASSERT(G_max .GT. G_min)
     134      112692 :       IF (potential_prv == eri_mme_coulomb .OR. potential_prv == eri_mme_longrange) THEN
     135      112684 :          minimax_Rc = (G_max/G_min)**2
     136           8 :       ELSEIF (potential_prv == eri_mme_yukawa) THEN
     137           8 :          minimax_Rc = (G_max**2 + pot_par**2)/(G_min**2 + pot_par**2)
     138             :       ENDIF
     139             : 
     140      112692 :       CALL get_exp_minimax_coeff(n_minimax, minimax_Rc, minimax_aw, err_minimax)
     141             : 
     142      563460 :       ALLOCATE (a(n_minimax)); ALLOCATE (w(n_minimax))
     143     1601966 :       a(:) = minimax_aw(:n_minimax)
     144     1601966 :       w(:) = minimax_aw(n_minimax + 1:)
     145      107514 :       SELECT CASE (potential_prv)
     146             :          ! Scale minimax coefficients to incorporate different Fourier transforms
     147             :       CASE (eri_mme_coulomb)
     148             :          ! FT = 1/G**2
     149     1493308 :          a(:) = a/G_min**2
     150     1493308 :          w(:) = w/G_min**2
     151             :       CASE (eri_mme_yukawa)
     152             :          ! FT = 1/(G**2 + pot_par**2)
     153         128 :          w(:) = w*EXP((-a*pot_par**2)/(G_min**2 + pot_par**2))/(G_min**2 + pot_par**2)
     154         128 :          a(:) = a/(G_min**2 + pot_par**2)
     155             :       CASE (eri_mme_longrange)
     156             :          ! FT = exp(-(G/pot_par)**2)/G**2
     157             :          ! approximating 1/G**2 as for Coulomb:
     158      108530 :          a(:) = a/G_min**2
     159      108530 :          w(:) = w/G_min**2
     160             :          ! incorporate exponential factor:
     161      221222 :          a(:) = a + 1.0_dp/pot_par**2
     162             :       END SELECT
     163     6069788 :       minimax_aw = [a(:), w(:)]
     164             : 
     165      112692 :       IF (PRESENT(err_minimax)) THEN
     166      112692 :          IF (potential_prv == eri_mme_coulomb) THEN
     167      107514 :             err_minimax = err_minimax/G_min**2
     168        5178 :          ELSEIF (potential_prv == eri_mme_yukawa) THEN
     169           8 :             err_minimax = err_minimax/(G_min**2 + pot_par**2)
     170        5170 :          ELSEIF (potential_prv == eri_mme_longrange) THEN
     171        5170 :             err_minimax = err_minimax/G_min**2 ! approx. of Coulomb
     172        5170 :             err_minimax = err_minimax*EXP(-G_min**2/pot_par**2) ! exponential factor
     173             :          ENDIF
     174             :       ENDIF
     175             : 
     176      112692 :    END SUBROUTINE get_minimax_coeff_v_gspace
     177             : 
     178             : ! **************************************************************************************************
     179             : !> \brief Expand 1d product of cartesian (or hermite) gaussians into single hermite gaussians:
     180             : !>        Find E_t^{lm} s.t.
     181             : !>        F(l, a, r-R1) * F(m, b, r-R2) = sum_{t=0}^{l+m} E_t^{lm} H(t, p, r-R_P)
     182             : !>        with p = a + b, R_P = (a*R1 + b*R2)/p. The function F can be either Cartesian
     183             : !>        Gaussian or Hermite Gaussian.
     184             : !> \param l ...
     185             : !> \param m ...
     186             : !> \param a ...
     187             : !> \param b ...
     188             : !> \param R1 ...
     189             : !> \param R2 ...
     190             : !> \param H_or_C_product 1: cartesian product, 2: hermite product
     191             : !> \param E ...
     192             : ! **************************************************************************************************
     193     2232345 :    PURE SUBROUTINE create_gaussian_overlap_dist_to_hermite(l, m, a, b, R1, R2, H_or_C_product, E)
     194             :       INTEGER, INTENT(IN)                                :: l, m
     195             :       REAL(KIND=dp), INTENT(IN)                          :: a, b, R1, R2
     196             :       INTEGER, INTENT(IN)                                :: H_or_C_product
     197             :       REAL(KIND=dp), DIMENSION(-1:l+m+1, -1:l, -1:m), &
     198             :          INTENT(OUT)                                     :: E
     199             : 
     200             :       INTEGER                                            :: ll, mm, t
     201             :       REAL(KIND=dp)                                      :: c1, c2, c3
     202             : 
     203    81759456 :       E(:, :, :) = 0.0_dp
     204     2232345 :       E(0, 0, 0) = EXP(-a*b/(a + b)*(R1 - R2)**2) ! cost: exp_w flops
     205             : 
     206     2232345 :       c1 = 0.5_dp/(a + b)
     207     2232345 :       c2 = (b/(a + b))*(R2 - R1)
     208     2232345 :       c3 = (a/(a + b))*(R1 - R2)
     209             : 
     210     2232345 :       IF (H_or_C_product .EQ. 1) THEN ! Cartesian overlap dist
     211     2627658 :          DO mm = 0, m
     212     5013384 :             DO ll = 0, l
     213    10674180 :                DO t = 0, ll + mm + 1
     214     6673956 :                   IF (ll .LT. l) THEN
     215             :                      E(t, ll + 1, mm) = c1*E(t - 1, ll, mm) + & ! cost: 8 flops
     216             :                                         c2*E(t, ll, mm) + &
     217     1988118 :                                         (t + 1)*E(t + 1, ll, mm)
     218             :                   ENDIF
     219     9059682 :                   IF (mm .LT. m) THEN
     220             :                      E(t, ll, mm + 1) = c1*E(t - 1, ll, mm) + & ! cost: 8 flops
     221             :                                         c3*E(t, ll, mm) + &
     222     2342994 :                                         (t + 1)*E(t + 1, ll, mm)
     223             :                   ENDIF
     224             :                ENDDO
     225             :             ENDDO
     226             :          ENDDO
     227             :       ELSE ! Hermite overlap dist
     228     2960643 :          DO mm = 0, m
     229     5620305 :             DO ll = 0, l
     230    11931531 :                DO t = 0, ll + mm + 1
     231     7530411 :                   IF (ll .LT. l) THEN
     232             :                      E(t, ll + 1, mm) = a*(2*c1*E(t - 1, ll, mm) + & ! cost: 16 flops
     233             :                                            2*c2*E(t, ll, mm) + &
     234             :                                            2*(t + 1)*E(t + 1, ll, mm) - &
     235     2526669 :                                            2*ll*E(t, ll - 1, mm))
     236             :                   ENDIF
     237    10190073 :                   IF (mm .LT. m) THEN
     238             :                      E(t, ll, mm + 1) = b*(2*c1*E(t - 1, ll, mm) + & ! cost: 16 flops
     239             :                                            2*c3*E(t, ll, mm) + &
     240             :                                            2*(t + 1)*E(t + 1, ll, mm) - &
     241     2529750 :                                            2*mm*E(t, ll, mm - 1))
     242             : 
     243             :                   ENDIF
     244             :                ENDDO
     245             :             ENDDO
     246             :          ENDDO
     247             :       ENDIF
     248             : 
     249     2232345 :    END SUBROUTINE create_gaussian_overlap_dist_to_hermite
     250             : END MODULE eri_mme_gaussian

Generated by: LCOV version 1.15