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