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 Interface to direct methods for electron repulsion integrals for MP2.
10 : ! **************************************************************************************************
11 : #:def conditional(n)
12 : $:'' if n else '.NOT.'
13 : #:enddef
14 :
15 : MODULE mp2_eri
16 : USE ai_contraction_sphi, ONLY: ab_contract, &
17 : abc_contract
18 : USE atomic_kind_types, ONLY: atomic_kind_type, &
19 : get_atomic_kind_set
20 : USE basis_set_types, ONLY: gto_basis_set_p_type, &
21 : gto_basis_set_type
22 : USE cell_types, ONLY: cell_type, &
23 : pbc
24 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_finalize, &
25 : cp_eri_mme_init_read_input, &
26 : cp_eri_mme_param, &
27 : cp_eri_mme_set_params
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p, &
30 : dbcsr_p_type
31 : USE eri_mme_integrate, ONLY: eri_mme_2c_integrate, &
32 : eri_mme_3c_integrate
33 : USE eri_mme_test, ONLY: eri_mme_2c_perf_acc_test, &
34 : eri_mme_3c_perf_acc_test
35 : USE eri_mme_types, ONLY: eri_mme_param, &
36 : eri_mme_set_potential, eri_mme_coulomb, eri_mme_longrange
37 : USE input_constants, ONLY: do_eri_gpw, &
38 : do_eri_mme, &
39 : do_eri_os, &
40 : do_potential_coulomb, &
41 : do_potential_long
42 : USE input_section_types, ONLY: section_vals_get_subs_vals, &
43 : section_vals_type, &
44 : section_vals_val_get
45 : USE kinds, ONLY: dp
46 : USE coulomb_integral_interface, ONLY: compute_coulomb_2c, &
47 : compute_coulomb_3c, &
48 : coulomb_integral_cleanup, &
49 : coulomb_integral_context_type, &
50 : coulomb_integral_init
51 : USE libint_2c_3c, ONLY: libint_potential_type
52 : USE orbital_pointers, ONLY: coset, &
53 : init_orbital_pointers, &
54 : ncoset
55 : USE particle_types, ONLY: particle_type
56 : USE qs_environment_types, ONLY: get_qs_env, &
57 : qs_environment_type
58 : USE qs_integral_utils, ONLY: basis_set_list_setup
59 : USE qs_kind_types, ONLY: get_qs_kind, &
60 : get_qs_kind_set, &
61 : qs_kind_type
62 : USE qs_neighbor_list_types, ONLY: get_iterator_info, &
63 : get_neighbor_list_set_p, &
64 : neighbor_list_iterate, &
65 : neighbor_list_iterator_create, &
66 : neighbor_list_iterator_p_type, &
67 : neighbor_list_iterator_release, &
68 : neighbor_list_set_p_type
69 : USE util, ONLY: get_limit
70 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_update_local_counts
71 : #include "./base/base_uses.f90"
72 :
73 : IMPLICIT NONE
74 :
75 : PRIVATE
76 :
77 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
78 :
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_eri'
80 :
81 : PUBLIC :: &
82 : mp2_eri_2c_integrate, &
83 : mp2_eri_3c_integrate, &
84 : mp2_eri_allocate_forces, &
85 : mp2_eri_deallocate_forces, &
86 : mp2_eri_force, &
87 : integrate_set_2c, &
88 : convert_potential_type
89 :
90 : TYPE mp2_eri_force
91 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: forces
92 : END TYPE mp2_eri_force
93 :
94 : CONTAINS
95 :
96 : ! **************************************************************************************************
97 : !> \brief high-level integration routine for 2c integrals over CP2K basis sets.
98 : !> Contiguous column-wise distribution and parallelization over pairs of sets.
99 : !> \param param ...
100 : !> \param para_env mpi environment for local columns
101 : !> \param potential_parameter ...
102 : !> \param qs_env ...
103 : !> \param basis_type_a ...
104 : !> \param basis_type_b ...
105 : !> \param hab columns of ERI matrix
106 : !> \param first_b first column of hab
107 : !> \param last_b last column of hab
108 : !> \param eri_method ...
109 : !> \param pab ...
110 : !> \param force_a ...
111 : !> \param force_b ...
112 : !> \param hdab ...
113 : !> \param hadb ...
114 : !> \param reflection_z_a ...
115 : !> \param reflection_z_b ...
116 : !> \param do_reflection_a ...
117 : !> \param do_reflection_b ...
118 : ! **************************************************************************************************
119 342 : SUBROUTINE mp2_eri_2c_integrate(param, potential_parameter, para_env, qs_env, basis_type_a, basis_type_b, hab, first_b, &
120 342 : last_b, eri_method, pab, force_a, force_b, hdab, hadb, &
121 : reflection_z_a, reflection_z_b, do_reflection_a, do_reflection_b)
122 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
123 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
124 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
125 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
126 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type_a, basis_type_b
127 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: hab
128 : INTEGER, INTENT(IN) :: first_b, last_b
129 : INTEGER, INTENT(IN), OPTIONAL :: eri_method
130 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
131 : OPTIONAL :: pab
132 : TYPE(mp2_eri_force), ALLOCATABLE, &
133 : DIMENSION(:), INTENT(OUT), OPTIONAL :: force_a, force_b
134 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
135 : OPTIONAL :: hdab, hadb
136 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: reflection_z_a, reflection_z_b
137 : LOGICAL, INTENT(IN), OPTIONAL :: do_reflection_a, do_reflection_b
138 :
139 : CHARACTER(len=*), PARAMETER :: routineN = 'mp2_eri_2c_integrate'
140 :
141 : INTEGER :: atom_a, atom_b, atom_end, atom_start, first_set, G_count, handle, iatom, ikind, &
142 : iset, jatom, jkind, jset, jset_end, jset_start, last_set, max_am_a, max_am_b, &
143 : my_eri_method, my_setpair, n_setpair, natom, nkind, nseta, nseta_total, &
144 : nsetb, nsetb_total, offset_a_end, &
145 : offset_a_start, offset_b_end, offset_b_start, R_count, set_end, set_offset_end, &
146 : set_offset_start, set_start, sgfa, sgfb
147 342 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
148 342 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eri_offsets
149 342 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
150 342 : npgfb, nsgfa, nsgfb
151 342 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
152 : LOGICAL :: map_it_here, my_do_reflection_a, &
153 : my_do_reflection_b
154 : REAL(KIND=dp) :: dab
155 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb
156 342 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
157 342 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
158 : TYPE(cell_type), POINTER :: cell
159 : TYPE(coulomb_integral_context_type) :: integral_context
160 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
161 342 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
162 342 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
163 :
164 342 : CALL timeset(routineN, handle)
165 :
166 342 : my_eri_method = do_eri_mme
167 342 : IF (PRESENT(eri_method)) my_eri_method = eri_method
168 :
169 342 : my_do_reflection_a = .FALSE.
170 342 : IF (PRESENT(do_reflection_a) .AND. PRESENT(reflection_z_a)) my_do_reflection_a = do_reflection_a
171 :
172 342 : my_do_reflection_b = .FALSE.
173 342 : IF (PRESENT(do_reflection_b) .AND. PRESENT(reflection_z_b)) my_do_reflection_b = do_reflection_b
174 :
175 342 : G_count = 0; R_count = 0
176 : ! get mapping between ERIs and atoms, sets, set offsets
177 342 : CALL get_eri_offsets(qs_env, basis_type_b, eri_offsets)
178 :
179 342 : atom_start = eri_offsets(first_b, 1)
180 342 : set_start = eri_offsets(first_b, 2)
181 342 : set_offset_start = eri_offsets(first_b, 3)
182 :
183 342 : atom_end = eri_offsets(last_b, 1)
184 342 : set_end = eri_offsets(last_b, 2)
185 342 : set_offset_end = eri_offsets(last_b, 3)
186 :
187 : ! get QS stuff
188 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
189 342 : cell=cell, particle_set=particle_set, natom=natom, nkind=nkind)
190 :
191 342 : IF (my_eri_method == do_eri_os) THEN
192 186 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=max_am_a, basis_type=basis_type_a)
193 186 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=max_am_b, basis_type=basis_type_b)
194 186 : CALL coulomb_integral_init(integral_context, max_am_2c=MAX(max_am_a, max_am_b))
195 : END IF
196 :
197 342 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, natom_of_kind=natom_of_kind, atom_of_kind=atom_of_kind)
198 :
199 342 : IF (PRESENT(force_a)) CALL mp2_eri_allocate_forces(force_a, natom_of_kind)
200 342 : IF (PRESENT(force_b)) CALL mp2_eri_allocate_forces(force_b, natom_of_kind)
201 :
202 : ! get total number of local set pairs to integrate
203 342 : nseta_total = 0
204 1304 : DO iatom = 1, natom
205 962 : ikind = kind_of(iatom)
206 962 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type_a)
207 1304 : nseta_total = nseta_total + basis_set_a%nset
208 : END DO
209 :
210 : nsetb_total = 0
211 954 : DO jatom = atom_start, atom_end
212 612 : jkind = kind_of(jatom)
213 612 : CALL get_qs_kind(qs_kind=qs_kind_set(jkind), basis_set=basis_set_b, basis_type=basis_type_b)
214 954 : nsetb_total = nsetb_total + basis_set_b%nset
215 : END DO
216 :
217 : n_setpair = nseta_total*nsetb_total
218 :
219 342 : my_setpair = 0
220 :
221 342 : offset_a_end = 0
222 1304 : DO iatom = 1, natom
223 962 : ikind = kind_of(iatom)
224 962 : atom_a = atom_of_kind(iatom)
225 962 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type_a)
226 :
227 962 : first_sgfa => basis_set_a%first_sgf
228 962 : la_max => basis_set_a%lmax
229 962 : la_min => basis_set_a%lmin
230 962 : nseta = basis_set_a%nset
231 962 : nsgfa => basis_set_a%nsgf_set
232 962 : sphi_a => basis_set_a%sphi
233 962 : rpgfa => basis_set_a%pgf_radius
234 962 : zeta => basis_set_a%zet
235 962 : npgfa => basis_set_a%npgf
236 :
237 962 : ra(:) = pbc(particle_set(iatom)%r, cell)
238 :
239 962 : IF (my_do_reflection_a) THEN
240 0 : ra(3) = 2.0_dp*reflection_z_a - ra(3)
241 : END IF
242 :
243 10114 : DO iset = 1, nseta
244 8810 : offset_a_start = offset_a_end
245 8810 : offset_a_end = offset_a_end + nsgfa(iset)
246 8810 : sgfa = first_sgfa(1, iset)
247 :
248 8810 : offset_b_end = 0
249 26943 : DO jatom = atom_start, atom_end
250 17171 : jkind = kind_of(jatom)
251 17171 : atom_b = atom_of_kind(jatom)
252 17171 : CALL get_qs_kind(qs_kind=qs_kind_set(jkind), basis_set=basis_set_b, basis_type=basis_type_b)
253 :
254 17171 : first_sgfb => basis_set_b%first_sgf
255 17171 : lb_max => basis_set_b%lmax
256 17171 : lb_min => basis_set_b%lmin
257 17171 : nsetb = basis_set_b%nset
258 17171 : nsgfb => basis_set_b%nsgf_set
259 17171 : sphi_b => basis_set_b%sphi
260 17171 : rpgfb => basis_set_b%pgf_radius
261 17171 : zetb => basis_set_b%zet
262 17171 : npgfb => basis_set_b%npgf
263 :
264 17171 : rb(:) = pbc(particle_set(jatom)%r, cell)
265 :
266 17171 : IF (my_do_reflection_b) THEN
267 0 : rb(3) = 2.0_dp*reflection_z_b - rb(3)
268 : END IF
269 :
270 68684 : rab(:) = ra(:) - rb(:) ! pbc not needed?
271 : dab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
272 :
273 17171 : jset_start = 1; jset_end = nsetb
274 17171 : IF (jatom == atom_start) jset_start = set_start
275 17171 : IF (jatom == atom_end) jset_end = set_end
276 :
277 160294 : DO jset = jset_start, jset_end
278 134313 : first_set = 1; last_set = nsgfb(jset)
279 134313 : IF (jset == jset_start .AND. jatom == atom_start) first_set = set_offset_start
280 134313 : IF (jset == jset_end .AND. jatom == atom_end) last_set = set_offset_end
281 :
282 134313 : offset_b_start = offset_b_end
283 134313 : offset_b_end = offset_b_end + last_set + 1 - first_set
284 134313 : sgfb = first_sgfb(1, jset)
285 134313 : my_setpair = my_setpair + 1
286 134313 : map_it_here = MODULO(my_setpair, para_env%num_pe) == para_env%mepos
287 :
288 151484 : IF (map_it_here) THEN
289 : #!some fypp magic to deal with combinations of optional arguments
290 : #:for doforce_1 in [0, 1]
291 : #:for doforce_2 in [0, 1]
292 262320 : IF (${conditional(doforce_1)}$PRESENT(force_a) .AND. &
293 : ${conditional(doforce_2)}$PRESENT(force_b)) THEN
294 :
295 : CALL integrate_set_2c( &
296 : param%par, potential_parameter, &
297 : la_min(iset), la_max(iset), &
298 : lb_min(jset), lb_max(jset), &
299 : npgfa(iset), npgfb(jset), &
300 : zeta(:, iset), zetb(:, jset), &
301 : ra, rb, &
302 : hab, nsgfa(iset), last_set - first_set + 1, &
303 : offset_a_start, offset_b_start, &
304 : 0, first_set - 1, &
305 : sphi_a, sphi_b, &
306 : sgfa, sgfb, nsgfa(iset), nsgfb(jset), &
307 : my_eri_method, &
308 : pab=pab, &
309 : $: 'force_a=force_a(ikind)%forces(:, atom_a), &'*doforce_1
310 : $: 'force_b=force_b(jkind)%forces(:, atom_b), &'*doforce_2
311 : hdab=hdab, hadb=hadb, &
312 : G_count=G_count, R_count=R_count, &
313 : do_reflection_a=do_reflection_a, do_reflection_b=do_reflection_b, &
314 : coulomb_context=integral_context, &
315 514496 : rpgfa=rpgfa(:, iset), rpgfb=rpgfb(:, jset))
316 : END IF
317 : #:endfor
318 : #:endfor
319 : END IF
320 : END DO
321 : END DO
322 : END DO
323 : END DO
324 :
325 342 : IF (my_eri_method == do_eri_os) CALL coulomb_integral_cleanup(integral_context)
326 :
327 342 : IF (my_eri_method == do_eri_mme) THEN
328 :
329 156 : CALL cp_eri_mme_update_local_counts(param, para_env, G_count_2c=G_count, R_count_2c=R_count)
330 :
331 : END IF
332 :
333 2196142 : CALL para_env%sum(hab)
334 342 : IF (PRESENT(hdab)) CALL para_env%sum(hdab)
335 342 : IF (PRESENT(hadb)) CALL para_env%sum(hadb)
336 :
337 342 : CALL timestop(handle)
338 684 : END SUBROUTINE mp2_eri_2c_integrate
339 :
340 : ! **************************************************************************************************
341 : !> \brief Integrate set pair and contract with sphi matrix.
342 : !> \param param ...
343 : !> \param potential_parameter ...
344 : !> \param la_min ...
345 : !> \param la_max ...
346 : !> \param lb_min ...
347 : !> \param lb_max ...
348 : !> \param npgfa ...
349 : !> \param npgfb ...
350 : !> \param zeta ...
351 : !> \param zetb ...
352 : !> \param ra ...
353 : !> \param rb ...
354 : !> \param hab ...
355 : !> \param n_hab_a ...
356 : !> \param n_hab_b ...
357 : !> \param offset_hab_a ...
358 : !> \param offset_hab_b ...
359 : !> \param offset_set_a ...
360 : !> \param offset_set_b ...
361 : !> \param sphi_a ...
362 : !> \param sphi_b ...
363 : !> \param sgfa ...
364 : !> \param sgfb ...
365 : !> \param nsgfa ...
366 : !> \param nsgfb ...
367 : !> \param eri_method ...
368 : !> \param pab ...
369 : !> \param force_a ...
370 : !> \param force_b ...
371 : !> \param hdab ...
372 : !> \param hadb ...
373 : !> \param G_count ...
374 : !> \param R_count ...
375 : !> \param do_reflection_a ...
376 : !> \param do_reflection_b ...
377 : !> \param rpgfa primitive Gaussian radii for basis A (required by the Libint library)
378 : !> \param rpgfb primitive Gaussian radii for basis B (required by the Libint library)
379 : !> \param coulomb_context context for the selected Coulomb integral library
380 : ! **************************************************************************************************
381 273616 : SUBROUTINE integrate_set_2c(param, potential_parameter, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, &
382 136808 : ra, rb, hab, n_hab_a, n_hab_b, offset_hab_a, offset_hab_b, &
383 136808 : offset_set_a, offset_set_b, sphi_a, sphi_b, sgfa, sgfb, nsgfa, nsgfb, &
384 136808 : eri_method, pab, force_a, force_b, hdab, hadb, G_count, R_count, &
385 136808 : do_reflection_a, do_reflection_b, rpgfa, rpgfb, coulomb_context)
386 : TYPE(eri_mme_param), INTENT(INOUT) :: param
387 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
388 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, npgfa
389 : REAL(KIND=dp), DIMENSION(npgfa), INTENT(IN) :: zeta
390 : INTEGER, INTENT(IN) :: npgfb
391 : REAL(KIND=dp), DIMENSION(npgfb), INTENT(IN) :: zetb
392 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
393 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: hab
394 : INTEGER, INTENT(IN) :: n_hab_a, n_hab_b, offset_hab_a, &
395 : offset_hab_b, offset_set_a, &
396 : offset_set_b
397 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sphi_a
398 : INTEGER, INTENT(IN) :: sgfa, nsgfa
399 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sphi_b
400 : INTEGER, INTENT(IN) :: sgfb, nsgfb, eri_method
401 : TYPE(coulomb_integral_context_type), INTENT(INOUT), OPTIONAL :: coulomb_context
402 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
403 : OPTIONAL :: pab
404 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
405 : OPTIONAL :: force_a, force_b
406 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT), &
407 : OPTIONAL :: hdab, hadb
408 : INTEGER, INTENT(INOUT), OPTIONAL :: G_count, R_count
409 : LOGICAL, INTENT(IN), OPTIONAL :: do_reflection_a, do_reflection_b
410 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: rpgfa, rpgfb
411 :
412 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_set_2c'
413 :
414 : INTEGER :: ax, ay, az, bx, by, bz, hab_a_end, hab_a_start, hab_b_end, hab_b_start, handle, &
415 : i_xyz, ico, icox, icoy, icoz, ipgf, jco, jcox, jcoy, jcoz, jpgf, la, la_max_d, lb, &
416 : lb_max_d, na, nb, ncoa, ncob, set_a_end, set_a_start, set_b_end, set_b_start, &
417 : sphi_a_start, sphi_b_start
418 : INTEGER, DIMENSION(3) :: la_xyz, lb_xyz
419 : LOGICAL :: calculate_forces, my_do_reflection_a, &
420 : my_do_reflection_b, do_force_a, do_force_b
421 136808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: hab_contr, hab_uncontr, &
422 136808 : hab_uncontr_d, pab_hh, pab_hs, &
423 136808 : pab_ss
424 136808 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hadb_contr, hadb_uncontr, hdab_contr, &
425 136808 : hdab_uncontr
426 :
427 : ! note: tested only for one exponent per pair (npgfa = npgfb = 1)
428 136808 : CALL timeset(routineN, handle)
429 :
430 136808 : my_do_reflection_a = .FALSE.
431 136808 : IF (PRESENT(do_reflection_a)) my_do_reflection_a = do_reflection_a
432 :
433 136808 : my_do_reflection_b = .FALSE.
434 136808 : IF (PRESENT(do_reflection_b)) my_do_reflection_b = do_reflection_b
435 :
436 136808 : do_force_a = PRESENT(force_a) .OR. PRESENT(hdab)
437 136808 : do_force_b = PRESENT(force_b) .OR. PRESENT(hadb)
438 136808 : calculate_forces = do_force_a .OR. do_force_b
439 :
440 136808 : IF (PRESENT(force_a) .OR. PRESENT(force_b)) THEN
441 10144 : CPASSERT(PRESENT(pab))
442 30432 : CPASSERT(ALL(SHAPE(pab) == SHAPE(hab)))
443 : END IF
444 :
445 136808 : la_max_d = la_max
446 136808 : lb_max_d = lb_max
447 :
448 136808 : IF (calculate_forces) THEN
449 15792 : IF (do_force_a) la_max_d = la_max + 1
450 15792 : IF (do_force_b) lb_max_d = lb_max + 1
451 : END IF
452 :
453 136808 : ncoa = npgfa*ncoset(la_max)
454 136808 : ncob = npgfb*ncoset(lb_max)
455 :
456 547232 : ALLOCATE (hab_uncontr_d(ncoset(la_max_d), ncoset(lb_max_d))); hab_uncontr_d(:, :) = 0.0_dp
457 547232 : ALLOCATE (hab_uncontr(ncoa, ncob)); hab_uncontr(:, :) = 0.0_dp
458 136808 : IF (PRESENT(hdab)) THEN
459 22592 : ALLOCATE (hdab_uncontr(3, ncoa, ncob)); hdab_uncontr(:, :, :) = 0.0_dp
460 : END IF
461 136808 : IF (PRESENT(hadb)) THEN
462 0 : ALLOCATE (hadb_uncontr(3, ncoa, ncob)); hadb_uncontr(:, :, :) = 0.0_dp
463 : END IF
464 :
465 136808 : hab_a_start = offset_hab_a + 1; hab_a_end = offset_hab_a + n_hab_a
466 136808 : hab_b_start = offset_hab_b + 1; hab_b_end = offset_hab_b + n_hab_b
467 136808 : set_a_start = offset_set_a + 1; set_a_end = offset_set_a + n_hab_a
468 136808 : set_b_start = offset_set_b + 1; set_b_end = offset_set_b + n_hab_b
469 :
470 136808 : IF (eri_method == do_eri_mme) THEN
471 48982 : CALL eri_mme_set_potential(param, convert_potential_type(potential_parameter%potential_type), potential_parameter%omega)
472 :
473 48982 : IF (calculate_forces .AND. PRESENT(pab)) THEN
474 : ! uncontracted hermite-gaussian representation of density matrix
475 10144 : sphi_a_start = sgfa - 1 + set_a_start
476 10144 : sphi_b_start = sgfb - 1 + set_b_start
477 :
478 40576 : ALLOCATE (pab_ss(n_hab_a, n_hab_b))
479 121184 : pab_ss(:, :) = pab(hab_a_start:hab_a_end, hab_b_start:hab_b_end)
480 60864 : ALLOCATE (pab_hs(ncoa, n_hab_b)); ALLOCATE (pab_hh(ncoa, ncob))
481 : CALL dgemm("N", "N", ncoa, n_hab_b, n_hab_a, 1.0_dp, &
482 10144 : sphi_a(:, sphi_a_start), SIZE(sphi_a, 1), pab_ss, n_hab_a, 0.0_dp, pab_hs, ncoa)
483 : CALL dgemm("N", "T", ncoa, ncob, n_hab_b, 1.0_dp, &
484 10144 : pab_hs, ncoa, sphi_b(:, sphi_b_start), SIZE(sphi_b, 1), 0.0_dp, pab_hh, ncoa)
485 : END IF
486 :
487 97964 : DO ipgf = 1, npgfa
488 48982 : na = (ipgf - 1)*ncoset(la_max)
489 146946 : DO jpgf = 1, npgfb
490 48982 : nb = (jpgf - 1)*ncoset(lb_max)
491 48982 : hab_uncontr_d(:, :) = 0.0_dp
492 : CALL eri_mme_2c_integrate(param, &
493 : la_min, la_max_d, lb_min, lb_max_d, &
494 195928 : zeta(ipgf), zetb(jpgf), ra - rb, hab_uncontr_d, 0, 0, G_count, R_count)
495 :
496 : hab_uncontr(na + 1:na + ncoset(la_max), nb + 1:nb + ncoset(lb_max)) = &
497 1599874 : hab_uncontr_d(:ncoset(la_max), :ncoset(lb_max))
498 :
499 97964 : IF (calculate_forces) THEN
500 31584 : DO lb = lb_min, lb_max
501 62368 : DO bx = 0, lb
502 99346 : DO by = 0, lb - bx
503 52770 : bz = lb - bx - by
504 52770 : jco = coset(bx, by, bz)
505 52770 : jcox = coset(bx + 1, by, bz)
506 52770 : jcoy = coset(bx, by + 1, bz)
507 52770 : jcoz = coset(bx, by, bz + 1)
508 136324 : DO la = la_min, la_max
509 209820 : DO ax = 0, la
510 336900 : DO ay = 0, la - ax
511 179850 : az = la - ax - ay
512 719400 : la_xyz = [ax, ay, az]
513 719400 : lb_xyz = [bx, by, bz]
514 179850 : ico = coset(ax, ay, az)
515 179850 : icox = coset(ax + 1, ay, az)
516 179850 : icoy = coset(ax, ay + 1, az)
517 179850 : icoz = coset(ax, ay, az + 1)
518 179850 : IF (PRESENT(force_a)) THEN
519 : force_a(:) = force_a(:) + 2.0_dp*zeta(ipgf)* &
520 : [pab_hh(na + ico, nb + jco)*hab_uncontr_d(icox, jco), &
521 : pab_hh(na + ico, nb + jco)*hab_uncontr_d(icoy, jco), &
522 473800 : pab_hh(na + ico, nb + jco)*hab_uncontr_d(icoz, jco)]
523 : END IF
524 179850 : IF (PRESENT(force_b)) THEN
525 : force_b(:) = force_b(:) + 2.0_dp*zetb(jpgf)* &
526 : [pab_hh(na + ico, nb + jco)*hab_uncontr_d(ico, jcox), &
527 : pab_hh(na + ico, nb + jco)*hab_uncontr_d(ico, jcoy), &
528 0 : pab_hh(na + ico, nb + jco)*hab_uncontr_d(ico, jcoz)]
529 : END IF
530 179850 : IF (PRESENT(hdab)) THEN
531 : hdab_uncontr(1:3, na + ico, nb + jco) = 2.0_dp*zeta(ipgf)* &
532 : [hab_uncontr_d(icox, jco), &
533 : hab_uncontr_d(icoy, jco), &
534 245600 : hab_uncontr_d(icoz, jco)]
535 : END IF
536 284130 : IF (PRESENT(hadb)) THEN
537 : hadb_uncontr(1:3, na + ico, nb + jco) = 2.0_dp*zetb(jpgf)* &
538 : [hab_uncontr_d(ico, jcox), &
539 : hab_uncontr_d(ico, jcoy), &
540 0 : hab_uncontr_d(ico, jcoz)]
541 : END IF
542 : END DO
543 : END DO
544 : END DO
545 : END DO
546 : END DO
547 : END DO
548 : END IF
549 :
550 : END DO
551 : END DO
552 :
553 87826 : ELSE IF (eri_method == do_eri_os) THEN
554 :
555 87826 : IF (calculate_forces) CPABORT("NYI")
556 87826 : IF (.NOT. PRESENT(rpgfa) .OR. .NOT. PRESENT(rpgfb)) THEN
557 0 : CPABORT("Two-center OS library requires primitive Gaussian radii")
558 : END IF
559 :
560 87826 : CPASSERT(PRESENT(coulomb_context))
561 : CALL compute_coulomb_2c(coulomb_context, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, &
562 87826 : rpgfa, rpgfb, ra, rb, hab_uncontr, potential_parameter)
563 :
564 0 : ELSE IF (eri_method == do_eri_gpw) THEN
565 :
566 0 : CPABORT("GPW not enabled in the ERI interface.")
567 :
568 : END IF
569 :
570 547232 : ALLOCATE (hab_contr(nsgfa, nsgfb))
571 136808 : IF (PRESENT(hdab)) THEN
572 22592 : ALLOCATE (hdab_contr(3, nsgfa, nsgfb))
573 : END IF
574 136808 : IF (PRESENT(hadb)) THEN
575 0 : ALLOCATE (hadb_contr(3, nsgfa, nsgfb))
576 : END IF
577 :
578 136808 : CALL ab_contract(hab_contr, hab_uncontr, sphi_a(:, sgfa:), sphi_b(:, sgfb:), ncoa, ncob, nsgfa, nsgfb)
579 :
580 136808 : IF (calculate_forces) THEN
581 63168 : DO i_xyz = 1, 3
582 47376 : IF (PRESENT(hdab)) THEN
583 : CALL ab_contract(hdab_contr(i_xyz, :, :), hdab_uncontr(i_xyz, :, :), &
584 16944 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), ncoa, ncob, nsgfa, nsgfb)
585 : END IF
586 63168 : IF (PRESENT(hadb)) THEN
587 : CALL ab_contract(hadb_contr(i_xyz, :, :), hadb_uncontr(i_xyz, :, :), &
588 0 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), ncoa, ncob, nsgfa, nsgfb)
589 : END IF
590 : END DO
591 : END IF
592 :
593 1620064 : hab(hab_a_start:hab_a_end, hab_b_start:hab_b_end) = hab_contr(set_a_start:set_a_end, set_b_start:set_b_end)
594 :
595 136808 : IF (calculate_forces) THEN
596 15792 : IF (PRESENT(hdab)) hdab(:, hab_a_start:hab_a_end, hab_b_start:hab_b_end) = &
597 207536 : hdab_contr(:, set_a_start:set_a_end, set_b_start:set_b_end)
598 15792 : IF (PRESENT(hadb)) hadb(:, hab_a_start:hab_a_end, hab_b_start:hab_b_end) = &
599 0 : hadb_contr(:, set_a_start:set_a_end, set_b_start:set_b_end)
600 : END IF
601 :
602 136808 : CALL timestop(handle)
603 :
604 273616 : END SUBROUTINE integrate_set_2c
605 :
606 : ! **************************************************************************************************
607 : !> \brief high-level integration routine for 3c integrals (ab|c) over CP2K basis sets.
608 : !> For each local function of c, (ab|c) is written to a DBCSR matrix mat_ab.
609 : !> \param param ...
610 : !> \param potential_parameter ...
611 : !> \param para_env ...
612 : !> \param qs_env ...
613 : !> \param first_c start index of local range of c
614 : !> \param last_c end index of local range of c
615 : !> \param mat_ab DBCSR matrices for each c
616 : !> \param basis_type_a ...
617 : !> \param basis_type_b ...
618 : !> \param basis_type_c ...
619 : !> \param sab_nl neighbor list for a, b
620 : !> \param eri_method ...
621 : !> \param pabc ...
622 : !> \param force_a ...
623 : !> \param force_b ...
624 : !> \param force_c ...
625 : !> \param mat_dabc ...
626 : !> \param mat_adbc ...
627 : !> \param mat_abdc ...
628 : ! **************************************************************************************************
629 210 : SUBROUTINE mp2_eri_3c_integrate(param, potential_parameter, para_env, qs_env, &
630 210 : first_c, last_c, mat_ab, &
631 : basis_type_a, basis_type_b, basis_type_c, &
632 : sab_nl, eri_method, &
633 210 : pabc, force_a, force_b, force_c, &
634 210 : mat_dabc, mat_adbc, mat_abdc)
635 : TYPE(cp_eri_mme_param), INTENT(INOUT) :: param
636 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
637 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
638 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
639 : INTEGER, INTENT(IN) :: first_c, last_c
640 : TYPE(dbcsr_p_type), DIMENSION(last_c - first_c + 1), &
641 : INTENT(INOUT) :: mat_ab
642 : CHARACTER(LEN=*), INTENT(IN) :: basis_type_a, basis_type_b, basis_type_c
643 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
644 : POINTER :: sab_nl
645 : INTEGER, INTENT(IN), OPTIONAL :: eri_method
646 : TYPE(dbcsr_p_type), DIMENSION(last_c - first_c + 1), &
647 : INTENT(INOUT), OPTIONAL :: pabc
648 : TYPE(mp2_eri_force), ALLOCATABLE, &
649 : DIMENSION(:), INTENT(OUT), OPTIONAL :: force_a, force_b, force_c
650 : TYPE(dbcsr_p_type), &
651 : DIMENSION(3, last_c - first_c + 1), INTENT(INOUT), &
652 : OPTIONAL :: mat_dabc, mat_adbc, mat_abdc
653 :
654 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_eri_3c_integrate'
655 :
656 : INTEGER :: atom_a, atom_b, atom_c, atom_end, atom_start, first_set, GG_count, GR_count, &
657 : handle, i_xyz, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, katom, &
658 : kkind, kset, kset_end, kset_start, last_jatom, last_set, max_am_a, max_am_b, max_am_c, &
659 : mepos, my_eri_method, na, natom, &
660 : nb, nc, nkind, nseta, nsetb, nsetc, nthread, offset_a_end, offset_a_start, offset_b_end, &
661 : offset_b_start, offset_c_end, offset_c_start, RR_count, set_end, set_offset_end, &
662 : set_offset_start, set_start, sgfa, sgfb, sgfc
663 210 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
664 210 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eri_offsets
665 210 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
666 210 : lc_min, npgfa, npgfb, npgfc, nsgfa, &
667 210 : nsgfb, nsgfc
668 210 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb, first_sgfc
669 : LOGICAL :: calculate_forces, do_symmetric, found, to_be_asserted
670 : REAL(KIND=dp) :: dab
671 210 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: habc, pabc_block
672 210 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: habdc, hadbc, hdabc
673 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb, rc
674 210 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
675 210 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: munu_block, pab_block, rpgfa, rpgfb, rpgfc, sphi_a, &
676 210 : sphi_b, sphi_c, zeta, zetb, zetc
677 210 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
678 : TYPE(cell_type), POINTER :: cell
679 210 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list_a, basis_set_list_b
680 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, basis_set_c
681 : TYPE(coulomb_integral_context_type) :: integral_context
682 : TYPE(neighbor_list_iterator_p_type), &
683 210 : DIMENSION(:), POINTER :: nl_iterator
684 210 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
685 210 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
686 :
687 210 : CALL timeset(routineN, handle)
688 :
689 : calculate_forces = PRESENT(force_a) .OR. PRESENT(force_b) .OR. PRESENT(force_c) .OR. &
690 210 : PRESENT(mat_dabc) .OR. PRESENT(mat_adbc) .OR. PRESENT(mat_abdc)
691 :
692 210 : my_eri_method = do_eri_mme
693 210 : IF (PRESENT(eri_method)) my_eri_method = eri_method
694 :
695 210 : IF (PRESENT(force_a) .OR. PRESENT(force_b) .OR. PRESENT(force_c)) THEN
696 26 : CPASSERT(PRESENT(pabc))
697 : END IF
698 :
699 210 : GG_count = 0; GR_count = 0; RR_count = 0
700 :
701 210 : nthread = 1
702 :
703 : ! get mapping between ERIs and atoms, sets, set offsets
704 210 : CALL get_eri_offsets(qs_env, basis_type_c, eri_offsets)
705 :
706 210 : atom_start = eri_offsets(first_c, 1)
707 210 : set_start = eri_offsets(first_c, 2)
708 210 : set_offset_start = eri_offsets(first_c, 3)
709 :
710 210 : atom_end = eri_offsets(last_c, 1)
711 210 : set_end = eri_offsets(last_c, 2)
712 210 : set_offset_end = eri_offsets(last_c, 3)
713 :
714 : ! get QS stuff
715 : CALL get_qs_env(qs_env, &
716 : atomic_kind_set=atomic_kind_set, &
717 : natom=natom, &
718 : qs_kind_set=qs_kind_set, &
719 : particle_set=particle_set, &
720 : cell=cell, &
721 210 : nkind=nkind)
722 :
723 210 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of, natom_of_kind=natom_of_kind)
724 :
725 210 : IF (PRESENT(force_a)) CALL mp2_eri_allocate_forces(force_a, natom_of_kind)
726 210 : IF (PRESENT(force_b)) CALL mp2_eri_allocate_forces(force_b, natom_of_kind)
727 210 : IF (PRESENT(force_c)) CALL mp2_eri_allocate_forces(force_c, natom_of_kind)
728 :
729 210 : nc = last_c - first_c + 1
730 :
731 : ! check for symmetry
732 210 : CPASSERT(SIZE(sab_nl) > 0)
733 210 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
734 :
735 210 : IF (do_symmetric) THEN
736 210 : CPASSERT(basis_type_a == basis_type_b)
737 : END IF
738 :
739 1512 : ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
740 210 : CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
741 210 : CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
742 :
743 210 : IF (my_eri_method == do_eri_os) THEN
744 74 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=max_am_a, basis_type=basis_type_a)
745 74 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=max_am_b, basis_type=basis_type_b)
746 74 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxlgto=max_am_c, basis_type=basis_type_c)
747 74 : CALL coulomb_integral_init(integral_context, max_am_3c=MAX(max_am_a, max_am_b, max_am_c))
748 : END IF
749 :
750 210 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
751 :
752 210 : mepos = 0
753 :
754 11966 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
755 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, inode=inode, &
756 11756 : iatom=iatom, jatom=jatom, r=rab)
757 :
758 : ! exclude periodic images because method is periodic intrinsically
759 11756 : IF (inode == 1) last_jatom = 0
760 :
761 11756 : IF (jatom /= last_jatom) THEN
762 984 : last_jatom = jatom
763 : ELSE
764 : CYCLE
765 : END IF
766 :
767 984 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
768 : ! When RI_AUX NONE is invoked, the pointers to basis_set_a and basis_set_b are created,
769 : ! but not filled. Therefore, we check for the association and the number of entries.
770 984 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
771 3618 : IF (SUM(basis_set_a%nsgf_set) <= 0) CYCLE
772 984 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
773 984 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
774 3618 : IF (SUM(basis_set_b%nsgf_set) <= 0) CYCLE
775 984 : atom_a = atom_of_kind(iatom)
776 984 : atom_b = atom_of_kind(jatom)
777 :
778 984 : first_sgfa => basis_set_a%first_sgf
779 984 : la_max => basis_set_a%lmax
780 984 : la_min => basis_set_a%lmin
781 984 : npgfa => basis_set_a%npgf
782 984 : nseta = basis_set_a%nset
783 984 : nsgfa => basis_set_a%nsgf_set
784 984 : rpgfa => basis_set_a%pgf_radius
785 984 : set_radius_a => basis_set_a%set_radius
786 984 : sphi_a => basis_set_a%sphi
787 984 : zeta => basis_set_a%zet
788 3618 : na = SUM(nsgfa)
789 :
790 984 : ra(:) = pbc(particle_set(iatom)%r, cell)
791 :
792 : ! basis jkind
793 984 : first_sgfb => basis_set_b%first_sgf
794 984 : lb_max => basis_set_b%lmax
795 984 : lb_min => basis_set_b%lmin
796 984 : npgfb => basis_set_b%npgf
797 984 : nsetb = basis_set_b%nset
798 984 : nsgfb => basis_set_b%nsgf_set
799 984 : rpgfb => basis_set_b%pgf_radius
800 984 : set_radius_b => basis_set_b%set_radius
801 984 : sphi_b => basis_set_b%sphi
802 984 : zetb => basis_set_b%zet
803 3618 : nb = SUM(nsgfb)
804 :
805 984 : rb(:) = pbc(particle_set(jatom)%r, cell)
806 :
807 984 : IF (do_symmetric) THEN
808 984 : IF (iatom <= jatom) THEN
809 656 : irow = iatom
810 656 : icol = jatom
811 : ELSE
812 328 : irow = jatom
813 328 : icol = iatom
814 : END IF
815 : ELSE
816 0 : irow = iatom
817 0 : icol = jatom
818 : END IF
819 :
820 4920 : ALLOCATE (habc(na, nb, nc))
821 984 : habc(:, :, :) = 0.0_dp ! needs to be initialized due to screening
822 984 : IF (PRESENT(mat_dabc)) THEN
823 0 : ALLOCATE (hdabc(3, na, nb, nc))
824 0 : hdabc(:, :, :, :) = 0.0_dp
825 : END IF
826 984 : IF (PRESENT(mat_adbc)) THEN
827 0 : ALLOCATE (hadbc(3, na, nb, nc))
828 0 : hadbc(:, :, :, :) = 0.0_dp
829 : END IF
830 984 : IF (PRESENT(mat_abdc)) THEN
831 0 : ALLOCATE (habdc(3, na, nb, nc))
832 0 : habdc(:, :, :, :) = 0.0_dp
833 : END IF
834 :
835 984 : IF (calculate_forces .AND. PRESENT(pabc)) THEN
836 540 : ALLOCATE (pabc_block(na, nb, nc))
837 5613 : DO ic = 1, nc
838 5478 : NULLIFY (pab_block)
839 : CALL dbcsr_get_block_p(matrix=pabc(ic)%matrix, &
840 5478 : row=irow, col=icol, block=pab_block, found=found)
841 5478 : CPASSERT(found)
842 11091 : IF (irow == iatom) THEN
843 3652 : to_be_asserted = SIZE(pab_block, 1) == SIZE(pabc_block, 1) .AND. SIZE(pab_block, 2) == SIZE(pabc_block, 2)
844 0 : CPASSERT(to_be_asserted)
845 188567 : pabc_block(:, :, ic) = pab_block(:, :)
846 : ELSE
847 1826 : to_be_asserted = SIZE(pab_block, 2) == SIZE(pabc_block, 1) .AND. SIZE(pab_block, 1) == SIZE(pabc_block, 2)
848 0 : CPASSERT(to_be_asserted)
849 72496 : pabc_block(:, :, ic) = TRANSPOSE(pab_block(:, :))
850 : END IF
851 : END DO
852 : END IF
853 :
854 3936 : rab(:) = pbc(rab, cell)
855 984 : dab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
856 :
857 984 : offset_a_end = 0
858 3618 : DO iset = 1, nseta
859 2634 : offset_a_start = offset_a_end
860 2634 : offset_a_end = offset_a_end + nsgfa(iset)
861 2634 : sgfa = first_sgfa(1, iset)
862 :
863 2634 : offset_b_end = 0
864 12450 : DO jset = 1, nsetb
865 8832 : offset_b_start = offset_b_end
866 8832 : offset_b_end = offset_b_end + nsgfb(jset)
867 :
868 8832 : sgfb = first_sgfb(1, jset)
869 :
870 : ! Screening
871 8832 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
872 :
873 : offset_c_end = 0
874 25552 : DO katom = atom_start, atom_end
875 :
876 14090 : atom_c = atom_of_kind(katom)
877 :
878 14090 : kkind = kind_of(katom)
879 14090 : CALL get_qs_kind(qs_kind=qs_kind_set(kkind), basis_set=basis_set_c, basis_type=basis_type_c)
880 14090 : first_sgfc => basis_set_c%first_sgf
881 14090 : lc_max => basis_set_c%lmax
882 14090 : lc_min => basis_set_c%lmin
883 14090 : nsetc = basis_set_c%nset
884 14090 : nsgfc => basis_set_c%nsgf_set
885 14090 : rpgfc => basis_set_c%pgf_radius
886 14090 : sphi_c => basis_set_c%sphi
887 14090 : zetc => basis_set_c%zet
888 14090 : npgfc => basis_set_c%npgf
889 :
890 14090 : rc(:) = pbc(particle_set(katom)%r, cell)
891 :
892 14090 : kset_start = 1; kset_end = nsetc
893 14090 : IF (katom == atom_start) kset_start = set_start
894 14090 : IF (katom == atom_end) kset_end = set_end
895 :
896 149818 : DO kset = kset_start, kset_end
897 126896 : first_set = 1; last_set = nsgfc(kset)
898 126896 : IF (kset == kset_start .AND. katom == atom_start) first_set = set_offset_start
899 126896 : IF (kset == kset_end .AND. katom == atom_end) last_set = set_offset_end
900 :
901 126896 : offset_c_start = offset_c_end
902 126896 : offset_c_end = offset_c_end + last_set + 1 - first_set
903 126896 : sgfc = first_sgfc(1, kset)
904 :
905 : #!some fypp magic to deal with combinations of optional arguments
906 : #:for pabc_present in [0, 1]
907 : #:for doforce_1 in [0, 1]
908 : #:for doforce_2 in [0, 1]
909 : #:for doforce_3 in [0, 1]
910 : #:for dabc in [0, 1]
911 : #:for adbc in [0, 1]
912 : #:for abdc in [0, 1]
913 : IF (${conditional(doforce_1)}$PRESENT(force_a) .AND. &
914 : ${conditional(doforce_2)}$PRESENT(force_b) .AND. &
915 : ${conditional(doforce_3)}$PRESENT(force_c) .AND. &
916 : ${conditional(pabc_present)}$PRESENT(pabc) .AND. &
917 : ${conditional(dabc)}$PRESENT(mat_dabc) .AND. &
918 253792 : ${conditional(adbc)}$PRESENT(mat_adbc) .AND. &
919 14090 : ${conditional(abdc)}$PRESENT(mat_abdc)) THEN
920 : CALL integrate_set_3c( &
921 : param%par, potential_parameter, &
922 : la_min(iset), la_max(iset), &
923 : lb_min(jset), lb_max(jset), &
924 : lc_min(kset), lc_max(kset), &
925 : npgfa(iset), npgfb(jset), npgfc(kset), &
926 : zeta(:, iset), zetb(:, jset), zetc(:, kset), &
927 : rpgfa(:, iset), rpgfb(:, jset), rpgfc(:, kset), &
928 : ra, rb, rc, &
929 : habc, &
930 : nsgfa(iset), nsgfb(jset), last_set - first_set + 1, &
931 : offset_a_start, offset_b_start, offset_c_start, &
932 : 0, 0, first_set - 1, &
933 : sphi_a, sphi_b, sphi_c, &
934 : sgfa, sgfb, sgfc, &
935 : nsgfa(iset), nsgfb(jset), nsgfc(kset), &
936 : my_eri_method, &
937 : coulomb_context=integral_context, &
938 : $: 'pabc=pabc_block, &'*pabc_present
939 : $: 'force_a=force_a(ikind)%forces(:, atom_a), &'*doforce_1
940 : $: 'force_b=force_b(jkind)%forces(:, atom_b), &'*doforce_2
941 : $: 'force_c=force_c(kkind)%forces(:, atom_c), &'*doforce_3
942 : do_symmetric=do_symmetric, &
943 : on_diagonal=iatom == jatom, &
944 : $: 'hdabc=hdabc, &'*dabc
945 : $: 'hadbc=hadbc, &'*adbc
946 : $: 'habdc=habdc, &'*abdc
947 126896 : GG_count=GG_count, GR_count=GR_count, RR_count=RR_count)
948 : END IF
949 : #:endfor
950 : #:endfor
951 : #:endfor
952 : #:endfor
953 : #:endfor
954 : #:endfor
955 : #:endfor
956 : END DO
957 : END DO
958 : END DO
959 : END DO
960 :
961 984 : IF (calculate_forces .AND. PRESENT(pabc)) DEALLOCATE (pabc_block)
962 36702 : DO ic = 1, nc
963 35718 : NULLIFY (munu_block)
964 : CALL dbcsr_get_block_p(matrix=mat_ab(ic)%matrix, &
965 35718 : row=irow, col=icol, block=munu_block, found=found)
966 35718 : CPASSERT(found)
967 2838622 : munu_block(:, :) = 0.0_dp
968 72420 : IF (irow == iatom) THEN
969 23812 : to_be_asserted = SIZE(munu_block, 1) == SIZE(habc, 1) .AND. SIZE(munu_block, 2) == SIZE(habc, 2)
970 0 : CPASSERT(to_be_asserted)
971 2035637 : munu_block(:, :) = habc(:, :, ic)
972 : ELSE
973 11906 : to_be_asserted = SIZE(munu_block, 2) == SIZE(habc, 1) .AND. SIZE(munu_block, 1) == SIZE(habc, 2)
974 0 : CPASSERT(to_be_asserted)
975 802985 : munu_block(:, :) = TRANSPOSE(habc(:, :, ic))
976 : END IF
977 : END DO
978 984 : DEALLOCATE (habc)
979 1194 : IF (calculate_forces) THEN
980 5613 : DO ic = 1, nc
981 22047 : DO i_xyz = 1, 3
982 16434 : IF (PRESENT(mat_dabc)) THEN
983 0 : NULLIFY (munu_block)
984 : CALL dbcsr_get_block_p(matrix=mat_dabc(i_xyz, ic)%matrix, &
985 0 : row=irow, col=icol, block=munu_block, found=found)
986 0 : CPASSERT(found)
987 0 : munu_block(:, :) = 0.0_dp
988 0 : IF (irow == iatom) THEN
989 0 : munu_block(:, :) = hdabc(i_xyz, :, :, ic)
990 : ELSE
991 0 : munu_block(:, :) = TRANSPOSE(hdabc(i_xyz, :, :, ic))
992 : END IF
993 : END IF
994 16434 : IF (PRESENT(mat_adbc)) THEN
995 0 : NULLIFY (munu_block)
996 : CALL dbcsr_get_block_p(matrix=mat_adbc(i_xyz, ic)%matrix, &
997 0 : row=irow, col=icol, block=munu_block, found=found)
998 0 : CPASSERT(found)
999 0 : munu_block(:, :) = 0.0_dp
1000 0 : IF (irow == iatom) THEN
1001 0 : munu_block(:, :) = hadbc(i_xyz, :, :, ic)
1002 : ELSE
1003 0 : munu_block(:, :) = TRANSPOSE(hadbc(i_xyz, :, :, ic))
1004 : END IF
1005 : END IF
1006 21912 : IF (PRESENT(mat_abdc)) THEN
1007 0 : NULLIFY (munu_block)
1008 : CALL dbcsr_get_block_p(matrix=mat_abdc(i_xyz, ic)%matrix, &
1009 0 : row=irow, col=icol, block=munu_block, found=found)
1010 0 : CPASSERT(found)
1011 0 : munu_block(:, :) = 0.0_dp
1012 0 : IF (irow == iatom) THEN
1013 0 : munu_block(:, :) = habdc(i_xyz, :, :, ic)
1014 : ELSE
1015 0 : munu_block(:, :) = TRANSPOSE(habdc(i_xyz, :, :, ic))
1016 : END IF
1017 : END IF
1018 : END DO
1019 : END DO
1020 135 : IF (PRESENT(mat_dabc)) DEALLOCATE (hdabc)
1021 135 : IF (PRESENT(mat_adbc)) DEALLOCATE (hadbc)
1022 135 : IF (PRESENT(mat_abdc)) DEALLOCATE (habdc)
1023 : END IF
1024 : END DO
1025 :
1026 210 : DEALLOCATE (basis_set_list_a, basis_set_list_b)
1027 210 : CALL neighbor_list_iterator_release(nl_iterator)
1028 :
1029 210 : IF (my_eri_method == do_eri_os) CALL coulomb_integral_cleanup(integral_context)
1030 :
1031 210 : CALL cp_eri_mme_update_local_counts(param, para_env, GG_count_3c=GG_count, GR_count_3c=GR_count, RR_count_3c=RR_count)
1032 :
1033 210 : CALL timestop(handle)
1034 630 : END SUBROUTINE mp2_eri_3c_integrate
1035 :
1036 : ! **************************************************************************************************
1037 : !> \brief Integrate set triple and contract with sphi matrix
1038 : !> \param param ...
1039 : !> \param potential_parameter ...
1040 : !> \param la_min ...
1041 : !> \param la_max ...
1042 : !> \param lb_min ...
1043 : !> \param lb_max ...
1044 : !> \param lc_min ...
1045 : !> \param lc_max ...
1046 : !> \param npgfa ...
1047 : !> \param npgfb ...
1048 : !> \param npgfc ...
1049 : !> \param zeta ...
1050 : !> \param zetb ...
1051 : !> \param zetc ...
1052 : !> \param ra ...
1053 : !> \param rb ...
1054 : !> \param rc ...
1055 : !> \param habc ...
1056 : !> \param n_habc_a ...
1057 : !> \param n_habc_b ...
1058 : !> \param n_habc_c ...
1059 : !> \param offset_habc_a ...
1060 : !> \param offset_habc_b ...
1061 : !> \param offset_habc_c ...
1062 : !> \param offset_set_a ...
1063 : !> \param offset_set_b ...
1064 : !> \param offset_set_c ...
1065 : !> \param sphi_a ...
1066 : !> \param sphi_b ...
1067 : !> \param sphi_c ...
1068 : !> \param sgfa ...
1069 : !> \param sgfb ...
1070 : !> \param sgfc ...
1071 : !> \param nsgfa ...
1072 : !> \param nsgfb ...
1073 : !> \param nsgfc ...
1074 : !> \param eri_method ...
1075 : !> \param pabc ...
1076 : !> \param force_a ...
1077 : !> \param force_b ...
1078 : !> \param force_c ...
1079 : !> \param do_symmetric ...
1080 : !> \param on_diagonal ...
1081 : !> \param hdabc ...
1082 : !> \param hadbc ...
1083 : !> \param habdc ...
1084 : !> \param GG_count ...
1085 : !> \param GR_count ...
1086 : !> \param RR_count ...
1087 : !> \param coulomb_context context for the selected Coulomb integral library
1088 : !> \note
1089 : ! **************************************************************************************************
1090 126896 : SUBROUTINE integrate_set_3c(param, potential_parameter, &
1091 : la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
1092 : npgfa, npgfb, npgfc, &
1093 126896 : zeta, zetb, zetc, &
1094 126896 : rpgfa, rpgfb, rpgfc, &
1095 : ra, rb, rc, &
1096 253792 : habc, &
1097 : n_habc_a, n_habc_b, n_habc_c, &
1098 : offset_habc_a, offset_habc_b, offset_habc_c, &
1099 : offset_set_a, offset_set_b, offset_set_c, &
1100 126896 : sphi_a, sphi_b, sphi_c, &
1101 : sgfa, sgfb, sgfc, &
1102 : nsgfa, nsgfb, nsgfc, &
1103 : eri_method, coulomb_context, &
1104 126896 : pabc, &
1105 : force_a, force_b, force_c, &
1106 : do_symmetric, on_diagonal, &
1107 126896 : hdabc, hadbc, habdc, &
1108 : GG_count, GR_count, RR_count)
1109 :
1110 : TYPE(eri_mme_param), INTENT(INOUT) :: param
1111 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1112 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
1113 : lc_max, npgfa, npgfb, npgfc
1114 : REAL(KIND=dp), DIMENSION(npgfa), INTENT(IN) :: zeta
1115 : REAL(KIND=dp), DIMENSION(npgfb), INTENT(IN) :: zetb
1116 : REAL(KIND=dp), DIMENSION(npgfc), INTENT(IN) :: zetc
1117 : REAL(KIND=dp), DIMENSION(npgfa), INTENT(IN) :: rpgfa
1118 : REAL(KIND=dp), DIMENSION(npgfb), INTENT(IN) :: rpgfb
1119 : REAL(KIND=dp), DIMENSION(npgfc), INTENT(IN) :: rpgfc
1120 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb, rc
1121 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
1122 : INTEGER, INTENT(IN) :: n_habc_a, n_habc_b, n_habc_c, offset_habc_a, offset_habc_b, &
1123 : offset_habc_c, offset_set_a, offset_set_b, offset_set_c
1124 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sphi_a, sphi_b, sphi_c
1125 : INTEGER, INTENT(IN) :: sgfa, sgfb, sgfc, nsgfa, nsgfb, nsgfc, &
1126 : eri_method
1127 : TYPE(coulomb_integral_context_type), INTENT(INOUT), OPTIONAL :: coulomb_context
1128 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
1129 : OPTIONAL :: pabc
1130 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
1131 : OPTIONAL :: force_a, force_b, force_c
1132 : LOGICAL, INTENT(IN) :: do_symmetric
1133 : LOGICAL, INTENT(IN), OPTIONAL :: on_diagonal
1134 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
1135 : INTENT(OUT), OPTIONAL :: hdabc, hadbc, habdc
1136 : INTEGER, INTENT(INOUT), OPTIONAL :: GG_count, GR_count, RR_count
1137 :
1138 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_set_3c'
1139 :
1140 : INTEGER :: ax, ay, az, bx, by, bz, cx, cy, cz, habc_a_end, habc_a_start, habc_b_end, &
1141 : habc_b_start, habc_c_end, habc_c_start, handle, i_xyz, ico, icoc, icox, icoy, icoz, ipgf, &
1142 : jco, jcox, jcoy, jcoz, jpgf, kco, kcox, kcoy, kcoz, kpgf, la, la_max_d, lb, &
1143 : lb_max_d, lc, lc_max_d, na, nb, nc, ncoa, ncoa_d, ncob, ncob_d, ncoc, ncoc_d, &
1144 : set_a_end, set_a_start, set_b_end, set_b_start, &
1145 : set_c_end, set_c_start, sphi_a_start, sphi_b_start, sphi_c_start
1146 : INTEGER, DIMENSION(3) :: la_xyz, lb_xyz
1147 : LOGICAL :: calculate_forces, do_force_a, do_force_b, &
1148 : do_force_c
1149 : REAL(KIND=dp) :: w
1150 126896 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab_hh, pab_hs
1151 126896 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: habc_contr, habc_uncontr, &
1152 126896 : habc_uncontr_d, pabc_hhh, &
1153 126896 : pabc_hsh, pabc_hss, pabc_sss
1154 126896 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: habdc_contr, habdc_uncontr, hadbc_contr, &
1155 126896 : hadbc_uncontr, hdabc_contr, hdabc_uncontr
1156 :
1157 126896 : CALL timeset(routineN, handle)
1158 :
1159 126896 : do_force_a = PRESENT(force_a) .OR. PRESENT(hdabc)
1160 126896 : do_force_b = PRESENT(force_b) .OR. PRESENT(hadbc)
1161 126896 : do_force_c = PRESENT(force_c) .OR. PRESENT(habdc)
1162 126896 : calculate_forces = do_force_a .OR. do_force_b .OR. do_force_c
1163 :
1164 126896 : IF (do_symmetric) THEN
1165 126896 : CPASSERT(PRESENT(on_diagonal))
1166 : END IF
1167 :
1168 126896 : la_max_d = la_max
1169 126896 : lb_max_d = lb_max
1170 126896 : lc_max_d = lc_max
1171 :
1172 126896 : IF (calculate_forces) THEN
1173 10386 : IF (do_force_a) la_max_d = la_max + 1
1174 10386 : IF (do_force_b) lb_max_d = lb_max + 1
1175 10386 : IF (do_force_c) lc_max_d = lc_max + 1
1176 : END IF
1177 :
1178 126896 : ncoa = npgfa*ncoset(la_max)
1179 126896 : ncob = npgfb*ncoset(lb_max)
1180 126896 : ncoc = npgfc*ncoset(lc_max)
1181 :
1182 126896 : ncoa_d = npgfa*ncoset(la_max_d)
1183 126896 : ncob_d = npgfb*ncoset(lb_max_d)
1184 126896 : ncoc_d = npgfc*ncoset(lc_max_d)
1185 :
1186 634480 : ALLOCATE (habc_uncontr_d(ncoset(la_max_d), ncoset(lb_max_d), ncoset(lc_max_d)))
1187 126896 : habc_uncontr_d(:, :, :) = 0.0_dp
1188 634480 : ALLOCATE (habc_uncontr(ncoa, ncob, ncoc)); habc_uncontr(:, :, :) = 0.0_dp
1189 126896 : IF (PRESENT(hdabc)) THEN
1190 0 : ALLOCATE (hdabc_uncontr(3, ncoa, ncob, ncoc)); hdabc_uncontr(:, :, :, :) = 0.0_dp
1191 : END IF
1192 126896 : IF (PRESENT(hadbc)) THEN
1193 0 : ALLOCATE (hadbc_uncontr(3, ncoa, ncob, ncoc)); hadbc_uncontr(:, :, :, :) = 0.0_dp
1194 : END IF
1195 126896 : IF (PRESENT(habdc)) THEN
1196 0 : ALLOCATE (habdc_uncontr(3, ncoa, ncob, ncoc)); habdc_uncontr(:, :, :, :) = 0.0_dp
1197 : END IF
1198 :
1199 126896 : habc_a_start = offset_habc_a + 1; habc_a_end = offset_habc_a + n_habc_a
1200 126896 : habc_b_start = offset_habc_b + 1; habc_b_end = offset_habc_b + n_habc_b
1201 126896 : habc_c_start = offset_habc_c + 1; habc_c_end = offset_habc_c + n_habc_c
1202 126896 : set_a_start = offset_set_a + 1; set_a_end = offset_set_a + n_habc_a
1203 126896 : set_b_start = offset_set_b + 1; set_b_end = offset_set_b + n_habc_b
1204 126896 : set_c_start = offset_set_c + 1; set_c_end = offset_set_c + n_habc_c
1205 :
1206 126896 : IF (eri_method == do_eri_mme) THEN
1207 36342 : CALL eri_mme_set_potential(param, convert_potential_type(potential_parameter%potential_type), potential_parameter%omega)
1208 :
1209 36342 : IF (calculate_forces .AND. PRESENT(pabc)) THEN
1210 : ! uncontracted hermite-gaussian representation of density matrix
1211 10386 : sphi_a_start = sgfa - 1 + set_a_start
1212 10386 : sphi_b_start = sgfb - 1 + set_b_start
1213 10386 : sphi_c_start = sgfc - 1 + set_c_start
1214 :
1215 51930 : ALLOCATE (pabc_sss(n_habc_a, n_habc_b, n_habc_c))
1216 342953 : pabc_sss(:, :, :) = pabc(habc_a_start:habc_a_end, habc_b_start:habc_b_end, habc_c_start:habc_c_end)
1217 51930 : ALLOCATE (pabc_hss(ncoa, n_habc_b, n_habc_c))
1218 51930 : ALLOCATE (pabc_hsh(ncoa, n_habc_b, ncoc))
1219 41544 : ALLOCATE (pabc_hhh(ncoa, ncob, ncoc))
1220 41544 : ALLOCATE (pab_hs(ncoa, n_habc_b))
1221 41544 : ALLOCATE (pab_hh(ncoa, ncob))
1222 :
1223 : CALL dgemm("N", "N", ncoa, n_habc_b*n_habc_c, n_habc_a, 1.0_dp, &
1224 10386 : sphi_a(:, sphi_a_start), SIZE(sphi_a, 1), pabc_sss, n_habc_a, 0.0_dp, pabc_hss, ncoa)
1225 : CALL dgemm("N", "T", ncoa*n_habc_b, ncoc, n_habc_c, 1.0_dp, &
1226 10386 : pabc_hss, ncoa*n_habc_b, sphi_c(:, sphi_c_start), SIZE(sphi_c, 1), 0.0_dp, pabc_hsh, ncoa*n_habc_b)
1227 :
1228 68148 : DO icoc = 1, ncoc
1229 1108764 : pab_hs(:, :) = pabc_hsh(:, :, icoc)
1230 : CALL dgemm("N", "T", ncoa, ncob, n_habc_b, 1.0_dp, &
1231 57762 : pab_hs, ncoa, sphi_b(:, sphi_b_start), SIZE(sphi_b, 1), 0.0_dp, pab_hh, ncoa)
1232 2262300 : pabc_hhh(:, :, icoc) = pab_hh(:, :)
1233 : END DO
1234 : END IF
1235 :
1236 105732 : DO ipgf = 1, npgfa
1237 69390 : na = (ipgf - 1)*ncoset(la_max)
1238 254142 : DO jpgf = 1, npgfb
1239 148410 : nb = (jpgf - 1)*ncoset(lb_max)
1240 366210 : DO kpgf = 1, npgfc
1241 148410 : nc = (kpgf - 1)*ncoset(lc_max)
1242 148410 : habc_uncontr_d(:, :, :) = 0.0_dp
1243 : CALL eri_mme_3c_integrate(param, &
1244 : la_min, la_max_d, lb_min, lb_max_d, lc_min, lc_max_d, &
1245 : zeta(ipgf), zetb(jpgf), zetc(kpgf), ra, rb, rc, habc_uncontr_d, 0, 0, 0, &
1246 148410 : GG_count, GR_count, RR_count)
1247 :
1248 : habc_uncontr(na + 1:na + ncoset(la_max), nb + 1:nb + ncoset(lb_max), nc + 1:nc + ncoset(lc_max)) = &
1249 8154180 : habc_uncontr_d(:ncoset(la_max), :ncoset(lb_max), :ncoset(lc_max))
1250 :
1251 296820 : IF (calculate_forces) THEN
1252 82500 : DO lc = lc_min, lc_max
1253 164700 : DO cx = 0, lc
1254 266850 : DO cy = 0, lc - cx
1255 143400 : cz = lc - cx - cy
1256 143400 : kco = coset(cx, cy, cz)
1257 143400 : kcox = coset(cx + 1, cy, cz)
1258 143400 : kcoy = coset(cx, cy + 1, cz)
1259 143400 : kcoz = coset(cx, cy, cz + 1)
1260 414600 : DO lb = lb_min, lb_max
1261 592800 : DO bx = 0, lb
1262 788400 : DO by = 0, lb - bx
1263 339000 : bz = lb - bx - by
1264 339000 : jco = coset(bx, by, bz)
1265 339000 : jcox = coset(bx + 1, by, bz)
1266 339000 : jcoy = coset(bx, by + 1, bz)
1267 339000 : jcoz = coset(bx, by, bz + 1)
1268 1075620 : DO la = la_min, la_max
1269 1500105 : DO ax = 0, la
1270 2078460 : DO ay = 0, la - ax
1271 917355 : az = la - ax - ay
1272 3669420 : la_xyz = [ax, ay, az]
1273 3669420 : lb_xyz = [bx, by, bz]
1274 917355 : ico = coset(ax, ay, az)
1275 917355 : icox = coset(ax + 1, ay, az)
1276 917355 : icoy = coset(ax, ay + 1, az)
1277 917355 : icoz = coset(ax, ay, az + 1)
1278 :
1279 917355 : w = 1.0_dp
1280 917355 : IF (do_symmetric .AND. .NOT. on_diagonal) w = 2.0_dp
1281 :
1282 917355 : IF (PRESENT(force_a)) THEN
1283 : force_a = force_a + 2.0_dp*w*zeta(ipgf)* &
1284 : [pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(icox, jco, kco), &
1285 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(icoy, jco, kco), &
1286 3669420 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(icoz, jco, kco)]
1287 :
1288 : END IF
1289 917355 : IF (PRESENT(force_b)) THEN
1290 : force_b = force_b + 2.0_dp*w*zetb(jpgf)* &
1291 : [pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jcox, kco), &
1292 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jcoy, kco), &
1293 3669420 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jcoz, kco)]
1294 : END IF
1295 917355 : IF (PRESENT(force_c)) THEN
1296 : force_c = force_c + 2.0_dp*w*zetc(kpgf)* &
1297 : [pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jco, kcox), &
1298 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jco, kcoy), &
1299 3669420 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jco, kcoz)]
1300 : END IF
1301 :
1302 917355 : IF (PRESENT(hdabc)) THEN
1303 : hdabc_uncontr(1:3, na + ico, nb + jco, nc + kco) = 2.0_dp*zeta(ipgf)* &
1304 : [habc_uncontr_d(icox, jco, kco), &
1305 : habc_uncontr_d(icoy, jco, kco), &
1306 0 : habc_uncontr_d(icoz, jco, kco)]
1307 : END IF
1308 917355 : IF (PRESENT(hadbc)) THEN
1309 : hadbc_uncontr(1:3, na + ico, nb + jco, nc + kco) = 2.0_dp*zetb(jpgf)* &
1310 : [habc_uncontr_d(ico, jcox, kco), &
1311 : habc_uncontr_d(ico, jcoy, kco), &
1312 0 : habc_uncontr_d(ico, jcoz, kco)]
1313 : END IF
1314 1602240 : IF (PRESENT(habdc)) THEN
1315 : habdc_uncontr(1:3, na + ico, nb + jco, nc + kco) = 2.0_dp*zetc(kpgf)* &
1316 : [habc_uncontr_d(ico, jco, kcox), &
1317 : habc_uncontr_d(ico, jco, kcoy), &
1318 0 : habc_uncontr_d(ico, jco, kcoz)]
1319 : END IF
1320 : END DO
1321 : END DO
1322 : END DO
1323 : END DO
1324 : END DO
1325 : END DO
1326 : END DO
1327 : END DO
1328 : END DO
1329 : END IF
1330 :
1331 : END DO
1332 : END DO
1333 : END DO
1334 :
1335 90554 : ELSE IF (eri_method == do_eri_os) THEN
1336 :
1337 90554 : IF (calculate_forces) CPABORT("NYI")
1338 :
1339 90554 : CPASSERT(PRESENT(coulomb_context))
1340 : CALL compute_coulomb_3c(coulomb_context, la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
1341 : npgfa, npgfb, npgfc, zeta, zetb, zetc, rpgfa, rpgfb, rpgfc, &
1342 90554 : ra, rb, rc, habc_uncontr, potential_parameter)
1343 :
1344 0 : ELSE IF (eri_method == do_eri_gpw) THEN
1345 :
1346 0 : CPABORT("GPW not enabled in the ERI interface.")
1347 :
1348 : END IF
1349 :
1350 634480 : ALLOCATE (habc_contr(nsgfa, nsgfb, nsgfc))
1351 126896 : IF (PRESENT(hdabc)) THEN
1352 0 : ALLOCATE (hdabc_contr(3, nsgfa, nsgfb, nsgfc))
1353 : END IF
1354 126896 : IF (PRESENT(hadbc)) THEN
1355 0 : ALLOCATE (hadbc_contr(3, nsgfa, nsgfb, nsgfc))
1356 : END IF
1357 126896 : IF (PRESENT(habdc)) THEN
1358 0 : ALLOCATE (habdc_contr(3, nsgfa, nsgfb, nsgfc))
1359 : END IF
1360 :
1361 : CALL abc_contract(habc_contr, habc_uncontr, &
1362 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
1363 126896 : ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
1364 :
1365 126896 : IF (calculate_forces) THEN
1366 41544 : DO i_xyz = 1, 3
1367 31158 : IF (PRESENT(hdabc)) THEN
1368 : CALL abc_contract(hdabc_contr(i_xyz, :, :, :), hdabc_uncontr(i_xyz, :, :, :), &
1369 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
1370 0 : ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
1371 : END IF
1372 31158 : IF (PRESENT(hadbc)) THEN
1373 : CALL abc_contract(hadbc_contr(i_xyz, :, :, :), hadbc_uncontr(i_xyz, :, :, :), &
1374 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
1375 0 : ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
1376 : END IF
1377 41544 : IF (PRESENT(habdc)) THEN
1378 : CALL abc_contract(habdc_contr(i_xyz, :, :, :), habdc_uncontr(i_xyz, :, :, :), &
1379 : sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
1380 0 : ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
1381 : END IF
1382 : END DO
1383 : END IF
1384 :
1385 : habc(habc_a_start:habc_a_end, habc_b_start:habc_b_end, habc_c_start:habc_c_end) = &
1386 4036258 : habc_contr(set_a_start:set_a_end, set_b_start:set_b_end, set_c_start:set_c_end)
1387 :
1388 126896 : IF (calculate_forces) THEN
1389 10386 : IF (PRESENT(hdabc)) hdabc(:, habc_a_start:habc_a_end, habc_b_start:habc_b_end, habc_c_start:habc_c_end) = &
1390 0 : hdabc_contr(:, set_a_start:set_a_end, set_b_start:set_b_end, set_c_start:set_c_end)
1391 10386 : IF (PRESENT(hadbc)) hadbc(:, habc_a_start:habc_a_end, habc_b_start:habc_b_end, habc_c_start:habc_c_end) = &
1392 0 : hadbc_contr(:, set_a_start:set_a_end, set_b_start:set_b_end, set_c_start:set_c_end)
1393 10386 : IF (PRESENT(habdc)) habdc(:, habc_a_start:habc_a_end, habc_b_start:habc_b_end, habc_c_start:habc_c_end) = &
1394 0 : habdc_contr(:, set_a_start:set_a_end, set_b_start:set_b_end, set_c_start:set_c_end)
1395 : END IF
1396 :
1397 126896 : CALL timestop(handle)
1398 :
1399 253792 : END SUBROUTINE integrate_set_3c
1400 :
1401 : ! **************************************************************************************************
1402 : !> \brief get pointer to atom, pointer to set and offset in a set for each spherical orbital of a
1403 : !> basis.
1404 : !> \param qs_env ...
1405 : !> \param basis_type ...
1406 : !> \param eri_offsets (:,1) atom numbers
1407 : !> (:,2) set numbers
1408 : !> (:,3) set offsets
1409 : ! **************************************************************************************************
1410 552 : SUBROUTINE get_eri_offsets(qs_env, basis_type, eri_offsets)
1411 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1412 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: basis_type
1413 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: eri_offsets
1414 :
1415 : INTEGER :: dimen_basis, iatom, ikind, iset, isgf, &
1416 : natom, nkind, nset, nsgf, offset, &
1417 : set_offset
1418 552 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
1419 552 : INTEGER, DIMENSION(:), POINTER :: nsgf_set
1420 552 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1421 : TYPE(gto_basis_set_type), POINTER :: basis_set
1422 552 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1423 552 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1424 :
1425 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1426 552 : particle_set=particle_set, natom=natom, nkind=nkind)
1427 :
1428 552 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
1429 :
1430 552 : dimen_basis = 0
1431 2056 : DO iatom = 1, natom
1432 1504 : ikind = kind_of(iatom)
1433 1504 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type=basis_type)
1434 2056 : dimen_basis = dimen_basis + nsgf
1435 : END DO
1436 :
1437 1656 : ALLOCATE (eri_offsets(dimen_basis, 3))
1438 :
1439 552 : offset = 0
1440 2056 : DO iatom = 1, natom
1441 1504 : ikind = kind_of(iatom)
1442 1504 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
1443 1504 : nset = basis_set%nset
1444 1504 : nsgf_set => basis_set%nsgf_set
1445 15756 : DO iset = 1, nset
1446 13700 : set_offset = 0
1447 52180 : DO isgf = 1, nsgf_set(iset)
1448 38480 : set_offset = set_offset + 1
1449 167620 : eri_offsets(offset + set_offset, :) = [iatom, iset, set_offset]
1450 : END DO
1451 15204 : offset = offset + nsgf_set(iset)
1452 : END DO
1453 : END DO
1454 1104 : END SUBROUTINE get_eri_offsets
1455 :
1456 : ! **************************************************************************************************
1457 : !> \brief ...
1458 : !> \param force ...
1459 : !> \param natom_of_kind ...
1460 : ! **************************************************************************************************
1461 104 : PURE SUBROUTINE mp2_eri_allocate_forces(force, natom_of_kind)
1462 : TYPE(mp2_eri_force), ALLOCATABLE, &
1463 : DIMENSION(:), INTENT(OUT) :: force
1464 : INTEGER, DIMENSION(:), INTENT(IN) :: natom_of_kind
1465 :
1466 : INTEGER :: ikind, n, nkind
1467 :
1468 104 : nkind = SIZE(natom_of_kind)
1469 :
1470 496 : ALLOCATE (force(nkind))
1471 :
1472 288 : DO ikind = 1, nkind
1473 184 : n = natom_of_kind(ikind)
1474 552 : ALLOCATE (force(ikind)%forces(3, n))
1475 1440 : force(ikind)%forces(:, :) = 0.0_dp
1476 : END DO
1477 104 : END SUBROUTINE mp2_eri_allocate_forces
1478 :
1479 : ! **************************************************************************************************
1480 : !> \brief ...
1481 : !> \param force ...
1482 : ! **************************************************************************************************
1483 104 : PURE SUBROUTINE mp2_eri_deallocate_forces(force)
1484 : TYPE(mp2_eri_force), ALLOCATABLE, &
1485 : DIMENSION(:), INTENT(INOUT) :: force
1486 :
1487 : INTEGER :: ikind, nkind
1488 :
1489 104 : IF (ALLOCATED(force)) THEN
1490 104 : nkind = SIZE(force)
1491 288 : DO ikind = 1, nkind
1492 288 : IF (ALLOCATED(force(ikind)%forces)) DEALLOCATE (force(ikind)%forces)
1493 : END DO
1494 :
1495 288 : DEALLOCATE (force)
1496 : END IF
1497 104 : END SUBROUTINE mp2_eri_deallocate_forces
1498 :
1499 85324 : FUNCTION convert_potential_type(potential_type) RESULT(res)
1500 : INTEGER, INTENT(IN) :: potential_type
1501 : INTEGER :: res
1502 :
1503 85324 : IF (potential_type == do_potential_coulomb) THEN
1504 : res = eri_mme_coulomb
1505 12541 : ELSE IF (potential_type == do_potential_long) THEN
1506 : res = eri_mme_longrange
1507 : ELSE
1508 0 : CPABORT("MME potential not implemented!")
1509 : END IF
1510 :
1511 85324 : END FUNCTION convert_potential_type
1512 :
1513 0 : END MODULE mp2_eri
|
