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