18 USE dbcsr_api,
ONLY: dbcsr_get_block_p,&
29 #include "./base/base_uses.f90"
35 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_o3c_methods'
48 TYPE(o3c_container_type),
POINTER :: o3c
49 LOGICAL,
INTENT(IN),
OPTIONAL :: calculate_forces
50 TYPE(dbcsr_p_type),
DIMENSION(:),
OPTIONAL, &
53 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calculate_o3c_integrals'
55 INTEGER :: egfa, egfb, egfc, handle, i, iatom, icol, ikind, irow, iset, ispin, j, jatom, &
56 jkind, jset, katom, kkind, kset, mepos, ncoa, ncob, ncoc, ni, nj, nk, nseta, nsetb, &
57 nsetc, nspin, nthread, sgfa, sgfb, sgfc
58 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, lc_max, &
59 lc_min, npgfa, npgfb, npgfc, nsgfa, &
61 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb, first_sgfc
62 LOGICAL :: do_force, found, trans
63 REAL(kind=
dp) :: dij, dik, djk, fpre
64 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: pmat
65 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: sabc, sabc_contr
66 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: iabdc, iadbc, idabc, sabdc, sdabc
67 REAL(kind=
dp),
DIMENSION(3) :: rij, rik, rjk
68 REAL(kind=
dp),
DIMENSION(:),
POINTER :: set_radius_a, set_radius_b, set_radius_c
69 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: fi, fj, fk, pblock, rpgfa, rpgfb, rpgfc, &
70 sphi_a, sphi_b, sphi_c, tvec, zeta, &
72 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: iabc
73 TYPE(gto_basis_set_p_type),
DIMENSION(:),
POINTER :: basis_set_list_a, basis_set_list_b, &
75 TYPE(gto_basis_set_type),
POINTER :: basis_set_a, basis_set_b, basis_set_c
76 TYPE(o3c_iterator_type) :: o3c_iterator
78 CALL timeset(routinen, handle)
81 IF (
PRESENT(calculate_forces)) do_force = calculate_forces
86 basis_set_list_b=basis_set_list_b, basis_set_list_c=basis_set_list_c)
106 DO WHILE (
o3c_iterate(o3c_iterator, mepos=mepos) == 0)
108 ikind=ikind, jkind=jkind, kkind=kkind, rij=rij, rik=rik, &
109 integral=iabc, tvec=tvec, force_i=fi, force_j=fj, force_k=fk)
110 cpassert(.NOT.
ASSOCIATED(iabc))
111 cpassert(.NOT.
ASSOCIATED(tvec))
112 cpassert(.NOT.
ASSOCIATED(fi))
113 cpassert(.NOT.
ASSOCIATED(fj))
114 cpassert(.NOT.
ASSOCIATED(fk))
116 basis_set_a => basis_set_list_a(ikind)%gto_basis_set
117 basis_set_b => basis_set_list_b(jkind)%gto_basis_set
118 basis_set_c => basis_set_list_c(kkind)%gto_basis_set
120 first_sgfa => basis_set_a%first_sgf
121 la_max => basis_set_a%lmax
122 la_min => basis_set_a%lmin
123 npgfa => basis_set_a%npgf
124 nseta = basis_set_a%nset
125 nsgfa => basis_set_a%nsgf_set
126 rpgfa => basis_set_a%pgf_radius
127 set_radius_a => basis_set_a%set_radius
128 sphi_a => basis_set_a%sphi
129 zeta => basis_set_a%zet
131 first_sgfb => basis_set_b%first_sgf
132 lb_max => basis_set_b%lmax
133 lb_min => basis_set_b%lmin
134 npgfb => basis_set_b%npgf
135 nsetb = basis_set_b%nset
136 nsgfb => basis_set_b%nsgf_set
137 rpgfb => basis_set_b%pgf_radius
138 set_radius_b => basis_set_b%set_radius
139 sphi_b => basis_set_b%sphi
140 zetb => basis_set_b%zet
142 first_sgfc => basis_set_c%first_sgf
143 lc_max => basis_set_c%lmax
144 lc_min => basis_set_c%lmin
145 npgfc => basis_set_c%npgf
146 nsetc = basis_set_c%nset
147 nsgfc => basis_set_c%nsgf_set
148 rpgfc => basis_set_c%pgf_radius
149 set_radius_c => basis_set_c%set_radius
150 sphi_c => basis_set_c%sphi
151 zetc => basis_set_c%zet
157 ALLOCATE (iabc(ni, nj, nk))
158 iabc(:, :, :) = 0.0_dp
160 ALLOCATE (fi(nk, 3), fj(nk, 3), fk(nk, 3))
164 ALLOCATE (idabc(ni, nj, nk, 3))
165 idabc(:, :, :, :) = 0.0_dp
166 ALLOCATE (iadbc(ni, nj, nk, 3))
167 iadbc(:, :, :, :) = 0.0_dp
168 ALLOCATE (iabdc(ni, nj, nk, 3))
169 iabdc(:, :, :, :) = 0.0_dp
173 ALLOCATE (tvec(nk, nspin))
176 rjk(1:3) = rik(1:3) - rij(1:3)
183 IF (set_radius_a(iset) + set_radius_b(jset) < dij) cycle
185 IF (set_radius_a(iset) + set_radius_c(kset) < dik) cycle
186 IF (set_radius_b(jset) + set_radius_c(kset) < djk) cycle
188 ncoa = npgfa(iset)*
ncoset(la_max(iset))
189 ncob = npgfb(jset)*
ncoset(lb_max(jset))
190 ncoc = npgfc(kset)*
ncoset(lc_max(kset))
192 sgfa = first_sgfa(1, iset)
193 sgfb = first_sgfb(1, jset)
194 sgfc = first_sgfc(1, kset)
196 egfa = sgfa + nsgfa(iset) - 1
197 egfb = sgfb + nsgfb(jset) - 1
198 egfc = sgfc + nsgfc(kset) - 1
200 IF (ncoa*ncob*ncoc > 0)
THEN
201 ALLOCATE (sabc(ncoa, ncob, ncoc))
202 sabc(:, :, :) = 0.0_dp
204 ALLOCATE (sdabc(ncoa, ncob, ncoc, 3))
205 sdabc(:, :, :, :) = 0.0_dp
206 ALLOCATE (sabdc(ncoa, ncob, ncoc, 3))
207 sabdc(:, :, :, :) = 0.0_dp
208 CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
209 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
210 lc_max(kset), npgfc(kset), zetc(:, kset), rpgfc(:, kset), lc_min(kset), &
211 rij, dij, rik, dik, rjk, djk, sabc, sdabc, sabdc)
213 CALL overlap3(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
214 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
215 lc_max(kset), npgfc(kset), zetc(:, kset), rpgfc(:, kset), lc_min(kset), &
216 rij, dij, rik, dik, rjk, djk, sabc)
218 ALLOCATE (sabc_contr(nsgfa(iset), nsgfb(jset), nsgfc(kset)))
221 sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
222 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
223 iabc(sgfa:egfa, sgfb:egfb, sgfc:egfc) = &
224 sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
228 sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
229 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
230 idabc(sgfa:egfa, sgfb:egfb, sgfc:egfc, i) = &
231 sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
233 sphi_a(:, sgfa:), sphi_b(:, sgfb:), sphi_c(:, sgfc:), &
234 ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfc(kset))
235 iabdc(sgfa:egfa, sgfb:egfb, sgfc:egfc, i) = &
236 sabc_contr(1:nsgfa(iset), 1:nsgfb(jset), 1:nsgfc(kset))
240 DEALLOCATE (sabc_contr)
244 DEALLOCATE (sdabc, sabdc)
251 iadbc(:, :, :, :) = -idabc(:, :, :, :) - iabdc(:, :, :, :)
255 iatom=iatom, jatom=jatom, katom=katom)
258 IF (iatom <= jatom)
THEN
267 IF (iatom == jatom)
THEN
272 ALLOCATE (pmat(ni, nj))
275 CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
276 row=irow, col=icol, block=pblock, found=found)
279 pmat(:, :) = pmat(:, :) + transpose(pblock(:, :))
281 pmat(:, :) = pmat(:, :) + pblock(:, :)
287 fi(j, i) = fpre*sum(pmat(:, :)*idabc(:, :, j, i))
288 fj(j, i) = fpre*sum(pmat(:, :)*iadbc(:, :, j, i))
289 fk(j, i) = fpre*sum(pmat(:, :)*iabdc(:, :, j, i))
294 DEALLOCATE (idabc, iadbc, iabdc)
298 integral=iabc, tvec=tvec, force_i=fi, force_j=fj, force_k=fk)
304 CALL timestop(handle)
315 TYPE(o3c_container_type),
POINTER :: o3c
316 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_p
318 CHARACTER(LEN=*),
PARAMETER :: routinen =
'contract12_o3c'
320 INTEGER :: handle, iatom, icol, ik, irow, ispin, &
321 jatom, mepos, nk, nspin, nthread
322 LOGICAL :: found, ijsymmetric, trans
323 REAL(kind=
dp) :: fpre
324 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pblock, tvec
325 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: iabc
326 TYPE(o3c_iterator_type) :: o3c_iterator
328 CALL timeset(routinen, handle)
330 nspin =
SIZE(matrix_p, 1)
332 cpassert(ijsymmetric)
345 DO WHILE (
o3c_iterate(o3c_iterator, mepos=mepos) == 0)
347 integral=iabc, tvec=tvec)
350 IF (iatom <= jatom)
THEN
359 IF (iatom == jatom)
THEN
366 CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
367 row=irow, col=icol, block=pblock, found=found)
371 tvec(ik, ispin) = fpre*sum(transpose(pblock(:, :))*iabc(:, :, ik))
375 tvec(ik, ispin) = fpre*sum(pblock(:, :)*iabc(:, :, ik))
385 CALL timestop(handle)
397 TYPE(o3c_container_type),
POINTER :: o3c
398 TYPE(o3c_vec_type),
DIMENSION(:),
POINTER :: vec
399 TYPE(dbcsr_type) :: matrix
401 CHARACTER(LEN=*),
PARAMETER :: routinen =
'contract3_o3c'
403 INTEGER :: handle, iatom, icol, ik, irow, jatom, &
404 katom, mepos, nk, nthread, s1, s2
405 LOGICAL :: found, ijsymmetric, trans
406 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: work
407 REAL(kind=
dp),
DIMENSION(:),
POINTER :: v
408 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: pblock
409 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: iabc
410 TYPE(o3c_iterator_type) :: o3c_iterator
412 CALL timeset(routinen, handle)
415 cpassert(ijsymmetric)
428 DO WHILE (
o3c_iterate(o3c_iterator, mepos=mepos) == 0)
435 IF (iatom <= jatom)
THEN
445 CALL dbcsr_get_block_p(matrix=matrix, row=irow, col=icol, block=pblock, found=found)
448 s1 =
SIZE(pblock, 1); s2 =
SIZE(pblock, 2)
449 ALLOCATE (work(s1, s2))
454 CALL daxpy(s1*s2, v(ik), transpose(iabc(:, :, ik)), 1, work(:, :), 1)
458 CALL daxpy(s1*s2, v(ik), iabc(:, :, ik), 1, work(:, :), 1)
467 CALL dbcsr_get_block_p(matrix=matrix, row=irow, col=icol, block=pblock, found=found)
468 CALL daxpy(s1*s2, 1.0_dp, work(:, :), 1, pblock(:, :), 1)
478 CALL timestop(handle)
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
subroutine, public abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
contract three-center overlap integrals (a,b,c) and transfer to spherical Gaussians
subroutine, public overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max_set, npgfb, zetb, rpgfb, lb_min_set, lc_max_set, npgfc, zetc, rpgfc, lc_min_set, rab, dab, rac, dac, rbc, dbc, sabc, sdabc, sabdc, int_abc_ext)
Calculation of three-center overlap integrals [a|b|c] over primitive Cartesian Gaussian functions.
Defines the basic variable types.
integer, parameter, public dp
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
Methods used with 3-center overlap type integrals containers.
subroutine, public contract3_o3c(o3c, vec, matrix)
Contraction of 3-tensor over index 3 h(ij) = h(ij) + sum_k (ijk)*v(k)
subroutine, public calculate_o3c_integrals(o3c, calculate_forces, matrix_p)
...
subroutine, public contract12_o3c(o3c, matrix_p)
Contraction of 3-tensor over indices 1 and 2 (assuming symmetry) t(k) = sum_ij (ijk)*p(ij)
3-center overlap type integrals containers
subroutine, public set_o3c_container(o3c_iterator, mepos, integral, tvec, force_i, force_j, force_k)
...
subroutine, public get_o3c_iterator_info(o3c_iterator, mepos, iatom, jatom, katom, ikind, jkind, kkind, rij, rik, cellj, cellk, integral, tvec, force_i, force_j, force_k)
...
subroutine, public o3c_iterator_create(o3c, o3c_iterator, nthread)
...
subroutine, public o3c_iterator_release(o3c_iterator)
...
subroutine, public get_o3c_vec(o3c_vec, i, vec, n)
...
integer function, public o3c_iterate(o3c_iterator, mepos)
...
subroutine, public get_o3c_container(o3c, ijsymmetric, nspin, nijpairs, ijpair, basis_set_list_a, basis_set_list_b, basis_set_list_c, sab_nl, sac_nl)
...