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 Common interface for two- and three-center Coulomb integrals.
10 : ! **************************************************************************************************
11 : MODULE coulomb_integral_interface
12 : USE ai_coulomb, ONLY: coulomb2,&
13 : coulomb3
14 : USE integral_library_types, ONLY: active_integral_library,&
15 : coulomb_operator_type,&
16 : library_libint,&
17 : library_native
18 : USE kinds, ONLY: dp
19 : USE libint_2c_3c, ONLY: eri_2center,&
20 : eri_3center
21 : USE libint_wrapper, ONLY: cp_libint_cleanup_2eri,&
22 : cp_libint_cleanup_3eri,&
23 : cp_libint_init_2eri,&
24 : cp_libint_init_3eri,&
25 : cp_libint_set_contrdepth,&
26 : cp_libint_t
27 : USE orbital_pointers, ONLY: ncoset
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 :
32 : PRIVATE
33 :
34 : PUBLIC :: coulomb_integral_cleanup, &
35 : coulomb_integral_init, &
36 : coulomb_operator_type, &
37 : compute_coulomb_2c, &
38 : compute_coulomb_3c
39 :
40 : TYPE, PUBLIC :: coulomb_integral_context_type
41 : PRIVATE
42 : TYPE(cp_libint_t) :: lib_2c
43 : TYPE(cp_libint_t) :: lib_3c
44 : LOGICAL :: initialized_2c = .FALSE.
45 : LOGICAL :: initialized_3c = .FALSE.
46 : END TYPE coulomb_integral_context_type
47 :
48 : CONTAINS
49 :
50 : ! **************************************************************************************************
51 : !> \brief Initialize the engines needed by the requested Coulomb integral operations.
52 : !> \param context operation context containing the selected-library engines
53 : !> \param max_am_2c maximum angular momentum for the two-center operation
54 : !> \param max_am_3c maximum angular momentum for the three-center operation
55 : ! **************************************************************************************************
56 260 : SUBROUTINE coulomb_integral_init(context, max_am_2c, max_am_3c)
57 : TYPE(coulomb_integral_context_type), INTENT(INOUT) :: context
58 : INTEGER, INTENT(IN), OPTIONAL :: max_am_2c, max_am_3c
59 :
60 260 : IF (PRESENT(max_am_2c)) THEN
61 192 : SELECT CASE (active_integral_library%coulomb2_library)
62 : CASE (library_libint)
63 6 : CALL cp_libint_init_2eri(context%lib_2c, max_am_2c)
64 6 : CALL cp_libint_set_contrdepth(context%lib_2c, 1)
65 : CASE (library_native)
66 : ! Nothing to do.
67 : CASE DEFAULT
68 186 : CPABORT("Requested two-center Coulomb library is not implemented")
69 : END SELECT
70 186 : context%initialized_2c = .TRUE.
71 : END IF
72 :
73 260 : IF (PRESENT(max_am_3c)) THEN
74 74 : SELECT CASE (active_integral_library%coulomb3_library)
75 : CASE (library_libint)
76 0 : CALL cp_libint_init_3eri(context%lib_3c, max_am_3c)
77 0 : CALL cp_libint_set_contrdepth(context%lib_3c, 1)
78 : CASE (library_native)
79 : ! Nothing to do.
80 : CASE DEFAULT
81 74 : CPABORT("Requested three-center Coulomb library is not implemented")
82 : END SELECT
83 74 : context%initialized_3c = .TRUE.
84 : END IF
85 260 : END SUBROUTINE coulomb_integral_init
86 :
87 : ! **************************************************************************************************
88 : !> \brief Release the engines initialized in a Coulomb integral context.
89 : !> \param context operation context containing the selected-library engines
90 : ! **************************************************************************************************
91 260 : SUBROUTINE coulomb_integral_cleanup(context)
92 : TYPE(coulomb_integral_context_type), INTENT(INOUT) :: context
93 :
94 260 : IF (context%initialized_2c) THEN
95 192 : SELECT CASE (active_integral_library%coulomb2_library)
96 : CASE (library_libint)
97 6 : CALL cp_libint_cleanup_2eri(context%lib_2c)
98 : CASE (library_native)
99 : ! Nothing to do.
100 : CASE DEFAULT
101 186 : CPABORT("Requested two-center Coulomb library is not implemented")
102 : END SELECT
103 186 : context%initialized_2c = .FALSE.
104 : END IF
105 :
106 260 : IF (context%initialized_3c) THEN
107 74 : SELECT CASE (active_integral_library%coulomb3_library)
108 : CASE (library_libint)
109 0 : CALL cp_libint_cleanup_3eri(context%lib_3c)
110 : CASE (library_native)
111 : ! Nothing to do.
112 : CASE DEFAULT
113 74 : CPABORT("Requested three-center Coulomb library is not implemented")
114 : END SELECT
115 74 : context%initialized_3c = .FALSE.
116 : END IF
117 260 : END SUBROUTINE coulomb_integral_cleanup
118 :
119 : ! **************************************************************************************************
120 : !> \brief Evaluate one uncontracted two-center Coulomb integral block.
121 : !> \param context operation context containing the selected-library engine
122 : !> \param la_min minimum angular momentum on center A
123 : !> \param la_max maximum angular momentum on center A
124 : !> \param lb_min minimum angular momentum on center B
125 : !> \param lb_max maximum angular momentum on center B
126 : !> \param npgfa number of primitive functions on center A
127 : !> \param npgfb number of primitive functions on center B
128 : !> \param zeta exponents on center A
129 : !> \param zetb exponents on center B
130 : !> \param rpgfa primitive radii on center A
131 : !> \param rpgfb primitive radii on center B
132 : !> \param ra position of center A
133 : !> \param rb position of center B
134 : !> \param hab uncontracted integral block
135 : !> \param potential_parameter Coulomb operator parameters
136 : ! **************************************************************************************************
137 175652 : SUBROUTINE compute_coulomb_2c(context, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, rpgfa, rpgfb, &
138 87826 : ra, rb, hab, potential_parameter)
139 : TYPE(coulomb_integral_context_type), INTENT(INOUT) :: context
140 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, npgfa, &
141 : npgfb
142 : REAL(dp), DIMENSION(npgfa), INTENT(IN) :: zeta
143 : REAL(dp), DIMENSION(npgfb), INTENT(IN) :: zetb
144 : REAL(dp), DIMENSION(npgfa), INTENT(IN) :: rpgfa
145 : REAL(dp), DIMENSION(npgfb), INTENT(IN) :: rpgfb
146 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra, rb
147 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: hab
148 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
149 :
150 : INTEGER :: ncoa, ncob
151 : REAL(dp) :: rab2
152 87826 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: f_work
153 87826 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: v_work
154 :
155 87826 : ncoa = npgfa*ncoset(la_max)
156 87826 : ncob = npgfb*ncoset(lb_max)
157 87826 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
158 :
159 90797 : SELECT CASE (active_integral_library%coulomb2_library)
160 : CASE (library_libint)
161 2971 : CPASSERT(context%initialized_2c)
162 : CALL eri_2center(hab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
163 2971 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, SQRT(rab2), context%lib_2c, potential_parameter)
164 : CASE (library_native)
165 254565 : ALLOCATE (f_work(0:la_max + lb_max + 2))
166 424275 : ALLOCATE (v_work(ncoa, ncob, la_max + lb_max + 1))
167 :
168 : CALL coulomb2(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, &
169 339420 : rb - ra, rab2, hab, v_work, f_work, screening=.FALSE.)
170 :
171 84855 : DEALLOCATE (v_work, f_work)
172 : CASE DEFAULT
173 87826 : CPABORT("Requested two-center Coulomb library is not implemented")
174 : END SELECT
175 87826 : END SUBROUTINE compute_coulomb_2c
176 :
177 : ! **************************************************************************************************
178 : !> \brief Evaluate one uncontracted three-center Coulomb integral block.
179 : !> \param context operation context containing the selected-library engine
180 : !> \param la_min minimum angular momentum on center A
181 : !> \param la_max maximum angular momentum on center A
182 : !> \param lb_min minimum angular momentum on center B
183 : !> \param lb_max maximum angular momentum on center B
184 : !> \param lc_min minimum angular momentum on center C
185 : !> \param lc_max maximum angular momentum on center C
186 : !> \param npgfa number of primitive functions on center A
187 : !> \param npgfb number of primitive functions on center B
188 : !> \param npgfc number of primitive functions on center C
189 : !> \param zeta exponents on center A
190 : !> \param zetb exponents on center B
191 : !> \param zetc exponents on center C
192 : !> \param rpgfa primitive radii on center A
193 : !> \param rpgfb primitive radii on center B
194 : !> \param rpgfc primitive radii on center C
195 : !> \param ra position of center A
196 : !> \param rb position of center B
197 : !> \param rc position of center C
198 : !> \param habc uncontracted integral block
199 : !> \param potential_parameter Coulomb operator parameters
200 : ! **************************************************************************************************
201 90554 : SUBROUTINE compute_coulomb_3c(context, la_min, la_max, lb_min, lb_max, lc_min, lc_max, npgfa, npgfb, npgfc, &
202 90554 : zeta, zetb, zetc, rpgfa, rpgfb, rpgfc, ra, rb, rc, habc, potential_parameter)
203 : TYPE(coulomb_integral_context_type), INTENT(INOUT) :: context
204 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
205 : lc_max, npgfa, npgfb, npgfc
206 : REAL(dp), DIMENSION(npgfa), INTENT(IN) :: zeta
207 : REAL(dp), DIMENSION(npgfb), INTENT(IN) :: zetb
208 : REAL(dp), DIMENSION(npgfc), INTENT(IN) :: zetc
209 : REAL(dp), DIMENSION(npgfa), INTENT(IN) :: rpgfa
210 : REAL(dp), DIMENSION(npgfb), INTENT(IN) :: rpgfb
211 : REAL(dp), DIMENSION(npgfc), INTENT(IN) :: rpgfc
212 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra, rb, rc
213 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
214 : TYPE(coulomb_operator_type), INTENT(IN) :: potential_parameter
215 :
216 : INTEGER :: kpgf, nc_end, nc_start, ncoa, ncob, ncoc
217 : REAL(dp) :: rab2, rac2, rbc2
218 90554 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: f_work, gccc, rpgfa_work, rpgfb_work, &
219 90554 : rpgfc_work
220 90554 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vabc
221 90554 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: v_work
222 :
223 90554 : ncoa = npgfa*ncoset(la_max)
224 90554 : ncob = npgfb*ncoset(lb_max)
225 90554 : ncoc = npgfc*ncoset(lc_max)
226 90554 : rab2 = (rb(1) - ra(1))**2 + (rb(2) - ra(2))**2 + (rb(3) - ra(3))**2
227 90554 : rac2 = (rc(1) - ra(1))**2 + (rc(2) - ra(2))**2 + (rc(3) - ra(3))**2
228 90554 : rbc2 = (rc(1) - rb(1))**2 + (rc(2) - rb(2))**2 + (rc(3) - rb(3))**2
229 :
230 90554 : SELECT CASE (active_integral_library%coulomb3_library)
231 : CASE (library_libint)
232 0 : CPASSERT(context%initialized_3c)
233 : CALL eri_3center(habc, la_min, la_max, npgfa, zeta, rpgfa, ra, &
234 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, &
235 : lc_min, lc_max, npgfc, zetc, rpgfc, rc, &
236 0 : SQRT(rab2), SQRT(rac2), SQRT(rbc2), context%lib_3c, potential_parameter)
237 : CASE (library_native)
238 271662 : ALLOCATE (f_work(0:la_max + lb_max + lc_max + 2))
239 90554 : f_work(:) = 0.0_dp
240 543324 : ALLOCATE (v_work(ncoa, ncob, ncoc, la_max + lb_max + lc_max + 1))
241 90554 : v_work(:, :, :, :) = 0.0_dp
242 :
243 : ! no screening
244 271662 : ALLOCATE (rpgfa_work(npgfa))
245 271662 : ALLOCATE (rpgfb_work(npgfb))
246 271662 : ALLOCATE (rpgfc_work(npgfc))
247 240788 : rpgfa_work(:) = 1.0E10_dp
248 240788 : rpgfb_work(:) = 1.0E10_dp
249 182556 : rpgfc_work(:) = 1.0E10_dp
250 271662 : ALLOCATE (gccc(ncoc))
251 90554 : gccc(:) = 0.0_dp
252 362216 : ALLOCATE (vabc(ncoa, ncob))
253 90554 : vabc(:, :) = 0.0_dp
254 :
255 : ! in the RI basis, there is only a single primitive Gaussian
256 90554 : kpgf = 1
257 :
258 : ! ncoset is indexed by -1,..., ,5
259 : ! e.g. for pure s: lc_min = lc_max = 0, nc_start = 1
260 90554 : nc_start = ncoset(lc_min - 1) + 1
261 90554 : nc_end = ncoset(lc_max)
262 :
263 : CALL coulomb3(la_max, npgfa, zeta(:), rpgfa_work(:), la_min, &
264 : lb_max, npgfb, zetb(:), rpgfb_work(:), lb_min, &
265 : lc_max, zetc(kpgf), rpgfc_work(kpgf), lc_min, &
266 : gccc, rb - ra, rab2, rc - ra, rac2, rbc2, &
267 633878 : vabc, habc(:, :, nc_start:nc_end), v_work, f_work)
268 :
269 90554 : DEALLOCATE (v_work, f_work, rpgfa_work, rpgfb_work, rpgfc_work, gccc, vabc)
270 : CASE DEFAULT
271 90554 : CPABORT("Requested three-center Coulomb library is not implemented")
272 : END SELECT
273 90554 : END SUBROUTINE compute_coulomb_3c
274 :
275 0 : END MODULE coulomb_integral_interface
|
