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 Calculates integrals for LRIGPW method
10 : !> lri : local resolution of the identity
11 : !> \par History
12 : !> created JGH [08.2012]
13 : !> Dorothea Golze [02.2014] (1) extended, re-structured, cleaned
14 : !> (2) heavily debugged
15 : !> split off JGH [11.2017]
16 : !> \authors JGH
17 : !> Dorothea Golze
18 : ! **************************************************************************************************
19 : MODULE lri_integrals
20 : USE basis_set_types, ONLY: gto_basis_set_type
21 : USE generic_os_integrals, ONLY: int_overlap_ab_os,&
22 : int_overlap_aba_os,&
23 : int_overlap_abb_os
24 : USE generic_shg_integrals, ONLY: get_abb_same_kind,&
25 : int_overlap_ab_shg_low,&
26 : int_overlap_aba_shg_low,&
27 : int_overlap_abb_shg_low,&
28 : lri_precalc_angular_shg_part
29 : USE kinds, ONLY: dp
30 : USE lri_environment_types, ONLY: lri_environment_type,&
31 : lri_int_type
32 : #include "./base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : ! **************************************************************************************************
39 : TYPE int_type
40 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: sabint
41 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: sooint
42 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: abaint
43 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: abbint
44 : END TYPE int_type
45 : !
46 : TYPE dint_type
47 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dsabint
48 : REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dsooint
49 : REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: dabaint
50 : REAL(KIND=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: dabbint
51 : END TYPE dint_type
52 : ! **************************************************************************************************
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_integrals'
55 :
56 : PUBLIC :: lri_int, lri_int2, lri_dint, lri_dint2, int_type, dint_type, allocate_int_type, deallocate_int_type
57 :
58 : ! **************************************************************************************************
59 :
60 : CONTAINS
61 :
62 : ! **************************************************************************************************
63 : !> \brief calcuates the lri integrals using solid harmonic Gaussians
64 : !> \param lri_env ...
65 : !> \param lrii ...
66 : !> \param rab distance vector
67 : !> \param obasa orb basis on A
68 : !> \param obasb orb basis on B
69 : !> \param fbasa aux basis on A
70 : !> \param fbasb aux basis on B
71 : !> \param iatom index atom A
72 : !> \param jatom index atom B
73 : !> \param ikind kind atom A
74 : !> \param jkind kind atom B
75 : !> \param calculate_forces ...
76 : ! **************************************************************************************************
77 0 : SUBROUTINE lri_int(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
78 : iatom, jatom, ikind, jkind, calculate_forces)
79 :
80 : TYPE(lri_environment_type), POINTER :: lri_env
81 : TYPE(lri_int_type), POINTER :: lrii
82 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
83 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
84 : INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
85 : LOGICAL, INTENT(IN) :: calculate_forces
86 :
87 0 : IF (lri_env%use_shg_integrals) THEN
88 : CALL lri_int_shg(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
89 0 : iatom, jatom, ikind, jkind, calculate_forces)
90 : ELSE
91 : CALL lri_int_os(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
92 0 : iatom, jatom, ikind, calculate_forces)
93 : END IF
94 :
95 0 : END SUBROUTINE lri_int
96 :
97 : ! **************************************************************************************************
98 : !> \brief calcuates the lri integrals using solid harmonic Gaussians
99 : !> \param lri_env ...
100 : !> \param lrii ...
101 : !> \param lriint ...
102 : !> \param rab distance vector
103 : !> \param obasa orb basis on A
104 : !> \param obasb orb basis on B
105 : !> \param fbasa aux basis on A
106 : !> \param fbasb aux basis on B
107 : !> \param iatom index atom A
108 : !> \param jatom index atom B
109 : !> \param ikind kind atom A
110 : !> \param jkind kind atom B
111 : ! **************************************************************************************************
112 9414 : SUBROUTINE lri_int2(lri_env, lrii, lriint, rab, obasa, obasb, fbasa, fbasb, iatom, jatom, ikind, jkind)
113 :
114 : TYPE(lri_environment_type), POINTER :: lri_env
115 : TYPE(lri_int_type), POINTER :: lrii
116 : TYPE(int_type) :: lriint
117 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
118 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
119 : INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
120 :
121 : INTEGER :: nba, nfa
122 9414 : INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, &
123 9414 : obb_index
124 : REAL(KIND=dp) :: dab
125 9414 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dummy1, Waux_mat
126 9414 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dummy2, dWaux_mat
127 9414 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb
128 9414 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
129 :
130 37656 : dab = SQRT(SUM(rab*rab))
131 9414 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
132 135 : CPASSERT(ALLOCATED(lriint%sabint))
133 135 : CPASSERT(ALLOCATED(lriint%sooint))
134 135 : CPASSERT(ALLOCATED(lriint%abaint))
135 135 : nba = obasa%nsgf
136 135 : nfa = fbasa%nsgf
137 1494117 : lriint%sabint(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
138 10809 : lriint%sooint(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
139 1070984 : lriint%abaint(1:nba, 1:nba, 1:nfa) = lri_env%bas_prop(ikind)%ovlp3
140 : ELSE
141 9279 : IF (lri_env%use_shg_integrals) THEN
142 8593 : scon_oba => lri_env%bas_prop(ikind)%scon_orb
143 8593 : scon_obb => lri_env%bas_prop(jkind)%scon_orb
144 8593 : scon_fba => lri_env%bas_prop(ikind)%scon_ri
145 8593 : scon_fbb => lri_env%bas_prop(jkind)%scon_ri
146 8593 : scona_mix => lri_env%bas_prop(ikind)%scon_mix
147 8593 : sconb_mix => lri_env%bas_prop(jkind)%scon_mix
148 8593 : oba_index => lri_env%bas_prop(ikind)%orb_index
149 8593 : fba_index => lri_env%bas_prop(ikind)%ri_index
150 8593 : obb_index => lri_env%bas_prop(jkind)%orb_index
151 8593 : fbb_index => lri_env%bas_prop(jkind)%ri_index
152 : !
153 : CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, Waux_mat, dWaux_mat, &
154 8593 : .FALSE.)
155 : !*** (fa,fb)
156 8593 : IF (ALLOCATED(lriint%sabint)) THEN
157 : CALL int_overlap_ab_shg_low(lriint%sabint, dummy1, rab, fbasa, fbasb, scon_fba, scon_fbb, &
158 1305 : Waux_mat, dWaux_mat, .TRUE., .FALSE., contraction_high=.FALSE.)
159 : END IF
160 : !*** (a,b)
161 8593 : IF (ALLOCATED(lriint%sooint)) THEN
162 : CALL int_overlap_ab_shg_low(lriint%sooint, dummy1, rab, obasa, obasb, scon_oba, scon_obb, &
163 8593 : Waux_mat, dWaux_mat, .TRUE., .FALSE., contraction_high=.TRUE.)
164 : END IF
165 : !*** (a,b,fa)
166 8593 : IF (ALLOCATED(lriint%abaint)) THEN
167 : CALL int_overlap_aba_shg_low(lriint%abaint, dummy2, rab, obasa, obasb, fbasa, &
168 : scon_obb, scona_mix, oba_index, fba_index, &
169 : lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
170 : lri_env%cg_shg%ncg_none0, &
171 8593 : Waux_mat, dWaux_mat, .TRUE., .FALSE.)
172 : END IF
173 : !*** (a,b,fb)
174 8593 : IF (ALLOCATED(lriint%abbint)) THEN
175 8593 : IF (ikind == jkind) THEN
176 8043 : CPASSERT(ALLOCATED(lriint%abaint))
177 : CALL get_abb_same_kind(lriint%abbint, dummy2, lriint%abaint, dummy2, &
178 8043 : rab, obasa, fbasa, .TRUE., .FALSE.)
179 : ELSE
180 : CALL int_overlap_abb_shg_low(lriint%abbint, dummy2, rab, obasa, obasb, fbasb, &
181 : scon_oba, sconb_mix, obb_index, fbb_index, &
182 : lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
183 : lri_env%cg_shg%ncg_none0, &
184 550 : Waux_mat, dWaux_mat, .TRUE., .FALSE.)
185 : END IF
186 : END IF
187 8593 : DEALLOCATE (Waux_mat, dWaux_mat)
188 : ELSE
189 : !*** (fa,fb)
190 686 : IF (ALLOCATED(lriint%sabint)) THEN
191 : CALL int_overlap_ab_os(lriint%sabint, dummy1, rab, fbasa, fbasb, &
192 686 : .FALSE., lri_env%debug, lrii%dmax_ab)
193 : END IF
194 : !*** (a,b)
195 686 : IF (ALLOCATED(lriint%sooint)) THEN
196 : CALL int_overlap_ab_os(lriint%sooint, dummy1, rab, obasa, obasb, &
197 686 : .FALSE., lri_env%debug, lrii%dmax_oo)
198 : END IF
199 : !*** (a,b,fa)
200 686 : IF (ALLOCATED(lriint%abaint)) THEN
201 : CALL int_overlap_aba_os(lriint%abaint, dummy2, rab, obasa, obasb, fbasa, &
202 686 : .FALSE., lri_env%debug, lrii%dmax_aba)
203 : END IF
204 : !*** (a,b,fb)
205 686 : IF (ALLOCATED(lriint%abbint)) THEN
206 : CALL int_overlap_abb_os(lriint%abbint, dummy2, rab, obasa, obasb, fbasb, &
207 686 : .FALSE., lri_env%debug, lrii%dmax_abb)
208 : END IF
209 : END IF
210 : END IF
211 :
212 18828 : END SUBROUTINE lri_int2
213 :
214 : ! **************************************************************************************************
215 :
216 : ! **************************************************************************************************
217 : !> \brief ...
218 : !> \param lri_env ...
219 : !> \param lrii ...
220 : !> \param rab ...
221 : !> \param obasa ...
222 : !> \param obasb ...
223 : !> \param fbasa ...
224 : !> \param fbasb ...
225 : !> \param iatom ...
226 : !> \param jatom ...
227 : !> \param ikind ...
228 : !> \param jkind ...
229 : !> \param calculate_forces ...
230 : ! **************************************************************************************************
231 0 : SUBROUTINE lri_dint(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
232 : iatom, jatom, ikind, jkind, calculate_forces)
233 :
234 : TYPE(lri_environment_type), POINTER :: lri_env
235 : TYPE(lri_int_type), POINTER :: lrii
236 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
237 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
238 : INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
239 : LOGICAL, INTENT(IN) :: calculate_forces
240 :
241 0 : IF (lri_env%use_shg_integrals) THEN
242 : CALL lri_int_shg(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
243 0 : iatom, jatom, ikind, jkind, calculate_forces)
244 : ELSE
245 : CALL lri_int_os(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
246 0 : iatom, jatom, ikind, calculate_forces)
247 : END IF
248 :
249 0 : END SUBROUTINE lri_dint
250 :
251 : ! **************************************************************************************************
252 : !> \brief ...
253 : !> \param lri_env ...
254 : !> \param lrii ...
255 : !> \param lridint ...
256 : !> \param rab ...
257 : !> \param obasa ...
258 : !> \param obasb ...
259 : !> \param fbasa ...
260 : !> \param fbasb ...
261 : !> \param iatom ...
262 : !> \param jatom ...
263 : !> \param ikind ...
264 : !> \param jkind ...
265 : ! **************************************************************************************************
266 6586 : SUBROUTINE lri_dint2(lri_env, lrii, lridint, rab, obasa, obasb, fbasa, fbasb, &
267 : iatom, jatom, ikind, jkind)
268 :
269 : TYPE(lri_environment_type), POINTER :: lri_env
270 : TYPE(lri_int_type), POINTER :: lrii
271 : TYPE(dint_type) :: lridint
272 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
273 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
274 : INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
275 :
276 : INTEGER :: nba, nbb, nfa, nfb
277 6586 : INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, &
278 6586 : obb_index
279 : REAL(KIND=dp) :: dab
280 6586 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dummy1
281 6586 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dummy2, Waux_mat
282 6586 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
283 6586 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb
284 6586 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
285 :
286 26344 : dab = SQRT(SUM(rab*rab))
287 6586 : nba = obasa%nsgf
288 6586 : nbb = obasb%nsgf
289 6586 : nfa = fbasa%nsgf
290 6586 : nfb = fbasb%nsgf
291 :
292 6586 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
293 : ! nothing to do
294 : ELSE
295 6586 : IF (lrii%calc_force_pair) THEN
296 6586 : IF (lri_env%use_shg_integrals) THEN
297 6140 : scon_oba => lri_env%bas_prop(ikind)%scon_orb
298 6140 : scon_obb => lri_env%bas_prop(jkind)%scon_orb
299 6140 : scon_fba => lri_env%bas_prop(ikind)%scon_ri
300 6140 : scon_fbb => lri_env%bas_prop(jkind)%scon_ri
301 6140 : scona_mix => lri_env%bas_prop(ikind)%scon_mix
302 6140 : sconb_mix => lri_env%bas_prop(jkind)%scon_mix
303 6140 : oba_index => lri_env%bas_prop(ikind)%orb_index
304 6140 : fba_index => lri_env%bas_prop(ikind)%ri_index
305 6140 : obb_index => lri_env%bas_prop(jkind)%orb_index
306 6140 : fbb_index => lri_env%bas_prop(jkind)%ri_index
307 :
308 : CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, &
309 6140 : Waux_mat, dWaux_mat, .TRUE.)
310 : !*** (fa,fb)
311 6140 : IF (ALLOCATED(lridint%dsabint)) THEN
312 : CALL int_overlap_ab_shg_low(dummy1, lridint%dsabint, rab, fbasa, fbasb, scon_fba, scon_fbb, &
313 764 : Waux_mat, dWaux_mat, .FALSE., .TRUE., contraction_high=.FALSE.)
314 : END IF
315 : !*** (a,b)
316 6140 : IF (ALLOCATED(lridint%dsooint)) THEN
317 : CALL int_overlap_ab_shg_low(dummy1, lridint%dsooint, rab, obasa, obasb, scon_oba, scon_obb, &
318 6140 : Waux_mat, dWaux_mat, .FALSE., .TRUE., contraction_high=.TRUE.)
319 : END IF
320 : !*** (a,b,fa)
321 6140 : IF (ALLOCATED(lridint%dabaint)) THEN
322 : CALL int_overlap_aba_shg_low(dummy2, lridint%dabaint, rab, obasa, obasb, fbasa, &
323 : scon_obb, scona_mix, oba_index, fba_index, &
324 : lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
325 : lri_env%cg_shg%ncg_none0, &
326 6140 : Waux_mat, dWaux_mat, .FALSE., .TRUE.)
327 : END IF
328 : !*** (a,b,fb)
329 6140 : IF (ALLOCATED(lridint%dabbint)) THEN
330 6140 : IF (ikind == jkind) THEN
331 5962 : CPASSERT(ALLOCATED(lridint%dabaint))
332 : CALL get_abb_same_kind(dummy2, lridint%dabbint, dummy2, lridint%dabaint, &
333 5962 : rab, obasa, fbasa, .FALSE., .TRUE.)
334 : ELSE
335 : CALL int_overlap_abb_shg_low(dummy2, lridint%dabbint, rab, obasa, obasb, fbasb, &
336 : scon_oba, sconb_mix, obb_index, fbb_index, &
337 : lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
338 : lri_env%cg_shg%ncg_none0, &
339 178 : Waux_mat, dWaux_mat, .FALSE., .TRUE.)
340 : END IF
341 : END IF
342 6140 : DEALLOCATE (Waux_mat, dWaux_mat)
343 :
344 : ELSE
345 :
346 : !*** (fa,fb)
347 446 : IF (ALLOCATED(lridint%dsabint)) THEN
348 1784 : ALLOCATE (dummy1(nfa, nfb))
349 : CALL int_overlap_ab_os(dummy1, lridint%dsabint, rab, fbasa, fbasb, &
350 446 : .TRUE., lri_env%debug, lrii%dmax_ab)
351 446 : DEALLOCATE (dummy1)
352 : END IF
353 : !*** (a,b)
354 446 : IF (ALLOCATED(lridint%dsooint)) THEN
355 1784 : ALLOCATE (dummy1(nba, nbb))
356 : CALL int_overlap_ab_os(dummy1, lridint%dsooint, rab, obasa, obasb, &
357 446 : .TRUE., lri_env%debug, lrii%dmax_oo)
358 446 : DEALLOCATE (dummy1)
359 : END IF
360 : !*** (a,b,fa)
361 446 : IF (ALLOCATED(lridint%dabaint)) THEN
362 2230 : ALLOCATE (dummy2(nba, nbb, nfa))
363 : CALL int_overlap_aba_os(dummy2, lridint%dabaint, rab, obasa, obasb, fbasa, &
364 446 : .TRUE., lri_env%debug, lrii%dmax_aba)
365 446 : DEALLOCATE (dummy2)
366 : END IF
367 : !*** (a,b,fb)
368 446 : IF (ALLOCATED(lridint%dabbint)) THEN
369 2230 : ALLOCATE (dummy2(nba, nbb, nfb))
370 : CALL int_overlap_abb_os(dummy2, lridint%dabbint, rab, obasa, obasb, fbasb, &
371 446 : .TRUE., lri_env%debug, lrii%dmax_abb)
372 446 : DEALLOCATE (dummy2)
373 : END IF
374 : END IF
375 : END IF
376 : END IF
377 :
378 13172 : END SUBROUTINE lri_dint2
379 :
380 : ! **************************************************************************************************
381 : !> \brief calculates the lri intergrals according to the Obara-Saika (OS)
382 : !> scheme
383 : !> \param lri_env ...
384 : !> \param lrii ...
385 : !> \param rab distance vector
386 : !> \param obasa orb basis on center A
387 : !> \param obasb orb basis on center B
388 : !> \param fbasa aux basis on center A
389 : !> \param fbasb aux basis on center B
390 : !> \param iatom index atom A
391 : !> \param jatom index atom B
392 : !> \param ikind kind atom A
393 : !> \param calculate_forces ...
394 : ! **************************************************************************************************
395 0 : SUBROUTINE lri_int_os(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
396 : iatom, jatom, ikind, calculate_forces)
397 :
398 : TYPE(lri_environment_type), POINTER :: lri_env
399 : TYPE(lri_int_type), POINTER :: lrii
400 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
401 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
402 : INTEGER, INTENT(IN) :: iatom, jatom, ikind
403 : LOGICAL, INTENT(IN) :: calculate_forces
404 :
405 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_os'
406 :
407 : INTEGER :: handle, nba, nbb, nfa, nfb
408 : LOGICAL :: do_force
409 : REAL(KIND=dp) :: dab
410 :
411 0 : CALL timeset(routineN, handle)
412 :
413 0 : dab = SQRT(SUM(rab*rab))
414 0 : nba = obasa%nsgf
415 0 : nbb = obasb%nsgf
416 0 : nfa = fbasa%nsgf
417 0 : nfb = fbasb%nsgf
418 :
419 : !*** calculate overlap integrals; for iatom=jatom this is the self-overlap
420 0 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
421 : !*** (fa,fa)
422 0 : lrii%sab(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
423 : !*** (a,a)
424 0 : lrii%soo(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
425 : !*** (a,a,fa)
426 : CALL int_overlap_aba_os(lrii%abaint, rab=rab, oba=obasa, obb=obasb, &
427 : fba=fbasa, calculate_forces=.FALSE., debug=lri_env%debug, &
428 0 : dmax=lrii%dmax_aba)
429 0 : IF (calculate_forces) THEN
430 0 : lrii%dsab = 0._dp
431 0 : lrii%dsoo = 0._dp
432 0 : lrii%dabdaint = 0.0_dp
433 0 : lrii%dabbint = 0.0_dp
434 : END IF
435 : ELSE
436 0 : IF (calculate_forces) THEN
437 0 : do_force = lrii%calc_force_pair
438 : ELSE
439 0 : do_force = .FALSE.
440 : END IF
441 : !*** (fa,fb)
442 : CALL int_overlap_ab_os(lrii%sab, lrii%dsab, rab, fbasa, fbasb, &
443 0 : do_force, lri_env%debug, lrii%dmax_ab)
444 : !*** (a,b)
445 : CALL int_overlap_ab_os(lrii%soo, lrii%dsoo, rab, obasa, obasb, &
446 0 : do_force, lri_env%debug, lrii%dmax_oo)
447 : !*** (a,b,fa)
448 : CALL int_overlap_aba_os(lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, &
449 0 : do_force, lri_env%debug, lrii%dmax_aba)
450 : !*** (a,b,fb)
451 : CALL int_overlap_abb_os(lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, &
452 0 : do_force, lri_env%debug, lrii%dmax_abb)
453 : END IF
454 :
455 0 : CALL timestop(handle)
456 :
457 0 : END SUBROUTINE lri_int_os
458 :
459 : ! **************************************************************************************************
460 : !> \brief calcuates the lri integrals using solid harmonic Gaussians
461 : !> \param lri_env ...
462 : !> \param lrii ...
463 : !> \param rab distance vector
464 : !> \param obasa orb basis on A
465 : !> \param obasb orb basis on B
466 : !> \param fbasa aux basis on A
467 : !> \param fbasb aux basis on B
468 : !> \param iatom index atom A
469 : !> \param jatom index atom B
470 : !> \param ikind kind atom A
471 : !> \param jkind kind atom B
472 : !> \param calculate_forces ...
473 : ! **************************************************************************************************
474 0 : SUBROUTINE lri_int_shg(lri_env, lrii, rab, obasa, obasb, fbasa, fbasb, &
475 : iatom, jatom, ikind, jkind, calculate_forces)
476 :
477 : TYPE(lri_environment_type), POINTER :: lri_env
478 : TYPE(lri_int_type), POINTER :: lrii
479 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab
480 : TYPE(gto_basis_set_type), POINTER :: obasa, obasb, fbasa, fbasb
481 : INTEGER, INTENT(IN) :: iatom, jatom, ikind, jkind
482 : LOGICAL, INTENT(IN) :: calculate_forces
483 :
484 : CHARACTER(LEN=*), PARAMETER :: routineN = 'lri_int_shg'
485 :
486 : INTEGER :: handle, nba, nbb, nfa, nfb
487 0 : INTEGER, DIMENSION(:, :, :), POINTER :: fba_index, fbb_index, oba_index, &
488 0 : obb_index
489 : LOGICAL :: do_force
490 : REAL(KIND=dp) :: dab
491 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Waux_mat
492 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dWaux_mat
493 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_fba, scon_fbb, scon_oba, scon_obb
494 0 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scona_mix, sconb_mix
495 :
496 0 : CALL timeset(routineN, handle)
497 0 : NULLIFY (scon_oba, scon_obb, scon_fba, scon_fbb, scona_mix, sconb_mix, &
498 0 : oba_index, obb_index, fba_index, fbb_index)
499 0 : dab = SQRT(SUM(rab*rab))
500 0 : nba = obasa%nsgf
501 0 : nbb = obasb%nsgf
502 0 : nfa = fbasa%nsgf
503 0 : nfb = fbasb%nsgf
504 :
505 : !*** calculate overlap integrals; for iatom=jatom this is the self-overlap
506 0 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
507 : !*** (fa,fa)
508 0 : lrii%sab(1:nfa, 1:nfa) = lri_env%bas_prop(ikind)%ri_ovlp(1:nfa, 1:nfa)
509 : !*** (a,a)
510 0 : lrii%soo(1:nba, 1:nba) = lri_env%bas_prop(ikind)%orb_ovlp(1:nba, 1:nba)
511 : !*** (a,a,fa)
512 0 : lrii%abaint(1:nba, 1:nba, 1:nfa) = lri_env%bas_prop(ikind)%ovlp3
513 0 : IF (calculate_forces) THEN
514 0 : lrii%dsab = 0._dp
515 0 : lrii%dsoo = 0._dp
516 0 : lrii%dabdaint = 0.0_dp
517 0 : lrii%dabbint = 0.0_dp
518 : END IF
519 : ELSE
520 0 : IF (calculate_forces) THEN
521 0 : do_force = lrii%calc_force_pair
522 : ELSE
523 0 : do_force = .FALSE.
524 : END IF
525 0 : scon_oba => lri_env%bas_prop(ikind)%scon_orb
526 0 : scon_obb => lri_env%bas_prop(jkind)%scon_orb
527 0 : scon_fba => lri_env%bas_prop(ikind)%scon_ri
528 0 : scon_fbb => lri_env%bas_prop(jkind)%scon_ri
529 0 : scona_mix => lri_env%bas_prop(ikind)%scon_mix
530 0 : sconb_mix => lri_env%bas_prop(jkind)%scon_mix
531 0 : oba_index => lri_env%bas_prop(ikind)%orb_index
532 0 : fba_index => lri_env%bas_prop(ikind)%ri_index
533 0 : obb_index => lri_env%bas_prop(jkind)%orb_index
534 0 : fbb_index => lri_env%bas_prop(jkind)%ri_index
535 : CALL lri_precalc_angular_shg_part(obasa, obasb, fbasa, fbasb, rab, Waux_mat, dWaux_mat, &
536 : do_force)
537 : !*** (fa,fb)
538 : CALL int_overlap_ab_shg_low(lrii%sab, lrii%dsab, rab, fbasa, fbasb, scon_fba, scon_fbb, &
539 0 : Waux_mat, dWaux_mat, .TRUE., do_force, contraction_high=.FALSE.)
540 : !*** (a,b)
541 : CALL int_overlap_ab_shg_low(lrii%soo, lrii%dsoo, rab, obasa, obasb, scon_oba, scon_obb, &
542 0 : Waux_mat, dWaux_mat, .TRUE., do_force, contraction_high=.TRUE.)
543 : !*** (a,b,fa)
544 : CALL int_overlap_aba_shg_low(lrii%abaint, lrii%dabdaint, rab, obasa, obasb, fbasa, &
545 : scon_obb, scona_mix, oba_index, fba_index, &
546 : lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
547 : lri_env%cg_shg%ncg_none0, &
548 0 : Waux_mat, dWaux_mat, .TRUE., do_force)
549 : !*** (a,b,fb)
550 0 : IF (ikind == jkind) THEN
551 : CALL get_abb_same_kind(lrii%abbint, lrii%dabbint, lrii%abaint, lrii%dabdaint, &
552 0 : rab, obasa, fbasa, .TRUE., do_force)
553 : ELSE
554 : CALL int_overlap_abb_shg_low(lrii%abbint, lrii%dabbint, rab, obasa, obasb, fbasb, &
555 : scon_oba, sconb_mix, obb_index, fbb_index, &
556 : lri_env%cg_shg%cg_coeff, lri_env%cg_shg%cg_none0_list, &
557 : lri_env%cg_shg%ncg_none0, &
558 0 : Waux_mat, dWaux_mat, .TRUE., do_force)
559 : END IF
560 :
561 0 : DEALLOCATE (Waux_mat, dWaux_mat)
562 : END IF
563 :
564 0 : CALL timestop(handle)
565 :
566 0 : END SUBROUTINE lri_int_shg
567 :
568 : ! **************************************************************************************************
569 : !> \brief ...
570 : !> \param lriint ...
571 : !> \param lridint ...
572 : !> \param nba ...
573 : !> \param nbb ...
574 : !> \param nfa ...
575 : !> \param nfb ...
576 : !> \param skip_sab ...
577 : !> \param skip_soo ...
578 : !> \param skip_aba ...
579 : !> \param skip_abb ...
580 : !> \param skip_dsab ...
581 : !> \param skip_dsoo ...
582 : !> \param skip_daba ...
583 : !> \param skip_dabb ...
584 : ! **************************************************************************************************
585 16000 : SUBROUTINE allocate_int_type(lriint, lridint, nba, nbb, nfa, nfb, &
586 : skip_sab, skip_soo, skip_aba, skip_abb, &
587 : skip_dsab, skip_dsoo, skip_daba, skip_dabb)
588 :
589 : TYPE(int_type), INTENT(INOUT), OPTIONAL :: lriint
590 : TYPE(dint_type), INTENT(INOUT), OPTIONAL :: lridint
591 : INTEGER, INTENT(IN) :: nba, nbb, nfa, nfb
592 : LOGICAL, INTENT(IN), OPTIONAL :: skip_sab, skip_soo, skip_aba, skip_abb, &
593 : skip_dsab, skip_dsoo, skip_daba, &
594 : skip_dabb
595 :
596 : LOGICAL :: do_aba, do_abb, do_daba, do_dabb, &
597 : do_dsab, do_dsoo, do_sab, do_soo
598 :
599 16000 : IF (PRESENT(lriint)) THEN
600 9414 : do_sab = .TRUE.
601 9414 : IF (PRESENT(skip_sab)) do_sab = .NOT. skip_sab
602 9414 : do_soo = .TRUE.
603 9414 : IF (PRESENT(skip_soo)) do_soo = .NOT. skip_soo
604 9414 : do_aba = .TRUE.
605 9414 : IF (PRESENT(skip_aba)) do_aba = .NOT. skip_aba
606 9414 : do_abb = .TRUE.
607 9414 : IF (PRESENT(skip_abb)) do_abb = .NOT. skip_abb
608 : !
609 9414 : IF (do_sab) THEN
610 2126 : IF (ALLOCATED(lriint%sabint)) DEALLOCATE (lriint%sabint)
611 8504 : ALLOCATE (lriint%sabint(nfa, nfb))
612 11684512 : lriint%sabint = 0.0_dp
613 : END IF
614 9414 : IF (do_soo) THEN
615 9414 : IF (ALLOCATED(lriint%sooint)) DEALLOCATE (lriint%sooint)
616 37656 : ALLOCATE (lriint%sooint(nba, nbb))
617 1480901 : lriint%sooint = 0.0_dp
618 : END IF
619 9414 : IF (do_aba) THEN
620 9414 : IF (ALLOCATED(lriint%abaint)) DEALLOCATE (lriint%abaint)
621 47070 : ALLOCATE (lriint%abaint(nba, nbb, nfa))
622 114627218 : lriint%abaint = 0.0_dp
623 : END IF
624 9414 : IF (do_abb) THEN
625 9414 : IF (ALLOCATED(lriint%abbint)) DEALLOCATE (lriint%abbint)
626 47070 : ALLOCATE (lriint%abbint(nba, nbb, nfb))
627 114858301 : lriint%abbint = 0.0_dp
628 : END IF
629 : END IF
630 : !
631 16000 : IF (PRESENT(lridint)) THEN
632 6586 : do_dsab = .TRUE.
633 6586 : IF (PRESENT(skip_dsab)) do_dsab = .NOT. skip_dsab
634 6586 : do_dsoo = .TRUE.
635 6586 : IF (PRESENT(skip_dsoo)) do_dsoo = .NOT. skip_dsoo
636 6586 : do_daba = .TRUE.
637 6586 : IF (PRESENT(skip_daba)) do_daba = .NOT. skip_daba
638 6586 : do_dabb = .TRUE.
639 6586 : IF (PRESENT(skip_dabb)) do_dabb = .NOT. skip_dabb
640 : !
641 6586 : IF (do_dsab) THEN
642 1210 : IF (ALLOCATED(lridint%dsabint)) DEALLOCATE (lridint%dsabint)
643 6050 : ALLOCATE (lridint%dsabint(nfa, nfb, 3))
644 14103661 : lridint%dsabint = 0.0_dp
645 : END IF
646 6586 : IF (do_dsoo) THEN
647 6586 : IF (ALLOCATED(lridint%dsooint)) DEALLOCATE (lridint%dsooint)
648 32930 : ALLOCATE (lridint%dsooint(nba, nbb, 3))
649 3410512 : lridint%dsooint = 0.0_dp
650 : END IF
651 6586 : IF (do_daba) THEN
652 6586 : IF (ALLOCATED(lridint%dabaint)) DEALLOCATE (lridint%dabaint)
653 39516 : ALLOCATE (lridint%dabaint(nba, nbb, nfa, 3))
654 254242525 : lridint%dabaint = 0.0_dp
655 : END IF
656 6586 : IF (do_dabb) THEN
657 6586 : IF (ALLOCATED(lridint%dabbint)) DEALLOCATE (lridint%dabbint)
658 39516 : ALLOCATE (lridint%dabbint(nba, nbb, nfb, 3))
659 254272045 : lridint%dabbint = 0.0_dp
660 : END IF
661 : END IF
662 :
663 16000 : END SUBROUTINE allocate_int_type
664 :
665 : ! **************************************************************************************************
666 : !> \brief ...
667 : !> \param lriint ...
668 : !> \param lridint ...
669 : ! **************************************************************************************************
670 16000 : SUBROUTINE deallocate_int_type(lriint, lridint)
671 :
672 : TYPE(int_type), INTENT(INOUT), OPTIONAL :: lriint
673 : TYPE(dint_type), INTENT(INOUT), OPTIONAL :: lridint
674 :
675 16000 : IF (PRESENT(lriint)) THEN
676 9414 : IF (ALLOCATED(lriint%sabint)) DEALLOCATE (lriint%sabint)
677 9414 : IF (ALLOCATED(lriint%sooint)) DEALLOCATE (lriint%sooint)
678 9414 : IF (ALLOCATED(lriint%abaint)) DEALLOCATE (lriint%abaint)
679 9414 : IF (ALLOCATED(lriint%abbint)) DEALLOCATE (lriint%abbint)
680 : END IF
681 : !
682 16000 : IF (PRESENT(lridint)) THEN
683 6586 : IF (ALLOCATED(lridint%dsabint)) DEALLOCATE (lridint%dsabint)
684 6586 : IF (ALLOCATED(lridint%dsooint)) DEALLOCATE (lridint%dsooint)
685 6586 : IF (ALLOCATED(lridint%dabaint)) DEALLOCATE (lridint%dabaint)
686 6586 : IF (ALLOCATED(lridint%dabbint)) DEALLOCATE (lridint%dabbint)
687 : END IF
688 :
689 16000 : END SUBROUTINE deallocate_int_type
690 :
691 0 : END MODULE lri_integrals
|
