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 280 : 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 280 : IF (PRESENT(max_am_2c)) THEN
61 396 : SELECT CASE (active_integral_library%coulomb2_library)
62 : CASE (library_libint)
63 196 : CALL cp_libint_init_2eri(context%lib_2c, max_am_2c)
64 196 : CALL cp_libint_set_contrdepth(context%lib_2c, 1)
65 : CASE (library_native)
66 : ! Nothing to do.
67 : CASE DEFAULT
68 200 : CPABORT("Requested two-center Coulomb library is not implemented")
69 : END SELECT
70 200 : context%initialized_2c = .TRUE.
71 : END IF
72 :
73 280 : IF (PRESENT(max_am_3c)) THEN
74 174 : SELECT CASE (active_integral_library%coulomb3_library)
75 : CASE (library_libint)
76 86 : CALL cp_libint_init_3eri(context%lib_3c, max_am_3c)
77 86 : CALL cp_libint_set_contrdepth(context%lib_3c, 1)
78 : CASE (library_native)
79 : ! Nothing to do.
80 : CASE DEFAULT
81 88 : CPABORT("Requested three-center Coulomb library is not implemented")
82 : END SELECT
83 88 : context%initialized_3c = .TRUE.
84 : END IF
85 280 : 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 280 : SUBROUTINE coulomb_integral_cleanup(context)
92 : TYPE(coulomb_integral_context_type), INTENT(INOUT) :: context
93 :
94 280 : IF (context%initialized_2c) THEN
95 396 : SELECT CASE (active_integral_library%coulomb2_library)
96 : CASE (library_libint)
97 196 : CALL cp_libint_cleanup_2eri(context%lib_2c)
98 : CASE (library_native)
99 : ! Nothing to do.
100 : CASE DEFAULT
101 200 : CPABORT("Requested two-center Coulomb library is not implemented")
102 : END SELECT
103 200 : context%initialized_2c = .FALSE.
104 : END IF
105 :
106 280 : IF (context%initialized_3c) THEN
107 174 : SELECT CASE (active_integral_library%coulomb3_library)
108 : CASE (library_libint)
109 86 : CALL cp_libint_cleanup_3eri(context%lib_3c)
110 : CASE (library_native)
111 : ! Nothing to do.
112 : CASE DEFAULT
113 88 : CPABORT("Requested three-center Coulomb library is not implemented")
114 : END SELECT
115 88 : context%initialized_3c = .FALSE.
116 : END IF
117 280 : 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 188944 : SUBROUTINE compute_coulomb_2c(context, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, rpgfa, rpgfb, &
138 94472 : 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 94472 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: f_work
153 94472 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: v_work
154 :
155 94472 : ncoa = npgfa*ncoset(la_max)
156 94472 : ncob = npgfb*ncoset(lb_max)
157 94472 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
158 :
159 187540 : SELECT CASE (active_integral_library%coulomb2_library)
160 : CASE (library_libint)
161 93068 : CPASSERT(context%initialized_2c)
162 : CALL eri_2center(hab, la_min, la_max, npgfa, zeta, rpgfa, ra, &
163 93068 : lb_min, lb_max, npgfb, zetb, rpgfb, rb, SQRT(rab2), context%lib_2c, potential_parameter)
164 : CASE (library_native)
165 4212 : ALLOCATE (f_work(0:la_max + lb_max + 2))
166 7020 : 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 5616 : rb - ra, rab2, hab, v_work, f_work, screening=.FALSE.)
170 :
171 1404 : DEALLOCATE (v_work, f_work)
172 : CASE DEFAULT
173 94472 : CPABORT("Requested two-center Coulomb library is not implemented")
174 : END SELECT
175 94472 : 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 95922 : SUBROUTINE compute_coulomb_3c(context, la_min, la_max, lb_min, lb_max, lc_min, lc_max, npgfa, npgfb, npgfc, &
202 95922 : 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 95922 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: f_work, gccc, rpgfa_work, rpgfb_work, &
219 95922 : rpgfc_work
220 95922 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: vabc
221 95922 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: v_work
222 :
223 95922 : ncoa = npgfa*ncoset(la_max)
224 95922 : ncob = npgfb*ncoset(lb_max)
225 95922 : ncoc = ncoset(lc_max)
226 95922 : rab2 = (rb(1) - ra(1))**2 + (rb(2) - ra(2))**2 + (rb(3) - ra(3))**2
227 95922 : rac2 = (rc(1) - ra(1))**2 + (rc(2) - ra(2))**2 + (rc(3) - ra(3))**2
228 95922 : rbc2 = (rc(1) - rb(1))**2 + (rc(2) - rb(2))**2 + (rc(3) - rb(3))**2
229 :
230 190224 : SELECT CASE (active_integral_library%coulomb3_library)
231 : CASE (library_libint)
232 94302 : 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 94302 : SQRT(rab2), SQRT(rac2), SQRT(rbc2), context%lib_3c, potential_parameter)
237 : CASE (library_native)
238 4860 : ALLOCATE (f_work(0:la_max + lb_max + lc_max + 2))
239 1620 : f_work(:) = 0.0_dp
240 9720 : ALLOCATE (v_work(ncoa, ncob, ncoc, la_max + lb_max + lc_max + 1))
241 1620 : v_work(:, :, :, :) = 0.0_dp
242 :
243 : ! no screening
244 4860 : ALLOCATE (rpgfa_work(npgfa))
245 4860 : ALLOCATE (rpgfb_work(npgfb))
246 4860 : ALLOCATE (rpgfc_work(npgfc))
247 4320 : rpgfa_work(:) = 1.0E10_dp
248 4320 : rpgfb_work(:) = 1.0E10_dp
249 3240 : rpgfc_work(:) = 1.0E10_dp
250 4860 : ALLOCATE (gccc(ncoc))
251 1620 : gccc(:) = 0.0_dp
252 6480 : ALLOCATE (vabc(ncoa, ncob))
253 1620 : vabc(:, :) = 0.0_dp
254 :
255 : ! coulomb3 handles a single primitive Gaussian on center C per call.
256 3240 : DO kpgf = 1, npgfc
257 : ! Each primitive occupies an ncoc-wide block in the uncontracted array.
258 1620 : nc_start = (kpgf - 1)*ncoc + ncoset(lc_min - 1) + 1
259 1620 : nc_end = kpgf*ncoc
260 :
261 : CALL coulomb3(la_max, npgfa, zeta(:), rpgfa_work(:), la_min, &
262 : lb_max, npgfb, zetb(:), rpgfb_work(:), lb_min, &
263 : lc_max, zetc(kpgf), rpgfc_work(kpgf), lc_min, &
264 : gccc, rb - ra, rab2, rc - ra, rac2, rbc2, &
265 12960 : vabc, habc(:, :, nc_start:nc_end), v_work, f_work)
266 : END DO
267 :
268 1620 : DEALLOCATE (v_work, f_work, rpgfa_work, rpgfb_work, rpgfc_work, gccc, vabc)
269 : CASE DEFAULT
270 95922 : CPABORT("Requested three-center Coulomb library is not implemented")
271 : END SELECT
272 95922 : END SUBROUTINE compute_coulomb_3c
273 :
274 0 : END MODULE coulomb_integral_interface
|
