Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of Coulomb integrals over Cartesian Gaussian-type functions
10 : !> (electron repulsion integrals, ERIs).
11 : !> \par Literature
12 : !> S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
13 : !> \par History
14 : !> none
15 : !> \par Parameters
16 : !> - ax,ay,az : Angular momentum index numbers of orbital a.
17 : !> - bx,by,bz : Angular momentum index numbers of orbital b.
18 : !> - coset : Cartesian orbital set pointer.
19 : !> - dab : Distance between the atomic centers a and b.
20 : !> - dac : Distance between the atomic centers a and c.
21 : !> - dbc : Distance between the atomic centers b and c.
22 : !> - l{a,b,c} : Angular momentum quantum number of shell a, b or c.
23 : !> - l{a,b,c}_max: Maximum angular momentum quantum number of shell a, b or c.
24 : !> - l{a,b,c}_min: Minimum angular momentum quantum number of shell a, b or c.
25 : !> - ncoset : Number of orbitals in a Cartesian orbital set.
26 : !> - npgf{a,b} : Degree of contraction of shell a or b.
27 : !> - rab : Distance vector between the atomic centers a and b.
28 : !> - rab2 : Square of the distance between the atomic centers a and b.
29 : !> - rac : Distance vector between the atomic centers a and c.
30 : !> - rac2 : Square of the distance between the atomic centers a and c.
31 : !> - rbc : Distance vector between the atomic centers b and c.
32 : !> - rbc2 : Square of the distance between the atomic centers b and c.
33 : !> - rpgf{a,b,c} : Radius of the primitive Gaussian-type function a, b or c.
34 : !> - zet{a,b,c} : Exponents of the Gaussian-type functions a, b or c.
35 : !> - zetp : Reciprocal of the sum of the exponents of orbital a and b.
36 : !> \author VW
37 : ! **************************************************************************************************
38 : MODULE ai_elec_field
39 : USE ai_cartesian_operators, ONLY: apply_cartesian_operator_pair,&
40 : cartesian_ladder_map
41 : USE ai_os_rr, ONLY: os_rr_coul
42 : USE kinds, ONLY: dp
43 : USE mathconstants, ONLY: pi
44 : USE orbital_pointers, ONLY: ncoset
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_elec_field'
49 : PRIVATE
50 :
51 : ! *** Public subroutines ***
52 :
53 : PUBLIC :: efg
54 :
55 : CONTAINS
56 :
57 : ! **************************************************************************************************
58 : !> \brief Calculation of the primitive electric field integrals over
59 : !> Cartesian Gaussian-type functions.
60 : !> \param la_max ...
61 : !> \param la_min ...
62 : !> \param npgfa ...
63 : !> \param rpgfa ...
64 : !> \param zeta ...
65 : !> \param lb_max ...
66 : !> \param lb_min ...
67 : !> \param npgfb ...
68 : !> \param rpgfb ...
69 : !> \param zetb ...
70 : !> \param rac ...
71 : !> \param rbc ...
72 : !> \param rab ...
73 : !> \param vab ...
74 : !> \param ldrr1 ...
75 : !> \param ldrr2 ...
76 : !> \param rr ...
77 : !> \date 02.03.2009
78 : !> \author VW
79 : !> \version 1.0
80 : ! **************************************************************************************************
81 9222 : SUBROUTINE efg(la_max, la_min, npgfa, rpgfa, zeta, &
82 18444 : lb_max, lb_min, npgfb, rpgfb, zetb, &
83 9222 : rac, rbc, rab, vab, ldrr1, ldrr2, rr)
84 : INTEGER, INTENT(IN) :: la_max, la_min, npgfa
85 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfa, zeta
86 : INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb
87 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rpgfb, zetb
88 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac, rbc, rab
89 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: vab
90 : INTEGER, INTENT(IN) :: ldrr1, ldrr2
91 : REAL(KIND=dp), DIMENSION(0:ldrr1-1, ldrr2, *), &
92 : INTENT(INOUT) :: rr
93 :
94 : INTEGER, DIMENSION(2, 6), PARAMETER :: &
95 : symmetric_components = RESHAPE([1, 1, 1, 2, 1, 3, 2, 2, 2, 3, 3, 3], [2, 6])
96 : INTEGER, DIMENSION(2, 9), PARAMETER :: &
97 : mixed_components = RESHAPE([1, 1, 2, 1, 3, 1, 1, 2, 2, 2, 3, 2, 1, 3, 2, 3, 3, 3], [2, 9])
98 :
99 : INTEGER :: axis, first_a, first_b, i, icomp, idir, &
100 : ipgf, j, jdir, jpgf, na, nb, ncap, &
101 : ncbp, ncol, nrow
102 9222 : INTEGER, ALLOCATABLE :: map_a(:, :, :), map_ap(:, :, :), &
103 : map_b(:, :, :), map_bp(:, :, :)
104 : REAL(KIND=dp) :: dab, f0, rab2, trace_third, xhi, za, zb, &
105 : zet
106 : REAL(KIND=dp), ALLOCATABLE :: hessian(:, :, :), tensor(:, :, :), &
107 9222 : work(:, :)
108 : REAL(KIND=dp), DIMENSION(3) :: rap, rbp, rcp
109 :
110 9222 : first_a = ncoset(la_min - 1) + 1
111 9222 : first_b = ncoset(lb_min - 1) + 1
112 9222 : nrow = ncoset(la_max) - first_a + 1
113 9222 : ncol = ncoset(lb_max) - first_b + 1
114 9222 : ncap = ncoset(la_max + 1)
115 9222 : ncbp = ncoset(lb_max + 1)
116 119886 : ALLOCATE (map_a(4, nrow, 3), map_ap(4, ncap, 3), map_b(4, ncol, 3), map_bp(4, ncbp, 3))
117 110664 : ALLOCATE (hessian(nrow, ncol, 6), tensor(nrow, ncol, 9), work(ncap, ncbp))
118 36888 : DO axis = 1, 3
119 27666 : CALL cartesian_ladder_map(la_max, la_min, axis, map_a(:, :, axis), source_lmin=0)
120 27666 : CALL cartesian_ladder_map(la_max + 1, 0, axis, map_ap(:, :, axis))
121 27666 : CALL cartesian_ladder_map(lb_max, lb_min, axis, map_b(:, :, axis), source_lmin=0)
122 36888 : CALL cartesian_ladder_map(lb_max + 1, 0, axis, map_bp(:, :, axis))
123 : END DO
124 :
125 : ! *** Calculate the distance of the centers a and c ***
126 :
127 9222 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
128 9222 : dab = SQRT(rab2)
129 :
130 : ! *** Loop over all pairs of primitive Gaussian-type functions ***
131 :
132 9222 : na = 0
133 :
134 21050 : DO ipgf = 1, npgfa
135 :
136 11828 : nb = 0
137 :
138 27026 : DO jpgf = 1, npgfb
139 :
140 : ! *** Screening ***
141 :
142 15198 : IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN
143 6888 : DO j = nb + 1, nb + ncoset(lb_max)
144 12801 : DO i = na + 1, na + ncoset(la_max)
145 5913 : vab(i, j, 1) = 0.0_dp
146 5913 : vab(i, j, 2) = 0.0_dp
147 5913 : vab(i, j, 3) = 0.0_dp
148 5913 : vab(i, j, 4) = 0.0_dp
149 5913 : vab(i, j, 5) = 0.0_dp
150 9999 : vab(i, j, 6) = 0.0_dp
151 : END DO
152 : END DO
153 2802 : nb = nb + ncoset(lb_max)
154 2802 : CYCLE
155 : END IF
156 :
157 : ! *** Calculate some prefactors ***
158 :
159 12396 : za = zeta(ipgf)
160 12396 : zb = zetb(jpgf)
161 12396 : zet = za + zb
162 12396 : xhi = za*zb/zet
163 49584 : rap = zb*rab/zet
164 49584 : rbp = -za*rab/zet
165 49584 : rcp = -(za*rac + zb*rbc)/zet
166 12396 : f0 = 2.0_dp*SQRT(zet/pi)*(pi/zet)**(1.5_dp)*EXP(-xhi*rab2)
167 :
168 : ! *** Calculate the recurrence relation
169 :
170 12396 : CALL os_rr_coul(rap, la_max + 2, rbp, lb_max + 2, rcp, zet, ldrr1, ldrr2, rr)
171 :
172 : ! Translational invariance gives D_C = -(D_A + D_B).
173 : ! The A/B ladders differentiate the basis functions with the potential held fixed.
174 : ! Form its Hessian from AA, BB and both ordered AB derivatives.
175 : CALL apply_cartesian_operator_pair(rr(0, 1:ncoset(la_max + 2), first_b:ncoset(lb_max)), &
176 : hessian, map_ap, map_a, symmetric_components, [.TRUE., .TRUE.], [za, za], &
177 37188 : [.FALSE., .FALSE.], work(:, 1:ncol))
178 : CALL apply_cartesian_operator_pair(rr(0, first_a:ncoset(la_max), 1:ncoset(lb_max + 2)), &
179 : tensor(:, :, 1:6), map_bp, map_b, symmetric_components, [.TRUE., .TRUE.], [zb, zb], &
180 37188 : [.TRUE., .TRUE.], work(1:nrow, :))
181 442044 : hessian(:, :, :) = hessian + tensor(:, :, 1:6)
182 : CALL apply_cartesian_operator_pair(rr(0, 1:ncap, 1:ncbp), tensor, map_a, map_b, &
183 37188 : mixed_components, [.TRUE., .TRUE.], [za, zb], [.FALSE., .TRUE.], work(1:nrow, :))
184 86772 : DO icomp = 1, 6
185 74376 : idir = symmetric_components(1, icomp)
186 74376 : jdir = symmetric_components(2, icomp)
187 : hessian(:, :, icomp) = hessian(:, :, icomp) + &
188 442044 : tensor(:, :, idir + 3*(jdir - 1)) + tensor(:, :, jdir + 3*(idir - 1))
189 : END DO
190 :
191 : ! Preserve the symmetric, traceless EFG ordering: xx, xy, xz, yy, yz, zz.
192 33442 : DO j = 1, ncol
193 71608 : DO i = 1, nrow
194 38166 : trace_third = (hessian(i, j, 1) + hessian(i, j, 4) + hessian(i, j, 6))/3.0_dp
195 38166 : hessian(i, j, 1) = hessian(i, j, 1) - trace_third
196 38166 : hessian(i, j, 4) = hessian(i, j, 4) - trace_third
197 59212 : hessian(i, j, 6) = hessian(i, j, 6) - trace_third
198 : END DO
199 : END DO
200 442044 : vab(na + first_a:na + ncoset(la_max), nb + first_b:nb + ncoset(lb_max), 1:6) = f0*hessian
201 :
202 24224 : nb = nb + ncoset(lb_max)
203 :
204 : END DO
205 :
206 21050 : na = na + ncoset(la_max)
207 :
208 : END DO
209 :
210 9222 : END SUBROUTINE efg
211 :
212 : END MODULE ai_elec_field
|
