(git:1f285aa)
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 ! **************************************************************************************************
16 MODULE 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 
204 CONTAINS
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 
1057 END MODULE semi_empirical_par_utils
1058 
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