Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 3964 : 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 3964 : use_p_orbitals = .TRUE.
221 3964 : z = sep%z
222 3964 : CPASSERT(z >= 0)
223 : ! Special case for Hydrogen.. If requested allow p-orbitals on it..
224 1768 : SELECT CASE (z)
225 : CASE (0, 2)
226 1768 : use_p_orbitals = .FALSE.
227 : CASE (1)
228 3964 : 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 3964 : natorb = 0
234 3964 : IF (nqs(z) > 0) natorb = natorb + 1
235 3964 : IF ((nqp(z) > 0) .OR. use_p_orbitals) natorb = natorb + 3
236 3964 : IF (extended_basis_set .AND. element_has_d(sep)) natorb = natorb + 5
237 3964 : IF (extended_basis_set .AND. element_has_f(sep)) natorb = natorb + 7
238 : ! Check and assignment
239 3964 : check = (natorb <= 4) .OR. (extended_basis_set)
240 0 : CPASSERT(check)
241 3964 : sep%natorb = natorb
242 3964 : sep%extended_basis_set = extended_basis_set
243 : ! Determine the Z eff
244 3964 : zeff = REAL(zval(z), KIND=dp)
245 3964 : sep%zeff = zeff
246 3964 : 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 4750 : 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 4750 : IF (sep%z < 0 .OR. sep%z > nelem) THEN
260 0 : CPABORT("Invalid atomic number !")
261 : ELSE
262 4750 : IF (l == 0) THEN
263 2240 : n = nqs(sep%z)
264 2510 : ELSEIF (l == 1) THEN
265 : ! Special case for Hydrogen.. If requested allow p-orbitals on it..
266 1772 : IF ((sep%z == 1) .AND. sep%p_orbitals_on_h) THEN
267 : n = 1
268 : ELSE
269 1758 : n = nqp(sep%z)
270 : END IF
271 738 : ELSEIF (l == 2) THEN
272 738 : n = nqd(sep%z)
273 0 : ELSEIF (l == 3) THEN
274 0 : n = nqf(sep%z)
275 : ELSE
276 0 : CPABORT("Invalid l quantum number !")
277 : END IF
278 4736 : IF (n < 0) THEN
279 0 : CPABORT("Invalid n quantum number !")
280 : END IF
281 : END IF
282 4750 : 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 2240 : SUBROUTINE convert_param_to_cp2k(sep)
291 : TYPE(semi_empirical_type), POINTER :: sep
292 :
293 11200 : sep%beta = sep%beta/evolt
294 2240 : sep%uss = sep%uss/evolt
295 2240 : sep%upp = sep%upp/evolt
296 2240 : sep%udd = sep%udd/evolt
297 2240 : sep%alp = sep%alp/bohr
298 2240 : sep%eisol = sep%eisol/evolt
299 2240 : sep%gss = sep%gss/evolt
300 2240 : sep%gsp = sep%gsp/evolt
301 2240 : sep%gpp = sep%gpp/evolt
302 2240 : sep%gp2 = sep%gp2/evolt
303 2240 : sep%gsd = sep%gsd/evolt
304 2240 : sep%gpd = sep%gpd/evolt
305 2240 : sep%gdd = sep%gdd/evolt
306 2240 : sep%hsp = sep%hsp/evolt
307 11200 : sep%fn1 = sep%fn1*bohr/evolt
308 11200 : sep%fn2 = sep%fn2/bohr/bohr
309 11200 : sep%fn3 = sep%fn3*bohr
310 47040 : sep%bfn1 = sep%bfn1*bohr/evolt
311 47040 : sep%bfn2 = sep%bfn2/bohr/bohr
312 47040 : sep%bfn3 = sep%bfn3*bohr
313 2240 : sep%f0sd = sep%f0sd
314 : sep%g2sd = sep%g2sd
315 2240 : sep%a = sep%a*bohr/evolt
316 2240 : sep%b = sep%b/bohr/bohr
317 2240 : sep%c = sep%c*bohr
318 6720 : sep%pre = sep%pre/evolt
319 6720 : sep%d = sep%d/bohr
320 :
321 2240 : 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 2240 : 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 2240 : IF (.NOT. sep%defined) RETURN
340 2240 : uss = sep%uss
341 2240 : upp = sep%upp
342 2240 : udd = sep%udd
343 2240 : gss = sep%gss
344 2240 : gpp = sep%gpp
345 2240 : gsp = sep%gsp
346 2240 : gp2 = sep%gp2
347 2240 : hsp = sep%hsp
348 2240 : zs = sep%sto_exponents(0)
349 2240 : zp = sep%sto_exponents(1)
350 2240 : ios = Nos(z)
351 2240 : iop = Nop(z)
352 2240 : iod = Nod(z)
353 :
354 2240 : p = 2.0_dp
355 2240 : p4 = p**4
356 : ! GSSC is the number of two-electron terms of type <SS|SS>
357 2240 : gssc = REAL(MAX(ios - 1, 0), KIND=dp)
358 2240 : k = iop
359 : ! GSPC is the number of two-electron terms of type <SS|PP>
360 2240 : gspc = REAL(ios*k, KIND=dp)
361 2240 : 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 2240 : 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 2240 : 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 2240 : hspc = REAL(-k, KIND=dp)
372 : ! Constraint the value of the STO exponent
373 2240 : zp = MAX(0.3_dp, zp)
374 : ! Take into account constraints on the values of the integrals
375 2240 : hpp = 0.5_dp*(gpp - gp2)
376 2240 : hpp = MAX(0.1_dp, hpp)
377 2240 : hsp = MAX(1.E-7_dp, hsp)
378 :
379 : ! Evaluation of EISOL
380 2240 : eisol = uss*ios + upp*iop + udd*iod + gss*gssc + gpp*gppc + gsp*gspc + gp2*gp2c + hsp*hspc
381 :
382 : ! Principal quantum number
383 2240 : qn = REAL(nqs(z), KIND=dp)
384 2240 : CPASSERT(qn > 0)
385 :
386 : ! Charge separation evaluation
387 2240 : 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 2240 : 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 2240 : jmax = 5
392 2240 : gdd1 = (hsp/(evolt*dd**2))**(1.0_dp/3.0_dp)
393 2240 : d1 = gdd1
394 2240 : d2 = gdd1 + 0.04_dp
395 12126 : DO j = 1, jmax
396 10324 : df = d2 - d1
397 10324 : hsp1 = 0.5_dp*d1 - 0.5_dp/SQRT(4.0_dp*dd**2 + 1.0_dp/d1**2)
398 10324 : hsp2 = 0.5_dp*d2 - 0.5_dp/SQRT(4.0_dp*dd**2 + 1.0_dp/d2**2)
399 10324 : IF (ABS(hsp2 - hsp1) < EPSILON(0.0_dp)) EXIT
400 9886 : d3 = d1 + df*(hsp/evolt - hsp1)/(hsp2 - hsp1)
401 9886 : d1 = d2
402 12126 : d2 = d3
403 : END DO
404 2240 : gqq = (p4*hpp/(evolt*48.0_dp*qq**4))**0.2_dp
405 2240 : q1 = gqq
406 2240 : q2 = gqq + 0.04_dp
407 13440 : DO j = 1, jmax
408 11200 : qf = q2 - q1
409 11200 : 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 11200 : 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 11200 : IF (ABS(hpp2 - hpp1) < EPSILON(0.0_dp)) EXIT
412 11200 : q3 = q1 + qf*(hpp/evolt - hpp1)/(hpp2 - hpp1)
413 11200 : q1 = q2
414 13440 : q2 = q3
415 : END DO
416 2240 : am = gss/evolt
417 2240 : ad = d2
418 2240 : aq = q2
419 2240 : IF (z == 1) THEN
420 438 : ad = am
421 438 : aq = am
422 438 : dd = 0.0_dp
423 438 : qq = 0.0_dp
424 : END IF
425 : ! Overwrite these parameters if they were undefined.. otherwise keep the defined
426 : ! value
427 2240 : IF (ABS(sep%eisol) < EPSILON(0.0_dp)) sep%eisol = eisol
428 2240 : IF (ABS(sep%dd) < EPSILON(0.0_dp)) sep%dd = dd
429 2240 : IF (ABS(sep%qq) < EPSILON(0.0_dp)) sep%qq = qq
430 2240 : IF (ABS(sep%am) < EPSILON(0.0_dp)) sep%am = am
431 2240 : IF (ABS(sep%ad) < EPSILON(0.0_dp)) sep%ad = ad
432 2240 : 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 2240 : 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 2240 : 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 2240 : IF (sep%extended_basis_set) sep%dorb = element_has_d(sep)
455 2240 : IF (sep%dorb) THEN
456 738 : CALL amn_l(sep, amn)
457 738 : CALL eval_1c_2el_spd(sep)
458 738 : CALL eval_cs_ko(sep, amn)
459 : END IF
460 2240 : IF (.NOT. sep%dorb) THEN
461 : ! Use the old integral module
462 1502 : IF (ABS(sep%am) > EPSILON(0.0_dp)) THEN
463 1502 : sep%ko(1) = 0.5_dp/sep%am
464 : END IF
465 1502 : IF (ABS(sep%ad) > EPSILON(0.0_dp) .AND. (sep%z /= 1)) THEN
466 1064 : sep%ko(2) = 0.5_dp/sep%ad
467 : END IF
468 1502 : IF (ABS(sep%aq) > EPSILON(0.0_dp) .AND. (sep%z /= 1)) THEN
469 1064 : sep%ko(3) = 0.5_dp/sep%aq
470 : END IF
471 1502 : sep%ko(7) = sep%ko(1)
472 1502 : sep%ko(9) = sep%ko(1)
473 1502 : sep%cs(2) = sep%dd
474 1502 : sep%cs(3) = sep%qq*SQRT(2.0_dp)
475 : ELSE
476 : ! Use the new integral module
477 738 : sep%ko(9) = sep%ko(1)
478 738 : 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 2240 : IF (ABS(sep%rho) > EPSILON(0.0_dp)) THEN
483 130 : sep%ko(9) = sep%rho
484 : END IF
485 2240 : 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 3464 : FUNCTION element_has_d(sep) RESULT(res)
496 : TYPE(semi_empirical_type), POINTER :: sep
497 : LOGICAL :: res
498 :
499 3464 : res = (nqd(sep%z) > 0 .AND. sep%sto_exponents(2) > EPSILON(0.0_dp))
500 3464 : 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 1732 : FUNCTION element_has_f(sep) RESULT(res)
511 : TYPE(semi_empirical_type), POINTER :: sep
512 : LOGICAL :: res
513 :
514 1732 : res = (nqf(sep%z) > 0 .AND. sep%sto_exponents(3) > EPSILON(0.0_dp))
515 1732 : 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 738 : 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 738 : z1 = sep%sto_exponents(0)
536 738 : z2 = sep%sto_exponents(1)
537 738 : z3 = sep%sto_exponents(2)
538 738 : 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 0 : "Please check if your parameterization supports the usage of s orbitals! ")
542 738 : amn = 0.0_dp
543 738 : nsp = nqs(sep%z)
544 738 : IF (sep%natorb >= 4) THEN
545 738 : 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 0 : "Please check if your parameterization supports the usage of p orbitals! ")
549 738 : amn(2) = amn_l_low(z1, z2, nsp, nsp, 1)
550 738 : amn(3) = amn_l_low(z2, z2, nsp, nsp, 2)
551 738 : IF (sep%dorb) THEN
552 738 : 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 0 : "Please check if your parameterization supports the usage of d orbitals! ")
556 738 : nd = nqd(sep%z)
557 738 : amn(4) = amn_l_low(z1, z3, nsp, nd, 2)
558 738 : amn(5) = amn_l_low(z2, z3, nsp, nd, 1)
559 738 : amn(6) = amn_l_low(z3, z3, nd, nd, 2)
560 : END IF
561 : END IF
562 738 : 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 2240 : 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 2240 : z1 = sep%sto_exponents(0)
583 2240 : z2 = sep%sto_exponents(1)
584 2240 : z3 = sep%sto_exponents(2)
585 2240 : CPASSERT(z1 > 0.0_dp)
586 2240 : amn = 0.0_dp
587 2240 : nsp = nqs(sep%z)
588 2240 : amn(1, 0) = amn_l_low(z1, z1, nsp, nsp, 0)
589 2240 : IF (sep%natorb >= 4) THEN
590 1772 : CPASSERT(z2 > 0.0_dp)
591 1772 : amn(2, 1) = amn_l_low(z1, z2, nsp, nsp, 1)
592 1772 : amn(3, 0) = amn_l_low(z2, z2, nsp, nsp, 0)
593 1772 : amn(3, 2) = amn_l_low(z2, z2, nsp, nsp, 2)
594 1772 : IF (sep%dorb) THEN
595 738 : CPASSERT(z3 > 0.0_dp)
596 738 : nd = nqd(sep%z)
597 738 : amn(4, 2) = amn_l_low(z1, z3, nsp, nd, 2)
598 738 : amn(5, 1) = amn_l_low(z2, z3, nsp, nd, 1)
599 738 : amn(6, 0) = amn_l_low(z3, z3, nd, nd, 0)
600 738 : amn(6, 2) = amn_l_low(z3, z3, nd, nd, 2)
601 : END IF
602 : END IF
603 2240 : 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 14198 : 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 14198 : (2.0_dp*z2/(z1 + z2))**n2*2.0_dp*SQRT(z1*z2)/(z1 + z2)**(l + 1)
623 :
624 14198 : 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 738 : 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 738 : fg = sep%gss
649 738 : sep%ko(1) = ko_ij(0, 1.0_dp, fg)
650 738 : IF (sep%natorb >= 4) THEN
651 : ! Other terms for SP basis
652 : ! SP
653 738 : d = amn(2)/SQRT(3.0_dp)
654 738 : fg = sep%hsp
655 738 : sep%cs(2) = d
656 738 : sep%ko(2) = ko_ij(1, d, fg)
657 : ! PP
658 738 : sep%ko(7) = sep%ko(1)
659 738 : d = SQRT(amn(3)*2.0_dp/5.0_dp)
660 738 : fg = 0.5_dp*(sep%gpp - sep%gp2)
661 738 : sep%cs(3) = d
662 738 : sep%ko(3) = ko_ij(2, d, fg)
663 : ! Terms involving d-orbitals
664 738 : IF (sep%dorb) THEN
665 : ! SD
666 738 : d = SQRT(amn(4)*2.0_dp/SQRT(15.0_dp))
667 738 : fg = sep%onec2el(19)
668 738 : sep%cs(4) = d
669 738 : sep%ko(4) = ko_ij(2, d, fg)
670 : ! PD
671 738 : d = amn(5)/SQRT(5.0_dp)
672 738 : fg = sep%onec2el(23) - 1.8_dp*sep%onec2el(35)
673 738 : sep%cs(5) = d
674 738 : sep%ko(5) = ko_ij(1, d, fg)
675 : ! DD
676 738 : fg = 0.2_dp*(sep%onec2el(29) + 2.0_dp*sep%onec2el(30) + 2.0_dp*sep%onec2el(31))
677 738 : sep%ko(8) = ko_ij(0, 1.0_dp, fg)
678 738 : d = SQRT(amn(6)*2.0_dp/7.0_dp)
679 738 : fg = sep%onec2el(44) - (20.0_dp/35.0_dp)*sep%onec2el(52)
680 738 : sep%cs(6) = d
681 738 : sep%ko(6) = ko_ij(2, d, fg)
682 : END IF
683 : END IF
684 738 : 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 738 : 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 738 : IF (sep%dorb) THEN
700 738 : s3 = SQRT(3.0_dp)
701 738 : s5 = SQRT(5.0_dp)
702 738 : 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 738 : r234, r246)
707 :
708 738 : IF (ABS(sep%f0sd) > EPSILON(0.0_dp)) THEN
709 436 : r016 = sep%f0sd
710 : END IF
711 738 : IF (ABS(sep%g2sd) > EPSILON(0.0_dp)) THEN
712 436 : r244 = sep%g2sd
713 : END IF
714 738 : CALL eisol_corr(sep, r016, r066, r244, r266, r466)
715 738 : sep%onec2el(1) = r016
716 738 : sep%onec2el(2) = 2.0_dp/(3.0_dp*s5)*r125
717 738 : sep%onec2el(3) = 1.0_dp/s15*r125
718 738 : sep%onec2el(4) = 2.0_dp/(5.0_dp*s5)*r234
719 738 : sep%onec2el(5) = r036 + 4.0_dp/35.0_dp*r236
720 738 : sep%onec2el(6) = r036 + 2.0_dp/35.0_dp*r236
721 738 : sep%onec2el(7) = r036 - 4.0_dp/35.0_dp*r236
722 738 : sep%onec2el(8) = -1.0_dp/(3.0_dp*s5)*r125
723 738 : sep%onec2el(9) = SQRT(3.0_dp/125.0_dp)*r234
724 738 : sep%onec2el(10) = s3/35.0_dp*r236
725 738 : sep%onec2el(11) = 3.0_dp/35.0_dp*r236
726 738 : sep%onec2el(12) = -1.0_dp/(5.0_dp*s5)*r234
727 738 : sep%onec2el(13) = r036 - 2.0_dp/35.0_dp*r236
728 738 : sep%onec2el(14) = -2.0_dp*s3/35.0_dp*r236
729 738 : sep%onec2el(15) = -sep%onec2el(3)
730 738 : sep%onec2el(16) = -sep%onec2el(11)
731 738 : sep%onec2el(17) = -sep%onec2el(9)
732 738 : sep%onec2el(18) = -sep%onec2el(14)
733 738 : sep%onec2el(19) = 1.0_dp/5.0_dp*r244
734 738 : sep%onec2el(20) = 2.0_dp/(7.0_dp*s5)*r246
735 738 : sep%onec2el(21) = sep%onec2el(20)/2.0_dp
736 738 : sep%onec2el(22) = -sep%onec2el(20)
737 738 : sep%onec2el(23) = 4.0_dp/15.0_dp*r155 + 27.0_dp/245.0_dp*r355
738 738 : sep%onec2el(24) = 2.0_dp*s3/15.0_dp*r155 - 9.0_dp*s3/245.0_dp*r355
739 738 : sep%onec2el(25) = 1.0_dp/15.0_dp*r155 + 18.0_dp/245.0_dp*r355
740 738 : sep%onec2el(26) = -s3/15.0_dp*r155 + 12.0_dp*s3/245.0_dp*r355
741 738 : sep%onec2el(27) = -s3/15.0_dp*r155 - 3.0_dp*s3/245.0_dp*r355
742 738 : sep%onec2el(28) = -sep%onec2el(27)
743 738 : sep%onec2el(29) = r066 + 4.0_dp/49.0_dp*r266 + 4.0_dp/49.0_dp*r466
744 738 : sep%onec2el(30) = r066 + 2.0_dp/49.0_dp*r266 - 24.0_dp/441.0_dp*r466
745 738 : sep%onec2el(31) = r066 - 4.0_dp/49.0_dp*r266 + 6.0_dp/441.0_dp*r466
746 738 : sep%onec2el(32) = SQRT(3.0_dp/245.0_dp)*r246
747 738 : sep%onec2el(33) = 1.0_dp/5.0_dp*r155 + 24.0_dp/245.0_dp*r355
748 738 : sep%onec2el(34) = 1.0_dp/5.0_dp*r155 - 6.0_dp/245.0_dp*r355
749 738 : sep%onec2el(35) = 3.0_dp/49.0_dp*r355
750 738 : sep%onec2el(36) = 1.0_dp/49.0_dp*r266 + 30.0_dp/441.0_dp*r466
751 738 : sep%onec2el(37) = s3/49.0_dp*r266 - 5.0_dp*s3/441.0_dp*r466
752 738 : sep%onec2el(38) = r066 - 2.0_dp/49.0_dp*r266 - 4.0_dp/441.0_dp*r466
753 738 : sep%onec2el(39) = -2.0_dp*s3/49.0_dp*r266 + 10.0_dp*s3/441.0_dp*r466
754 738 : sep%onec2el(40) = -sep%onec2el(32)
755 738 : sep%onec2el(41) = -sep%onec2el(34)
756 738 : sep%onec2el(42) = -sep%onec2el(35)
757 738 : sep%onec2el(43) = -sep%onec2el(37)
758 738 : sep%onec2el(44) = 3.0_dp/49.0_dp*r266 + 20.0_dp/441.0_dp*r466
759 738 : sep%onec2el(45) = -sep%onec2el(39)
760 738 : sep%onec2el(46) = 1.0_dp/5.0_dp*r155 - 3.0_dp/35.0_dp*r355
761 738 : sep%onec2el(47) = -sep%onec2el(46)
762 738 : sep%onec2el(48) = 4.0_dp/49.0_dp*r266 + 15.0_dp/441.0_dp*r466
763 738 : sep%onec2el(49) = 3.0_dp/49.0_dp*r266 - 5.0_dp/147.0_dp*r466
764 738 : sep%onec2el(50) = -sep%onec2el(49)
765 738 : sep%onec2el(51) = r066 + 4.0_dp/49.0_dp*r266 - 34.0_dp/441.0_dp*r466
766 738 : sep%onec2el(52) = 35.0_dp/441.0_dp*r466
767 738 : sep%f0dd = r066
768 738 : sep%f2dd = r266
769 738 : sep%f4dd = r466
770 738 : sep%f0sd = r016
771 738 : sep%g2sd = r244
772 738 : sep%f0pd = r036
773 738 : sep%f2pd = r236
774 738 : sep%g1pd = r155
775 738 : sep%g3pd = r355
776 : END IF
777 738 : 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 738 : 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 738 : ns = nqs(sep%z)
807 738 : nd = nqd(sep%z)
808 738 : es = sep%zn(0)
809 738 : ep = sep%zn(1)
810 738 : ed = sep%zn(2)
811 738 : r016 = sc_param_low(0, ns, es, ns, es, nd, ed, nd, ed)
812 738 : r036 = sc_param_low(0, ns, ep, ns, ep, nd, ed, nd, ed)
813 738 : r066 = sc_param_low(0, nd, ed, nd, ed, nd, ed, nd, ed)
814 738 : r155 = sc_param_low(1, ns, ep, nd, ed, ns, ep, nd, ed)
815 738 : r125 = sc_param_low(1, ns, es, ns, ep, ns, ep, nd, ed)
816 738 : r244 = sc_param_low(2, ns, es, nd, ed, ns, es, nd, ed)
817 738 : r236 = sc_param_low(2, ns, ep, ns, ep, nd, ed, nd, ed)
818 738 : r266 = sc_param_low(2, nd, ed, nd, ed, nd, ed, nd, ed)
819 738 : r234 = sc_param_low(2, ns, ep, ns, ep, ns, es, nd, ed)
820 738 : r246 = sc_param_low(2, ns, es, nd, ed, nd, ed, nd, ed)
821 738 : r355 = sc_param_low(3, ns, ep, nd, ed, ns, ep, nd, ed)
822 738 : r466 = sc_param_low(4, nd, ed, nd, ed, nd, ed, nd, ed)
823 738 : 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 8856 : 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 8856 : CPASSERT(ea > 0.0_dp)
862 8856 : CPASSERT(eb > 0.0_dp)
863 8856 : CPASSERT(ec > 0.0_dp)
864 8856 : CPASSERT(ed > 0.0_dp)
865 8856 : aea = LOG(ea)
866 8856 : aeb = LOG(eb)
867 8856 : aec = LOG(ec)
868 8856 : aed = LOG(ed)
869 8856 : nab = na + nb
870 8856 : ncd = nc + nd
871 8856 : ecd = ec + ed
872 8856 : eab = ea + eb
873 8856 : e = ecd + eab
874 8856 : n = nab + ncd
875 8856 : ae = LOG(e)
876 8856 : a2 = LOG(2.0_dp)
877 8856 : acd = LOG(ecd)
878 8856 : aab = LOG(eab)
879 8856 : ff = fac(n - 1)/SQRT(fac(2*na)*fac(2*nb)*fac(2*nc)*fac(2*nd))
880 8856 : tmp = na*aea + nb*aeb + nc*aec + nd*aed + 0.5_dp*(aea + aeb + aec + aed) + a2*(n + 2) - ae*n
881 8856 : c = evolt*ff*EXP(tmp)
882 8856 : s0 = 1.0_dp/e
883 8856 : s1 = 0.0_dp
884 8856 : s2 = 0.0_dp
885 8856 : m = ncd - k
886 57302 : DO i = 1, m
887 48446 : s0 = s0*e/ecd
888 57302 : s1 = s1 + s0*(binomial(ncd - k - 1, i - 1) - binomial(ncd + k, i - 1))/binomial(n - 1, i - 1)
889 : END DO
890 8856 : m1 = m
891 8856 : m2 = ncd + k
892 45756 : DO i = m1, m2
893 36900 : s0 = s0*e/ecd
894 45756 : s2 = s2 + s0*binomial(m2, i)/binomial(n - 1, i)
895 : END DO
896 8856 : s3 = EXP(ae*n - acd*(m2 + 1) - aab*(nab - k))/binomial(n - 1, m2)
897 8856 : res = c*(s1 - s2 + s3)
898 8856 : 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 738 : 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 738 : ir466(sep%z)*r466/49.0_dp
928 738 : 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 5166 : 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 5166 : CPASSERT(fg /= 0.0_dp)
955 : ! Term for SS
956 5166 : IF (l == 0) THEN
957 1476 : res = 0.5_dp*evolt/fg
958 1476 : RETURN
959 : END IF
960 : ! Term for Higher angular momentum
961 3690 : dsq = d*d
962 3690 : ev4 = evolt*0.25_dp
963 3690 : ev8 = evolt/8.0_dp
964 3690 : a1 = 0.1_dp
965 3690 : a2 = 5.0_dp
966 158670 : DO i = 1, niter
967 158670 : delta = a2 - a1
968 158670 : IF (delta < epsil) EXIT
969 154980 : y1 = a1 + delta*g1
970 154980 : y2 = a1 + delta*g2
971 154980 : IF (l == 1) THEN
972 61992 : f1 = (ev4*(1/y1 - 1/SQRT(y1**2 + dsq)) - fg)**2
973 61992 : f2 = (ev4*(1/y2 - 1/SQRT(y2**2 + dsq)) - fg)**2
974 92988 : ELSE IF (l == 2) THEN
975 92988 : 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 92988 : 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 158670 : IF (f1 < f2) THEN
979 : a2 = y2
980 : ELSE
981 74844 : a1 = y1
982 : END IF
983 : END DO
984 : ! Convergence reached.. define additive terms
985 3690 : IF (f1 >= f2) THEN
986 : res = a2
987 : ELSE
988 1832 : 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 3964 : 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 3964 : gss=gss, gsp=gsp, gpp=gpp, gp2=gp2, hsp=hsp)
1009 3964 : CPASSERT(defined)
1010 :
1011 3964 : isize = natorb*(natorb + 1)/2
1012 12408 : ALLOCATE (sep%w(isize, isize))
1013 : ! Initialize array
1014 1646300 : sep%w = 0.0_dp
1015 : ! Fill the array
1016 3964 : IF (natorb > 0) THEN
1017 2240 : ip = 1
1018 2240 : sep%w(ip, ip) = gss
1019 2240 : IF (natorb > 2) THEN
1020 1772 : ipx = ip + 2
1021 1772 : ipy = ip + 5
1022 1772 : ipz = ip + 9
1023 1772 : sep%w(ipx, ip) = gsp
1024 1772 : sep%w(ipy, ip) = gsp
1025 1772 : sep%w(ipz, ip) = gsp
1026 1772 : sep%w(ip, ipx) = gsp
1027 1772 : sep%w(ip, ipy) = gsp
1028 1772 : sep%w(ip, ipz) = gsp
1029 1772 : sep%w(ipx, ipx) = gpp
1030 1772 : sep%w(ipy, ipy) = gpp
1031 1772 : sep%w(ipz, ipz) = gpp
1032 1772 : sep%w(ipy, ipx) = gp2
1033 1772 : sep%w(ipz, ipx) = gp2
1034 1772 : sep%w(ipz, ipy) = gp2
1035 1772 : sep%w(ipx, ipy) = gp2
1036 1772 : sep%w(ipx, ipz) = gp2
1037 1772 : sep%w(ipy, ipz) = gp2
1038 1772 : sep%w(ip + 1, ip + 1) = hsp
1039 1772 : sep%w(ip + 3, ip + 3) = hsp
1040 1772 : sep%w(ip + 6, ip + 6) = hsp
1041 1772 : sep%w(ip + 4, ip + 4) = 0.5_dp*(gpp - gp2)
1042 1772 : sep%w(ip + 7, ip + 7) = 0.5_dp*(gpp - gp2)
1043 1772 : sep%w(ip + 8, ip + 8) = 0.5_dp*(gpp - gp2)
1044 1772 : IF (sep%dorb) THEN
1045 : ij0 = ip - 1
1046 180072 : DO i = 1, 243
1047 179334 : ij = int_ij(i)
1048 179334 : kl = int_kl(i)
1049 179334 : ind = int_onec2el(i)
1050 180072 : sep%w(ij + ij0, kl + ij0) = sep%onec2el(ind)/evolt
1051 : END DO
1052 : END IF
1053 : END IF
1054 : END IF
1055 3964 : END SUBROUTINE setup_1c_2el_int
1056 :
1057 : END MODULE semi_empirical_par_utils
1058 :
|
