Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 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 286 : SUBROUTINE mp2_eri_2c_integrate(param, potential_parameter, para_env, qs_env, basis_type_a, basis_type_b, hab, first_b, &
116 286 : 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 286 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
143 286 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eri_offsets
144 286 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
145 286 : npgfb, nsgfa, nsgfb
146 286 : 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 286 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_a, sphi_b, zeta, zetb
152 286 : 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 286 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
156 286 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
157 :
158 286 : CALL timeset(routineN, handle)
159 :
160 286 : my_eri_method = do_eri_mme
161 286 : IF (PRESENT(eri_method)) my_eri_method = eri_method
162 :
163 286 : my_do_reflection_a = .FALSE.
164 286 : IF (PRESENT(do_reflection_a) .AND. PRESENT(reflection_z_a)) my_do_reflection_a = do_reflection_a
165 :
166 286 : my_do_reflection_b = .FALSE.
167 286 : IF (PRESENT(do_reflection_b) .AND. PRESENT(reflection_z_b)) my_do_reflection_b = do_reflection_b
168 :
169 286 : G_count = 0; R_count = 0;
170 : ! get mapping between ERIs and atoms, sets, set offsets
171 286 : CALL get_eri_offsets(qs_env, basis_type_b, eri_offsets)
172 :
173 286 : atom_start = eri_offsets(first_b, 1)
174 286 : set_start = eri_offsets(first_b, 2)
175 286 : set_offset_start = eri_offsets(first_b, 3)
176 :
177 286 : atom_end = eri_offsets(last_b, 1)
178 286 : set_end = eri_offsets(last_b, 2)
179 286 : 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 286 : cell=cell, particle_set=particle_set, natom=natom, nkind=nkind)
184 :
185 286 : 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 286 : IF (PRESENT(force_a)) CALL mp2_eri_allocate_forces(force_a, natom_of_kind)
188 286 : 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 286 : nseta_total = 0
192 1088 : DO iatom = 1, natom
193 802 : ikind = kind_of(iatom)
194 802 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type_a)
195 1088 : nseta_total = nseta_total + basis_set_a%nset
196 : END DO
197 :
198 : nsetb_total = 0
199 796 : DO jatom = atom_start, atom_end
200 510 : jkind = kind_of(jatom)
201 510 : CALL get_qs_kind(qs_kind=qs_kind_set(jkind), basis_set=basis_set_b, basis_type=basis_type_b)
202 796 : nsetb_total = nsetb_total + basis_set_b%nset
203 : END DO
204 :
205 : n_setpair = nseta_total*nsetb_total
206 :
207 286 : my_setpair = 0
208 :
209 286 : offset_a_end = 0
210 1088 : DO iatom = 1, natom
211 802 : ikind = kind_of(iatom)
212 802 : atom_a = atom_of_kind(iatom)
213 802 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type=basis_type_a)
214 :
215 802 : first_sgfa => basis_set_a%first_sgf
216 802 : la_max => basis_set_a%lmax
217 802 : la_min => basis_set_a%lmin
218 802 : nseta = basis_set_a%nset
219 802 : nsgfa => basis_set_a%nsgf_set
220 802 : sphi_a => basis_set_a%sphi
221 802 : zeta => basis_set_a%zet
222 802 : npgfa => basis_set_a%npgf
223 :
224 802 : ra(:) = pbc(particle_set(iatom)%r, cell)
225 :
226 802 : IF (my_do_reflection_a) THEN
227 0 : ra(3) = 2.0_dp*reflection_z_a - ra(3)
228 : END IF
229 :
230 8268 : DO iset = 1, nseta
231 7180 : offset_a_start = offset_a_end
232 7180 : offset_a_end = offset_a_end + nsgfa(iset)
233 7180 : sgfa = first_sgfa(1, iset)
234 :
235 7180 : offset_b_end = 0
236 22131 : DO jatom = atom_start, atom_end
237 14149 : jkind = kind_of(jatom)
238 14149 : atom_b = atom_of_kind(jatom)
239 14149 : CALL get_qs_kind(qs_kind=qs_kind_set(jkind), basis_set=basis_set_b, basis_type=basis_type_b)
240 :
241 14149 : first_sgfb => basis_set_b%first_sgf
242 14149 : lb_max => basis_set_b%lmax
243 14149 : lb_min => basis_set_b%lmin
244 14149 : nsetb = basis_set_b%nset
245 14149 : nsgfb => basis_set_b%nsgf_set
246 14149 : sphi_b => basis_set_b%sphi
247 14149 : zetb => basis_set_b%zet
248 14149 : npgfb => basis_set_b%npgf
249 :
250 14149 : rb(:) = pbc(particle_set(jatom)%r, cell)
251 :
252 14149 : IF (my_do_reflection_b) THEN
253 0 : rb(3) = 2.0_dp*reflection_z_b - rb(3)
254 : END IF
255 :
256 56596 : rab(:) = ra(:) - rb(:) ! pbc not needed?
257 : dab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
258 :
259 14149 : jset_start = 1; jset_end = nsetb
260 14149 : IF (jatom == atom_start) jset_start = set_start
261 14149 : IF (jatom == atom_end) jset_end = set_end
262 :
263 130264 : DO jset = jset_start, jset_end
264 108935 : first_set = 1; last_set = nsgfb(jset)
265 108935 : IF (jset == jset_start .AND. jatom == atom_start) first_set = set_offset_start
266 108935 : IF (jset == jset_end .AND. jatom == atom_end) last_set = set_offset_end
267 :
268 108935 : offset_b_start = offset_b_end
269 108935 : offset_b_end = offset_b_end + last_set + 1 - first_set
270 108935 : sgfb = first_sgfb(1, jset)
271 108935 : my_setpair = my_setpair + 1
272 108935 : map_it_here = MODULO(my_setpair, para_env%num_pe) == para_env%mepos
273 :
274 123084 : 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 211564 : 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 412984 : 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 286 : 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 1761230 : CALL para_env%sum(hab)
316 286 : IF (PRESENT(hdab)) CALL para_env%sum(hdab)
317 286 : IF (PRESENT(hadb)) CALL para_env%sum(hadb)
318 :
319 286 : CALL timestop(handle)
320 572 : 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 222860 : SUBROUTINE integrate_set_2c(param, potential_parameter, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, &
361 111430 : ra, rb, hab, n_hab_a, n_hab_b, offset_hab_a, offset_hab_b, &
362 111430 : offset_set_a, offset_set_b, sphi_a, sphi_b, sgfa, sgfb, nsgfa, nsgfb, &
363 111430 : 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 111430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f_work
400 111430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: hab_contr, hab_uncontr, &
401 111430 : hab_uncontr_d, pab_hh, pab_hs, &
402 111430 : pab_ss
403 111430 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hadb_contr, hadb_uncontr, hdab_contr, &
404 111430 : hdab_uncontr, v_work
405 :
406 : ! note: tested only for one exponent per pair (npgfa = npgfb = 1)
407 111430 : CALL timeset(routineN, handle)
408 :
409 111430 : my_do_reflection_a = .FALSE.
410 111430 : IF (PRESENT(do_reflection_a)) my_do_reflection_a = do_reflection_a
411 :
412 111430 : my_do_reflection_b = .FALSE.
413 111430 : IF (PRESENT(do_reflection_b)) my_do_reflection_b = do_reflection_b
414 :
415 111430 : do_force_a = PRESENT(force_a) .OR. PRESENT(hdab)
416 111430 : do_force_b = PRESENT(force_b) .OR. PRESENT(hadb)
417 111430 : calculate_forces = do_force_a .OR. do_force_b
418 :
419 111430 : 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 111430 : la_max_d = la_max
425 111430 : lb_max_d = lb_max
426 :
427 111430 : 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 111430 : ncoa = npgfa*ncoset(la_max)
433 111430 : ncob = npgfb*ncoset(lb_max)
434 :
435 111430 : rab2 = (ra(1) - rb(1))**2 + (ra(2) - rb(2))**2 + (ra(3) - rb(3))**2
436 :
437 4543524 : ALLOCATE (hab_uncontr_d(ncoset(la_max_d), ncoset(lb_max_d))); hab_uncontr_d(:, :) = 0.0_dp
438 4012824 : ALLOCATE (hab_uncontr(ncoa, ncob)); hab_uncontr(:, :) = 0.0_dp
439 111430 : IF (PRESENT(hdab)) THEN
440 655468 : ALLOCATE (hdab_uncontr(3, ncoa, ncob)); hdab_uncontr(:, :, :) = 0.0_dp
441 : END IF
442 111430 : IF (PRESENT(hadb)) THEN
443 0 : ALLOCATE (hadb_uncontr(3, ncoa, ncob)); hadb_uncontr(:, :, :) = 0.0_dp
444 : END IF
445 :
446 111430 : hab_a_start = offset_hab_a + 1; hab_a_end = offset_hab_a + n_hab_a
447 111430 : hab_b_start = offset_hab_b + 1; hab_b_end = offset_hab_b + n_hab_b
448 111430 : set_a_start = offset_set_a + 1; set_a_end = offset_set_a + n_hab_a
449 111430 : set_b_start = offset_set_b + 1; set_b_end = offset_set_b + n_hab_b
450 :
451 111430 : 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 71008 : 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 62448 : ELSE IF (eri_method == do_eri_os) THEN
535 :
536 62448 : IF (calculate_forces) CPABORT("NYI")
537 :
538 187344 : ALLOCATE (f_work(0:la_max + lb_max + 2))
539 312240 : ALLOCATE (v_work(ncoa, ncob, la_max + lb_max + 1))
540 9567511 : v_work(:, :, :) = 0.0_dp
541 367812 : f_work = 0.0_dp
542 :
543 : CALL coulomb2_new(la_max, npgfa, zeta, la_min, lb_max, npgfb, zetb, lb_min, &
544 249792 : rb - ra, rab2, hab_uncontr, v_work, f_work)
545 :
546 62448 : 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 445720 : ALLOCATE (hab_contr(nsgfa, nsgfb))
555 111430 : IF (PRESENT(hdab)) THEN
556 22592 : ALLOCATE (hdab_contr(3, nsgfa, nsgfb))
557 : END IF
558 111430 : IF (PRESENT(hadb)) THEN
559 0 : ALLOCATE (hadb_contr(3, nsgfa, nsgfb))
560 : END IF
561 :
562 111430 : CALL ab_contract(hab_contr, hab_uncontr, sphi_a(:, sgfa:), sphi_b(:, sgfb:), ncoa, ncob, nsgfa, nsgfb)
563 :
564 111430 : 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 1306846 : 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 111430 : 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 111430 : CALL timestop(handle)
587 :
588 222860 : 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 164 : SUBROUTINE mp2_eri_3c_integrate(param, potential_parameter, para_env, qs_env, &
614 164 : first_c, last_c, mat_ab, &
615 : basis_type_a, basis_type_b, basis_type_c, &
616 : sab_nl, eri_method, &
617 164 : pabc, force_a, force_b, force_c, &
618 164 : 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 164 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind
647 164 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eri_offsets
648 164 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
649 164 : lc_min, npgfa, npgfb, npgfc, nsgfa, &
650 164 : nsgfb, nsgfc
651 164 : 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 164 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: habc, pabc_block
656 164 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: habdc, hadbc, hdabc
657 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb, rc
658 164 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
659 164 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: munu_block, pab_block, rpgfb, sphi_a, &
660 164 : sphi_b, sphi_c, zeta, zetb, zetc
661 164 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
662 : TYPE(cell_type), POINTER :: cell
663 164 : 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 164 : DIMENSION(:), POINTER :: nl_iterator
667 164 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
668 164 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
669 :
670 164 : CALL timeset(routineN, handle)
671 :
672 : calculate_forces = PRESENT(force_a) .OR. PRESENT(force_b) .OR. PRESENT(force_c) .OR. &
673 164 : PRESENT(mat_dabc) .OR. PRESENT(mat_adbc) .OR. PRESENT(mat_abdc)
674 :
675 164 : my_eri_method = do_eri_mme
676 164 : IF (PRESENT(eri_method)) my_eri_method = eri_method
677 :
678 164 : IF (PRESENT(force_a) .OR. PRESENT(force_b) .OR. PRESENT(force_c)) THEN
679 26 : CPASSERT(PRESENT(pabc))
680 : END IF
681 :
682 164 : GG_count = 0; GR_count = 0; RR_count = 0
683 :
684 164 : 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 164 : atom_start = eri_offsets(first_c, 1)
690 164 : set_start = eri_offsets(first_c, 2)
691 164 : set_offset_start = eri_offsets(first_c, 3)
692 :
693 164 : atom_end = eri_offsets(last_c, 1)
694 164 : set_end = eri_offsets(last_c, 2)
695 164 : 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 164 : nkind=nkind)
705 :
706 164 : 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 164 : IF (PRESENT(force_a)) CALL mp2_eri_allocate_forces(force_a, natom_of_kind)
709 164 : IF (PRESENT(force_b)) CALL mp2_eri_allocate_forces(force_b, natom_of_kind)
710 164 : IF (PRESENT(force_c)) CALL mp2_eri_allocate_forces(force_c, natom_of_kind)
711 :
712 164 : nc = last_c - first_c + 1
713 :
714 : ! check for symmetry
715 164 : CPASSERT(SIZE(sab_nl) > 0)
716 164 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
717 :
718 164 : IF (do_symmetric) THEN
719 164 : CPASSERT(basis_type_a == basis_type_b)
720 : END IF
721 :
722 1328 : ALLOCATE (basis_set_list_a(nkind), basis_set_list_b(nkind))
723 164 : CALL basis_set_list_setup(basis_set_list_a, basis_type_a, qs_kind_set)
724 164 : CALL basis_set_list_setup(basis_set_list_b, basis_type_b, qs_kind_set)
725 :
726 164 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
727 :
728 164 : mepos = 0
729 :
730 11674 : 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 11510 : iatom=iatom, jatom=jatom, r=rab)
733 :
734 : ! exclude periodic images because method is periodic intrinsically
735 11510 : IF (inode == 1) last_jatom = 0
736 :
737 11510 : IF (jatom /= last_jatom) THEN
738 738 : last_jatom = jatom
739 : ELSE
740 : CYCLE
741 : END IF
742 :
743 738 : basis_set_a => basis_set_list_a(ikind)%gto_basis_set
744 738 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
745 738 : basis_set_b => basis_set_list_b(jkind)%gto_basis_set
746 738 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
747 738 : atom_a = atom_of_kind(iatom)
748 738 : atom_b = atom_of_kind(jatom)
749 :
750 738 : first_sgfa => basis_set_a%first_sgf
751 738 : la_max => basis_set_a%lmax
752 738 : la_min => basis_set_a%lmin
753 738 : npgfa => basis_set_a%npgf
754 738 : nseta = basis_set_a%nset
755 738 : nsgfa => basis_set_a%nsgf_set
756 738 : set_radius_a => basis_set_a%set_radius
757 738 : sphi_a => basis_set_a%sphi
758 738 : zeta => basis_set_a%zet
759 2484 : na = SUM(nsgfa)
760 :
761 738 : ra(:) = pbc(particle_set(iatom)%r, cell)
762 :
763 : ! basis jkind
764 738 : first_sgfb => basis_set_b%first_sgf
765 738 : lb_max => basis_set_b%lmax
766 738 : lb_min => basis_set_b%lmin
767 738 : npgfb => basis_set_b%npgf
768 738 : nsetb = basis_set_b%nset
769 738 : nsgfb => basis_set_b%nsgf_set
770 738 : rpgfb => basis_set_b%pgf_radius
771 738 : set_radius_b => basis_set_b%set_radius
772 738 : sphi_b => basis_set_b%sphi
773 738 : zetb => basis_set_b%zet
774 2484 : nb = SUM(nsgfb)
775 :
776 738 : rb(:) = pbc(particle_set(jatom)%r, cell)
777 :
778 738 : IF (do_symmetric) THEN
779 738 : IF (iatom <= jatom) THEN
780 492 : irow = iatom
781 492 : icol = jatom
782 : ELSE
783 246 : irow = jatom
784 246 : icol = iatom
785 : END IF
786 : ELSE
787 0 : irow = iatom
788 0 : icol = jatom
789 : END IF
790 :
791 3690 : ALLOCATE (habc(na, nb, nc))
792 1299180 : habc(:, :, :) = 0.0_dp ! needs to be initialized due to screening
793 738 : IF (PRESENT(mat_dabc)) THEN
794 0 : ALLOCATE (hdabc(3, na, nb, nc))
795 0 : hdabc(:, :, :, :) = 0.0_dp
796 : END IF
797 738 : IF (PRESENT(mat_adbc)) THEN
798 0 : ALLOCATE (hadbc(3, na, nb, nc))
799 0 : hadbc(:, :, :, :) = 0.0_dp
800 : END IF
801 738 : IF (PRESENT(mat_abdc)) THEN
802 0 : ALLOCATE (habdc(3, na, nb, nc))
803 0 : habdc(:, :, :, :) = 0.0_dp
804 : END IF
805 :
806 738 : IF (calculate_forces .AND. PRESENT(pabc)) THEN
807 675 : ALLOCATE (pabc_block(na, nb, nc))
808 5613 : DO ic = 1, nc
809 5478 : NULLIFY (pab_block)
810 : CALL dbcsr_get_block_p(matrix=pabc(ic)%matrix, &
811 5478 : row=irow, col=icol, block=pab_block, found=found)
812 5478 : CPASSERT(found)
813 11091 : IF (irow .EQ. iatom) THEN
814 3652 : to_be_asserted = SIZE(pab_block, 1) .EQ. SIZE(pabc_block, 1) .AND. SIZE(pab_block, 2) .EQ. SIZE(pabc_block, 2)
815 0 : CPASSERT(to_be_asserted)
816 188567 : pabc_block(:, :, ic) = pab_block(:, :)
817 : ELSE
818 1826 : to_be_asserted = SIZE(pab_block, 2) .EQ. SIZE(pabc_block, 1) .AND. SIZE(pab_block, 1) .EQ. SIZE(pabc_block, 2)
819 0 : CPASSERT(to_be_asserted)
820 72496 : pabc_block(:, :, ic) = TRANSPOSE(pab_block(:, :))
821 : END IF
822 : END DO
823 : END IF
824 :
825 2952 : rab(:) = pbc(rab, cell)
826 738 : dab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
827 :
828 738 : offset_a_end = 0
829 2484 : DO iset = 1, nseta
830 1746 : offset_a_start = offset_a_end
831 1746 : offset_a_end = offset_a_end + nsgfa(iset)
832 1746 : sgfa = first_sgfa(1, iset)
833 :
834 1746 : offset_b_end = 0
835 7272 : DO jset = 1, nsetb
836 4788 : offset_b_start = offset_b_end
837 4788 : offset_b_end = offset_b_end + nsgfb(jset)
838 :
839 4788 : sgfb = first_sgfb(1, jset)
840 :
841 : ! Screening
842 4788 : IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
843 :
844 : offset_c_end = 0
845 14640 : DO katom = atom_start, atom_end
846 :
847 8106 : atom_c = atom_of_kind(katom)
848 :
849 8106 : kkind = kind_of(katom)
850 8106 : CALL get_qs_kind(qs_kind=qs_kind_set(kkind), basis_set=basis_set_c, basis_type=basis_type_c)
851 8106 : first_sgfc => basis_set_c%first_sgf
852 8106 : lc_max => basis_set_c%lmax
853 8106 : lc_min => basis_set_c%lmin
854 8106 : nsetc = basis_set_c%nset
855 8106 : nsgfc => basis_set_c%nsgf_set
856 8106 : sphi_c => basis_set_c%sphi
857 8106 : zetc => basis_set_c%zet
858 8106 : npgfc => basis_set_c%npgf
859 :
860 8106 : rc(:) = pbc(particle_set(katom)%r, cell)
861 :
862 8106 : kset_start = 1; kset_end = nsetc
863 8106 : IF (katom == atom_start) kset_start = set_start
864 8106 : IF (katom == atom_end) kset_end = set_end
865 :
866 72834 : DO kset = kset_start, kset_end
867 59940 : first_set = 1; last_set = nsgfc(kset)
868 59940 : IF (kset == kset_start .AND. katom == atom_start) first_set = set_offset_start
869 59940 : IF (kset == kset_end .AND. katom == atom_end) last_set = set_offset_end
870 :
871 59940 : offset_c_start = offset_c_end
872 59940 : offset_c_end = offset_c_end + last_set + 1 - first_set
873 59940 : sgfc = first_sgfc(1, kset)
874 :
875 : #!some fypp magic to deal with combinations of optional arguments
876 : #:for pabc_present in [0, 1]
877 : #:for doforce_1 in [0, 1]
878 : #:for doforce_2 in [0, 1]
879 : #:for doforce_3 in [0, 1]
880 : #:for dabc in [0, 1]
881 : #:for adbc in [0, 1]
882 : #:for abdc in [0, 1]
883 : IF (${conditional(doforce_1)}$PRESENT(force_a) .AND. &
884 : ${conditional(doforce_2)}$PRESENT(force_b) .AND. &
885 : ${conditional(doforce_3)}$PRESENT(force_c) .AND. &
886 : ${conditional(pabc_present)}$PRESENT(pabc) .AND. &
887 : ${conditional(dabc)}$PRESENT(mat_dabc) .AND. &
888 119880 : ${conditional(adbc)}$PRESENT(mat_adbc) .AND. &
889 8106 : ${conditional(abdc)}$PRESENT(mat_abdc)) THEN
890 : CALL integrate_set_3c( &
891 : param%par, potential_parameter, &
892 : la_min(iset), la_max(iset), &
893 : lb_min(jset), lb_max(jset), &
894 : lc_min(kset), lc_max(kset), &
895 : npgfa(iset), npgfb(jset), npgfc(kset), &
896 : zeta(:, iset), zetb(:, jset), zetc(:, kset), &
897 : ra, rb, rc, &
898 : habc, &
899 : nsgfa(iset), nsgfb(jset), last_set - first_set + 1, &
900 : offset_a_start, offset_b_start, offset_c_start, &
901 : 0, 0, first_set - 1, &
902 : sphi_a, sphi_b, sphi_c, &
903 : sgfa, sgfb, sgfc, &
904 : nsgfa(iset), nsgfb(jset), nsgfc(kset), &
905 : my_eri_method, &
906 : $: 'pabc=pabc_block, &'*pabc_present
907 : $: 'force_a=force_a(ikind)%forces(:, atom_a), &'*doforce_1
908 : $: 'force_b=force_b(jkind)%forces(:, atom_b), &'*doforce_2
909 : $: 'force_c=force_c(kkind)%forces(:, atom_c), &'*doforce_3
910 : do_symmetric=do_symmetric, &
911 : on_diagonal=iatom .EQ. jatom, &
912 : $: 'hdabc=hdabc, &'*dabc
913 : $: 'hadbc=hadbc, &'*adbc
914 : $: 'habdc=habdc, &'*abdc
915 59940 : GG_count=GG_count, GR_count=GR_count, RR_count=RR_count)
916 : END IF
917 : #:endfor
918 : #:endfor
919 : #:endfor
920 : #:endfor
921 : #:endfor
922 : #:endfor
923 : #:endfor
924 : END DO
925 : END DO
926 : END DO
927 : END DO
928 :
929 738 : IF (calculate_forces .AND. PRESENT(pabc)) DEALLOCATE (pabc_block)
930 26142 : DO ic = 1, nc
931 25404 : NULLIFY (munu_block)
932 : CALL dbcsr_get_block_p(matrix=mat_ab(ic)%matrix, &
933 25404 : row=irow, col=icol, block=munu_block, found=found)
934 25404 : CPASSERT(found)
935 1276620 : munu_block(:, :) = 0.0_dp
936 51546 : IF (irow .EQ. iatom) THEN
937 16936 : to_be_asserted = SIZE(munu_block, 1) .EQ. SIZE(habc, 1) .AND. SIZE(munu_block, 2) .EQ. SIZE(habc, 2)
938 0 : CPASSERT(to_be_asserted)
939 938377 : munu_block(:, :) = habc(:, :, ic)
940 : ELSE
941 8468 : to_be_asserted = SIZE(munu_block, 2) .EQ. SIZE(habc, 1) .AND. SIZE(munu_block, 1) .EQ. SIZE(habc, 2)
942 0 : CPASSERT(to_be_asserted)
943 338243 : munu_block(:, :) = TRANSPOSE(habc(:, :, ic))
944 : END IF
945 : END DO
946 738 : DEALLOCATE (habc)
947 902 : IF (calculate_forces) THEN
948 5613 : DO ic = 1, nc
949 22047 : DO i_xyz = 1, 3
950 16434 : IF (PRESENT(mat_dabc)) THEN
951 0 : NULLIFY (munu_block)
952 : CALL dbcsr_get_block_p(matrix=mat_dabc(i_xyz, ic)%matrix, &
953 0 : row=irow, col=icol, block=munu_block, found=found)
954 0 : CPASSERT(found)
955 0 : munu_block(:, :) = 0.0_dp
956 0 : IF (irow .EQ. iatom) THEN
957 0 : munu_block(:, :) = hdabc(i_xyz, :, :, ic)
958 : ELSE
959 0 : munu_block(:, :) = TRANSPOSE(hdabc(i_xyz, :, :, ic))
960 : END IF
961 : END IF
962 16434 : IF (PRESENT(mat_adbc)) THEN
963 0 : NULLIFY (munu_block)
964 : CALL dbcsr_get_block_p(matrix=mat_adbc(i_xyz, ic)%matrix, &
965 0 : row=irow, col=icol, block=munu_block, found=found)
966 0 : CPASSERT(found)
967 0 : munu_block(:, :) = 0.0_dp
968 0 : IF (irow .EQ. iatom) THEN
969 0 : munu_block(:, :) = hadbc(i_xyz, :, :, ic)
970 : ELSE
971 0 : munu_block(:, :) = TRANSPOSE(hadbc(i_xyz, :, :, ic))
972 : END IF
973 : END IF
974 21912 : IF (PRESENT(mat_abdc)) THEN
975 0 : NULLIFY (munu_block)
976 : CALL dbcsr_get_block_p(matrix=mat_abdc(i_xyz, ic)%matrix, &
977 0 : row=irow, col=icol, block=munu_block, found=found)
978 0 : CPASSERT(found)
979 0 : munu_block(:, :) = 0.0_dp
980 0 : IF (irow .EQ. iatom) THEN
981 0 : munu_block(:, :) = habdc(i_xyz, :, :, ic)
982 : ELSE
983 0 : munu_block(:, :) = TRANSPOSE(habdc(i_xyz, :, :, ic))
984 : END IF
985 : END IF
986 : END DO
987 : END DO
988 135 : IF (PRESENT(mat_dabc)) DEALLOCATE (hdabc)
989 135 : IF (PRESENT(mat_adbc)) DEALLOCATE (hadbc)
990 135 : IF (PRESENT(mat_abdc)) DEALLOCATE (habdc)
991 : END IF
992 : END DO
993 :
994 164 : DEALLOCATE (basis_set_list_a, basis_set_list_b)
995 164 : CALL neighbor_list_iterator_release(nl_iterator)
996 :
997 164 : 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)
998 :
999 164 : CALL timestop(handle)
1000 656 : END SUBROUTINE mp2_eri_3c_integrate
1001 :
1002 : ! **************************************************************************************************
1003 : !> \brief Integrate set triple and contract with sphi matrix
1004 : !> \param param ...
1005 : !> \param potential_parameter ...
1006 : !> \param la_min ...
1007 : !> \param la_max ...
1008 : !> \param lb_min ...
1009 : !> \param lb_max ...
1010 : !> \param lc_min ...
1011 : !> \param lc_max ...
1012 : !> \param npgfa ...
1013 : !> \param npgfb ...
1014 : !> \param npgfc ...
1015 : !> \param zeta ...
1016 : !> \param zetb ...
1017 : !> \param zetc ...
1018 : !> \param ra ...
1019 : !> \param rb ...
1020 : !> \param rc ...
1021 : !> \param habc ...
1022 : !> \param n_habc_a ...
1023 : !> \param n_habc_b ...
1024 : !> \param n_habc_c ...
1025 : !> \param offset_habc_a ...
1026 : !> \param offset_habc_b ...
1027 : !> \param offset_habc_c ...
1028 : !> \param offset_set_a ...
1029 : !> \param offset_set_b ...
1030 : !> \param offset_set_c ...
1031 : !> \param sphi_a ...
1032 : !> \param sphi_b ...
1033 : !> \param sphi_c ...
1034 : !> \param sgfa ...
1035 : !> \param sgfb ...
1036 : !> \param sgfc ...
1037 : !> \param nsgfa ...
1038 : !> \param nsgfb ...
1039 : !> \param nsgfc ...
1040 : !> \param eri_method ...
1041 : !> \param pabc ...
1042 : !> \param force_a ...
1043 : !> \param force_b ...
1044 : !> \param force_c ...
1045 : !> \param do_symmetric ...
1046 : !> \param on_diagonal ...
1047 : !> \param hdabc ...
1048 : !> \param hadbc ...
1049 : !> \param habdc ...
1050 : !> \param GG_count ...
1051 : !> \param GR_count ...
1052 : !> \param RR_count ...
1053 : !> \note
1054 : ! **************************************************************************************************
1055 59940 : SUBROUTINE integrate_set_3c(param, potential_parameter, &
1056 : la_min, la_max, lb_min, lb_max, lc_min, lc_max, &
1057 : npgfa, npgfb, npgfc, &
1058 59940 : zeta, zetb, zetc, &
1059 : ra, rb, rc, &
1060 119880 : habc, &
1061 : n_habc_a, n_habc_b, n_habc_c, &
1062 : offset_habc_a, offset_habc_b, offset_habc_c, &
1063 : offset_set_a, offset_set_b, offset_set_c, &
1064 59940 : sphi_a, sphi_b, sphi_c, &
1065 : sgfa, sgfb, sgfc, &
1066 : nsgfa, nsgfb, nsgfc, &
1067 : eri_method, &
1068 59940 : pabc, &
1069 : force_a, force_b, force_c, &
1070 : do_symmetric, on_diagonal, &
1071 59940 : hdabc, hadbc, habdc, &
1072 : GG_count, GR_count, RR_count)
1073 :
1074 : TYPE(eri_mme_param), INTENT(INOUT) :: param
1075 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1076 : INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, lc_min, &
1077 : lc_max, npgfa, npgfb, npgfc
1078 : REAL(KIND=dp), DIMENSION(npgfa), INTENT(IN) :: zeta
1079 : REAL(KIND=dp), DIMENSION(npgfb), INTENT(IN) :: zetb
1080 : REAL(KIND=dp), DIMENSION(npgfc), INTENT(IN) :: zetc
1081 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb, rc
1082 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: habc
1083 : INTEGER, INTENT(IN) :: n_habc_a, n_habc_b, n_habc_c, offset_habc_a, offset_habc_b, &
1084 : offset_habc_c, offset_set_a, offset_set_b, offset_set_c
1085 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sphi_a, sphi_b, sphi_c
1086 : INTEGER, INTENT(IN) :: sgfa, sgfb, sgfc, nsgfa, nsgfb, nsgfc, &
1087 : eri_method
1088 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
1089 : OPTIONAL :: pabc
1090 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT), &
1091 : OPTIONAL :: force_a, force_b, force_c
1092 : LOGICAL, INTENT(IN) :: do_symmetric
1093 : LOGICAL, INTENT(IN), OPTIONAL :: on_diagonal
1094 : REAL(KIND=dp), DIMENSION(:, :, :, :), &
1095 : INTENT(OUT), OPTIONAL :: hdabc, hadbc, habdc
1096 : INTEGER, INTENT(INOUT), OPTIONAL :: GG_count, GR_count, RR_count
1097 :
1098 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_set_3c'
1099 :
1100 : INTEGER :: ax, ay, az, bx, by, bz, cx, cy, cz, habc_a_end, habc_a_start, habc_b_end, &
1101 : habc_b_start, habc_c_end, habc_c_start, handle, i_xyz, ico, icoc, icox, icoy, icoz, ipgf, &
1102 : jco, jcox, jcoy, jcoz, jpgf, kco, kcox, kcoy, kcoz, kpgf, la, la_max_d, lb, &
1103 : lb_max_d, lc, lc_max_d, na, nb, nc, nc_end, nc_start, ncoa, ncoa_d, ncob, &
1104 : ncob_d, ncoc, ncoc_d, set_a_end, set_a_start, set_b_end, set_b_start, &
1105 : set_c_end, set_c_start, sphi_a_start, sphi_b_start, sphi_c_start
1106 : INTEGER, DIMENSION(3) :: la_xyz, lb_xyz
1107 : LOGICAL :: calculate_forces, do_force_a, do_force_b, &
1108 : do_force_c
1109 : REAL(KIND=dp) :: rab2, rac2, rbc2, w
1110 59940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: f_work, gccc, rpgfa, rpgfb, rpgfc
1111 59940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pab_hh, pab_hs, vabc
1112 59940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: habc_contr, habc_uncontr, &
1113 59940 : habc_uncontr_d, pabc_hhh, &
1114 59940 : pabc_hsh, pabc_hss, pabc_sss
1115 59940 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: habdc_contr, habdc_uncontr, hadbc_contr, &
1116 59940 : hadbc_uncontr, hdabc_contr, &
1117 59940 : hdabc_uncontr, v_work
1118 :
1119 59940 : CALL timeset(routineN, handle)
1120 :
1121 59940 : do_force_a = PRESENT(force_a) .OR. PRESENT(hdabc)
1122 59940 : do_force_b = PRESENT(force_b) .OR. PRESENT(hadbc)
1123 59940 : do_force_c = PRESENT(force_c) .OR. PRESENT(habdc)
1124 59940 : calculate_forces = do_force_a .OR. do_force_b .OR. do_force_c
1125 :
1126 59940 : IF (do_symmetric) THEN
1127 59940 : CPASSERT(PRESENT(on_diagonal))
1128 : END IF
1129 :
1130 59940 : la_max_d = la_max
1131 59940 : lb_max_d = lb_max
1132 59940 : lc_max_d = lc_max
1133 :
1134 59940 : IF (calculate_forces) THEN
1135 10386 : IF (do_force_a) la_max_d = la_max + 1
1136 10386 : IF (do_force_b) lb_max_d = lb_max + 1
1137 10386 : IF (do_force_c) lc_max_d = lc_max + 1
1138 : END IF
1139 :
1140 59940 : ncoa = npgfa*ncoset(la_max)
1141 59940 : ncob = npgfb*ncoset(lb_max)
1142 59940 : ncoc = npgfc*ncoset(lc_max)
1143 :
1144 59940 : ncoa_d = npgfa*ncoset(la_max_d)
1145 59940 : ncob_d = npgfb*ncoset(lb_max_d)
1146 59940 : ncoc_d = npgfc*ncoset(lc_max_d)
1147 :
1148 299700 : ALLOCATE (habc_uncontr_d(ncoset(la_max_d), ncoset(lb_max_d), ncoset(lc_max_d)))
1149 15444036 : habc_uncontr_d(:, :, :) = 0.0_dp
1150 10689432 : ALLOCATE (habc_uncontr(ncoa, ncob, ncoc)); habc_uncontr(:, :, :) = 0.0_dp
1151 59940 : IF (PRESENT(hdabc)) THEN
1152 0 : ALLOCATE (hdabc_uncontr(3, ncoa, ncob, ncoc)); hdabc_uncontr(:, :, :, :) = 0.0_dp
1153 : END IF
1154 59940 : IF (PRESENT(hadbc)) THEN
1155 0 : ALLOCATE (hadbc_uncontr(3, ncoa, ncob, ncoc)); hadbc_uncontr(:, :, :, :) = 0.0_dp
1156 : END IF
1157 59940 : IF (PRESENT(habdc)) THEN
1158 0 : ALLOCATE (habdc_uncontr(3, ncoa, ncob, ncoc)); habdc_uncontr(:, :, :, :) = 0.0_dp
1159 : END IF
1160 :
1161 59940 : habc_a_start = offset_habc_a + 1; habc_a_end = offset_habc_a + n_habc_a
1162 59940 : habc_b_start = offset_habc_b + 1; habc_b_end = offset_habc_b + n_habc_b
1163 59940 : habc_c_start = offset_habc_c + 1; habc_c_end = offset_habc_c + n_habc_c
1164 59940 : set_a_start = offset_set_a + 1; set_a_end = offset_set_a + n_habc_a
1165 59940 : set_b_start = offset_set_b + 1; set_b_end = offset_set_b + n_habc_b
1166 59940 : set_c_start = offset_set_c + 1; set_c_end = offset_set_c + n_habc_c
1167 :
1168 59940 : IF (eri_method == do_eri_mme) THEN
1169 36234 : CALL eri_mme_set_potential(param, convert_potential_type(potential_parameter%potential_type), potential_parameter%omega)
1170 :
1171 36234 : IF (calculate_forces .AND. PRESENT(pabc)) THEN
1172 : ! uncontracted hermite-gaussian representation of density matrix
1173 10386 : sphi_a_start = sgfa - 1 + set_a_start
1174 10386 : sphi_b_start = sgfb - 1 + set_b_start
1175 10386 : sphi_c_start = sgfc - 1 + set_c_start
1176 :
1177 51930 : ALLOCATE (pabc_sss(n_habc_a, n_habc_b, n_habc_c))
1178 342953 : pabc_sss(:, :, :) = pabc(habc_a_start:habc_a_end, habc_b_start:habc_b_end, habc_c_start:habc_c_end)
1179 51930 : ALLOCATE (pabc_hss(ncoa, n_habc_b, n_habc_c))
1180 51930 : ALLOCATE (pabc_hsh(ncoa, n_habc_b, ncoc))
1181 51930 : ALLOCATE (pabc_hhh(ncoa, ncob, ncoc))
1182 41544 : ALLOCATE (pab_hs(ncoa, n_habc_b))
1183 41544 : ALLOCATE (pab_hh(ncoa, ncob))
1184 :
1185 : CALL dgemm("N", "N", ncoa, n_habc_b*n_habc_c, n_habc_a, 1.0_dp, &
1186 10386 : sphi_a(:, sphi_a_start), SIZE(sphi_a, 1), pabc_sss, n_habc_a, 0.0_dp, pabc_hss, ncoa)
1187 : CALL dgemm("N", "T", ncoa*n_habc_b, ncoc, n_habc_c, 1.0_dp, &
1188 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)
1189 :
1190 68148 : DO icoc = 1, ncoc
1191 1108764 : pab_hs(:, :) = pabc_hsh(:, :, icoc)
1192 : CALL dgemm("N", "T", ncoa, ncob, n_habc_b, 1.0_dp, &
1193 57762 : pab_hs, ncoa, sphi_b(:, sphi_b_start), SIZE(sphi_b, 1), 0.0_dp, pab_hh, ncoa)
1194 2262300 : pabc_hhh(:, :, icoc) = pab_hh(:, :)
1195 : END DO
1196 : END IF
1197 :
1198 105444 : DO ipgf = 1, npgfa
1199 69210 : na = (ipgf - 1)*ncoset(la_max)
1200 253554 : DO jpgf = 1, npgfb
1201 148110 : nb = (jpgf - 1)*ncoset(lb_max)
1202 365430 : DO kpgf = 1, npgfc
1203 148110 : nc = (kpgf - 1)*ncoset(lc_max)
1204 37173750 : habc_uncontr_d(:, :, :) = 0.0_dp
1205 : CALL eri_mme_3c_integrate(param, &
1206 : la_min, la_max_d, lb_min, lb_max_d, lc_min, lc_max_d, &
1207 : zeta(ipgf), zetb(jpgf), zetc(kpgf), ra, rb, rc, habc_uncontr_d, 0, 0, 0, &
1208 148110 : GG_count, GR_count, RR_count)
1209 :
1210 : habc_uncontr(na + 1:na + ncoset(la_max), nb + 1:nb + ncoset(lb_max), nc + 1:nc + ncoset(lc_max)) = &
1211 8079720 : habc_uncontr_d(:ncoset(la_max), :ncoset(lb_max), :ncoset(lc_max))
1212 :
1213 296220 : IF (calculate_forces) THEN
1214 82500 : DO lc = lc_min, lc_max
1215 164700 : DO cx = 0, lc
1216 266850 : DO cy = 0, lc - cx
1217 143400 : cz = lc - cx - cy
1218 143400 : kco = coset(cx, cy, cz)
1219 143400 : kcox = coset(cx + 1, cy, cz)
1220 143400 : kcoy = coset(cx, cy + 1, cz)
1221 143400 : kcoz = coset(cx, cy, cz + 1)
1222 414600 : DO lb = lb_min, lb_max
1223 592800 : DO bx = 0, lb
1224 788400 : DO by = 0, lb - bx
1225 339000 : bz = lb - bx - by
1226 339000 : jco = coset(bx, by, bz)
1227 339000 : jcox = coset(bx + 1, by, bz)
1228 339000 : jcoy = coset(bx, by + 1, bz)
1229 339000 : jcoz = coset(bx, by, bz + 1)
1230 1075620 : DO la = la_min, la_max
1231 1500105 : DO ax = 0, la
1232 2078460 : DO ay = 0, la - ax
1233 917355 : az = la - ax - ay
1234 3669420 : la_xyz = [ax, ay, az]
1235 3669420 : lb_xyz = [bx, by, bz]
1236 917355 : ico = coset(ax, ay, az)
1237 917355 : icox = coset(ax + 1, ay, az)
1238 917355 : icoy = coset(ax, ay + 1, az)
1239 917355 : icoz = coset(ax, ay, az + 1)
1240 :
1241 917355 : w = 1.0_dp
1242 917355 : IF (do_symmetric .AND. .NOT. on_diagonal) w = 2.0_dp
1243 :
1244 917355 : IF (PRESENT(force_a)) THEN
1245 : force_a = force_a + 2.0_dp*w*zeta(ipgf)* &
1246 : [pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(icox, jco, kco), &
1247 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(icoy, jco, kco), &
1248 3669420 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(icoz, jco, kco)]
1249 :
1250 : END IF
1251 917355 : IF (PRESENT(force_b)) THEN
1252 : force_b = force_b + 2.0_dp*w*zetb(jpgf)* &
1253 : [pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jcox, kco), &
1254 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jcoy, kco), &
1255 3669420 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jcoz, kco)]
1256 : END IF
1257 917355 : IF (PRESENT(force_c)) THEN
1258 : force_c = force_c + 2.0_dp*w*zetc(kpgf)* &
1259 : [pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jco, kcox), &
1260 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jco, kcoy), &
1261 3669420 : pabc_hhh(na + ico, nb + jco, nc + kco)*habc_uncontr_d(ico, jco, kcoz)]
1262 : END IF
1263 :
1264 917355 : IF (PRESENT(hdabc)) THEN
1265 : hdabc_uncontr(1:3, na + ico, nb + jco, nc + kco) = 2.0_dp*zeta(ipgf)* &
1266 : [habc_uncontr_d(icox, jco, kco), &
1267 : habc_uncontr_d(icoy, jco, kco), &
1268 0 : habc_uncontr_d(icoz, jco, kco)]
1269 : END IF
1270 917355 : IF (PRESENT(hadbc)) THEN
1271 : hadbc_uncontr(1:3, na + ico, nb + jco, nc + kco) = 2.0_dp*zetb(jpgf)* &
1272 : [habc_uncontr_d(ico, jcox, kco), &
1273 : habc_uncontr_d(ico, jcoy, kco), &
1274 0 : habc_uncontr_d(ico, jcoz, kco)]
1275 : END IF
1276 1602240 : IF (PRESENT(habdc)) THEN
1277 : habdc_uncontr(1:3, na + ico, nb + jco, nc + kco) = 2.0_dp*zetc(kpgf)* &
1278 : [habc_uncontr_d(ico, jco, kcox), &
1279 : habc_uncontr_d(ico, jco, kcoy), &
1280 0 : habc_uncontr_d(ico, jco, kcoz)]
1281 : END IF
1282 : END DO
1283 : END DO
1284 : END DO
1285 : END DO
1286 : END DO
1287 : END DO
1288 : END DO
1289 : END DO
1290 : END DO
1291 : END IF
1292 :
1293 : END DO
1294 : END DO
1295 : END DO
1296 :
1297 23706 : ELSE IF (eri_method == do_eri_os) THEN
1298 :
1299 23706 : IF (calculate_forces) CPABORT("NYI")
1300 :
1301 71118 : ALLOCATE (f_work(0:la_max + lb_max + lc_max + 2))
1302 151704 : f_work(:) = 0.0_dp
1303 142236 : ALLOCATE (v_work(ncoa, ncob, ncoc, la_max + lb_max + lc_max + 1))
1304 20791632 : v_work(:, :, :, :) = 0.0_dp
1305 : ! no screening
1306 71118 : ALLOCATE (rpgfa(npgfa))
1307 71118 : ALLOCATE (rpgfb(npgfb))
1308 71118 : ALLOCATE (rpgfc(npgfc))
1309 61956 : rpgfa(:) = 1.0E10_dp
1310 61956 : rpgfb(:) = 1.0E10_dp
1311 47412 : rpgfc(:) = 1.0E10_dp
1312 71118 : ALLOCATE (gccc(ncoc))
1313 156600 : gccc(:) = 0.0_dp
1314 94824 : ALLOCATE (vabc(ncoa, ncob))
1315 710010 : vabc(:, :) = 0.0_dp
1316 23706 : rab2 = (rb(1) - ra(1))**2 + (rb(2) - ra(2))**2 + (rb(3) - ra(3))**2
1317 23706 : rac2 = (rc(1) - ra(1))**2 + (rc(2) - ra(2))**2 + (rc(3) - ra(3))**2
1318 23706 : rbc2 = (rc(1) - rb(1))**2 + (rc(2) - rb(2))**2 + (rc(3) - rb(3))**2
1319 :
1320 : ! in the RI basis, there is only a single primitive Gaussian
1321 23706 : kpgf = 1
1322 :
1323 23706 : IF (lc_max == 0) THEN
1324 : nc_start = 1
1325 : ELSE
1326 14256 : nc_start = ncoset(lc_max - 1) + 1
1327 : END IF
1328 23706 : 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 189648 : vabc, habc_uncontr(:, :, nc_start:nc_end), v_work, f_work)
1335 :
1336 23706 : 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 299700 : ALLOCATE (habc_contr(nsgfa, nsgfb, nsgfc))
1345 59940 : IF (PRESENT(hdabc)) THEN
1346 0 : ALLOCATE (hdabc_contr(3, nsgfa, nsgfb, nsgfc))
1347 : END IF
1348 59940 : IF (PRESENT(hadbc)) THEN
1349 0 : ALLOCATE (hadbc_contr(3, nsgfa, nsgfb, nsgfc))
1350 : END IF
1351 59940 : 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 59940 : ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
1358 :
1359 59940 : 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 1763298 : habc_contr(set_a_start:set_a_end, set_b_start:set_b_end, set_c_start:set_c_end)
1381 :
1382 59940 : 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 59940 : CALL timestop(handle)
1392 :
1393 119880 : 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 450 : 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 450 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
1413 450 : INTEGER, DIMENSION(:), POINTER :: nsgf_set
1414 450 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1415 : TYPE(gto_basis_set_type), POINTER :: basis_set
1416 450 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1417 450 : 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 450 : particle_set=particle_set, natom=natom, nkind=nkind)
1421 :
1422 450 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
1423 :
1424 450 : dimen_basis = 0
1425 1666 : DO iatom = 1, natom
1426 1216 : ikind = kind_of(iatom)
1427 1216 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type=basis_type)
1428 1666 : dimen_basis = dimen_basis + nsgf
1429 : END DO
1430 :
1431 1350 : ALLOCATE (eri_offsets(dimen_basis, 3))
1432 :
1433 450 : offset = 0
1434 1666 : DO iatom = 1, natom
1435 1216 : ikind = kind_of(iatom)
1436 1216 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set, basis_type=basis_type)
1437 1216 : nset = basis_set%nset
1438 1216 : nsgf_set => basis_set%nsgf_set
1439 12408 : DO iset = 1, nset
1440 10742 : set_offset = 0
1441 40620 : DO isgf = 1, nsgf_set(iset)
1442 29878 : set_offset = set_offset + 1
1443 130254 : eri_offsets(offset + set_offset, :) = [iatom, iset, set_offset]
1444 : END DO
1445 11958 : offset = offset + nsgf_set(iset)
1446 : END DO
1447 : END DO
1448 900 : 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 85216 : FUNCTION convert_potential_type(potential_type) RESULT(res)
1494 : INTEGER, INTENT(IN) :: potential_type
1495 : INTEGER :: res
1496 :
1497 85216 : 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 85216 : END FUNCTION convert_potential_type
1506 :
1507 0 : END MODULE mp2_eri
|
