(git:374b731)
Loading...
Searching...
No Matches
semi_empirical_par_utils.F
Go to the documentation of this file.
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 Utilities to post-process semi-empirical parameters
10!> \par History
11!> [tlaino] 03.2008 - Splitting from semi_empirical_parameters and
12!> keeping there only the setting of the SE params
13!> \author Teodoro Laino [tlaino] - University of Zurich
14!> \date 03.2008 [tlaino]
15! **************************************************************************************************
16MODULE semi_empirical_par_utils
17
18 USE kinds, ONLY: dp
19 USE mathconstants, ONLY: fac
20 USE mathlib, ONLY: binomial
21 USE physcon, ONLY: bohr,&
22 evolt
23 USE semi_empirical_int_arrays, ONLY: int_ij,&
24 int_kl,&
25 int_onec2el
26 USE semi_empirical_types, ONLY: get_se_param,&
27 semi_empirical_type
28#include "./base/base_uses.f90"
29
30 IMPLICIT NONE
31
32 PRIVATE
33
34 INTEGER, PARAMETER, PRIVATE :: nelem = 106
35
36 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_par_utils'
37
38! STANDARD MOPAC PARAMETERS USED FOR AM1, RM1, MNDO, PM3, PM6,
39! PM6-FM
40!
41! H He
42! Li Be B C N O F Ne
43! Na Mg Al Si P S Cl Ar
44! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
45! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe
46! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
47! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Ha 106
48
49! "s" shell
50 INTEGER, DIMENSION(0:nelem), PRIVATE :: Nos = (/-1, & ! 0
51 1, 2, & ! 2
52 1, 2, 2, 2, 2, 2, 2, 0, & ! 10
53 1, 2, 2, 2, 2, 2, 2, 0, & ! 18
54 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 0, & ! 36
55 1, 2, 2, 2, 1, 1, 2, 1, 1, 0, 1, 2, 2, 2, 2, 2, 2, 0, & ! 54
56 1, 2, 2, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, &
57 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 0, & ! 86
58 1, 1, 2, 4, 2, 2, 2, 2, 2, 2, 2, 2, 0, 3, -3, 1, 2, 1, -2, -1/)
59
60! "p" shell
61 INTEGER, DIMENSION(0:nelem), PRIVATE :: Nop = (/-1, & ! 0
62 0, 0, & ! 2
63 0, 0, 1, 2, 3, 4, 5, 6, & ! 10
64 0, 0, 1, 2, 3, 4, 5, 6, & ! 18
65 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & ! 36
66 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & ! 54
67 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
68 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 6, & ! 86
69 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
70
71! "d" shell
72 INTEGER, DIMENSION(0:nelem), PRIVATE :: Nod = (/-1, & ! 0
73 0, 0, & ! 2
74 0, 0, 0, 0, 0, 0, 0, 0, & ! 10
75 0, 0, 0, 0, 0, 0, 0, 0, & ! 18
76 0, 0, 1, 2, 3, 5, 5, 6, 7, 8, 10, 0, 0, 0, 0, 0, 0, 0, & ! 36
77 0, 0, 1, 2, 4, 5, 5, 7, 8, 10, 10, 0, 0, 0, 0, 0, 0, 0, & ! 54
78 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, &
79 2, 3, 5, 5, 6, 7, 9, 10, 0, 0, 0, 0, 0, 0, 0, & ! 86
80 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
81
82! H <Quantum Numbers for s, p, d and f orbitals> He
83! Li Be B C N O F Ne
84! Na Mg Al Si P S Cl Ar
85! K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
86! Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe
87! Cs Ba La Ce-Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
88! Fr Ra Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No Lr Rf Ha 106
89
90 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqs = (/-1, & ! 0
91 1, 1, & ! 2
92 2, 2, 2, 2, 2, 2, 2, 2, & ! 10
93 3, 3, 3, 3, 3, 3, 3, 3, & ! 18
94 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 36
95 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & ! 54
96 6, 6, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, &
97 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, -1, -1, -1, & ! 86
98 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
99
100 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqp = (/-1, & ! 0
101 -1, -1, & ! 2
102 2, 2, 2, 2, 2, 2, 2, 2, & ! 10
103 3, 3, 3, 3, 3, 3, 3, 3, & ! 18
104 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, & ! 36
105 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, & ! 54
106 6, 6, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 6, &
107 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, -1, -1, -1, & ! 86
108 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
109
110 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqd = (/-1, & ! 0
111 -1, -1, & ! 2
112 -1, -1, -1, -1, -1, -1, -1, -1, & ! 10
113 -1, -1, 3, 3, 3, 3, 3, -1, & ! 18
114 -1, -1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, -1, & ! 36
115 -1, -1, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, -1, & ! 54
116 -1, -1, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 5, &
117 5, 5, 5, 5, 5, 5, 5, 5, 5, -1, -1, -1, -1, -1, -1, & ! 86
118 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
119
120 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: nqf = (/-1, & ! 0
121 -1, -1, & ! 2
122 -1, -1, -1, -1, -1, -1, -1, -1, & ! 10
123 -1, -1, -1, -1, -1, -1, -1, -1, & ! 18
124 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! 36
125 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! 54
126 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, &
127 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, & ! 86
128 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
129
130 ! Element Valence
131 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: zval = (/-1, & ! 0
132 1, 2, & ! 2
133 1, 2, 3, 4, 5, 6, 7, 8, & ! 10
134 1, 2, 3, 4, 5, 6, 7, 8, & ! 18
135 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 36
136 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, 8, & ! 54
137 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 3, &
138 4, 5, 6, 7, 8, 9, 10, 11, 2, 3, 4, 5, 6, 7, -1, & ! 86
139 -1, -1, -1, 4, -1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1/)
140
141 ! Number of 1 center 2 electron integrals involving partially filled d shells
142 ! r016: <SS|DD>
143 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir016 = (/0, & ! 0
144 0, 0, & ! 2
145 0, 0, 0, 0, 0, 0, 0, 0, & ! 10
146 0, 0, 0, 0, 0, 0, 0, 0, & ! 18
147 0, 0, 2, 4, 6, 5, 10, 12, 14, 16, 10, 0, 0, 0, 0, 0, 0, 0, & ! 36
148 0, 0, 4, 4, 4, 5, 10, 7, 8, 0, 10, 0, 0, 0, 0, 0, 0, 0, & ! 54
149 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
150 4, 6, 8, 10, 12, 14, 9, 10, 0, 0, 0, 0, 0, 0, 0, & ! 86
151 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
152
153 ! r066: <DD|DD> "0" term
154 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir066 = (/0, & ! 0
155 0, 0, & ! 2
156 0, 0, 0, 0, 0, 0, 0, 0, & ! 10
157 0, 0, 0, 0, 0, 0, 0, 0, & ! 18
158 0, 0, 0, 1, 3, 10, 10, 15, 21, 28, 45, 0, 0, 0, 0, 0, 0, 0, & ! 36
159 0, 0, 0, 1, 6, 10, 10, 21, 28, 45, 45, 0, 0, 0, 0, 0, 0, 0, & ! 54
160 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
161 1, 3, 6, 10, 15, 21, 36, 45, 0, 0, 0, 0, 0, 0, 0, & ! 86
162 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
163
164 ! r244: <SD|SD>
165 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir244 = (/0, & ! 0
166 0, 0, & ! 2
167 0, 0, 0, 0, 0, 0, 0, 0, & ! 10
168 0, 0, 0, 0, 0, 0, 0, 0, & ! 18
169 0, 0, 1, 2, 3, 5, 5, 6, 7, 8, 5, 0, 0, 0, 0, 0, 0, 0, & ! 36
170 0, 0, 1, 2, 4, 5, 5, 5, 5, 0, 5, 0, 0, 0, 0, 0, 0, 0, & ! 54
171 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
172 2, 3, 4, 5, 6, 7, 5, 5, 0, 0, 0, 0, 0, 0, 0, & ! 86
173 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
174
175 ! r266: <DD|DD> "2" term
176 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir266 = (/0, & ! 0
177 0, 0, & ! 2
178 0, 0, 0, 0, 0, 0, 0, 0, & ! 10
179 0, 0, 0, 0, 0, 0, 0, 0, & ! 18
180 0, 0, 0, 8, 15, 35, 35, 35, 43, 50, 70, 0, 0, 0, 0, 0, 0, 0, & ! 36
181 0, 0, 0, 8, 21, 35, 35, 43, 50, 70, 70, 0, 0, 0, 0, 0, 0, 0, & ! 54
182 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
183 8, 15, 21, 35, 35, 43, 56, 70, 0, 0, 0, 0, 0, 0, 0, & ! 86
184 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
185
186 ! r466: <DD|DD> "4" term
187 INTEGER, DIMENSION(0:nelem), PARAMETER, PRIVATE :: ir466 = (/0, & ! 0
188 0, 0, & ! 2
189 0, 0, 0, 0, 0, 0, 0, 0, & ! 10
190 0, 0, 0, 0, 0, 0, 0, 0, & ! 18
191 0, 0, 0, 1, 8, 35, 35, 35, 36, 43, 70, 0, 0, 0, 0, 0, 0, 0, & ! 36
192 0, 0, 0, 1, 21, 35, 35, 36, 43, 70, 70, 0, 0, 0, 0, 0, 0, 0, & ! 54
193 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
194 1, 8, 21, 35, 35, 36, 56, 70, 0, 0, 0, 0, 0, 0, 0, & ! 86
195 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0/)
196
197 INTERFACE amn_l
198 MODULE PROCEDURE amn_l1, amn_l2
199 END INTERFACE
200
201 PUBLIC :: convert_param_to_cp2k, calpar, valence_electrons, get_se_basis, &
202 setup_1c_2el_int, amn_l
203
204CONTAINS
205
206! **************************************************************************************************
207!> \brief Gives back the number of valence electrons for element z and also the
208!> number of atomic orbitals for that specific element
209!> \param sep ...
210!> \param extended_basis_set ...
211! **************************************************************************************************
212 SUBROUTINE valence_electrons(sep, extended_basis_set)
213 TYPE(semi_empirical_type), POINTER :: sep
214 LOGICAL, INTENT(IN) :: extended_basis_set
215
216 INTEGER :: natorb, z
217 LOGICAL :: check, use_p_orbitals
218 REAL(kind=dp) :: zeff
219
220 use_p_orbitals = .true.
221 z = sep%z
222 cpassert(z >= 0)
223 ! Special case for Hydrogen.. If requested allow p-orbitals on it..
224 SELECT CASE (z)
225 CASE (0, 2)
226 use_p_orbitals = .false.
227 CASE (1)
228 use_p_orbitals = sep%p_orbitals_on_h
229 CASE DEFAULT
230 ! Nothing to do..
231 END SELECT
232 ! Determine the number of atomic orbitals
233 natorb = 0
234 IF (nqs(z) > 0) natorb = natorb + 1
235 IF ((nqp(z) > 0) .OR. use_p_orbitals) natorb = natorb + 3
236 IF (extended_basis_set .AND. element_has_d(sep)) natorb = natorb + 5
237 IF (extended_basis_set .AND. element_has_f(sep)) natorb = natorb + 7
238 ! Check and assignment
239 check = (natorb <= 4) .OR. (extended_basis_set)
240 cpassert(check)
241 sep%natorb = natorb
242 sep%extended_basis_set = extended_basis_set
243 ! Determine the Z eff
244 zeff = real(zval(z), kind=dp)
245 sep%zeff = zeff
246 END SUBROUTINE valence_electrons
247
248! **************************************************************************************************
249!> \brief Gives back the number of basis function for each l
250!> \param sep ...
251!> \param l ...
252!> \return ...
253! **************************************************************************************************
254 FUNCTION get_se_basis(sep, l) RESULT(n)
255 TYPE(semi_empirical_type), POINTER :: sep
256 INTEGER, INTENT(IN) :: l
257 INTEGER :: n
258
259 IF (sep%z < 0 .OR. sep%z > nelem) THEN
260 cpabort("Invalid atomic number !")
261 ELSE
262 IF (l == 0) THEN
263 n = nqs(sep%z)
264 ELSEIF (l == 1) THEN
265 ! Special case for Hydrogen.. If requested allow p-orbitals on it..
266 IF ((sep%z == 1) .AND. sep%p_orbitals_on_h) THEN
267 n = 1
268 ELSE
269 n = nqp(sep%z)
270 END IF
271 ELSEIF (l == 2) THEN
272 n = nqd(sep%z)
273 ELSEIF (l == 3) THEN
274 n = nqf(sep%z)
275 ELSE
276 cpabort("Invalid l quantum number !")
277 END IF
278 IF (n < 0) THEN
279 cpabort("Invalid n quantum number !")
280 END IF
281 END IF
282 END FUNCTION get_se_basis
283
284! **************************************************************************************************
285!> \brief Converts parameter units to internal
286!> \param sep ...
287!> \date 03.2008 [tlaino]
288!> \author Teodoro Laino [tlaino] - University of Zurich
289! **************************************************************************************************
290 SUBROUTINE convert_param_to_cp2k(sep)
291 TYPE(semi_empirical_type), POINTER :: sep
292
293 sep%beta = sep%beta/evolt
294 sep%uss = sep%uss/evolt
295 sep%upp = sep%upp/evolt
296 sep%udd = sep%udd/evolt
297 sep%alp = sep%alp/bohr
298 sep%eisol = sep%eisol/evolt
299 sep%gss = sep%gss/evolt
300 sep%gsp = sep%gsp/evolt
301 sep%gpp = sep%gpp/evolt
302 sep%gp2 = sep%gp2/evolt
303 sep%gsd = sep%gsd/evolt
304 sep%gpd = sep%gpd/evolt
305 sep%gdd = sep%gdd/evolt
306 sep%hsp = sep%hsp/evolt
307 sep%fn1 = sep%fn1*bohr/evolt
308 sep%fn2 = sep%fn2/bohr/bohr
309 sep%fn3 = sep%fn3*bohr
310 sep%bfn1 = sep%bfn1*bohr/evolt
311 sep%bfn2 = sep%bfn2/bohr/bohr
312 sep%bfn3 = sep%bfn3*bohr
313 sep%f0sd = sep%f0sd
314 sep%g2sd = sep%g2sd
315 sep%a = sep%a*bohr/evolt
316 sep%b = sep%b/bohr/bohr
317 sep%c = sep%c*bohr
318 sep%pre = sep%pre/evolt
319 sep%d = sep%d/bohr
320
321 END SUBROUTINE convert_param_to_cp2k
322
323! **************************************************************************************************
324!> \brief Calculates missing parameters
325!> \param z ...
326!> \param sep ...
327!> \date 03.2008 [tlaino]
328!> \author Teodoro Laino [tlaino] - University of Zurich
329! **************************************************************************************************
330 SUBROUTINE calpar(z, sep)
331 INTEGER :: z
332 TYPE(semi_empirical_type), POINTER :: sep
333
334 INTEGER :: iod, iop, ios, j, jmax, k, l
335 REAL(kind=dp) :: ad, am, aq, d1, d2, d3, dd, df, eisol, gdd1, gp2, gp2c, gpp, gppc, gqq, &
336 gsp, gspc, gss, gssc, hpp, hpp1, hpp2, hsp, hsp1, hsp2, hspc, p, p4, q1, q2, q3, qf, qn, &
337 qq, udd, upp, uss, zp, zs
338
339 IF (.NOT. sep%defined) RETURN
340 uss = sep%uss
341 upp = sep%upp
342 udd = sep%udd
343 gss = sep%gss
344 gpp = sep%gpp
345 gsp = sep%gsp
346 gp2 = sep%gp2
347 hsp = sep%hsp
348 zs = sep%sto_exponents(0)
349 zp = sep%sto_exponents(1)
350 ios = nos(z)
351 iop = nop(z)
352 iod = nod(z)
353
354 p = 2.0_dp
355 p4 = p**4
356 ! GSSC is the number of two-electron terms of type <SS|SS>
357 gssc = real(max(ios - 1, 0), kind=dp)
358 k = iop
359 ! GSPC is the number of two-electron terms of type <SS|PP>
360 gspc = real(ios*k, kind=dp)
361 l = min(k, 6 - k)
362 ! GP2C is the number of two-electron terms of type <PP|PP>
363 ! plus 0.5 of the number of HPP integrals.
364 ! (HPP is not used; instead it is replaced by 0.5(GPP-GP2))
365 gp2c = real((k*(k - 1))/2, kind=dp) + 0.5_dp*real((l*(l - 1))/2, kind=dp)
366 ! GPPC is minus 0.5 times the number of HPP integrals.
367 gppc = -0.5_dp*real((l*(l - 1))/2, kind=dp)
368 ! HSPC is the number of two-electron terms of type <SP|SP>.
369 ! (S and P must have the same spin. In all cases, if
370 ! P is non-zero, there are two S electrons)
371 hspc = real(-k, kind=dp)
372 ! Constraint the value of the STO exponent
373 zp = max(0.3_dp, zp)
374 ! Take into account constraints on the values of the integrals
375 hpp = 0.5_dp*(gpp - gp2)
376 hpp = max(0.1_dp, hpp)
377 hsp = max(1.e-7_dp, hsp)
378
379 ! Evaluation of EISOL
380 eisol = uss*ios + upp*iop + udd*iod + gss*gssc + gpp*gppc + gsp*gspc + gp2*gp2c + hsp*hspc
381
382 ! Principal quantum number
383 qn = real(nqs(z), kind=dp)
384 cpassert(qn > 0)
385
386 ! Charge separation evaluation
387 dd = (2.0_dp*qn + 1)*(4.0_dp*zs*zp)**(qn + 0.5_dp)/(zs + zp)**(2.0_dp*qn + 2)/sqrt(3.0_dp)
388 qq = sqrt((4.0_dp*qn*qn + 6.0_dp*qn + 2.0_dp)/20.0_dp)/zp
389
390 ! Calculation of the additive terms in atomic units
391 jmax = 5
392 gdd1 = (hsp/(evolt*dd**2))**(1.0_dp/3.0_dp)
393 d1 = gdd1
394 d2 = gdd1 + 0.04_dp
395 DO j = 1, jmax
396 df = d2 - d1
397 hsp1 = 0.5_dp*d1 - 0.5_dp/sqrt(4.0_dp*dd**2 + 1.0_dp/d1**2)
398 hsp2 = 0.5_dp*d2 - 0.5_dp/sqrt(4.0_dp*dd**2 + 1.0_dp/d2**2)
399 IF (abs(hsp2 - hsp1) < epsilon(0.0_dp)) EXIT
400 d3 = d1 + df*(hsp/evolt - hsp1)/(hsp2 - hsp1)
401 d1 = d2
402 d2 = d3
403 END DO
404 gqq = (p4*hpp/(evolt*48.0_dp*qq**4))**0.2_dp
405 q1 = gqq
406 q2 = gqq + 0.04_dp
407 DO j = 1, jmax
408 qf = q2 - q1
409 hpp1 = 0.25_dp*q1 - 0.5_dp/sqrt(4.0_dp*qq**2 + 1.0_dp/q1**2) + 0.25_dp/sqrt(8.0_dp*qq**2 + 1.0_dp/q1**2)
410 hpp2 = 0.25_dp*q2 - 0.5_dp/sqrt(4.0_dp*qq**2 + 1.0_dp/q2**2) + 0.25_dp/sqrt(8.0_dp*qq**2 + 1.0_dp/q2**2)
411 IF (abs(hpp2 - hpp1) < epsilon(0.0_dp)) EXIT
412 q3 = q1 + qf*(hpp/evolt - hpp1)/(hpp2 - hpp1)
413 q1 = q2
414 q2 = q3
415 END DO
416 am = gss/evolt
417 ad = d2
418 aq = q2
419 IF (z == 1) THEN
420 ad = am
421 aq = am
422 dd = 0.0_dp
423 qq = 0.0_dp
424 END IF
425 ! Overwrite these parameters if they were undefined.. otherwise keep the defined
426 ! value
427 IF (abs(sep%eisol) < epsilon(0.0_dp)) sep%eisol = eisol
428 IF (abs(sep%dd) < epsilon(0.0_dp)) sep%dd = dd
429 IF (abs(sep%qq) < epsilon(0.0_dp)) sep%qq = qq
430 IF (abs(sep%am) < epsilon(0.0_dp)) sep%am = am
431 IF (abs(sep%ad) < epsilon(0.0_dp)) sep%ad = ad
432 IF (abs(sep%aq) < epsilon(0.0_dp)) sep%aq = aq
433 ! Proceed with d-orbitals and fill the Kolpman-Ohno and Charge Separation
434 ! arrays
435 CALL calpar_d(sep)
436 END SUBROUTINE calpar
437
438! **************************************************************************************************
439!> \brief Finalize the initialization of parameters, defining additional
440!> parameters for d-orbitals
441!>
442!> \param sep ...
443!> \date 03.2008 [tlaino]
444!> \author Teodoro Laino [tlaino] - University of Zurich
445! **************************************************************************************************
446 SUBROUTINE calpar_d(sep)
447 TYPE(semi_empirical_type), POINTER :: sep
448
449 REAL(kind=dp), DIMENSION(6) :: amn
450
451! Determine if this element owns d-orbitals (only if the parametrization
452! supports the d-orbitals)
453
454 IF (sep%extended_basis_set) sep%dorb = element_has_d(sep)
455 IF (sep%dorb) THEN
456 CALL amn_l(sep, amn)
457 CALL eval_1c_2el_spd(sep)
458 CALL eval_cs_ko(sep, amn)
459 END IF
460 IF (.NOT. sep%dorb) THEN
461 ! Use the old integral module
462 IF (abs(sep%am) > epsilon(0.0_dp)) THEN
463 sep%ko(1) = 0.5_dp/sep%am
464 END IF
465 IF (abs(sep%ad) > epsilon(0.0_dp) .AND. (sep%z /= 1)) THEN
466 sep%ko(2) = 0.5_dp/sep%ad
467 END IF
468 IF (abs(sep%aq) > epsilon(0.0_dp) .AND. (sep%z /= 1)) THEN
469 sep%ko(3) = 0.5_dp/sep%aq
470 END IF
471 sep%ko(7) = sep%ko(1)
472 sep%ko(9) = sep%ko(1)
473 sep%cs(2) = sep%dd
474 sep%cs(3) = sep%qq*sqrt(2.0_dp)
475 ELSE
476 ! Use the new integral module
477 sep%ko(9) = sep%ko(1)
478 sep%aq = 0.5_dp/sep%ko(3)
479 END IF
480 ! In case the Klopman-Ohno CORE therm is provided let's overwrite the
481 ! computed one
482 IF (abs(sep%rho) > epsilon(0.0_dp)) THEN
483 sep%ko(9) = sep%rho
484 END IF
485 END SUBROUTINE calpar_d
486
487! **************************************************************************************************
488!> \brief Determines if the elements has d-orbitals
489!>
490!> \param sep ...
491!> \return ...
492!> \date 05.2008 [tlaino]
493!> \author Teodoro Laino [tlaino] - University of Zurich
494! **************************************************************************************************
495 FUNCTION element_has_d(sep) RESULT(res)
496 TYPE(semi_empirical_type), POINTER :: sep
497 LOGICAL :: res
498
499 res = (nqd(sep%z) > 0 .AND. sep%sto_exponents(2) > epsilon(0.0_dp))
500 END FUNCTION element_has_d
501
502! **************************************************************************************************
503!> \brief Determines if the elements has f-orbitals
504!>
505!> \param sep ...
506!> \return ...
507!> \date 05.2008 [tlaino]
508!> \author Teodoro Laino [tlaino] - University of Zurich
509! **************************************************************************************************
510 FUNCTION element_has_f(sep) RESULT(res)
511 TYPE(semi_empirical_type), POINTER :: sep
512 LOGICAL :: res
513
514 res = (nqf(sep%z) > 0 .AND. sep%sto_exponents(3) > epsilon(0.0_dp))
515 END FUNCTION element_has_f
516
517! **************************************************************************************************
518!> \brief Computes the A^{\mu \nu}_l values for the evaluation of the two-center
519!> two-electron integrals. The term is the one reported in Eq.(7) of TCA
520!>
521!> \param sep ...
522!> \param amn ...
523!> \date 03.2008 [tlaino]
524!> \par Notation Index: 1 (SS), 2 (SP), 3 (SD), 4 (PP), 5 (PD), 6 (DD)
525!>
526!> \author Teodoro Laino [tlaino] - University of Zurich
527! **************************************************************************************************
528 SUBROUTINE amn_l1(sep, amn)
529 TYPE(semi_empirical_type), POINTER :: sep
530 REAL(KIND=dp), DIMENSION(6), INTENT(OUT) :: amn
531
532 INTEGER :: nd, nsp
533 REAL(KIND=dp) :: z1, z2, z3
534
535 z1 = sep%sto_exponents(0)
536 z2 = sep%sto_exponents(1)
537 z3 = sep%sto_exponents(2)
538 IF (z1 <= 0.0_dp) &
539 CALL cp_abort(__location__, &
540 "Trying to use s-orbitals, but the STO exponents is set to 0. "// &
541 "Please check if your parameterization supports the usage of s orbitals! ")
542 amn = 0.0_dp
543 nsp = nqs(sep%z)
544 IF (sep%natorb >= 4) THEN
545 IF (z2 <= 0.0_dp) &
546 CALL cp_abort(__location__, &
547 "Trying to use p-orbitals, but the STO exponents is set to 0. "// &
548 "Please check if your parameterization supports the usage of p orbitals! ")
549 amn(2) = amn_l_low(z1, z2, nsp, nsp, 1)
550 amn(3) = amn_l_low(z2, z2, nsp, nsp, 2)
551 IF (sep%dorb) THEN
552 IF (z3 <= 0.0_dp) &
553 CALL cp_abort(__location__, &
554 "Trying to use d-orbitals, but the STO exponents is set to 0. "// &
555 "Please check if your parameterization supports the usage of d orbitals! ")
556 nd = nqd(sep%z)
557 amn(4) = amn_l_low(z1, z3, nsp, nd, 2)
558 amn(5) = amn_l_low(z2, z3, nsp, nd, 1)
559 amn(6) = amn_l_low(z3, z3, nd, nd, 2)
560 END IF
561 END IF
562 END SUBROUTINE amn_l1
563
564! **************************************************************************************************
565!> \brief Computes the A^{\mu \nu}_l values for the evaluation of the two-center
566!> two-electron integrals. The term is the one reported in Eq.(7) of TCA
567!>
568!> \param sep ...
569!> \param amn ...
570!> \date 09.2008 [tlaino]
571!> \par
572!>
573!> \author Teodoro Laino [tlaino] - University of Zurich
574! **************************************************************************************************
575 SUBROUTINE amn_l2(sep, amn)
576 TYPE(semi_empirical_type), POINTER :: sep
577 REAL(kind=dp), DIMENSION(6, 0:2), INTENT(OUT) :: amn
578
579 INTEGER :: nd, nsp
580 REAL(kind=dp) :: z1, z2, z3
581
582 z1 = sep%sto_exponents(0)
583 z2 = sep%sto_exponents(1)
584 z3 = sep%sto_exponents(2)
585 cpassert(z1 > 0.0_dp)
586 amn = 0.0_dp
587 nsp = nqs(sep%z)
588 amn(1, 0) = amn_l_low(z1, z1, nsp, nsp, 0)
589 IF (sep%natorb >= 4) THEN
590 cpassert(z2 > 0.0_dp)
591 amn(2, 1) = amn_l_low(z1, z2, nsp, nsp, 1)
592 amn(3, 0) = amn_l_low(z2, z2, nsp, nsp, 0)
593 amn(3, 2) = amn_l_low(z2, z2, nsp, nsp, 2)
594 IF (sep%dorb) THEN
595 cpassert(z3 > 0.0_dp)
596 nd = nqd(sep%z)
597 amn(4, 2) = amn_l_low(z1, z3, nsp, nd, 2)
598 amn(5, 1) = amn_l_low(z2, z3, nsp, nd, 1)
599 amn(6, 0) = amn_l_low(z3, z3, nd, nd, 0)
600 amn(6, 2) = amn_l_low(z3, z3, nd, nd, 2)
601 END IF
602 END IF
603 END SUBROUTINE amn_l2
604
605! **************************************************************************************************
606!> \brief Low level for computing Eq.(7) of TCA
607!> \param z1 ...
608!> \param z2 ...
609!> \param n1 ...
610!> \param n2 ...
611!> \param l ...
612!> \return ...
613!> \date 03.2008 [tlaino]
614!> \author Teodoro Laino [tlaino] - University of Zurich
615! **************************************************************************************************
616 FUNCTION amn_l_low(z1, z2, n1, n2, l) RESULT(amnl)
617 REAL(kind=dp), INTENT(IN) :: z1, z2
618 INTEGER, INTENT(IN) :: n1, n2, l
619 REAL(kind=dp) :: amnl
620
621 amnl = fac(n1 + n2 + l)/sqrt(fac(2*n1)*fac(2*n2))*(2.0_dp*z1/(z1 + z2))**n1* &
622 (2.0_dp*z2/(z1 + z2))**n2*2.0_dp*sqrt(z1*z2)/(z1 + z2)**(l + 1)
623
624 END FUNCTION amn_l_low
625
626! **************************************************************************************************
627!> \brief Calculation of chare separations and additive terms used for computing
628!> the two-center two-electron integrals with d-orbitals
629!> \param sep ...
630!> \param amn ...
631!> \date 03.2008 [tlaino]
632!> \par Notation
633!> -) Charge separations [sep%cs(1:6)] [see equations (12)-(16) of TCA]
634!> -) Additive terms of Klopman-Ohno terms [sep%ko(1:9)] [see equations
635!> (19)-(26) of TCA]
636!> -) Atomic core additive term stored in sep%ko(9): used in the calculation
637!> of the core-electron attractions and core-core repulsions
638!> \author Teodoro Laino [tlaino] - University of Zurich
639! **************************************************************************************************
640 SUBROUTINE eval_cs_ko(sep, amn)
641 TYPE(semi_empirical_type), POINTER :: sep
642 REAL(kind=dp), DIMENSION(6), INTENT(IN) :: amn
643
644 REAL(kind=dp) :: d, fg
645
646! SS term
647
648 fg = sep%gss
649 sep%ko(1) = ko_ij(0, 1.0_dp, fg)
650 IF (sep%natorb >= 4) THEN
651 ! Other terms for SP basis
652 ! SP
653 d = amn(2)/sqrt(3.0_dp)
654 fg = sep%hsp
655 sep%cs(2) = d
656 sep%ko(2) = ko_ij(1, d, fg)
657 ! PP
658 sep%ko(7) = sep%ko(1)
659 d = sqrt(amn(3)*2.0_dp/5.0_dp)
660 fg = 0.5_dp*(sep%gpp - sep%gp2)
661 sep%cs(3) = d
662 sep%ko(3) = ko_ij(2, d, fg)
663 ! Terms involving d-orbitals
664 IF (sep%dorb) THEN
665 ! SD
666 d = sqrt(amn(4)*2.0_dp/sqrt(15.0_dp))
667 fg = sep%onec2el(19)
668 sep%cs(4) = d
669 sep%ko(4) = ko_ij(2, d, fg)
670 ! PD
671 d = amn(5)/sqrt(5.0_dp)
672 fg = sep%onec2el(23) - 1.8_dp*sep%onec2el(35)
673 sep%cs(5) = d
674 sep%ko(5) = ko_ij(1, d, fg)
675 ! DD
676 fg = 0.2_dp*(sep%onec2el(29) + 2.0_dp*sep%onec2el(30) + 2.0_dp*sep%onec2el(31))
677 sep%ko(8) = ko_ij(0, 1.0_dp, fg)
678 d = sqrt(amn(6)*2.0_dp/7.0_dp)
679 fg = sep%onec2el(44) - (20.0_dp/35.0_dp)*sep%onec2el(52)
680 sep%cs(6) = d
681 sep%ko(6) = ko_ij(2, d, fg)
682 END IF
683 END IF
684 END SUBROUTINE eval_cs_ko
685
686! **************************************************************************************************
687!> \brief Computes the 1 center two-electrons integrals for a SPD basis
688!> \param sep ...
689!> \date 03.2008 [tlaino]
690!> \author Teodoro Laino [tlaino] - University of Zurich
691! **************************************************************************************************
692 SUBROUTINE eval_1c_2el_spd(sep)
693 TYPE(semi_empirical_type), POINTER :: sep
694
695 REAL(kind=dp) :: r016, r036, r066, r125, r155, r234, &
696 r236, r244, r246, r266, r355, r466, &
697 s15, s3, s5
698
699 IF (sep%dorb) THEN
700 s3 = sqrt(3.0_dp)
701 s5 = sqrt(5.0_dp)
702 s15 = sqrt(15.0_dp)
703
704 ! We evaluate now the Slater-Condon parameters (Rlij)
705 CALL sc_param(sep, r066, r266, r466, r016, r244, r036, r236, r155, r355, r125, &
706 r234, r246)
707
708 IF (abs(sep%f0sd) > epsilon(0.0_dp)) THEN
709 r016 = sep%f0sd
710 END IF
711 IF (abs(sep%g2sd) > epsilon(0.0_dp)) THEN
712 r244 = sep%g2sd
713 END IF
714 CALL eisol_corr(sep, r016, r066, r244, r266, r466)
715 sep%onec2el(1) = r016
716 sep%onec2el(2) = 2.0_dp/(3.0_dp*s5)*r125
717 sep%onec2el(3) = 1.0_dp/s15*r125
718 sep%onec2el(4) = 2.0_dp/(5.0_dp*s5)*r234
719 sep%onec2el(5) = r036 + 4.0_dp/35.0_dp*r236
720 sep%onec2el(6) = r036 + 2.0_dp/35.0_dp*r236
721 sep%onec2el(7) = r036 - 4.0_dp/35.0_dp*r236
722 sep%onec2el(8) = -1.0_dp/(3.0_dp*s5)*r125
723 sep%onec2el(9) = sqrt(3.0_dp/125.0_dp)*r234
724 sep%onec2el(10) = s3/35.0_dp*r236
725 sep%onec2el(11) = 3.0_dp/35.0_dp*r236
726 sep%onec2el(12) = -1.0_dp/(5.0_dp*s5)*r234
727 sep%onec2el(13) = r036 - 2.0_dp/35.0_dp*r236
728 sep%onec2el(14) = -2.0_dp*s3/35.0_dp*r236
729 sep%onec2el(15) = -sep%onec2el(3)
730 sep%onec2el(16) = -sep%onec2el(11)
731 sep%onec2el(17) = -sep%onec2el(9)
732 sep%onec2el(18) = -sep%onec2el(14)
733 sep%onec2el(19) = 1.0_dp/5.0_dp*r244
734 sep%onec2el(20) = 2.0_dp/(7.0_dp*s5)*r246
735 sep%onec2el(21) = sep%onec2el(20)/2.0_dp
736 sep%onec2el(22) = -sep%onec2el(20)
737 sep%onec2el(23) = 4.0_dp/15.0_dp*r155 + 27.0_dp/245.0_dp*r355
738 sep%onec2el(24) = 2.0_dp*s3/15.0_dp*r155 - 9.0_dp*s3/245.0_dp*r355
739 sep%onec2el(25) = 1.0_dp/15.0_dp*r155 + 18.0_dp/245.0_dp*r355
740 sep%onec2el(26) = -s3/15.0_dp*r155 + 12.0_dp*s3/245.0_dp*r355
741 sep%onec2el(27) = -s3/15.0_dp*r155 - 3.0_dp*s3/245.0_dp*r355
742 sep%onec2el(28) = -sep%onec2el(27)
743 sep%onec2el(29) = r066 + 4.0_dp/49.0_dp*r266 + 4.0_dp/49.0_dp*r466
744 sep%onec2el(30) = r066 + 2.0_dp/49.0_dp*r266 - 24.0_dp/441.0_dp*r466
745 sep%onec2el(31) = r066 - 4.0_dp/49.0_dp*r266 + 6.0_dp/441.0_dp*r466
746 sep%onec2el(32) = sqrt(3.0_dp/245.0_dp)*r246
747 sep%onec2el(33) = 1.0_dp/5.0_dp*r155 + 24.0_dp/245.0_dp*r355
748 sep%onec2el(34) = 1.0_dp/5.0_dp*r155 - 6.0_dp/245.0_dp*r355
749 sep%onec2el(35) = 3.0_dp/49.0_dp*r355
750 sep%onec2el(36) = 1.0_dp/49.0_dp*r266 + 30.0_dp/441.0_dp*r466
751 sep%onec2el(37) = s3/49.0_dp*r266 - 5.0_dp*s3/441.0_dp*r466
752 sep%onec2el(38) = r066 - 2.0_dp/49.0_dp*r266 - 4.0_dp/441.0_dp*r466
753 sep%onec2el(39) = -2.0_dp*s3/49.0_dp*r266 + 10.0_dp*s3/441.0_dp*r466
754 sep%onec2el(40) = -sep%onec2el(32)
755 sep%onec2el(41) = -sep%onec2el(34)
756 sep%onec2el(42) = -sep%onec2el(35)
757 sep%onec2el(43) = -sep%onec2el(37)
758 sep%onec2el(44) = 3.0_dp/49.0_dp*r266 + 20.0_dp/441.0_dp*r466
759 sep%onec2el(45) = -sep%onec2el(39)
760 sep%onec2el(46) = 1.0_dp/5.0_dp*r155 - 3.0_dp/35.0_dp*r355
761 sep%onec2el(47) = -sep%onec2el(46)
762 sep%onec2el(48) = 4.0_dp/49.0_dp*r266 + 15.0_dp/441.0_dp*r466
763 sep%onec2el(49) = 3.0_dp/49.0_dp*r266 - 5.0_dp/147.0_dp*r466
764 sep%onec2el(50) = -sep%onec2el(49)
765 sep%onec2el(51) = r066 + 4.0_dp/49.0_dp*r266 - 34.0_dp/441.0_dp*r466
766 sep%onec2el(52) = 35.0_dp/441.0_dp*r466
767 sep%f0dd = r066
768 sep%f2dd = r266
769 sep%f4dd = r466
770 sep%f0sd = r016
771 sep%g2sd = r244
772 sep%f0pd = r036
773 sep%f2pd = r236
774 sep%g1pd = r155
775 sep%g3pd = r355
776 END IF
777 END SUBROUTINE eval_1c_2el_spd
778
779! **************************************************************************************************
780!> \brief Slater-Condon parameters for 1 center 2 electrons integrals
781!> \param sep ...
782!> \param r066 ...
783!> \param r266 ...
784!> \param r466 ...
785!> \param r016 ...
786!> \param r244 ...
787!> \param r036 ...
788!> \param r236 ...
789!> \param r155 ...
790!> \param r355 ...
791!> \param r125 ...
792!> \param r234 ...
793!> \param r246 ...
794!> \date 03.2008 [tlaino]
795!> \author Teodoro Laino [tlaino] - University of Zurich
796! **************************************************************************************************
797 SUBROUTINE sc_param(sep, r066, r266, r466, r016, r244, r036, r236, r155, r355, &
798 r125, r234, r246)
799 TYPE(semi_empirical_type), POINTER :: sep
800 REAL(kind=dp), INTENT(out) :: r066, r266, r466, r016, r244, r036, &
801 r236, r155, r355, r125, r234, r246
802
803 INTEGER :: nd, ns
804 REAL(kind=dp) :: ed, ep, es
805
806 ns = nqs(sep%z)
807 nd = nqd(sep%z)
808 es = sep%zn(0)
809 ep = sep%zn(1)
810 ed = sep%zn(2)
811 r016 = sc_param_low(0, ns, es, ns, es, nd, ed, nd, ed)
812 r036 = sc_param_low(0, ns, ep, ns, ep, nd, ed, nd, ed)
813 r066 = sc_param_low(0, nd, ed, nd, ed, nd, ed, nd, ed)
814 r155 = sc_param_low(1, ns, ep, nd, ed, ns, ep, nd, ed)
815 r125 = sc_param_low(1, ns, es, ns, ep, ns, ep, nd, ed)
816 r244 = sc_param_low(2, ns, es, nd, ed, ns, es, nd, ed)
817 r236 = sc_param_low(2, ns, ep, ns, ep, nd, ed, nd, ed)
818 r266 = sc_param_low(2, nd, ed, nd, ed, nd, ed, nd, ed)
819 r234 = sc_param_low(2, ns, ep, ns, ep, ns, es, nd, ed)
820 r246 = sc_param_low(2, ns, es, nd, ed, nd, ed, nd, ed)
821 r355 = sc_param_low(3, ns, ep, nd, ed, ns, ep, nd, ed)
822 r466 = sc_param_low(4, nd, ed, nd, ed, nd, ed, nd, ed)
823 END SUBROUTINE sc_param
824
825! **************************************************************************************************
826!> \brief Slater-Condon parameters for 1 center 2 electrons integrals - Low level
827!> \param k ...
828!> \param na ...
829!> \param ea ...
830!> \param nb ...
831!> \param eb ...
832!> \param nc ...
833!> \param ec ...
834!> \param nd ...
835!> \param ed ...
836!> \return ...
837!> \date 03.2008 [tlaino]
838!> \par Notation
839!> -) k: Type of integral
840!> -) na,na: Principle Quantum Number of AO,corresponding to electron 1
841!> -) ea,eb: Exponents of AO,corresponding to electron 1
842!> -) nb,nc: Principle Quantum Number of AO,corresponding to electron 2
843!> -) ec,ed: Exponents of AO,corresponding to electron 2
844!> \author Teodoro Laino [tlaino] - University of Zurich
845! **************************************************************************************************
846 FUNCTION sc_param_low(k, na, ea, nb, eb, nc, ec, nd, ed) RESULT(res)
847 INTEGER, INTENT(in) :: k, na
848 REAL(kind=dp), INTENT(in) :: ea
849 INTEGER, INTENT(in) :: nb
850 REAL(kind=dp), INTENT(in) :: eb
851 INTEGER, INTENT(in) :: nc
852 REAL(kind=dp), INTENT(in) :: ec
853 INTEGER, INTENT(in) :: nd
854 REAL(kind=dp), INTENT(in) :: ed
855 REAL(kind=dp) :: res
856
857 INTEGER :: i, m, m1, m2, n, nab, ncd
858 REAL(kind=dp) :: a2, aab, acd, ae, aea, aeb, aec, aed, c, &
859 e, eab, ecd, ff, s0, s1, s2, s3, tmp
860
861 cpassert(ea > 0.0_dp)
862 cpassert(eb > 0.0_dp)
863 cpassert(ec > 0.0_dp)
864 cpassert(ed > 0.0_dp)
865 aea = log(ea)
866 aeb = log(eb)
867 aec = log(ec)
868 aed = log(ed)
869 nab = na + nb
870 ncd = nc + nd
871 ecd = ec + ed
872 eab = ea + eb
873 e = ecd + eab
874 n = nab + ncd
875 ae = log(e)
876 a2 = log(2.0_dp)
877 acd = log(ecd)
878 aab = log(eab)
879 ff = fac(n - 1)/sqrt(fac(2*na)*fac(2*nb)*fac(2*nc)*fac(2*nd))
880 tmp = na*aea + nb*aeb + nc*aec + nd*aed + 0.5_dp*(aea + aeb + aec + aed) + a2*(n + 2) - ae*n
881 c = evolt*ff*exp(tmp)
882 s0 = 1.0_dp/e
883 s1 = 0.0_dp
884 s2 = 0.0_dp
885 m = ncd - k
886 DO i = 1, m
887 s0 = s0*e/ecd
888 s1 = s1 + s0*(binomial(ncd - k - 1, i - 1) - binomial(ncd + k, i - 1))/binomial(n - 1, i - 1)
889 END DO
890 m1 = m
891 m2 = ncd + k
892 DO i = m1, m2
893 s0 = s0*e/ecd
894 s2 = s2 + s0*binomial(m2, i)/binomial(n - 1, i)
895 END DO
896 s3 = exp(ae*n - acd*(m2 + 1) - aab*(nab - k))/binomial(n - 1, m2)
897 res = c*(s1 - s2 + s3)
898 END FUNCTION sc_param_low
899
900! **************************************************************************************************
901!> \brief Corrects the EISOL fo the one-center terms coming from those atoms
902!> that have partially filled "d" shells
903!> \param sep ...
904!> \param r016 ...
905!> \param r066 ...
906!> \param r244 ...
907!> \param r266 ...
908!> \param r466 ...
909!> \date 03.2008 [tlaino]
910!> \par Notation
911!> r016: <SS|DD>
912!> r066: <DD|DD> "0" term
913!> r244: <SD|SD>
914!> r266: <DD|DD> "2" term
915!> r466: <DD|DD> "4" term
916!>
917!> \author Teodoro Laino [tlaino] - University of Zurich
918! **************************************************************************************************
919 SUBROUTINE eisol_corr(sep, r016, r066, r244, r266, r466)
920 TYPE(semi_empirical_type), POINTER :: sep
921 REAL(kind=dp), INTENT(in) :: r016, r066, r244, r266, r466
922
923 sep%eisol = sep%eisol + ir016(sep%z)*r016 + &
924 ir066(sep%z)*r066 - &
925 ir244(sep%z)*r244/5.0_dp - &
926 ir266(sep%z)*r266/49.0_dp - &
927 ir466(sep%z)*r466/49.0_dp
928 END SUBROUTINE eisol_corr
929
930! **************************************************************************************************
931!> \brief Computes the Klopman-Ohno additive terms for 2-center 2-electron
932!> integrals requiring that the corresponding 1-center 2-electron integral
933!> is reproduced from the 2-center one for r->0
934!> \param l ...
935!> \param d ...
936!> \param fg ...
937!> \return ...
938!> \date 03.2008 [tlaino]
939!> \author Teodoro Laino [tlaino] - University of Zurich
940! **************************************************************************************************
941 FUNCTION ko_ij(l, d, fg) RESULT(res)
942 INTEGER, INTENT(in) :: l
943 REAL(kind=dp), INTENT(in) :: d, fg
944 REAL(kind=dp) :: res
945
946 INTEGER, PARAMETER :: niter = 100
947 REAL(kind=dp), PARAMETER :: epsil = 1.0e-08_dp, g1 = 0.382_dp, &
948 g2 = 0.618_dp
949
950 INTEGER :: i
951 REAL(kind=dp) :: a1, a2, delta, dsq, ev4, ev8, f1, f2, &
952 y1, y2
953
954 cpassert(fg /= 0.0_dp)
955 ! Term for SS
956 IF (l == 0) THEN
957 res = 0.5_dp*evolt/fg
958 RETURN
959 END IF
960 ! Term for Higher angular momentum
961 dsq = d*d
962 ev4 = evolt*0.25_dp
963 ev8 = evolt/8.0_dp
964 a1 = 0.1_dp
965 a2 = 5.0_dp
966 DO i = 1, niter
967 delta = a2 - a1
968 IF (delta < epsil) EXIT
969 y1 = a1 + delta*g1
970 y2 = a1 + delta*g2
971 IF (l == 1) THEN
972 f1 = (ev4*(1/y1 - 1/sqrt(y1**2 + dsq)) - fg)**2
973 f2 = (ev4*(1/y2 - 1/sqrt(y2**2 + dsq)) - fg)**2
974 ELSE IF (l == 2) THEN
975 f1 = (ev8*(1.0_dp/y1 - 2.0_dp/sqrt(y1**2 + dsq*0.5_dp) + 1.0_dp/sqrt(y1**2 + dsq)) - fg)**2
976 f2 = (ev8*(1/y2 - 2.0_dp/sqrt(y2**2 + dsq*0.5_dp) + 1.0_dp/sqrt(y2**2 + dsq)) - fg)**2
977 END IF
978 IF (f1 < f2) THEN
979 a2 = y2
980 ELSE
981 a1 = y1
982 END IF
983 END DO
984 ! Convergence reached.. define additive terms
985 IF (f1 >= f2) THEN
986 res = a2
987 ELSE
988 res = a1
989 END IF
990 END FUNCTION ko_ij
991
992! **************************************************************************************************
993!> \brief Fills the 1 center 2 electron integrals for the construction of the
994!> one-electron fock matrix
995!> \param sep ...
996!> \date 04.2008 [tlaino]
997!> \author Teodoro Laino [tlaino] - University of Zurich
998! **************************************************************************************************
999 SUBROUTINE setup_1c_2el_int(sep)
1000 TYPE(semi_empirical_type), POINTER :: sep
1001
1002 INTEGER :: i, ij, ij0, ind, ip, ipx, ipy, ipz, &
1003 isize, kl, natorb
1004 LOGICAL :: defined
1005 REAL(kind=dp) :: gp2, gpp, gsp, gss, hsp
1006
1007 CALL get_se_param(sep, defined=defined, natorb=natorb, &
1008 gss=gss, gsp=gsp, gpp=gpp, gp2=gp2, hsp=hsp)
1009 cpassert(defined)
1010
1011 isize = natorb*(natorb + 1)/2
1012 ALLOCATE (sep%w(isize, isize))
1013 ! Initialize array
1014 sep%w = 0.0_dp
1015 ! Fill the array
1016 IF (natorb > 0) THEN
1017 ip = 1
1018 sep%w(ip, ip) = gss
1019 IF (natorb > 2) THEN
1020 ipx = ip + 2
1021 ipy = ip + 5
1022 ipz = ip + 9
1023 sep%w(ipx, ip) = gsp
1024 sep%w(ipy, ip) = gsp
1025 sep%w(ipz, ip) = gsp
1026 sep%w(ip, ipx) = gsp
1027 sep%w(ip, ipy) = gsp
1028 sep%w(ip, ipz) = gsp
1029 sep%w(ipx, ipx) = gpp
1030 sep%w(ipy, ipy) = gpp
1031 sep%w(ipz, ipz) = gpp
1032 sep%w(ipy, ipx) = gp2
1033 sep%w(ipz, ipx) = gp2
1034 sep%w(ipz, ipy) = gp2
1035 sep%w(ipx, ipy) = gp2
1036 sep%w(ipx, ipz) = gp2
1037 sep%w(ipy, ipz) = gp2
1038 sep%w(ip + 1, ip + 1) = hsp
1039 sep%w(ip + 3, ip + 3) = hsp
1040 sep%w(ip + 6, ip + 6) = hsp
1041 sep%w(ip + 4, ip + 4) = 0.5_dp*(gpp - gp2)
1042 sep%w(ip + 7, ip + 7) = 0.5_dp*(gpp - gp2)
1043 sep%w(ip + 8, ip + 8) = 0.5_dp*(gpp - gp2)
1044 IF (sep%dorb) THEN
1045 ij0 = ip - 1
1046 DO i = 1, 243
1047 ij = int_ij(i)
1048 kl = int_kl(i)
1049 ind = int_onec2el(i)
1050 sep%w(ij + ij0, kl + ij0) = sep%onec2el(ind)/evolt
1051 END DO
1052 END IF
1053 END IF
1054 END IF
1055 END SUBROUTINE setup_1c_2el_int
1056
1057END MODULE semi_empirical_par_utils
1058
semi_empirical_par_utils::amn_l
Definition semi_empirical_par_utils.F:197
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::binomial
elemental real(kind=dp) function, public binomial(n, k)
The binomial coefficient n over k for 0 <= k <= n is calculated, otherwise zero is returned.
Definition mathlib.F:205
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
semi_empirical_int_arrays
Arrays of parameters used in the semi-empirical calculations \References Everywhere in this module TC...
Definition semi_empirical_int_arrays.F:17
semi_empirical_int_arrays::int_kl
integer, dimension(243), public int_kl
Definition semi_empirical_int_arrays.F:107
semi_empirical_int_arrays::int_ij
integer, dimension(243), public int_ij
Definition semi_empirical_int_arrays.F:95
semi_empirical_int_arrays::int_onec2el
integer, dimension(243), public int_onec2el
Definition semi_empirical_int_arrays.F:119
semi_empirical_par_utils
Utilities to post-process semi-empirical parameters.
Definition semi_empirical_par_utils.F:16
semi_empirical_par_utils::setup_1c_2el_int
subroutine, public setup_1c_2el_int(sep)
Fills the 1 center 2 electron integrals for the construction of the one-electron fock matrix.
Definition semi_empirical_par_utils.F:1000
semi_empirical_par_utils::valence_electrons
subroutine, public valence_electrons(sep, extended_basis_set)
Gives back the number of valence electrons for element z and also the number of atomic orbitals for t...
Definition semi_empirical_par_utils.F:213
semi_empirical_par_utils::calpar
subroutine, public calpar(z, sep)
Calculates missing parameters.
Definition semi_empirical_par_utils.F:331
semi_empirical_par_utils::get_se_basis
integer function, public get_se_basis(sep, l)
Gives back the number of basis function for each l.
Definition semi_empirical_par_utils.F:255
semi_empirical_par_utils::convert_param_to_cp2k
subroutine, public convert_param_to_cp2k(sep)
Converts parameter units to internal.
Definition semi_empirical_par_utils.F:291
semi_empirical_types
Definition of the semi empirical parameter types.
Definition semi_empirical_types.F:12
semi_empirical_types::get_se_param
subroutine, public get_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
Get info from the semi-empirical type.
Definition semi_empirical_types.F:363
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57