1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation
8!>        of contracted, spherical Gaussian integrals using the solid harmonics. Initialization of
9!>        the contraction matrices
10!> \par History
11!>      created [08.2016]
12!> \author Dorothea Golze
13! **************************************************************************************************
14MODULE generic_shg_integrals_init
15   USE basis_set_types,                 ONLY: gto_basis_set_type
16   USE kinds,                           ONLY: dp
17   USE mathconstants,                   ONLY: fac,&
18                                              ifac,&
19                                              pi
20   USE memory_utilities,                ONLY: reallocate
21   USE orbital_pointers,                ONLY: indso,&
22                                              nsoset
23   USE spherical_harmonics,             ONLY: clebsch_gordon,&
24                                              clebsch_gordon_deallocate,&
25                                              clebsch_gordon_init
26#include "../base/base_uses.f90"
27
28   IMPLICIT NONE
29
30   PRIVATE
31
32   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals_init'
33
34   PUBLIC :: contraction_matrix_shg, contraction_matrix_shg_mix, contraction_matrix_shg_rx2m, &
35             get_clebsch_gordon_coefficients
36
37! **************************************************************************************************
38
39CONTAINS
40
41! **************************************************************************************************
42!> \brief contraction matrix for SHG integrals
43!> \param basis ...
44!> \param scon_shg contraction matrix
45! **************************************************************************************************
46   SUBROUTINE contraction_matrix_shg(basis, scon_shg)
47
48      TYPE(gto_basis_set_type), POINTER                  :: basis
49      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: scon_shg
50
51      CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg', &
52         routineP = moduleN//':'//routineN
53
54      INTEGER                                            :: ipgf, iset, ishell, l, maxpgf, maxshell, &
55                                                            nset
56      INTEGER, DIMENSION(:), POINTER                     :: npgf, nshell
57      REAL(KIND=dp)                                      :: aif, gcc
58      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: norm
59      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
60
61      nset = basis%nset
62      npgf => basis%npgf
63      nshell => basis%nshell
64      zet => basis%zet
65
66      maxpgf = SIZE(basis%gcc, 1)
67      maxshell = SIZE(basis%gcc, 2)
68      ALLOCATE (norm(basis%nset, maxshell))
69      ALLOCATE (scon_shg(maxpgf, maxshell, nset))
70      scon_shg = 0.0_dp
71
72      CALL basis_norm_shg(basis, norm)
73
74      DO iset = 1, nset
75         DO ishell = 1, nshell(iset)
76            l = basis%l(ishell, iset)
77            DO ipgf = 1, npgf(iset)
78               aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l)
79               gcc = basis%gcc(ipgf, ishell, iset)
80               scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif
81            END DO
82         END DO
83      END DO
84
85      DEALLOCATE (norm)
86
87   END SUBROUTINE contraction_matrix_shg
88
89!***************************************************************************************************
90!> \brief normalization solid harmonic Gaussians (SHG)
91!> \param basis ...
92!> \param norm ...
93! **************************************************************************************************
94   SUBROUTINE basis_norm_shg(basis, norm)
95
96      TYPE(gto_basis_set_type), POINTER                  :: basis
97      REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: norm
98
99      CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_shg', routineP = moduleN//':'//routineN
100
101      INTEGER                                            :: ipgf, iset, ishell, jpgf, l
102      REAL(KIND=dp)                                      :: aai, aaj, cci, ccj, expa, ppl
103
104      norm = 0._dp
105
106      DO iset = 1, basis%nset
107         DO ishell = 1, basis%nshell(iset)
108            l = basis%l(ishell, iset)
109            expa = 0.5_dp*REAL(2*l + 3, dp)
110            ppl = fac(2*l + 2)*pi**(1.5_dp)/fac(l + 1)
111            ppl = ppl/(2._dp**REAL(2*l + 1, dp))
112            ppl = ppl/REAL(2*l + 1, dp)
113            DO ipgf = 1, basis%npgf(iset)
114               cci = basis%gcc(ipgf, ishell, iset)
115               aai = basis%zet(ipgf, iset)
116               DO jpgf = 1, basis%npgf(iset)
117                  ccj = basis%gcc(jpgf, ishell, iset)
118                  aaj = basis%zet(jpgf, iset)
119                  norm(iset, ishell) = norm(iset, ishell) + cci*ccj*ppl/(aai + aaj)**expa
120               END DO
121            END DO
122            norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell))
123         END DO
124      END DO
125
126   END SUBROUTINE basis_norm_shg
127
128! **************************************************************************************************
129!> \brief mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis
130!>        at the same atom
131!> \param orb_basis orbital basis
132!> \param ri_basis   ...
133!> \param orb_index index for orbital basis
134!> \param ri_index index for ri basis
135!> \param scon_mix mixed contraction matrix
136! **************************************************************************************************
137   SUBROUTINE contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix)
138
139      TYPE(gto_basis_set_type), POINTER                  :: orb_basis, ri_basis
140      INTEGER, DIMENSION(:, :, :), POINTER               :: orb_index, ri_index
141      REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: scon_mix
142
143      CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg_mix', &
144         routineP = moduleN//':'//routineN
145
146      INTEGER :: forb, fri, iil, il, ipgf, iset, ishell, jpgf, jset, jshell, l, l1, l2, lmax_orb, &
147         lmax_ri, maxpgf_orb, maxpgf_ri, maxshell_orb, maxshell_ri, nf_orb, nf_ri, nl, nl_max, &
148         nset_orb, nset_ri
149      INTEGER, DIMENSION(:), POINTER                     :: npgf_orb, npgf_ri, nshell_orb, nshell_ri
150      REAL(KIND=dp)                                      :: cjf, const, const1, const2, gcc_orb, &
151                                                            gcc_ri, prefac, scon1, scon2, zet
152      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: shg_fac
153      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: norm_orb, norm_ri, zet_orb, zet_ri
154
155      nset_orb = orb_basis%nset
156      npgf_orb => orb_basis%npgf
157      nshell_orb => orb_basis%nshell
158      zet_orb => orb_basis%zet
159      nset_ri = ri_basis%nset
160      npgf_ri => ri_basis%npgf
161      nshell_ri => ri_basis%nshell
162      zet_ri => ri_basis%zet
163
164      maxpgf_orb = SIZE(orb_basis%gcc, 1)
165      maxshell_orb = SIZE(orb_basis%gcc, 2)
166      ALLOCATE (norm_orb(nset_orb, maxshell_orb))
167      maxpgf_ri = SIZE(ri_basis%gcc, 1)
168      maxshell_ri = SIZE(ri_basis%gcc, 2)
169      ALLOCATE (norm_ri(nset_ri, maxshell_ri))
170
171      CALL basis_norm_shg(orb_basis, norm_orb)
172      CALL basis_norm_shg(ri_basis, norm_ri)
173
174      ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb))
175      ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri))
176
177      !** index orbital basis set
178      nf_orb = 0
179      DO iset = 1, nset_orb
180         DO ishell = 1, nshell_orb(iset)
181            DO ipgf = 1, npgf_orb(iset)
182               nf_orb = nf_orb + 1
183               orb_index(ipgf, ishell, iset) = nf_orb
184            END DO
185         END DO
186      END DO
187
188      !** index ri basis set
189      nf_ri = 0
190      DO iset = 1, nset_ri
191         DO ishell = 1, nshell_ri(iset)
192            DO ipgf = 1, npgf_ri(iset)
193               nf_ri = nf_ri + 1
194               ri_index(ipgf, ishell, iset) = nf_ri
195            END DO
196         END DO
197      END DO
198
199      lmax_orb = MAXVAL(orb_basis%lmax)
200      lmax_ri = MAXVAL(ri_basis%lmax)
201      nl_max = INT((lmax_orb + lmax_ri)/2) + 1
202      ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max))
203      scon_mix = 0.0_dp
204
205      ALLOCATE (shg_fac(0:nl_max - 1))
206      shg_fac(0) = 1.0_dp
207
208      DO iset = 1, nset_orb
209         DO ishell = 1, nshell_orb(iset)
210            l1 = orb_basis%l(ishell, iset)
211            const1 = SQRT(1.0_dp/REAL(2*l1 + 1, dp))
212            DO jset = 1, nset_ri
213               DO jshell = 1, nshell_ri(jset)
214                  l2 = ri_basis%l(jshell, jset)
215                  const2 = SQRT(1.0_dp/REAL(2*l2 + 1, dp))
216                  nl = INT((l1 + l2)/2)
217                  IF (l1 == 0 .OR. l2 == 0) nl = 0
218                  DO il = 0, nl
219                     l = l1 + l2 - 2*il
220                     const = const1*const2*2.0_dp*SQRT(pi*REAL(2*l + 1, dp))
221                     DO iil = 1, il
222                        shg_fac(iil) = fac(l + iil - 1)*ifac(l)*REAL(l, dp) &
223                                       *fac(il)/fac(il - iil)/fac(iil)
224                     ENDDO
225                     DO ipgf = 1, npgf_orb(iset)
226                        forb = orb_index(ipgf, ishell, iset)
227                        gcc_orb = orb_basis%gcc(ipgf, ishell, iset)
228                        scon1 = norm_orb(iset, ishell)*gcc_orb
229                        DO jpgf = 1, npgf_ri(jset)
230                           fri = ri_index(jpgf, jshell, jset)
231                           gcc_ri = ri_basis%gcc(jpgf, jshell, jset)
232                           scon2 = norm_ri(jset, jshell)*gcc_ri
233                           zet = zet_orb(ipgf, iset) + zet_ri(jpgf, jset)
234                           cjf = 1.0_dp/((2._dp*zet)**l)
235                           prefac = const*cjf*scon1*scon2
236                           DO iil = 0, il
237                              scon_mix(iil + 1, fri, forb, il + 1) = prefac*shg_fac(iil)/zet**iil
238                           ENDDO
239                        ENDDO
240                     ENDDO
241                  ENDDO
242               ENDDO
243            ENDDO
244         ENDDO
245      ENDDO
246
247      DEALLOCATE (norm_orb, norm_ri, shg_fac)
248
249   END SUBROUTINE contraction_matrix_shg_mix
250
251! **************************************************************************************************
252!> \brief ...
253!> \param basis ...
254!> \param m ...
255!> \param scon_shg ...
256!> \param scon_rx2m ...
257! **************************************************************************************************
258   SUBROUTINE contraction_matrix_shg_rx2m(basis, m, scon_shg, scon_rx2m)
259
260      TYPE(gto_basis_set_type), POINTER                  :: basis
261      INTEGER, INTENT(IN)                                :: m
262      REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN)      :: scon_shg
263      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: scon_rx2m
264
265      CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg_rx2m', &
266         routineP = moduleN//':'//routineN
267
268      INTEGER                                            :: ipgf, iset, ishell, j, l, maxpgf, &
269                                                            maxshell, nset
270      INTEGER, DIMENSION(:), POINTER                     :: npgf, nshell
271      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: shg_fac
272      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
273
274      npgf => basis%npgf
275      nshell => basis%nshell
276      zet => basis%zet
277      nset = basis%nset
278
279      maxpgf = SIZE(basis%gcc, 1)
280      maxshell = SIZE(basis%gcc, 2)
281      ALLOCATE (scon_rx2m(maxpgf, m + 1, maxshell, nset))
282      scon_rx2m = 0.0_dp
283      ALLOCATE (shg_fac(0:m))
284      shg_fac(0) = 1.0_dp
285
286      DO iset = 1, nset
287         DO ishell = 1, nshell(iset)
288            l = basis%l(ishell, iset)
289            DO j = 1, m
290               shg_fac(j) = fac(l + j - 1)*ifac(l)*REAL(l, dp) &
291                            *fac(m)/fac(m - j)/fac(j)
292            ENDDO
293            DO ipgf = 1, npgf(iset)
294               DO j = 0, m
295                  scon_rx2m(ipgf, j + 1, ishell, iset) = scon_shg(ipgf, ishell, iset) &
296                                                         *shg_fac(j)/zet(ipgf, iset)**j
297               ENDDO
298            ENDDO
299         ENDDO
300      ENDDO
301
302      DEALLOCATE (shg_fac)
303
304   END SUBROUTINE contraction_matrix_shg_rx2m
305
306! **************************************************************************************************
307!> \brief calculate the Clebsch-Gordon (CG) coefficients for expansion of the
308!>        product of two spherical harmonic Gaussians
309!> \param my_cg matrix storing CG coefficients
310!> \param cg_none0_list list of none-zero CG coefficients
311!> \param ncg_none0 number of none-zero CG coefficients
312!> \param maxl1 maximal l quantum number of 1st spherical function
313!> \param maxl2 maximal l quantum number of 2nd spherical function
314! **************************************************************************************************
315   SUBROUTINE get_clebsch_gordon_coefficients(my_cg, cg_none0_list, ncg_none0, &
316                                              maxl1, maxl2)
317
318      REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: my_cg
319      INTEGER, DIMENSION(:, :, :), POINTER               :: cg_none0_list
320      INTEGER, DIMENSION(:, :), POINTER                  :: ncg_none0
321      INTEGER, INTENT(IN)                                :: maxl1, maxl2
322
323      CHARACTER(len=*), PARAMETER :: routineN = 'get_clebsch_gordon_coefficients', &
324         routineP = moduleN//':'//routineN
325
326      INTEGER                                            :: il, ilist, iso, iso1, iso2, l1, l1l2, &
327                                                            l2, lc1, lc2, lp, m1, m2, maxl, mm, &
328                                                            mp, nlist, nlist_max, nsfunc, nsfunc1, &
329                                                            nsfunc2
330      REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: rga
331
332      nlist_max = 6
333      nsfunc1 = nsoset(maxl1)
334      nsfunc2 = nsoset(maxl2)
335      maxl = maxl1 + maxl2
336      nsfunc = nsoset(maxl)
337
338      CALL clebsch_gordon_init(maxl)
339
340      ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc))
341      my_cg = 0.0_dp
342      ALLOCATE (ncg_none0(nsfunc1, nsfunc2))
343      ncg_none0 = 0
344      ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max))
345      cg_none0_list = 0
346
347      ALLOCATE (rga(maxl, 2))
348      rga = 0.0_dp
349      DO lc1 = 0, maxl1
350         DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
351            l1 = indso(1, iso1)
352            m1 = indso(2, iso1)
353            DO lc2 = 0, maxl2
354               DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
355                  l2 = indso(1, iso2)
356                  m2 = indso(2, iso2)
357                  CALL clebsch_gordon(l1, m1, l2, m2, rga)
358                  l1l2 = l1 + l2
359                  mp = m1 + m2
360                  mm = m1 - m2
361                  IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
362                     mp = -ABS(mp)
363                     mm = -ABS(mm)
364                  ELSE
365                     mp = ABS(mp)
366                     mm = ABS(mm)
367                  END IF
368                  DO lp = MOD(l1 + l2, 2), l1l2, 2
369                     il = lp/2 + 1
370                     IF (ABS(mp) <= lp) THEN
371                        IF (mp >= 0) THEN
372                           iso = nsoset(lp - 1) + lp + 1 + mp
373                        ELSE
374                           iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
375                        END IF
376                        my_cg(iso1, iso2, iso) = rga(il, 1)
377                     ENDIF
378                     IF (mp /= mm .AND. ABS(mm) <= lp) THEN
379                        IF (mm >= 0) THEN
380                           iso = nsoset(lp - 1) + lp + 1 + mm
381                        ELSE
382                           iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
383                        END IF
384                        my_cg(iso1, iso2, iso) = rga(il, 2)
385                     ENDIF
386                  END DO
387                  nlist = 0
388                  DO ilist = 1, nsfunc
389                     IF (ABS(my_cg(iso1, iso2, ilist)) > 1.E-8_dp) THEN
390                        nlist = nlist + 1
391                        IF (nlist > nlist_max) THEN
392                           CALL reallocate(cg_none0_list, 1, nsfunc1, 1, nsfunc2, 1, nlist)
393                           nlist_max = nlist
394                        ENDIF
395                        cg_none0_list(iso1, iso2, nlist) = ilist
396                     ENDIF
397                  ENDDO
398                  ncg_none0(iso1, iso2) = nlist
399               ENDDO ! iso2
400            ENDDO ! lc2
401         ENDDO ! iso1
402      ENDDO ! lc1
403
404      DEALLOCATE (rga)
405      CALL clebsch_gordon_deallocate()
406
407   END SUBROUTINE get_clebsch_gordon_coefficients
408
409END MODULE generic_shg_integrals_init
410