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