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 Analytical derivatives of Integrals for semi-empirical methods
10 : !> \author Teodoro Laino - Zurich University 04.2007 [tlaino]
11 : !> \par History
12 : !> 23.11.2007 jhu short range version of integrals
13 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
14 : !> for computing integrals
15 : !> Teodoro Laino (05.2008) [tlaino] - University of Zurich : analytical
16 : !> derivatives for d-orbitals
17 : ! **************************************************************************************************
18 : MODULE semi_empirical_int_ana
19 :
20 : USE input_constants, ONLY: do_method_am1, &
21 : do_method_pchg, &
22 : do_method_pdg, &
23 : do_method_pm3, &
24 : do_method_pm6, &
25 : do_method_pm6fm, &
26 : do_method_undef, &
27 : do_se_IS_kdso_d
28 : USE kinds, ONLY: dp
29 : USE multipole_types, ONLY: do_multipole_none
30 : USE physcon, ONLY: angstrom, &
31 : evolt
32 : USE semi_empirical_int_arrays, ONLY: &
33 : fac_x_to_z, ijkl_ind, ijkl_sym, inddd, inddp, indexa, indexb, indpp, int2c_type, l_index, &
34 : map_x_to_z, rij_threshold
35 : USE semi_empirical_int_num, ONLY: nucint_d_num, &
36 : nucint_sp_num, &
37 : terep_d_num, &
38 : terep_sp_num
39 : USE semi_empirical_int_utils, ONLY: d_ijkl_d, &
40 : d_ijkl_sp, &
41 : rot_2el_2c_first, &
42 : rotmat, &
43 : store_2el_2c_diag
44 : USE semi_empirical_types, ONLY: rotmat_create, &
45 : rotmat_release, &
46 : rotmat_type, &
47 : se_int_control_type, &
48 : se_int_screen_type, &
49 : se_taper_type, &
50 : semi_empirical_type, &
51 : setup_se_int_control_type
52 : USE taper_types, ONLY: dtaper_eval, &
53 : taper_eval
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 : PRIVATE
58 :
59 : #:include 'semi_empirical_int_debug.fypp'
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_ana'
62 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
63 : PUBLIC :: rotnuc_ana, rotint_ana, corecore_ana, corecore_el_ana
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Computes analytical gradients for semiempirical integrals
69 : !> \param sepi Atomic parameters of first atom
70 : !> \param sepj Atomic parameters of second atom
71 : !> \param rijv Coordinate vector i -> j
72 : !> \param itype ...
73 : !> \param e1b Array of electron-nuclear attraction integrals, Electron on atom ni attracting nucleus of nj.
74 : !> \param e2a Array of electron-nuclear attraction integrals, Electron on atom nj attracting nucleus of ni.
75 : !> \param de1b derivative of e1b term
76 : !> \param de2a derivative of e2a term
77 : !> \param se_int_control input parameters that control the calculation of SE
78 : !> integrals (shortrange, R3 residual, screening type)
79 : !> \param se_taper ...
80 : !> \par History
81 : !> 04.2007 created [tlaino]
82 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
83 : !> for computing integrals
84 : !> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the core-core part
85 : !> \author Teodoro Laino [tlaino] - Zurich University
86 : !> \note
87 : !> Analytical version of the MOPAC rotnuc routine
88 : ! **************************************************************************************************
89 16142300 : RECURSIVE SUBROUTINE rotnuc_ana(sepi, sepj, rijv, itype, e1b, e2a, de1b, de2a, &
90 : se_int_control, se_taper)
91 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
92 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
93 : INTEGER, INTENT(IN) :: itype
94 : REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a
95 : REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a
96 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
97 : TYPE(se_taper_type), POINTER :: se_taper
98 :
99 : INTEGER :: i, idd, idp, ind1, ind2, ipp, j, &
100 : last_orbital(2), m, n
101 : LOGICAL :: invert, l_de1b, l_de2a, l_e1b, l_e2a, &
102 : lgrad, task(2)
103 : REAL(KIND=dp) :: rij, xtmp
104 : REAL(KIND=dp), DIMENSION(10, 2) :: core, dcore
105 : REAL(KIND=dp), DIMENSION(3) :: drij
106 : REAL(KIND=dp), DIMENSION(3, 45) :: tmp_d
107 : REAL(KIND=dp), DIMENSION(45) :: tmp
108 : TYPE(rotmat_type), POINTER :: ij_matrix
109 :
110 16142300 : NULLIFY (ij_matrix)
111 64569200 : rij = DOT_PRODUCT(rijv, rijv)
112 : ! Initialization
113 16142300 : l_e1b = PRESENT(e1b)
114 16142300 : l_e2a = PRESENT(e2a)
115 16142300 : l_de1b = PRESENT(de1b)
116 16142300 : l_de2a = PRESENT(de2a)
117 16142300 : lgrad = l_de1b .OR. l_de2a
118 :
119 16142300 : IF (rij > rij_threshold) THEN
120 : ! Compute Integrals in diatomic frame opportunely inverted
121 16142300 : rij = SQRT(rij)
122 : ! Create the rotation matrix
123 16142300 : CALL rotmat_create(ij_matrix)
124 16142300 : CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert)
125 :
126 16142300 : IF (lgrad) THEN
127 8014152 : drij(1) = rijv(1)/rij
128 8014152 : drij(2) = rijv(2)/rij
129 8014152 : drij(3) = rijv(3)/rij
130 : ! Possibly Invert Frame
131 8014152 : IF (invert) THEN
132 879 : xtmp = drij(3)
133 879 : drij(3) = drij(1)
134 879 : drij(1) = xtmp
135 : END IF
136 : END IF
137 :
138 : CALL dcore_nucint_ana(sepi, sepj, rij, core=core, dcore=dcore, itype=itype, se_taper=se_taper, &
139 16142300 : se_int_control=se_int_control, lgrad=lgrad)
140 :
141 : ! Copy parameters over to arrays for do loop.
142 16142300 : last_orbital(1) = sepi%natorb
143 16142300 : last_orbital(2) = sepj%natorb
144 16142300 : task(1) = l_e1b
145 16142300 : task(2) = l_e2a
146 48426900 : DO n = 1, 2
147 32284600 : IF (.NOT. task(n)) CYCLE
148 28902913 : DO i = 1, last_orbital(n)
149 20318894 : ind1 = i - 1
150 73155897 : DO j = 1, i
151 44252984 : ind2 = j - 1
152 44252984 : m = (i*(i - 1))/2 + j
153 : ! Perform Rotations ...
154 64571878 : IF (ind2 == 0) THEN
155 20318894 : IF (ind1 == 0) THEN
156 : ! Type of Integral (SS/)
157 8584019 : tmp(m) = core(1, n)
158 11734875 : ELSE IF (ind1 < 4) THEN
159 : ! Type of Integral (SP/)
160 11618790 : tmp(m) = ij_matrix%sp(1, ind1)*core(2, n)
161 : ELSE
162 : ! Type of Integral (SD/)
163 116085 : tmp(m) = ij_matrix%sd(1, ind1 - 3)*core(5, n)
164 : END IF
165 23934090 : ELSE IF (ind2 < 4) THEN
166 23585835 : IF (ind1 < 4) THEN
167 : ! Type of Integral (PP/)
168 23237580 : ipp = indpp(ind1, ind2)
169 : tmp(m) = core(3, n)*ij_matrix%pp(ipp, 1, 1) + &
170 23237580 : core(4, n)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3))
171 : ELSE
172 : ! Type of Integral (PD/)
173 348255 : idp = inddp(ind1 - 3, ind2)
174 : tmp(m) = core(6, n)*ij_matrix%pd(idp, 1, 1) + &
175 348255 : core(8, n)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3))
176 : END IF
177 : ELSE
178 : ! Type of Integral (DD/)
179 348255 : idd = inddd(ind1 - 3, ind2 - 3)
180 : tmp(m) = core(7, n)*ij_matrix%dd(idd, 1, 1) + &
181 : core(9, n)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
182 348255 : core(10, n)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5))
183 : END IF
184 : END DO
185 : END DO
186 8584019 : IF (n == 1) THEN
187 49462846 : DO i = 1, sepi%atm_int_size
188 49462846 : e1b(i) = -tmp(i)
189 : END DO
190 : END IF
191 24726319 : IF (n == 2) THEN
192 3374157 : DO i = 1, sepj%atm_int_size
193 3374157 : e2a(i) = -tmp(i)
194 : END DO
195 : END IF
196 : END DO
197 16142300 : IF (invert .AND. l_e1b) CALL invert_integral(sepi, sepi, int1el=e1b)
198 16142300 : IF (invert .AND. l_e2a) CALL invert_integral(sepj, sepj, int1el=e2a)
199 :
200 : ! Possibly compute derivatives
201 16142300 : task(1) = l_de1b
202 16142300 : task(2) = l_de2a
203 48426900 : DO n = 1, 2
204 32284600 : IF (.NOT. task(n)) CYCLE
205 28141050 : DO i = 1, last_orbital(n)
206 19784366 : ind1 = i - 1
207 71181920 : DO j = 1, i
208 43040870 : ind2 = j - 1
209 43040870 : m = (i*(i - 1))/2 + j
210 : ! Perform Rotations ...
211 62825236 : IF (ind2 == 0) THEN
212 19784366 : IF (ind1 == 0) THEN
213 : ! Type of Integral (SS/)
214 8356684 : tmp_d(1, m) = dcore(1, n)*drij(1)
215 8356684 : tmp_d(2, m) = dcore(1, n)*drij(2)
216 8356684 : tmp_d(3, m) = dcore(1, n)*drij(3)
217 11427682 : ELSE IF (ind1 < 4) THEN
218 : ! Type of Integral (SP/)
219 : tmp_d(1, m) = ij_matrix%sp_d(1, 1, ind1)*core(2, n) + &
220 11327397 : ij_matrix%sp(1, ind1)*dcore(2, n)*drij(1)
221 :
222 : tmp_d(2, m) = ij_matrix%sp_d(2, 1, ind1)*core(2, n) + &
223 11327397 : ij_matrix%sp(1, ind1)*dcore(2, n)*drij(2)
224 :
225 : tmp_d(3, m) = ij_matrix%sp_d(3, 1, ind1)*core(2, n) + &
226 11327397 : ij_matrix%sp(1, ind1)*dcore(2, n)*drij(3)
227 : ELSE
228 : ! Type of Integral (SD/)
229 : tmp_d(1, m) = ij_matrix%sd_d(1, 1, ind1 - 3)*core(5, n) + &
230 100285 : ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(1)
231 :
232 : tmp_d(2, m) = ij_matrix%sd_d(2, 1, ind1 - 3)*core(5, n) + &
233 100285 : ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(2)
234 :
235 : tmp_d(3, m) = ij_matrix%sd_d(3, 1, ind1 - 3)*core(5, n) + &
236 100285 : ij_matrix%sd(1, ind1 - 3)*dcore(5, n)*drij(3)
237 : END IF
238 23256504 : ELSE IF (ind2 < 4) THEN
239 22955649 : IF (ind1 < 4) THEN
240 : ! Type of Integral (PP/)
241 22654794 : ipp = indpp(ind1, ind2)
242 : tmp_d(1, m) = dcore(3, n)*drij(1)*ij_matrix%pp(ipp, 1, 1) + &
243 : core(3, n)*ij_matrix%pp_d(1, ipp, 1, 1) + &
244 : dcore(4, n)*drij(1)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
245 22654794 : core(4, n)*(ij_matrix%pp_d(1, ipp, 2, 2) + ij_matrix%pp_d(1, ipp, 3, 3))
246 :
247 : tmp_d(2, m) = dcore(3, n)*drij(2)*ij_matrix%pp(ipp, 1, 1) + &
248 : core(3, n)*ij_matrix%pp_d(2, ipp, 1, 1) + &
249 : dcore(4, n)*drij(2)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
250 22654794 : core(4, n)*(ij_matrix%pp_d(2, ipp, 2, 2) + ij_matrix%pp_d(2, ipp, 3, 3))
251 :
252 : tmp_d(3, m) = dcore(3, n)*drij(3)*ij_matrix%pp(ipp, 1, 1) + &
253 : core(3, n)*ij_matrix%pp_d(3, ipp, 1, 1) + &
254 : dcore(4, n)*drij(3)*(ij_matrix%pp(ipp, 2, 2) + ij_matrix%pp(ipp, 3, 3)) + &
255 22654794 : core(4, n)*(ij_matrix%pp_d(3, ipp, 2, 2) + ij_matrix%pp_d(3, ipp, 3, 3))
256 : ELSE
257 : ! Type of Integral (PD/)
258 300855 : idp = inddp(ind1 - 3, ind2)
259 : tmp_d(1, m) = dcore(6, n)*drij(1)*ij_matrix%pd(idp, 1, 1) + &
260 : core(6, n)*ij_matrix%pd_d(1, idp, 1, 1) + &
261 : dcore(8, n)*drij(1)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
262 300855 : core(8, n)*(ij_matrix%pd_d(1, idp, 2, 2) + ij_matrix%pd_d(1, idp, 3, 3))
263 :
264 : tmp_d(2, m) = dcore(6, n)*drij(2)*ij_matrix%pd(idp, 1, 1) + &
265 : core(6, n)*ij_matrix%pd_d(2, idp, 1, 1) + &
266 : dcore(8, n)*drij(2)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
267 300855 : core(8, n)*(ij_matrix%pd_d(2, idp, 2, 2) + ij_matrix%pd_d(2, idp, 3, 3))
268 :
269 : tmp_d(3, m) = dcore(6, n)*drij(3)*ij_matrix%pd(idp, 1, 1) + &
270 : core(6, n)*ij_matrix%pd_d(3, idp, 1, 1) + &
271 : dcore(8, n)*drij(3)*(ij_matrix%pd(idp, 2, 2) + ij_matrix%pd(idp, 3, 3)) + &
272 300855 : core(8, n)*(ij_matrix%pd_d(3, idp, 2, 2) + ij_matrix%pd_d(3, idp, 3, 3))
273 : END IF
274 : ELSE
275 : ! Type of Integral (DD/)
276 300855 : idd = inddd(ind1 - 3, ind2 - 3)
277 : tmp_d(1, m) = dcore(7, n)*drij(1)*ij_matrix%dd(idd, 1, 1) + &
278 : core(7, n)*ij_matrix%dd_d(1, idd, 1, 1) + &
279 : dcore(9, n)*drij(1)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
280 : core(9, n)*(ij_matrix%dd_d(1, idd, 2, 2) + ij_matrix%dd_d(1, idd, 3, 3)) + &
281 : dcore(10, n)*drij(1)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
282 300855 : core(10, n)*(ij_matrix%dd_d(1, idd, 4, 4) + ij_matrix%dd_d(1, idd, 5, 5))
283 :
284 : tmp_d(2, m) = dcore(7, n)*drij(2)*ij_matrix%dd(idd, 1, 1) + &
285 : core(7, n)*ij_matrix%dd_d(2, idd, 1, 1) + &
286 : dcore(9, n)*drij(2)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
287 : core(9, n)*(ij_matrix%dd_d(2, idd, 2, 2) + ij_matrix%dd_d(2, idd, 3, 3)) + &
288 : dcore(10, n)*drij(2)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
289 300855 : core(10, n)*(ij_matrix%dd_d(2, idd, 4, 4) + ij_matrix%dd_d(2, idd, 5, 5))
290 :
291 : tmp_d(3, m) = dcore(7, n)*drij(3)*ij_matrix%dd(idd, 1, 1) + &
292 : core(7, n)*ij_matrix%dd_d(3, idd, 1, 1) + &
293 : dcore(9, n)*drij(3)*(ij_matrix%dd(idd, 2, 2) + ij_matrix%dd(idd, 3, 3)) + &
294 : core(9, n)*(ij_matrix%dd_d(3, idd, 2, 2) + ij_matrix%dd_d(3, idd, 3, 3)) + &
295 : dcore(10, n)*drij(3)*(ij_matrix%dd(idd, 4, 4) + ij_matrix%dd(idd, 5, 5)) + &
296 300855 : core(10, n)*(ij_matrix%dd_d(3, idd, 4, 4) + ij_matrix%dd_d(3, idd, 5, 5))
297 : END IF
298 : END DO
299 : END DO
300 8356684 : IF (n == 1) THEN
301 48741983 : DO i = 1, sepi%atm_int_size
302 40727831 : de1b(1, i) = -tmp_d(1, i)
303 40727831 : de1b(2, i) = -tmp_d(2, i)
304 48741983 : de1b(3, i) = -tmp_d(3, i)
305 : END DO
306 : END IF
307 24498984 : IF (n == 2) THEN
308 2655571 : DO i = 1, sepj%atm_int_size
309 2313039 : de2a(1, i) = -tmp_d(1, i)
310 2313039 : de2a(2, i) = -tmp_d(2, i)
311 2655571 : de2a(3, i) = -tmp_d(3, i)
312 : END DO
313 : END IF
314 : END DO
315 16142300 : IF (invert .AND. l_de1b) CALL invert_derivative(sepi, sepi, dint1el=de1b)
316 16142300 : IF (invert .AND. l_de2a) CALL invert_derivative(sepj, sepj, dint1el=de2a)
317 16142300 : CALL rotmat_release(ij_matrix)
318 :
319 : ! Possibly debug the analytical values versus the numerical ones
320 : IF (debug_this_module) THEN
321 : CALL check_drotnuc_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, e1b, e2a, de1b, de2a)
322 : END IF
323 : END IF
324 16142300 : END SUBROUTINE rotnuc_ana
325 :
326 : ! **************************************************************************************************
327 : !> \brief Computes analytical gradients for semiempirical core-core interaction.
328 : !> \param sepi Atomic parameters of first atom
329 : !> \param sepj Atomic parameters of second atom
330 : !> \param rijv Coordinate vector i -> j
331 : !> \param itype ...
332 : !> \param enuc nuclear-nuclear repulsion term.
333 : !> \param denuc derivative of nuclear-nuclear repulsion term.
334 : !> \param se_int_control input parameters that control the calculation of SE
335 : !> integrals (shortrange, R3 residual, screening type)
336 : !> \param se_taper ...
337 : !> \par History
338 : !> 04.2007 created [tlaino]
339 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
340 : !> for computing integrals
341 : !> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the
342 : !> core-core part
343 : !> \author Teodoro Laino [tlaino] - Zurich University
344 : !> \note
345 : !> Analytical version of the MOPAC rotnuc routine
346 : ! **************************************************************************************************
347 23294450 : RECURSIVE SUBROUTINE corecore_ana(sepi, sepj, rijv, itype, enuc, denuc, se_int_control, se_taper)
348 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
349 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
350 : INTEGER, INTENT(IN) :: itype
351 : REAL(dp), INTENT(OUT), OPTIONAL :: enuc
352 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc
353 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
354 : TYPE(se_taper_type), POINTER :: se_taper
355 :
356 : INTEGER :: ig, nt
357 : LOGICAL :: l_denuc, l_enuc
358 : REAL(dp) :: aab, alpi, alpj, apdg, ax, dai, daj, dax, dbi, dbj, denuc_loc, dqcorr, drija, &
359 : dscale, dssss, dssss_sr, dtmp, dzz, enuc_loc, pai, paj, pbi, pbj, qcorr, rij, rija, &
360 : scale, ssss, ssss_sr, tmp, xab, xtmp, zaf, zbf, zz
361 : REAL(dp), DIMENSION(3) :: drij
362 : REAL(dp), DIMENSION(4) :: fni1, fni2, fni3, fnj1, fnj2, fnj3
363 : TYPE(se_int_control_type) :: se_int_control_off
364 :
365 93177800 : rij = DOT_PRODUCT(rijv, rijv)
366 : ! Initialization
367 23294450 : l_enuc = PRESENT(enuc)
368 23294450 : l_denuc = PRESENT(denuc)
369 23294450 : IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
370 : ! Compute Integrals in diatomic frame
371 23294450 : rij = SQRT(rij)
372 : CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., &
373 : do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, &
374 23294450 : max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
375 : CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, &
376 23294450 : se_int_control=se_int_control_off, lgrad=l_denuc)
377 : ! In case let's compute the short-range part of the (ss|ss) integral
378 23294450 : IF (se_int_control%shortrange) THEN
379 : CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, &
380 248393 : se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc)
381 : ELSE
382 23046057 : ssss_sr = ssss
383 23046057 : dssss_sr = dssss
384 : END IF
385 : ! Zeroing local method dependent core-core corrections
386 23294450 : enuc_loc = 0.0_dp
387 23294450 : denuc_loc = 0.0_dp
388 23294450 : qcorr = 0.0_dp
389 23294450 : scale = 0.0_dp
390 23294450 : dscale = 0.0_dp
391 23294450 : dqcorr = 0.0_dp
392 23294450 : zz = sepi%zeff*sepj%zeff
393 : ! Core Core electrostatic contribution
394 23294450 : IF (l_enuc) enuc_loc = zz*ssss_sr
395 23294450 : IF (l_denuc) denuc_loc = zz*dssss_sr
396 : ! Method dependent code
397 23294450 : tmp = zz*ssss
398 23294450 : IF (l_denuc) dtmp = zz*dssss
399 23294450 : IF (itype /= do_method_pm6 .AND. itype /= do_method_pm6fm) THEN
400 15726037 : alpi = sepi%alp
401 15726037 : alpj = sepj%alp
402 15726037 : scale = EXP(-alpi*rij) + EXP(-alpj*rij)
403 15726037 : IF (l_denuc) THEN
404 7796641 : dscale = -alpi*EXP(-alpi*rij) - alpj*EXP(-alpj*rij)
405 : END IF
406 15726037 : nt = sepi%z + sepj%z
407 15726037 : IF (nt == 8 .OR. nt == 9) THEN
408 2993026 : IF (sepi%z == 7 .OR. sepi%z == 8) THEN
409 2961656 : scale = scale + (angstrom*rij - 1._dp)*EXP(-alpi*rij)
410 2961656 : IF (l_denuc) THEN
411 1475382 : dscale = dscale + angstrom*EXP(-alpi*rij) - (angstrom*rij - 1._dp)*alpi*EXP(-alpi*rij)
412 : END IF
413 : END IF
414 2993026 : IF (sepj%z == 7 .OR. sepj%z == 8) THEN
415 31370 : scale = scale + (angstrom*rij - 1._dp)*EXP(-alpj*rij)
416 31370 : IF (l_denuc) THEN
417 10395 : dscale = dscale + angstrom*EXP(-alpj*rij) - (angstrom*rij - 1._dp)*alpj*EXP(-alpj*rij)
418 : END IF
419 : END IF
420 : END IF
421 15726037 : IF (l_denuc) THEN
422 7796641 : dscale = SIGN(1.0_dp, scale*tmp)*(dscale*tmp + scale*dtmp)
423 7796641 : dzz = -zz/rij**2
424 : END IF
425 15726037 : scale = ABS(scale*tmp)
426 15726037 : zz = zz/rij
427 15726037 : IF (itype == do_method_am1 .OR. itype == do_method_pm3 .OR. itype == do_method_pdg) THEN
428 302787 : IF (itype == do_method_am1 .AND. sepi%z == 5) THEN
429 : !special case AM1 Boron
430 : SELECT CASE (sepj%z)
431 : CASE DEFAULT
432 0 : nt = 1
433 : CASE (1)
434 0 : nt = 2
435 : CASE (6)
436 0 : nt = 3
437 : CASE (9, 17, 35, 53)
438 0 : nt = 4
439 : END SELECT
440 0 : fni1(1) = sepi%bfn1(1, nt)
441 0 : fni1(2) = sepi%bfn1(2, nt)
442 0 : fni1(3) = sepi%bfn1(3, nt)
443 0 : fni1(4) = sepi%bfn1(4, nt)
444 0 : fni2(1) = sepi%bfn2(1, nt)
445 0 : fni2(2) = sepi%bfn2(2, nt)
446 0 : fni2(3) = sepi%bfn2(3, nt)
447 0 : fni2(4) = sepi%bfn2(4, nt)
448 0 : fni3(1) = sepi%bfn3(1, nt)
449 0 : fni3(2) = sepi%bfn3(2, nt)
450 0 : fni3(3) = sepi%bfn3(3, nt)
451 0 : fni3(4) = sepi%bfn3(4, nt)
452 : ELSE
453 302787 : fni1(1) = sepi%fn1(1)
454 302787 : fni1(2) = sepi%fn1(2)
455 302787 : fni1(3) = sepi%fn1(3)
456 302787 : fni1(4) = sepi%fn1(4)
457 302787 : fni2(1) = sepi%fn2(1)
458 302787 : fni2(2) = sepi%fn2(2)
459 302787 : fni2(3) = sepi%fn2(3)
460 302787 : fni2(4) = sepi%fn2(4)
461 302787 : fni3(1) = sepi%fn3(1)
462 302787 : fni3(2) = sepi%fn3(2)
463 302787 : fni3(3) = sepi%fn3(3)
464 302787 : fni3(4) = sepi%fn3(4)
465 : END IF
466 302787 : IF (itype == do_method_am1 .AND. sepj%z == 5) THEN
467 : !special case AM1 Boron
468 : SELECT CASE (sepi%z)
469 : CASE DEFAULT
470 0 : nt = 1
471 : CASE (1)
472 0 : nt = 2
473 : CASE (6)
474 0 : nt = 3
475 : CASE (9, 17, 35, 53)
476 0 : nt = 4
477 : END SELECT
478 0 : fnj1(1) = sepj%bfn1(1, nt)
479 0 : fnj1(2) = sepj%bfn1(2, nt)
480 0 : fnj1(3) = sepj%bfn1(3, nt)
481 0 : fnj1(4) = sepj%bfn1(4, nt)
482 0 : fnj2(1) = sepj%bfn2(1, nt)
483 0 : fnj2(2) = sepj%bfn2(2, nt)
484 0 : fnj2(3) = sepj%bfn2(3, nt)
485 0 : fnj2(4) = sepj%bfn2(4, nt)
486 0 : fnj3(1) = sepj%bfn3(1, nt)
487 0 : fnj3(2) = sepj%bfn3(2, nt)
488 0 : fnj3(3) = sepj%bfn3(3, nt)
489 0 : fnj3(4) = sepj%bfn3(4, nt)
490 : ELSE
491 302787 : fnj1(1) = sepj%fn1(1)
492 302787 : fnj1(2) = sepj%fn1(2)
493 302787 : fnj1(3) = sepj%fn1(3)
494 302787 : fnj1(4) = sepj%fn1(4)
495 302787 : fnj2(1) = sepj%fn2(1)
496 302787 : fnj2(2) = sepj%fn2(2)
497 302787 : fnj2(3) = sepj%fn2(3)
498 302787 : fnj2(4) = sepj%fn2(4)
499 302787 : fnj3(1) = sepj%fn3(1)
500 302787 : fnj3(2) = sepj%fn3(2)
501 302787 : fnj3(3) = sepj%fn3(3)
502 302787 : fnj3(4) = sepj%fn3(4)
503 : END IF
504 : ! AM1/PM3/PDG correction to nuclear repulsion
505 1513935 : DO ig = 1, SIZE(fni1)
506 1211148 : IF (ABS(fni1(ig)) > 0._dp) THEN
507 836741 : ax = fni2(ig)*(rij - fni3(ig))**2
508 836741 : IF (ax <= 25._dp) THEN
509 238319 : scale = scale + zz*fni1(ig)*EXP(-ax)
510 238319 : IF (l_denuc) THEN
511 79200 : dax = fni2(ig)*2.0_dp*(rij - fni3(ig))
512 79200 : dscale = dscale + dzz*fni1(ig)*EXP(-ax) - dax*zz*fni1(ig)*EXP(-ax)
513 : END IF
514 : END IF
515 : END IF
516 1513935 : IF (ABS(fnj1(ig)) > 0._dp) THEN
517 836729 : ax = fnj2(ig)*(rij - fnj3(ig))**2
518 836729 : IF (ax <= 25._dp) THEN
519 237956 : scale = scale + zz*fnj1(ig)*EXP(-ax)
520 237956 : IF (l_denuc) THEN
521 79084 : dax = fnj2(ig)*2.0_dp*(rij - fnj3(ig))
522 79084 : dscale = dscale + dzz*fnj1(ig)*EXP(-ax) - dax*zz*fnj1(ig)*EXP(-ax)
523 : END IF
524 : END IF
525 : END IF
526 : END DO
527 : END IF
528 15726037 : IF (itype == do_method_pdg) THEN
529 : ! PDDG function
530 9 : zaf = sepi%zeff/nt
531 9 : zbf = sepj%zeff/nt
532 9 : pai = sepi%pre(1)
533 9 : pbi = sepi%pre(2)
534 9 : paj = sepj%pre(1)
535 9 : pbj = sepj%pre(2)
536 9 : dai = sepi%d(1)
537 9 : dbi = sepi%d(2)
538 9 : daj = sepj%d(1)
539 9 : dbj = sepj%d(2)
540 9 : apdg = 10._dp*angstrom**2
541 : qcorr = (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2) + &
542 : (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2) + &
543 : (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2) + &
544 9 : (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2)
545 9 : IF (l_denuc) THEN
546 : dqcorr = (zaf*pai + zbf*paj)*EXP(-apdg*(rij - dai - daj)**2)*(-2.0_dp*apdg*(rij - dai - daj)) + &
547 : (zaf*pai + zbf*pbj)*EXP(-apdg*(rij - dai - dbj)**2)*(-2.0_dp*apdg*(rij - dai - dbj)) + &
548 : (zaf*pbi + zbf*paj)*EXP(-apdg*(rij - dbi - daj)**2)*(-2.0_dp*apdg*(rij - dbi - daj)) + &
549 3 : (zaf*pbi + zbf*pbj)*EXP(-apdg*(rij - dbi - dbj)**2)*(-2.0_dp*apdg*(rij - dbi - dbj))
550 : END IF
551 15726028 : ELSEIF (itype == do_method_pchg) THEN
552 : qcorr = 0.0_dp
553 : scale = 0.0_dp
554 : dscale = 0.0_dp
555 : dqcorr = 0.0_dp
556 : ELSE
557 382131 : qcorr = 0.0_dp
558 382131 : dqcorr = 0.0_dp
559 : END IF
560 : ELSE
561 : ! PM6 core-core terms
562 7568413 : scale = tmp
563 7568413 : IF (l_denuc) dscale = dtmp
564 7568413 : drija = angstrom
565 7568413 : rija = rij*drija
566 7568413 : xab = sepi%xab(sepj%z)
567 7568413 : aab = sepi%aab(sepj%z)
568 7568413 : IF ((sepi%z == 1 .AND. (sepj%z == 6 .OR. sepj%z == 7 .OR. sepj%z == 8)) .OR. &
569 : (sepj%z == 1 .AND. (sepi%z == 6 .OR. sepi%z == 7 .OR. sepi%z == 8))) THEN
570 : ! Special Case O-H or N-H or C-H
571 253301 : IF (l_denuc) dscale = dscale*(2._dp*xab*EXP(-aab*rija*rija)) - &
572 55301 : scale*2._dp*xab*EXP(-aab*rija*rija)*(2.0_dp*aab*rija)*drija
573 253301 : IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*rija*rija))
574 7315112 : ELSEIF (sepi%z == 6 .AND. sepj%z == 6) THEN
575 : ! Special Case C-C
576 352747 : IF (l_denuc) dscale = &
577 : dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija)) &
578 : - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija &
579 12485 : - scale*9.28_dp*EXP(-5.98_dp*rija)*5.98_dp*drija
580 352747 : IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) + 9.28_dp*EXP(-5.98_dp*rija))
581 6962365 : ELSEIF ((sepi%z == 8 .AND. sepj%z == 14) .OR. &
582 : (sepj%z == 8 .AND. sepi%z == 14)) THEN
583 : ! Special Case Si-O
584 97448 : IF (l_denuc) dscale = &
585 : dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2)) &
586 : - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija + &
587 2 : scale*0.0007_dp*EXP(-(rija - 2.9_dp)**2)*(2.0_dp*(rija - 2.9_dp)*drija)
588 97448 : IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)) - 0.0007_dp*EXP(-(rija - 2.9_dp)**2))
589 : ELSE
590 : ! General Case
591 : ! Factor of 2 found by experiment
592 6864917 : IF (l_denuc) dscale = dscale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))) &
593 105775 : - scale*2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6))*aab*(1.0_dp + 6.0_dp*0.0003_dp*rija**5)*drija
594 6864917 : IF (l_enuc) scale = scale*(2._dp*xab*EXP(-aab*(rija + 0.0003_dp*rija**6)))
595 : END IF
596 : ! General correction term a*exp(-b*(rij-c)^2)
597 7568413 : xtmp = 1.e-8_dp/evolt*((REAL(sepi%z, dp)**(1._dp/3._dp) + REAL(sepj%z, dp)**(1._dp/3._dp))/rija)**12
598 7568413 : IF (l_enuc) THEN
599 : qcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*zz/rij + &
600 : (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*zz/rij + &
601 : ! Hard core repulsion
602 7394850 : xtmp
603 : END IF
604 7568413 : IF (l_denuc) THEN
605 : dqcorr = (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2)*(-2.0_dp*sepi%b*(rij - sepi%c)))*zz/rij - &
606 : (sepi%a*EXP(-sepi%b*(rij - sepi%c)**2))*zz/rij**2 + &
607 : (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2)*(-2.0_dp*sepj%b*(rij - sepj%c)))*zz/rij - &
608 : (sepj%a*EXP(-sepj%b*(rij - sepj%c)**2))*zz/rij**2 + &
609 : ! Hard core repulsion
610 173563 : (-12.0_dp*xtmp/rija*drija)
611 : END IF
612 : END IF
613 :
614 : ! Only at the very end let's sum-up the several contributions energy/derivatives
615 : ! This assignment should be method independent
616 23294450 : IF (l_enuc) THEN
617 15324246 : enuc = enuc_loc + scale + qcorr
618 : END IF
619 23294450 : IF (l_denuc) THEN
620 7970204 : drij(1) = rijv(1)/rij
621 7970204 : drij(2) = rijv(2)/rij
622 7970204 : drij(3) = rijv(3)/rij
623 31880816 : denuc = (denuc_loc + dscale + dqcorr)*drij
624 : END IF
625 : ! Debug statement
626 : IF (debug_this_module) THEN
627 : CALL check_dcorecore_ana(sepi, sepj, rijv, itype, se_int_control, se_taper, enuc, denuc)
628 : END IF
629 : END IF
630 23294450 : END SUBROUTINE corecore_ana
631 :
632 : ! **************************************************************************************************
633 : !> \brief Computes analytical gradients for semiempirical core-core electrostatic
634 : !> interaction only.
635 : !> \param sepi Atomic parameters of first atom
636 : !> \param sepj Atomic parameters of second atom
637 : !> \param rijv Coordinate vector i -> j
638 : !> \param itype ...
639 : !> \param enuc nuclear-nuclear electrostatic repulsion term.
640 : !> \param denuc derivative of nuclear-nuclear electrostatic
641 : !> repulsion term.
642 : !> \param se_int_control input parameters that control the calculation of SE
643 : !> integrals (shortrange, R3 residual, screening type)
644 : !> \param se_taper ...
645 : !> \par History
646 : !> 04.2007 created [tlaino]
647 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
648 : !> for computing integrals
649 : !> Teodoro Laino [tlaino] - University of Zurich 04.2008 : removed the
650 : !> core-core part
651 : !> \author Teodoro Laino [tlaino] - Zurich University
652 : !> \note
653 : !> Analytical version of the MOPAC rotnuc routine
654 : ! **************************************************************************************************
655 1466621 : RECURSIVE SUBROUTINE corecore_el_ana(sepi, sepj, rijv, itype, enuc, denuc, &
656 : se_int_control, se_taper)
657 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
658 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
659 : INTEGER, INTENT(IN) :: itype
660 : REAL(dp), INTENT(OUT), OPTIONAL :: enuc
661 : REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc
662 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
663 : TYPE(se_taper_type), POINTER :: se_taper
664 :
665 : LOGICAL :: l_denuc, l_enuc
666 : REAL(dp) :: drij(3), dssss, dssss_sr, rij, ssss, &
667 : ssss_sr, tmp, zz
668 : TYPE(se_int_control_type) :: se_int_control_off
669 :
670 5866484 : rij = DOT_PRODUCT(rijv, rijv)
671 : ! Initialization
672 1466621 : l_enuc = PRESENT(enuc)
673 1466621 : l_denuc = PRESENT(denuc)
674 1466621 : IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN
675 : ! Compute Integrals in diatomic frame
676 1466621 : rij = SQRT(rij)
677 : CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., &
678 : do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, &
679 1466621 : max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
680 : CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss, dssss=dssss, itype=itype, se_taper=se_taper, &
681 1466621 : se_int_control=se_int_control_off, lgrad=l_denuc)
682 : ! In case let's compute the short-range part of the (ss|ss) integral
683 1466621 : IF (se_int_control%shortrange .OR. se_int_control%pc_coulomb_int) THEN
684 : CALL dssss_nucint_ana(sepi, sepj, rij, ssss=ssss_sr, dssss=dssss_sr, itype=itype, &
685 1466621 : se_taper=se_taper, se_int_control=se_int_control, lgrad=l_denuc)
686 : ELSE
687 0 : ssss_sr = ssss
688 0 : dssss_sr = dssss
689 : END IF
690 1466621 : zz = sepi%zeff*sepj%zeff
691 : ! Core Core electrostatic contribution
692 1466621 : IF (l_enuc) enuc = zz*ssss_sr
693 1466621 : IF (l_denuc) THEN
694 43876 : drij(1) = rijv(1)/rij
695 43876 : drij(2) = rijv(2)/rij
696 43876 : drij(3) = rijv(3)/rij
697 43876 : tmp = zz*dssss_sr
698 175504 : denuc = tmp*drij
699 : END IF
700 : END IF
701 1466621 : END SUBROUTINE corecore_el_ana
702 :
703 : ! **************************************************************************************************
704 : !> \brief Exploits inversion symmetry to avoid divergence
705 : !> \param sepi ...
706 : !> \param sepj ...
707 : !> \param int1el ...
708 : !> \param int2el ...
709 : !> \par History
710 : !> 04.2007 created [tlaino]
711 : !> 05.2008 New driver for integral invertion (supports d-orbitals)
712 : !> \author Teodoro Laino - Zurich University
713 : ! **************************************************************************************************
714 13358 : SUBROUTINE invert_integral(sepi, sepj, int1el, int2el)
715 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
716 : REAL(dp), DIMENSION(:), INTENT(INOUT), OPTIONAL :: int1el, int2el
717 :
718 : INTEGER :: fdim, gind, gknd, i, imap, ind, j, jmap, &
719 : jnd, k, kmap, knd, l, lmap, lnd, ndim, &
720 : sdim, tdim, tind
721 : REAL(KIND=dp) :: ifac, jfac, kfac, lfac
722 : REAL(KIND=dp), DIMENSION(2025) :: tmp2el
723 : REAL(KIND=dp), DIMENSION(45) :: tmp1el
724 :
725 : ! One-electron integral
726 :
727 13358 : IF (PRESENT(int1el)) THEN
728 8726 : fdim = sepi%atm_int_size
729 8726 : ndim = 0
730 108031 : DO i = 1, fdim
731 108031 : tmp1el(i) = 0.0_dp
732 : END DO
733 44845 : DO i = 1, sepi%natorb
734 144150 : DO j = 1, i
735 99305 : ndim = ndim + 1
736 :
737 : ! Get the integral in the original frame (along z)
738 334034 : DO ind = 1, 2
739 198610 : imap = map_x_to_z(ind, i)
740 198610 : IF (imap == 0) CYCLE
741 104345 : ifac = fac_x_to_z(ind, i)
742 412340 : DO jnd = 1, 2
743 208690 : jmap = map_x_to_z(jnd, j)
744 208690 : IF (jmap == 0) CYCLE
745 108965 : jfac = fac_x_to_z(jnd, j)
746 108965 : gind = indexb(imap, jmap)
747 :
748 407300 : tmp1el(ndim) = tmp1el(ndim) + ifac*jfac*int1el(gind)
749 : END DO
750 : END DO
751 : END DO
752 : END DO
753 108031 : DO i = 1, fdim
754 108031 : int1el(i) = tmp1el(i)
755 : END DO
756 : END IF
757 :
758 : ! Two electron integrals
759 13358 : IF (PRESENT(int2el)) THEN
760 4632 : sdim = sepi%atm_int_size
761 4632 : tdim = sepj%atm_int_size
762 4632 : fdim = sdim*tdim
763 4632 : ndim = 0
764 976572 : DO i = 1, fdim
765 976572 : tmp2el(i) = 0.0_dp
766 : END DO
767 24516 : DO i = 1, sepi%natorb
768 80424 : DO j = 1, i
769 353124 : DO k = 1, sepj%natorb
770 1305180 : DO l = 1, k
771 971940 : ndim = ndim + 1
772 :
773 : ! Get the integral in the original frame (along z)
774 3193152 : DO ind = 1, 2
775 1943880 : imap = map_x_to_z(ind, i)
776 1943880 : IF (imap == 0) CYCLE
777 1120980 : ifac = fac_x_to_z(ind, i)
778 4334880 : DO jnd = 1, 2
779 2241960 : jmap = map_x_to_z(jnd, j)
780 2241960 : IF (jmap == 0) CYCLE
781 1257600 : jfac = fac_x_to_z(jnd, j)
782 1257600 : gind = indexb(imap, jmap)
783 :
784 : ! Get the integral in the original frame (along z)
785 5716680 : DO knd = 1, 2
786 2515200 : kmap = map_x_to_z(knd, k)
787 2515200 : IF (kmap == 0) CYCLE
788 1484832 : kfac = fac_x_to_z(knd, k)
789 6696456 : DO lnd = 1, 2
790 2969664 : lmap = map_x_to_z(lnd, l)
791 2969664 : IF (lmap == 0) CYCLE
792 1693128 : lfac = fac_x_to_z(lnd, l)
793 1693128 : gknd = indexb(kmap, lmap)
794 :
795 1693128 : tind = (gind - 1)*tdim + gknd
796 5484864 : tmp2el(ndim) = tmp2el(ndim) + ifac*jfac*lfac*kfac*int2el(tind)
797 : END DO
798 : END DO
799 :
800 : END DO
801 : END DO
802 :
803 : END DO
804 : END DO
805 : END DO
806 : END DO
807 976572 : DO i = 1, fdim
808 976572 : int2el(i) = tmp2el(i)
809 : END DO
810 : END IF
811 13358 : END SUBROUTINE invert_integral
812 :
813 : ! **************************************************************************************************
814 : !> \brief Exploits inversion symmetry to avoid divergence
815 : !> \param sepi ...
816 : !> \param sepj ...
817 : !> \param dint1el ...
818 : !> \param dint2el ...
819 : !> \par History
820 : !> 04.2007 created [tlaino]
821 : !> \author Teodoro Laino - Zurich University
822 : ! **************************************************************************************************
823 2700 : SUBROUTINE invert_derivative(sepi, sepj, dint1el, dint2el)
824 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
825 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
826 : OPTIONAL :: dint1el, dint2el
827 :
828 : INTEGER :: i, m
829 : REAL(KIND=dp) :: tmp
830 :
831 : ! Integral part
832 :
833 10800 : DO i = 1, 3
834 8100 : IF (PRESENT(dint1el)) THEN
835 5259 : CALL invert_integral(sepi, sepj, int1el=dint1el(i, :))
836 : END IF
837 10800 : IF (PRESENT(dint2el)) THEN
838 2841 : CALL invert_integral(sepi, sepj, int2el=dint2el(i, :))
839 : END IF
840 : END DO
841 :
842 : ! Derivatives part
843 2700 : IF (PRESENT(dint1el)) THEN
844 80638 : DO m = 1, SIZE(dint1el, 2)
845 78885 : tmp = dint1el(3, m)
846 78885 : dint1el(3, m) = dint1el(1, m)
847 80638 : dint1el(1, m) = tmp
848 : END DO
849 : END IF
850 2700 : IF (PRESENT(dint2el)) THEN
851 1918622 : DO m = 1, SIZE(dint2el, 2)
852 1917675 : tmp = dint2el(3, m)
853 1917675 : dint2el(3, m) = dint2el(1, m)
854 1918622 : dint2el(1, m) = tmp
855 : END DO
856 : END IF
857 2700 : END SUBROUTINE invert_derivative
858 :
859 : ! **************************************************************************************************
860 : !> \brief Calculates the ssss integral and analytical derivatives (main driver)
861 : !> \param sepi parameters of atom i
862 : !> \param sepj parameters of atom j
863 : !> \param rij interatomic distance
864 : !> \param ssss ...
865 : !> \param dssss derivative of (ssss) integral
866 : !> derivatives are intended w.r.t. rij
867 : !> \param itype ...
868 : !> \param se_taper ...
869 : !> \param se_int_control input parameters that control the calculation of SE
870 : !> integrals (shortrange, R3 residual, screening type)
871 : !> \param lgrad ...
872 : !> \par History
873 : !> 03.2007 created [tlaino]
874 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
875 : !> for computing integrals
876 : !> \author Teodoro Laino - Zurich University
877 : !> \note
878 : !> Analytical version - Analytical evaluation of gradients
879 : !> Teodoro Laino - Zurich University 04.2007
880 : !>
881 : ! **************************************************************************************************
882 26476085 : SUBROUTINE dssss_nucint_ana(sepi, sepj, rij, ssss, dssss, itype, se_taper, se_int_control, &
883 : lgrad)
884 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
885 : REAL(dp), INTENT(IN) :: rij
886 : REAL(dp), INTENT(OUT) :: ssss, dssss
887 : INTEGER, INTENT(IN) :: itype
888 : TYPE(se_taper_type), POINTER :: se_taper
889 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
890 : LOGICAL, INTENT(IN) :: lgrad
891 :
892 : REAL(KIND=dp) :: dft, ft
893 : TYPE(se_int_screen_type) :: se_int_screen
894 :
895 : ! Compute the Tapering function
896 :
897 26476085 : ft = 1.0_dp
898 26476085 : dft = 0.0_dp
899 26476085 : IF (itype /= do_method_pchg) THEN
900 11132188 : ft = taper_eval(se_taper%taper, rij)
901 11132188 : dft = dtaper_eval(se_taper%taper, rij)
902 : END IF
903 : ! Evaluate additional taper function for dumped integrals
904 26476085 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
905 3429047 : se_int_screen%ft = 1.0_dp
906 3429047 : se_int_screen%dft = 0.0_dp
907 3429047 : IF (itype /= do_method_pchg) THEN
908 3429047 : se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
909 3429047 : se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
910 : END IF
911 : END IF
912 :
913 : ! Value of the integrals for sp shell
914 : CALL nucint_sp_num(sepi, sepj, rij, ssss=ssss, itype=itype, se_int_control=se_int_control, &
915 26476085 : se_int_screen=se_int_screen)
916 :
917 26476085 : IF (lgrad) THEN
918 : ! Integrals derivatives for sp shell
919 : CALL dnucint_sp_ana(sepi, sepj, rij, dssss=dssss, itype=itype, se_int_control=se_int_control, &
920 8139105 : se_int_screen=se_int_screen)
921 : END IF
922 :
923 : ! Tapering the value of the integrals
924 26476085 : IF (lgrad) THEN
925 8139105 : dssss = ft*dssss + dft*ssss
926 : END IF
927 26476085 : ssss = ft*ssss
928 :
929 : ! Debug Procedure.. Check valifity of analytical gradients of nucint
930 : IF (debug_this_module .AND. lgrad) THEN
931 : CALL check_dssss_nucint_ana(sepi, sepj, rij, dssss, itype, se_int_control, se_taper=se_taper)
932 : END IF
933 26476085 : END SUBROUTINE dssss_nucint_ana
934 :
935 : ! **************************************************************************************************
936 : !> \brief Calculates the nuclear attraction integrals and analytical integrals (main driver)
937 : !> \param sepi parameters of atom i
938 : !> \param sepj parameters of atom j
939 : !> \param rij interatomic distance
940 : !> \param core ...
941 : !> \param dcore derivative of 4 X 2 array of electron-core attraction integrals
942 : !> derivatives are intended w.r.t. rij
943 : !> The storage of the nuclear attraction integrals core(kl/ij) iS
944 : !> (SS/)=1, (SO/)=2, (OO/)=3, (PP/)=4
945 : !> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
946 : !> \param itype type of semi_empirical model
947 : !> extension to the original routine to compute qm/mm
948 : !> integrals
949 : !> \param se_taper ...
950 : !> \param se_int_control input parameters that control the calculation of SE
951 : !> integrals (shortrange, R3 residual, screening type)
952 : !> \param lgrad ...
953 : !> \par History
954 : !> 03.2007 created [tlaino]
955 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
956 : !> for computing integrals
957 : !> \author Teodoro Laino - Zurich University
958 : !> \note
959 : !> Analytical version - Analytical evaluation of gradients
960 : !> Teodoro Laino - Zurich University 04.2007
961 : !>
962 : ! **************************************************************************************************
963 16142300 : SUBROUTINE dcore_nucint_ana(sepi, sepj, rij, core, dcore, itype, se_taper, &
964 : se_int_control, lgrad)
965 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
966 : REAL(dp), INTENT(IN) :: rij
967 : REAL(dp), DIMENSION(10, 2), INTENT(OUT) :: core, dcore
968 : INTEGER, INTENT(IN) :: itype
969 : TYPE(se_taper_type), POINTER :: se_taper
970 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
971 : LOGICAL, INTENT(IN) :: lgrad
972 :
973 : INTEGER :: i
974 : REAL(KIND=dp) :: dft, ft
975 : TYPE(se_int_screen_type) :: se_int_screen
976 :
977 : ! Compute the Tapering function
978 :
979 16142300 : ft = 1.0_dp
980 16142300 : dft = 0.0_dp
981 16142300 : IF (itype /= do_method_pchg) THEN
982 798403 : ft = taper_eval(se_taper%taper, rij)
983 798403 : dft = dtaper_eval(se_taper%taper, rij)
984 : END IF
985 : ! Evaluate additional taper function for dumped integrals
986 16142300 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
987 260198 : se_int_screen%ft = 1.0_dp
988 260198 : se_int_screen%dft = 0.0_dp
989 260198 : IF (itype /= do_method_pchg) THEN
990 260198 : se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
991 260198 : se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
992 : END IF
993 : END IF
994 :
995 : ! Value of the integrals for sp shell
996 : CALL nucint_sp_num(sepi, sepj, rij, core=core, itype=itype, &
997 16142300 : se_int_control=se_int_control, se_int_screen=se_int_screen)
998 :
999 16142300 : IF (sepi%dorb .OR. sepj%dorb) THEN
1000 : ! Compute the contribution from d-orbitals
1001 : CALL nucint_d_num(sepi, sepj, rij, core, itype, &
1002 37922 : se_int_control=se_int_control, se_int_screen=se_int_screen)
1003 : END IF
1004 :
1005 16142300 : IF (lgrad) THEN
1006 : ! Integrals derivatives for sp shell
1007 : CALL dnucint_sp_ana(sepi, sepj, rij, dcore=dcore, itype=itype, &
1008 8014152 : se_int_control=se_int_control, se_int_screen=se_int_screen)
1009 :
1010 8014152 : IF (sepi%dorb .OR. sepj%dorb) THEN
1011 : ! Integral derivatives involving d-orbitals
1012 : CALL dnucint_d_ana(sepi, sepj, rij, dcore=dcore, itype=itype, &
1013 17442 : se_int_control=se_int_control, se_int_screen=se_int_screen)
1014 : END IF
1015 : END IF
1016 :
1017 : ! Tapering the value of the integrals
1018 16142300 : IF (lgrad) THEN
1019 26875194 : DO i = 1, sepi%core_size
1020 26875194 : dcore(i, 1) = ft*dcore(i, 1) + dft*core(i, 1)
1021 : END DO
1022 8957533 : DO i = 1, sepj%core_size
1023 8957533 : dcore(i, 2) = ft*dcore(i, 2) + dft*core(i, 2)
1024 : END DO
1025 : END IF
1026 54133615 : DO i = 1, sepi%core_size
1027 54133615 : core(i, 1) = ft*core(i, 1)
1028 : END DO
1029 18297519 : DO i = 1, sepj%core_size
1030 18297519 : core(i, 2) = ft*core(i, 2)
1031 : END DO
1032 :
1033 : ! Debug Procedure.. Check valifity of analytical gradients of nucint
1034 : IF (debug_this_module .AND. lgrad) THEN
1035 : CALL check_dcore_nucint_ana(sepi, sepj, rij, dcore, itype, se_int_control, se_taper=se_taper)
1036 : END IF
1037 16142300 : END SUBROUTINE dcore_nucint_ana
1038 :
1039 : ! **************************************************************************************************
1040 : !> \brief Calculates the nuclear attraction integrals and derivatives for sp basis
1041 : !> \param sepi parameters of atom i
1042 : !> \param sepj parameters of atom j
1043 : !> \param rij interatomic distance
1044 : !> \param dssss derivative of (ssss) integral
1045 : !> derivatives are intended w.r.t. rij
1046 : !> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
1047 : !> \param dcore derivative of 4 X 2 array of electron-core attraction integrals
1048 : !> The storage of the nuclear attraction integrals core(kl/ij) iS
1049 : !> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5,
1050 : !> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10
1051 : !> \param itype type of semi_empirical model
1052 : !> extension to the original routine to compute qm/mm
1053 : !> integrals
1054 : !> \param se_int_control input parameters that control the calculation of SE
1055 : !> integrals (shortrange, R3 residual, screening type)
1056 : !> \param se_int_screen ...
1057 : !> \par History
1058 : !> 04.2007 created [tlaino]
1059 : !> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
1060 : !> for computing integrals
1061 : !> 05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting
1062 : !> \author Teodoro Laino - Zurich University
1063 : !> \note
1064 : !> Analytical version - Analytical evaluation of gradients
1065 : !> Teodoro Laino - Zurich University 04.2007
1066 : !> routine adapted from mopac7 (repp)
1067 : !> vector version written by Ernest R. Davidson, Indiana University
1068 : ! **************************************************************************************************
1069 16153257 : SUBROUTINE dnucint_sp_ana(sepi, sepj, rij, dssss, dcore, itype, se_int_control, &
1070 : se_int_screen)
1071 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1072 : REAL(dp), INTENT(IN) :: rij
1073 : REAL(dp), INTENT(INOUT), OPTIONAL :: dssss
1074 : REAL(dp), DIMENSION(10, 2), INTENT(INOUT), &
1075 : OPTIONAL :: dcore
1076 : INTEGER, INTENT(IN) :: itype
1077 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1078 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1079 :
1080 : INTEGER :: ij, kl
1081 : LOGICAL :: l_core, l_ssss, si, sj
1082 :
1083 16153257 : l_core = PRESENT(dcore)
1084 16153257 : l_ssss = PRESENT(dssss)
1085 16153257 : IF (.NOT. (l_core .OR. l_ssss)) RETURN
1086 :
1087 16153257 : si = (sepi%natorb > 1)
1088 16153257 : sj = (sepj%natorb > 1)
1089 :
1090 16153257 : ij = indexa(1, 1)
1091 16153257 : IF (l_ssss) THEN
1092 : ! Store the value for the derivative of <S S | S S > (Used for computing the core-core interactions)
1093 8139105 : dssss = d_ijkl_sp(sepi, sepj, ij, ij, 0, 0, 0, 0, -1, rij, se_int_control, se_int_screen, itype)
1094 : END IF
1095 :
1096 16153257 : IF (l_core) THEN
1097 : ! <S S | S S >
1098 8014152 : kl = indexa(1, 1)
1099 8014152 : dcore(1, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1100 8014152 : IF (si) THEN
1101 : ! <S P | S S >
1102 3595276 : kl = indexa(2, 1)
1103 3595276 : dcore(2, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 0, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1104 : ! <P P | S S >
1105 3595276 : kl = indexa(2, 2)
1106 3595276 : dcore(3, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1107 : ! <P+ P+ | S S >
1108 3595276 : kl = indexa(3, 3)
1109 3595276 : dcore(4, 1) = d_ijkl_sp(sepi, sepj, kl, ij, 1, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1110 : END IF
1111 :
1112 : ! <S S | S S >
1113 8014152 : kl = indexa(1, 1)
1114 8014152 : dcore(1, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1115 8014152 : IF (sj) THEN
1116 : ! <S S | S P >
1117 180523 : kl = indexa(2, 1)
1118 180523 : dcore(2, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 0, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1119 : ! <S S | P P >
1120 180523 : kl = indexa(2, 2)
1121 180523 : dcore(3, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1122 : ! <S S | P+ P+ >
1123 180523 : kl = indexa(3, 3)
1124 180523 : dcore(4, 2) = d_ijkl_sp(sepi, sepj, ij, kl, 0, 0, 1, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1125 : END IF
1126 : END IF
1127 : END SUBROUTINE dnucint_sp_ana
1128 :
1129 : ! **************************************************************************************************
1130 : !> \brief Calculates the analytical derivative of the nuclear attraction
1131 : !> integrals involving d orbitals
1132 : !> \param sepi parameters of atom i
1133 : !> \param sepj parameters of atom j
1134 : !> \param rij interatomic distance
1135 : !> \param dcore 4 X 2 array of electron-core attraction integrals
1136 : !> The storage of the nuclear attraction integrals core(kl/ij) iS
1137 : !> (SS/)=1, (SP/)=2, (PP/)=3, (P+P+/)=4, (SD/)=5,
1138 : !> (DP/)=6, (DD/)=7, (D+P+)=8, (D+D+/)=9, (D#D#)=10
1139 : !>
1140 : !> where ij=1 if the orbitals centred on atom i, =2 if on atom j.
1141 : !> \param itype type of semi_empirical model
1142 : !> extension to the original routine to compute qm/mm
1143 : !> integrals
1144 : !> \param se_int_control input parameters that control the calculation of SE
1145 : !> integrals (shortrange, R3 residual, screening type)
1146 : !> \param se_int_screen ...
1147 : !> \author
1148 : !> Teodoro Laino (05.2008) [tlaino] - University of Zurich: created
1149 : ! **************************************************************************************************
1150 17442 : SUBROUTINE dnucint_d_ana(sepi, sepj, rij, dcore, itype, se_int_control, &
1151 : se_int_screen)
1152 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1153 : REAL(dp), INTENT(IN) :: rij
1154 : REAL(dp), DIMENSION(10, 2), INTENT(INOUT) :: dcore
1155 : INTEGER, INTENT(IN) :: itype
1156 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1157 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1158 :
1159 : INTEGER :: ij, kl
1160 :
1161 : ! Check if d-orbitals are present
1162 :
1163 17442 : IF (sepi%dorb .OR. sepj%dorb) THEN
1164 17442 : ij = indexa(1, 1)
1165 17442 : IF (sepj%dorb) THEN
1166 : ! <S S | D S>
1167 9880 : kl = indexa(5, 1)
1168 9880 : dcore(5, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 0, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1169 : ! <S S | D P >
1170 9880 : kl = indexa(5, 2)
1171 9880 : dcore(6, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1172 : ! <S S | D D >
1173 9880 : kl = indexa(5, 5)
1174 9880 : dcore(7, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1175 : ! <S S | D+P+>
1176 9880 : kl = indexa(6, 3)
1177 9880 : dcore(8, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 1, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1178 : ! <S S | D+D+>
1179 9880 : kl = indexa(6, 6)
1180 9880 : dcore(9, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1181 : ! <S S | D#D#>
1182 9880 : kl = indexa(8, 8)
1183 9880 : dcore(10, 2) = d_ijkl_d(sepi, sepj, ij, kl, 0, 0, 2, 2, 1, rij, se_int_control, se_int_screen, itype)*sepi%zeff
1184 : END IF
1185 17442 : IF (sepi%dorb) THEN
1186 : ! <D S | S S>
1187 10177 : kl = indexa(5, 1)
1188 10177 : dcore(5, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 0, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1189 : ! <D P | S S >
1190 10177 : kl = indexa(5, 2)
1191 10177 : dcore(6, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1192 : ! <D D | S S >
1193 10177 : kl = indexa(5, 5)
1194 10177 : dcore(7, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1195 : ! <D+P+| S S >
1196 10177 : kl = indexa(6, 3)
1197 10177 : dcore(8, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 1, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1198 : ! <D+D+| S S >
1199 10177 : kl = indexa(6, 6)
1200 10177 : dcore(9, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1201 : ! <D#D#| S S >
1202 10177 : kl = indexa(8, 8)
1203 10177 : dcore(10, 1) = d_ijkl_d(sepi, sepj, kl, ij, 2, 2, 0, 0, 2, rij, se_int_control, se_int_screen, itype)*sepj%zeff
1204 : END IF
1205 : END IF
1206 17442 : END SUBROUTINE dnucint_d_ana
1207 :
1208 : ! **************************************************************************************************
1209 : !> \brief calculates the derivative of the two-particle interactions
1210 : !> \param sepi Atomic parameters of first atom
1211 : !> \param sepj Atomic parameters of second atom
1212 : !> \param rijv Coordinate vector i -> j
1213 : !> \param w Array of two-electron repulsion integrals.
1214 : !> \param dw ...
1215 : !> \param se_int_control ...
1216 : !> \param se_taper ...
1217 : !> \par History
1218 : !> 04.2007 created [tlaino]
1219 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : new driver
1220 : !> for computing integrals
1221 : !> \author Teodoro Laino - Zurich University
1222 : !> \note
1223 : !> Analytical version - Analytical evaluation of gradients
1224 : !> Teodoro Laino - Zurich University 04.2007
1225 : !> routine adapted from mopac7 (repp)
1226 : !> vector version written by Ernest R. Davidson, Indiana University
1227 : ! **************************************************************************************************
1228 1144109 : RECURSIVE SUBROUTINE rotint_ana(sepi, sepj, rijv, w, dw, se_int_control, se_taper)
1229 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1230 : REAL(dp), DIMENSION(3), INTENT(IN) :: rijv
1231 : REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL :: w
1232 : REAL(dp), DIMENSION(3, 2025), INTENT(OUT), &
1233 : OPTIONAL :: dw
1234 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1235 : TYPE(se_taper_type), POINTER :: se_taper
1236 :
1237 : INTEGER :: i, i1, ii, ij, ij1, iminus, istep, &
1238 : iw_loc, j, j1, jj, k, kk, kl, l, &
1239 : limij, limkl, mm
1240 : LOGICAL :: invert, l_w, lgrad
1241 : LOGICAL, DIMENSION(45, 45) :: logv, logv_d
1242 : REAL(dp) :: rij, xtmp
1243 : REAL(dp), DIMENSION(3) :: drij
1244 : REAL(KIND=dp) :: cc, cc_d(3), wrepp, wrepp_d(3)
1245 : REAL(KIND=dp), DIMENSION(2025) :: ww
1246 : REAL(KIND=dp), DIMENSION(3, 2025) :: ww_d
1247 : REAL(KIND=dp), DIMENSION(3, 45, 45) :: v_d
1248 : REAL(KIND=dp), DIMENSION(45, 45) :: v
1249 : REAL(KIND=dp), DIMENSION(491) :: rep, rep_d
1250 : TYPE(rotmat_type), POINTER :: ij_matrix
1251 :
1252 1144109 : NULLIFY (ij_matrix)
1253 1144109 : l_w = PRESENT(w)
1254 1144109 : lgrad = PRESENT(dw)
1255 1144109 : IF (.NOT. (l_w .OR. lgrad)) RETURN
1256 :
1257 4576436 : rij = DOT_PRODUCT(rijv, rijv)
1258 1144109 : IF (rij > rij_threshold) THEN
1259 : ! The repulsion integrals over molecular frame (w) are stored in the
1260 : ! order in which they will later be used. ie. (i,j/k,l) where
1261 : ! j.le.i and l.le.k and l varies most rapidly and i least
1262 : ! rapidly. (anti-normal computer storage)
1263 1144109 : rij = SQRT(rij)
1264 :
1265 : ! Create the rotation matrix
1266 1144109 : CALL rotmat_create(ij_matrix)
1267 1144109 : CALL rotmat(sepi, sepj, rijv, rij, ij_matrix, do_derivatives=lgrad, do_invert=invert)
1268 :
1269 : ! Compute integrals in diatomic frame as well their derivatives (if requested)
1270 1144109 : CALL dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, se_int_control, lgrad=lgrad)
1271 :
1272 1144109 : IF (lgrad) THEN
1273 511306 : drij(1) = rijv(1)/rij
1274 511306 : drij(2) = rijv(2)/rij
1275 511306 : drij(3) = rijv(3)/rij
1276 : ! Possibly Invert Frame
1277 511306 : IF (invert) THEN
1278 947 : xtmp = drij(3)
1279 947 : drij(3) = drij(1)
1280 947 : drij(1) = xtmp
1281 : END IF
1282 : END IF
1283 :
1284 1144109 : ii = sepi%natorb
1285 1144109 : kk = sepj%natorb
1286 : ! First step in rotation of integrals
1287 : CALL rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, rep, logv, ij_matrix, &
1288 1144109 : v, lgrad, rep_d, v_d, logv_d, drij)
1289 :
1290 : ! Integrals if requested
1291 1144109 : IF (l_w) THEN
1292 : ! Rotate Integrals
1293 632803 : IF (ii*kk > 0) THEN
1294 632803 : limij = sepi%atm_int_size
1295 632803 : limkl = sepj%atm_int_size
1296 632803 : istep = limkl*limij
1297 35148717 : DO i1 = 1, istep
1298 35148717 : ww(i1) = 0.0_dp
1299 : END DO
1300 : ! Second step in rotation of integrals
1301 2456506 : DO i1 = 1, ii
1302 7017769 : DO j1 = 1, i1
1303 4561263 : ij = indexa(i1, j1)
1304 4561263 : jj = indexb(i1, j1)
1305 4561263 : mm = int2c_type(jj)
1306 19448726 : DO k = 1, kk
1307 52140937 : DO l = 1, k
1308 34515914 : kl = indexb(k, l)
1309 47579674 : IF (logv(ij, kl)) THEN
1310 21417646 : wrepp = v(ij, kl)
1311 3259225 : SELECT CASE (mm)
1312 : CASE (1)
1313 : ! (SS/)
1314 3259225 : i = 1
1315 3259225 : j = 1
1316 3259225 : iw_loc = (indexb(i, j) - 1)*limkl + kl
1317 3259225 : ww(iw_loc) = wrepp
1318 : CASE (2)
1319 : ! (SP/)
1320 20426976 : j = 1
1321 20426976 : DO i = 1, 3
1322 15320232 : iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
1323 20426976 : ww(iw_loc) = ww(iw_loc) + ij_matrix%sp(i1 - 1, i)*wrepp
1324 : END DO
1325 : CASE (3)
1326 : ! (PP/)
1327 40700512 : DO i = 1, 3
1328 30525384 : cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
1329 30525384 : iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
1330 30525384 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
1331 30525384 : iminus = i - 1
1332 40700512 : IF (iminus /= 0) THEN
1333 50875640 : DO j = 1, iminus
1334 30525384 : cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
1335 30525384 : iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
1336 50875640 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
1337 : END DO
1338 : END IF
1339 : END DO
1340 : CASE (4)
1341 : ! (SD/)
1342 2577996 : j = 1
1343 2577996 : DO i = 1, 5
1344 2148330 : iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
1345 2577996 : ww(iw_loc) = ww(iw_loc) + ij_matrix%sd(i1 - 4, i)*wrepp
1346 : END DO
1347 : CASE (5)
1348 : ! (DP/)
1349 6637362 : DO i = 1, 5
1350 23230767 : DO j = 1, 3
1351 16593405 : iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
1352 16593405 : ij1 = 3*(i - 1) + j
1353 22124540 : ww(iw_loc) = ww(iw_loc) + ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp
1354 : END DO
1355 : END DO
1356 : CASE (6)
1357 : ! (DD/)
1358 28120926 : DO i = 1, 5
1359 6703280 : cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
1360 6703280 : iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
1361 6703280 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
1362 6703280 : iminus = i - 1
1363 8043936 : IF (iminus /= 0) THEN
1364 18769184 : DO j = 1, iminus
1365 13406560 : ij1 = inddd(i, j)
1366 13406560 : cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
1367 13406560 : iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
1368 18769184 : ww(iw_loc) = ww(iw_loc) + cc*wrepp
1369 : END DO
1370 : END IF
1371 : END DO
1372 : END SELECT
1373 : END IF
1374 : END DO
1375 : END DO
1376 : END DO
1377 : END DO
1378 : ! Store two electron integrals in the triangular format
1379 632803 : CALL store_2el_2c_diag(limij, limkl, ww(1:istep), w)
1380 632803 : IF (invert) CALL invert_integral(sepi, sepj, int2el=w)
1381 : END IF
1382 :
1383 : IF (debug_this_module) THEN
1384 : ! Check value of integrals
1385 : CALL check_rotint_ana(sepi, sepj, rijv, w, se_int_control=se_int_control, se_taper=se_taper)
1386 : END IF
1387 : END IF
1388 :
1389 : ! Gradients if requested
1390 1144109 : IF (lgrad) THEN
1391 : ! Rotate Integrals derivatives
1392 511306 : IF (ii*kk > 0) THEN
1393 511306 : limij = sepi%atm_int_size
1394 511306 : limkl = sepj%atm_int_size
1395 511306 : istep = limkl*limij
1396 30915998 : DO i1 = 1, istep
1397 30404692 : ww_d(1, i1) = 0.0_dp
1398 30404692 : ww_d(2, i1) = 0.0_dp
1399 30915998 : ww_d(3, i1) = 0.0_dp
1400 : END DO
1401 :
1402 : ! Second step in rotation of integrals
1403 2018237 : DO i1 = 1, ii
1404 5821698 : DO j1 = 1, i1
1405 3803461 : ij = indexa(i1, j1)
1406 3803461 : jj = indexb(i1, j1)
1407 3803461 : mm = int2c_type(jj)
1408 16580730 : DO k = 1, kk
1409 45478491 : DO l = 1, k
1410 30404692 : kl = indexb(k, l)
1411 41675030 : IF (logv_d(ij, kl)) THEN
1412 19007788 : wrepp_d(1) = v_d(1, ij, kl)
1413 19007788 : wrepp_d(2) = v_d(2, ij, kl)
1414 19007788 : wrepp_d(3) = v_d(3, ij, kl)
1415 19007788 : wrepp = v(ij, kl)
1416 2769357 : SELECT CASE (mm)
1417 : CASE (1)
1418 : ! (SS/)
1419 2769357 : i = 1
1420 2769357 : j = 1
1421 2769357 : iw_loc = (indexb(i, j) - 1)*limkl + kl
1422 2769357 : ww_d(1, iw_loc) = wrepp_d(1)
1423 2769357 : ww_d(2, iw_loc) = wrepp_d(2)
1424 2769357 : ww_d(3, iw_loc) = wrepp_d(3)
1425 : CASE (2)
1426 : ! (SP/)
1427 18275868 : j = 1
1428 18275868 : DO i = 1, 3
1429 13706901 : iw_loc = (indexb(i + 1, j) - 1)*limkl + kl
1430 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sp_d(1, i1 - 1, i)*wrepp + &
1431 13706901 : ij_matrix%sp(i1 - 1, i)*wrepp_d(1)
1432 :
1433 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sp_d(2, i1 - 1, i)*wrepp + &
1434 13706901 : ij_matrix%sp(i1 - 1, i)*wrepp_d(2)
1435 :
1436 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sp_d(3, i1 - 1, i)*wrepp + &
1437 18275868 : ij_matrix%sp(i1 - 1, i)*wrepp_d(3)
1438 : END DO
1439 : CASE (3)
1440 : ! (PP/)
1441 36203184 : DO i = 1, 3
1442 27152388 : cc = ij_matrix%pp(i, i1 - 1, j1 - 1)
1443 27152388 : cc_d(1) = ij_matrix%pp_d(1, i, i1 - 1, j1 - 1)
1444 27152388 : cc_d(2) = ij_matrix%pp_d(2, i, i1 - 1, j1 - 1)
1445 27152388 : cc_d(3) = ij_matrix%pp_d(3, i, i1 - 1, j1 - 1)
1446 27152388 : iw_loc = (indexb(i + 1, i + 1) - 1)*limkl + kl
1447 27152388 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1448 27152388 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1449 27152388 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1450 27152388 : iminus = i - 1
1451 36203184 : IF (iminus /= 0) THEN
1452 45253980 : DO j = 1, iminus
1453 27152388 : cc = ij_matrix%pp(1 + i + j, i1 - 1, j1 - 1)
1454 27152388 : cc_d(1) = ij_matrix%pp_d(1, 1 + i + j, i1 - 1, j1 - 1)
1455 27152388 : cc_d(2) = ij_matrix%pp_d(2, 1 + i + j, i1 - 1, j1 - 1)
1456 27152388 : cc_d(3) = ij_matrix%pp_d(3, 1 + i + j, i1 - 1, j1 - 1)
1457 27152388 : iw_loc = (indexb(i + 1, j + 1) - 1)*limkl + kl
1458 27152388 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1459 27152388 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1460 45253980 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1461 : END DO
1462 : END IF
1463 : END DO
1464 : CASE (4)
1465 : ! (SD/)
1466 2351124 : j = 1
1467 2351124 : DO i = 1, 5
1468 1959270 : iw_loc = (indexb(i + 4, j) - 1)*limkl + kl
1469 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%sd_d(1, i1 - 4, i)*wrepp + &
1470 1959270 : ij_matrix%sd(i1 - 4, i)*wrepp_d(1)
1471 :
1472 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%sd_d(2, i1 - 4, i)*wrepp + &
1473 1959270 : ij_matrix%sd(i1 - 4, i)*wrepp_d(2)
1474 :
1475 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%sd_d(3, i1 - 4, i)*wrepp + &
1476 2351124 : ij_matrix%sd(i1 - 4, i)*wrepp_d(3)
1477 : END DO
1478 : CASE (5)
1479 : ! (DP/)
1480 6077760 : DO i = 1, 5
1481 21272160 : DO j = 1, 3
1482 15194400 : iw_loc = (indexb(i + 4, j + 1) - 1)*limkl + kl
1483 15194400 : ij1 = 3*(i - 1) + j
1484 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + ij_matrix%pd_d(1, ij1, i1 - 4, j1 - 1)*wrepp + &
1485 15194400 : ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(1)
1486 :
1487 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + ij_matrix%pd_d(2, ij1, i1 - 4, j1 - 1)*wrepp + &
1488 15194400 : ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(2)
1489 :
1490 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + ij_matrix%pd_d(3, ij1, i1 - 4, j1 - 1)*wrepp + &
1491 20259200 : ij_matrix%pd(ij1, i1 - 4, j1 - 1)*wrepp_d(3)
1492 : END DO
1493 : END DO
1494 : CASE (6)
1495 : ! (DD/)
1496 25077058 : DO i = 1, 5
1497 6069270 : cc = ij_matrix%dd(i, i1 - 4, j1 - 4)
1498 24277080 : cc_d = ij_matrix%dd_d(:, i, i1 - 4, j1 - 4)
1499 6069270 : iw_loc = (indexb(i + 4, i + 4) - 1)*limkl + kl
1500 6069270 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1501 6069270 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1502 6069270 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1503 6069270 : iminus = i - 1
1504 7283124 : IF (iminus /= 0) THEN
1505 16993956 : DO j = 1, iminus
1506 12138540 : ij1 = inddd(i, j)
1507 12138540 : cc = ij_matrix%dd(ij1, i1 - 4, j1 - 4)
1508 12138540 : cc_d(1) = ij_matrix%dd_d(1, ij1, i1 - 4, j1 - 4)
1509 12138540 : cc_d(2) = ij_matrix%dd_d(2, ij1, i1 - 4, j1 - 4)
1510 12138540 : cc_d(3) = ij_matrix%dd_d(3, ij1, i1 - 4, j1 - 4)
1511 12138540 : iw_loc = (indexb(i + 4, j + 4) - 1)*limkl + kl
1512 12138540 : ww_d(1, iw_loc) = ww_d(1, iw_loc) + cc_d(1)*wrepp + cc*wrepp_d(1)
1513 12138540 : ww_d(2, iw_loc) = ww_d(2, iw_loc) + cc_d(2)*wrepp + cc*wrepp_d(2)
1514 16993956 : ww_d(3, iw_loc) = ww_d(3, iw_loc) + cc_d(3)*wrepp + cc*wrepp_d(3)
1515 : END DO
1516 : END IF
1517 : END DO
1518 : END SELECT
1519 : END IF
1520 : END DO
1521 : END DO
1522 : END DO
1523 : END DO
1524 : ! Store two electron integrals in the triangular format
1525 : CALL store_2el_2c_diag(limij, limkl, ww_dx=ww_d(1, 1:istep), ww_dy=ww_d(2, 1:istep), ww_dz=ww_d(3, 1:istep), &
1526 91725382 : dw=dw)
1527 511306 : IF (invert) CALL invert_derivative(sepi, sepj, dint2el=dw)
1528 : END IF
1529 :
1530 : IF (debug_this_module) THEN
1531 : ! Check derivatives
1532 : CALL check_rotint_ana(sepi, sepj, rijv, dw=dw, se_int_control=se_int_control, se_taper=se_taper)
1533 : END IF
1534 : END IF
1535 1144109 : CALL rotmat_release(ij_matrix)
1536 : END IF
1537 : END SUBROUTINE rotint_ana
1538 :
1539 : ! **************************************************************************************************
1540 : !> \brief Calculates the derivative and the value of two-electron repulsion
1541 : !> integrals and the nuclear attraction integrals w.r.t. |r|
1542 : !> \param sepi parameters of atom i
1543 : !> \param sepj parameters of atom j
1544 : !> \param rij interatomic distance
1545 : !> \param rep rray of two-electron repulsion integrals
1546 : !> \param rep_d array of two-electron repulsion integrals derivatives
1547 : !> \param se_taper ...
1548 : !> \param se_int_control input parameters that control the calculation of SE
1549 : !> integrals (shortrange, R3 residual, screening type)
1550 : !> \param lgrad ...
1551 : !> \par History
1552 : !> 03.2008 created [tlaino]
1553 : !> \author Teodoro Laino [tlaino] - Zurich University
1554 : ! **************************************************************************************************
1555 1144109 : RECURSIVE SUBROUTINE dterep_ana(sepi, sepj, rij, rep, rep_d, se_taper, &
1556 : se_int_control, lgrad)
1557 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1558 : REAL(KIND=dp), INTENT(IN) :: rij
1559 : REAL(KIND=dp), DIMENSION(491), INTENT(OUT) :: rep, rep_d
1560 : TYPE(se_taper_type), POINTER :: se_taper
1561 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1562 : LOGICAL, INTENT(IN) :: lgrad
1563 :
1564 : INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
1565 : lj, lk, ll, numb
1566 : REAL(KIND=dp) :: dft, ft, ft1
1567 : TYPE(se_int_screen_type) :: se_int_screen
1568 :
1569 : ! Compute the tapering function and its derivatives
1570 :
1571 1144109 : ft = taper_eval(se_taper%taper, rij)
1572 1144109 : dft = 0.0_dp
1573 1144109 : ft1 = ft
1574 1144109 : IF (lgrad) THEN
1575 511306 : ft1 = 1.0_dp
1576 511306 : dft = dtaper_eval(se_taper%taper, rij)
1577 : END IF
1578 : ! Evaluate additional taper function for dumped integrals
1579 1144109 : IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN
1580 260198 : se_int_screen%ft = taper_eval(se_taper%taper_add, rij)
1581 260198 : IF (lgrad) &
1582 130072 : se_int_screen%dft = dtaper_eval(se_taper%taper_add, rij)
1583 : END IF
1584 :
1585 : ! Integral Values for sp shells only
1586 : CALL terep_sp_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
1587 1144109 : se_int_screen=se_int_screen, ft=ft1)
1588 :
1589 1144109 : IF (sepi%dorb .OR. sepj%dorb) THEN
1590 : ! Compute the contribution from d-orbitals
1591 : CALL terep_d_num(sepi, sepj, rij, rep, se_int_control=se_int_control, &
1592 54500 : se_int_screen=se_int_screen, ft=ft1)
1593 : END IF
1594 :
1595 1144109 : IF (lgrad) THEN
1596 : ! Integral Derivatives
1597 : CALL dterep_sp_ana(sepi, sepj, rij, rep_d, rep, se_int_control, &
1598 511306 : se_int_screen, ft, dft)
1599 :
1600 511306 : IF (sepi%dorb .OR. sepj%dorb) THEN
1601 : ! Compute the derivatives from d-orbitals
1602 : CALL dterep_d_ana(sepi, sepj, rij, rep_d, rep, se_int_control, &
1603 25368 : se_int_screen, ft, dft)
1604 : END IF
1605 :
1606 : ! Tapering Integral values
1607 511306 : lasti = sepi%natorb
1608 511306 : lastj = sepj%natorb
1609 2018237 : DO i = 1, lasti
1610 1506931 : li = l_index(i)
1611 5821698 : DO j = 1, i
1612 3803461 : lj = l_index(j)
1613 3803461 : ij = indexa(i, j)
1614 16580730 : DO k = 1, lastj
1615 11270338 : lk = l_index(k)
1616 45478491 : DO l = 1, k
1617 30404692 : ll = l_index(l)
1618 30404692 : kl = indexa(k, l)
1619 30404692 : numb = ijkl_ind(ij, kl)
1620 41675030 : IF (numb > 0) rep(numb) = rep(numb)*ft
1621 : END DO
1622 : END DO
1623 : END DO
1624 : END DO
1625 : END IF
1626 :
1627 : ! Possibly debug 2el 2cent integrals and derivatives
1628 : IF (debug_this_module) THEN
1629 : CALL check_dterep_ana(sepi, sepj, rij, rep, rep_d, se_int_control, se_taper=se_taper, &
1630 : lgrad=lgrad)
1631 : END IF
1632 1144109 : END SUBROUTINE dterep_ana
1633 :
1634 : ! **************************************************************************************************
1635 : !> \brief Calculates the derivative and the value of two-electron repulsion
1636 : !> integrals and the nuclear attraction integrals w.r.t. |r| - sp shells only
1637 : !> \param sepi parameters of atom i
1638 : !> \param sepj parameters of atom j
1639 : !> \param rij interatomic distance
1640 : !> \param drep array of derivatives of two-electron repulsion integrals
1641 : !> \param rep array of two-electron repulsion integrals
1642 : !> \param se_int_control input parameters that control the calculation of SE
1643 : !> integrals (shortrange, R3 residual, screening type)
1644 : !> \param se_int_screen ...
1645 : !> \param ft ...
1646 : !> \param dft ...
1647 : !> \par History
1648 : !> 04.2007 created [tlaino]
1649 : !> Teodoro Laino (03.2008) [tlaino] - University of Zurich : new driver
1650 : !> for computing integrals
1651 : !> 05.2008 Teodoro Laino [tlaino] - University of Zurich: major rewriting
1652 : !> \author Teodoro Laino - Zurich University
1653 : !> \note
1654 : !> Analytical version - Analytical evaluation of gradients
1655 : !> Teodoro Laino - Zurich University 04.2007
1656 : !> routine adapted from mopac7 (repp)
1657 : !> vector version written by Ernest R. Davidson, Indiana University
1658 : ! **************************************************************************************************
1659 511306 : SUBROUTINE dterep_sp_ana(sepi, sepj, rij, drep, rep, se_int_control, &
1660 : se_int_screen, ft, dft)
1661 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1662 : REAL(dp), INTENT(IN) :: rij
1663 : REAL(dp), DIMENSION(491), INTENT(OUT) :: drep
1664 : REAL(dp), DIMENSION(491), INTENT(IN) :: rep
1665 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1666 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1667 : REAL(dp), INTENT(IN) :: ft, dft
1668 :
1669 : INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
1670 : lj, lk, ll, nold, numb
1671 : REAL(KIND=dp) :: tmp
1672 :
1673 511306 : lasti = sepi%natorb
1674 511306 : lastj = sepj%natorb
1675 1941917 : DO i = 1, MIN(lasti, 4)
1676 1430611 : li = l_index(i)
1677 5211138 : DO j = 1, i
1678 3269221 : lj = l_index(j)
1679 3269221 : ij = indexa(i, j)
1680 13175340 : DO k = 1, MIN(lastj, 4)
1681 8475508 : lk = l_index(k)
1682 30632811 : DO l = 1, k
1683 18888082 : ll = l_index(l)
1684 18888082 : kl = indexa(k, l)
1685 :
1686 18888082 : numb = ijkl_ind(ij, kl)
1687 27363590 : IF (numb > 0) THEN
1688 6851959 : nold = ijkl_sym(numb)
1689 6851959 : IF (nold > 0) THEN
1690 2207937 : drep(numb) = drep(nold)
1691 4644022 : ELSE IF (nold < 0) THEN
1692 0 : drep(numb) = -drep(-nold)
1693 4644022 : ELSE IF (nold == 0) THEN
1694 : tmp = d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
1695 4644022 : se_int_control, se_int_screen, do_method_undef)
1696 4644022 : drep(numb) = dft*rep(numb) + ft*tmp
1697 : END IF
1698 : END IF
1699 : END DO
1700 : END DO
1701 : END DO
1702 : END DO
1703 511306 : END SUBROUTINE dterep_sp_ana
1704 :
1705 : ! **************************************************************************************************
1706 : !> \brief Calculates the derivatives of the two-electron repulsion integrals - d shell only
1707 : !> \param sepi parameters of atom i
1708 : !> \param sepj parameters of atom j
1709 : !> \param rij interatomic distance
1710 : !> \param drep ...
1711 : !> \param rep array of two-electron repulsion integrals
1712 : !> \param se_int_control input parameters that control the calculation of
1713 : !> integrals (shortrange, R3 residual, screening type)
1714 : !> \param se_int_screen ...
1715 : !> \param ft ...
1716 : !> \param dft ...
1717 : !> \par History
1718 : !> Teodoro Laino (05.2008) [tlaino] - University of Zurich : new driver
1719 : !> for computing integral derivatives for d-orbitals
1720 : ! **************************************************************************************************
1721 25368 : SUBROUTINE dterep_d_ana(sepi, sepj, rij, drep, rep, se_int_control, &
1722 : se_int_screen, ft, dft)
1723 : TYPE(semi_empirical_type), POINTER :: sepi, sepj
1724 : REAL(dp), INTENT(IN) :: rij
1725 : REAL(dp), DIMENSION(491), INTENT(INOUT) :: drep
1726 : REAL(dp), DIMENSION(491), INTENT(IN) :: rep
1727 : TYPE(se_int_control_type), INTENT(IN) :: se_int_control
1728 : TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen
1729 : REAL(dp), INTENT(IN) :: ft, dft
1730 :
1731 : INTEGER :: i, ij, j, k, kl, l, lasti, lastj, li, &
1732 : lj, lk, ll, nold, numb
1733 : REAL(KIND=dp) :: tmp
1734 :
1735 25368 : lasti = sepi%natorb
1736 25368 : lastj = sepj%natorb
1737 196725 : DO i = 1, lasti
1738 171357 : li = l_index(i)
1739 965340 : DO j = 1, i
1740 768615 : lj = l_index(j)
1741 768615 : ij = indexa(i, j)
1742 4560222 : DO k = 1, lastj
1743 3620250 : lk = l_index(k)
1744 17912985 : DO l = 1, k
1745 13524120 : ll = l_index(l)
1746 13524120 : kl = indexa(k, l)
1747 :
1748 13524120 : numb = ijkl_ind(ij, kl)
1749 : ! From 1 to 34 we store integrals involving sp shells
1750 17144370 : IF (numb > 34) THEN
1751 2836440 : nold = ijkl_sym(numb)
1752 2836440 : IF (nold > 34) THEN
1753 1324558 : drep(numb) = drep(nold)
1754 1511882 : ELSE IF (nold < -34) THEN
1755 259002 : drep(numb) = -drep(-nold)
1756 1252880 : ELSE IF (nold == 0) THEN
1757 : tmp = d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, 0, rij, &
1758 1252880 : se_int_control, se_int_screen, do_method_undef)
1759 1252880 : drep(numb) = dft*rep(numb) + ft*tmp
1760 : END IF
1761 : END IF
1762 : END DO
1763 : END DO
1764 : END DO
1765 : END DO
1766 25368 : END SUBROUTINE dterep_d_ana
1767 :
1768 : END MODULE semi_empirical_int_ana
|
