49 LOGICAL,
INTENT(IN),
OPTIONAL :: calculate_forces
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
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))
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)