Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 2483557 : 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 9934228 : ALLOCATE (h_to_c(-1:l_max + 1, 0:l_max))
59 58912129 : h_to_c(:, :) = 0.0_dp
60 2483557 : h_to_c(0, 0) = 1.0_dp
61 8074188 : DO l = 0, l_max - 1
62 25730729 : DO k = 0, l + 1
63 23247172 : 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 : END DO
65 : END DO
66 :
67 2483557 : 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 5198884 : 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 20795536 : 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 5198884 : 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 172380 : 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 172380 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: a, w
109 :
110 172380 : IF (PRESENT(potential)) THEN
111 170684 : potential_prv = potential
112 : ELSE
113 : potential_prv = eri_mme_coulomb
114 : END IF
115 :
116 170684 : IF (potential_prv > 3) THEN
117 0 : CPABORT("unknown potential")
118 : END IF
119 :
120 172380 : 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 : END IF
123 :
124 172380 : 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 172380 : G_max = SQRT(3.0_dp*2.0_dp*cutoff)
133 172380 : CPASSERT(G_max > G_min)
134 172380 : IF (potential_prv == eri_mme_coulomb .OR. potential_prv == eri_mme_longrange) THEN
135 172372 : 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 : END IF
139 :
140 172380 : CALL get_exp_minimax_coeff(n_minimax, minimax_Rc, minimax_aw, err_minimax)
141 :
142 689520 : ALLOCATE (a(n_minimax)); ALLOCATE (w(n_minimax))
143 2855414 : a(:) = minimax_aw(:n_minimax)
144 2855414 : w(:) = minimax_aw(n_minimax + 1:)
145 147282 : SELECT CASE (potential_prv)
146 : ! Scale minimax coefficients to incorporate different Fourier transforms
147 : CASE (eri_mme_coulomb)
148 : ! FT = 1/G**2
149 2328436 : a(:) = a/G_min**2
150 2328436 : 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 526850 : a(:) = a/G_min**2
159 526850 : w(:) = w/G_min**2
160 : ! incorporate exponential factor:
161 699230 : a(:) = a + 1.0_dp/pot_par**2
162 : END SELECT
163 10904516 : minimax_aw = [a(:), w(:)]
164 :
165 172380 : IF (PRESENT(err_minimax)) THEN
166 172380 : IF (potential_prv == eri_mme_coulomb) THEN
167 147282 : err_minimax = err_minimax/G_min**2
168 25098 : ELSEIF (potential_prv == eri_mme_yukawa) THEN
169 8 : err_minimax = err_minimax/(G_min**2 + pot_par**2)
170 25090 : ELSEIF (potential_prv == eri_mme_longrange) THEN
171 25090 : err_minimax = err_minimax/G_min**2 ! approx. of Coulomb
172 25090 : err_minimax = err_minimax*EXP(-G_min**2/pot_par**2) ! exponential factor
173 : END IF
174 : END IF
175 :
176 172380 : 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 == 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 < 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 : END IF
219 9059682 : IF (mm < 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 : END IF
224 : END DO
225 : END DO
226 : END DO
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 < 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 : END IF
237 10190073 : IF (mm < 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 : END IF
244 : END DO
245 : END DO
246 : END DO
247 : END IF
248 :
249 2232345 : END SUBROUTINE create_gaussian_overlap_dist_to_hermite
250 : END MODULE eri_mme_gaussian
|
