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 Calculation of EHT matrix elements in xTB
10 : !> Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
11 : !> JCTC 13, 1989-2009, (2017)
12 : !> DOI: 10.1021/acs.jctc.7b00118
13 : !> \author JGH
14 : ! **************************************************************************************************
15 : MODULE xtb_hcore
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind,&
18 : get_atomic_kind_set
19 : USE cp_control_types, ONLY: dft_control_type,&
20 : xtb_control_type
21 : USE kinds, ONLY: dp
22 : USE physcon, ONLY: evolt
23 : USE qs_environment_types, ONLY: get_qs_env,&
24 : qs_environment_type
25 : USE qs_kind_types, ONLY: get_qs_kind,&
26 : get_qs_kind_set,&
27 : qs_kind_type
28 : USE xtb_parameters, ONLY: early3d,&
29 : metal,&
30 : pp_gfn0,&
31 : xtb_set_kab
32 : USE xtb_types, ONLY: get_xtb_atom_param,&
33 : xtb_atom_type
34 : #include "./base/base_uses.f90"
35 :
36 : IMPLICIT NONE
37 :
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_hcore'
41 :
42 : PUBLIC :: gfn0_huckel, gfn1_huckel, gfn0_kpair, gfn1_kpair
43 :
44 : CONTAINS
45 :
46 : ! **************************************************************************************************
47 : !> \brief ...
48 : !> \param qs_env ...
49 : !> \param cnumbers ...
50 : !> \param charges ...
51 : !> \param huckel ...
52 : !> \param dhuckel ...
53 : !> \param dqhuckel ...
54 : !> \param calculate_forces ...
55 : ! **************************************************************************************************
56 1702 : SUBROUTINE gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
57 : TYPE(qs_environment_type), POINTER :: qs_env
58 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cnumbers, charges
59 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: huckel, dhuckel, dqhuckel
60 : LOGICAL, INTENT(IN) :: calculate_forces
61 :
62 : INTEGER :: i, iatom, ikind, l, natom, nshell
63 1702 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
64 : INTEGER, DIMENSION(25) :: lval
65 : REAL(KIND=dp) :: kqat2
66 : REAL(KIND=dp), DIMENSION(5) :: hena, kcn, kq
67 1702 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
68 : TYPE(dft_control_type), POINTER :: dft_control
69 1702 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
70 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a
71 : TYPE(xtb_control_type), POINTER :: xtb_control
72 :
73 : CALL get_qs_env(qs_env=qs_env, &
74 : atomic_kind_set=atomic_kind_set, &
75 : qs_kind_set=qs_kind_set, &
76 1702 : dft_control=dft_control)
77 1702 : xtb_control => dft_control%qs_control%xtb_control
78 :
79 1702 : CALL get_qs_env(qs_env=qs_env, natom=natom)
80 :
81 5106 : ALLOCATE (huckel(5, natom))
82 1702 : IF (calculate_forces) THEN
83 490 : ALLOCATE (dhuckel(5, natom), dqhuckel(5, natom))
84 : END IF
85 1702 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
86 9676 : DO iatom = 1, natom
87 7974 : ikind = kind_of(iatom)
88 7974 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
89 : CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, &
90 7974 : kcn=kcn, kq=kq, kqat2=kqat2, hen=hena)
91 47844 : kcn = kcn/evolt
92 47844 : kq = kq/evolt
93 7974 : kqat2 = kqat2/evolt
94 47844 : huckel(:, iatom) = 0.0_dp
95 26368 : DO i = 1, nshell
96 18394 : l = lval(i) + 1
97 : huckel(i, iatom) = hena(i) - kcn(l)*cnumbers(iatom) &
98 26368 : - kq(l)*charges(iatom) - kqat2*charges(iatom)**2
99 : END DO
100 17650 : IF (calculate_forces) THEN
101 2544 : dhuckel(:, iatom) = 0.0_dp
102 2544 : dqhuckel(:, iatom) = 0.0_dp
103 1316 : DO i = 1, nshell
104 892 : l = lval(i) + 1
105 892 : dhuckel(i, iatom) = -kcn(l)
106 1316 : dqhuckel(i, iatom) = -kq(l) - 2.0_dp*kqat2*charges(iatom)
107 : END DO
108 : END IF
109 : END DO
110 :
111 3404 : END SUBROUTINE gfn0_huckel
112 :
113 : ! **************************************************************************************************
114 : !> \brief ...
115 : !> \param qs_env ...
116 : !> \param cnumbers ...
117 : !> \param huckel ...
118 : !> \param dhuckel ...
119 : !> \param calculate_forces ...
120 : ! **************************************************************************************************
121 3454 : SUBROUTINE gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
122 : TYPE(qs_environment_type), POINTER :: qs_env
123 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cnumbers
124 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: huckel, dhuckel
125 : LOGICAL, INTENT(IN) :: calculate_forces
126 :
127 : INTEGER :: i, iatom, ikind, natom, nkind, nshell, za
128 3454 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
129 : INTEGER, DIMENSION(25) :: lval
130 : REAL(KIND=dp) :: kcnd, kcnp, kcns
131 3454 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: kcnlk
132 : REAL(KIND=dp), DIMENSION(5) :: hena
133 3454 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
134 : TYPE(dft_control_type), POINTER :: dft_control
135 3454 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
136 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a
137 : TYPE(xtb_control_type), POINTER :: xtb_control
138 :
139 : CALL get_qs_env(qs_env=qs_env, &
140 : atomic_kind_set=atomic_kind_set, &
141 : qs_kind_set=qs_kind_set, &
142 3454 : dft_control=dft_control)
143 3454 : xtb_control => dft_control%qs_control%xtb_control
144 :
145 3454 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom)
146 :
147 3454 : kcns = xtb_control%kcns
148 3454 : kcnp = xtb_control%kcnp
149 3454 : kcnd = xtb_control%kcnd
150 :
151 : ! Calculate Huckel parameters
152 : ! Eq 12
153 : ! huckel(nshell,natom)
154 10362 : ALLOCATE (kcnlk(0:3, nkind))
155 12082 : DO ikind = 1, nkind
156 8628 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
157 12082 : IF (metal(za)) THEN
158 250 : kcnlk(0:3, ikind) = 0.0_dp
159 8578 : ELSEIF (early3d(za)) THEN
160 0 : kcnlk(0, ikind) = kcns
161 0 : kcnlk(1, ikind) = kcnp
162 0 : kcnlk(2, ikind) = 0.005_dp
163 0 : kcnlk(3, ikind) = 0.0_dp
164 : ELSE
165 8578 : kcnlk(0, ikind) = kcns
166 8578 : kcnlk(1, ikind) = kcnp
167 8578 : kcnlk(2, ikind) = kcnd
168 8578 : kcnlk(3, ikind) = 0.0_dp
169 : END IF
170 : END DO
171 :
172 10362 : ALLOCATE (huckel(5, natom))
173 3454 : IF (calculate_forces) THEN
174 920 : ALLOCATE (dhuckel(5, natom))
175 : END IF
176 3454 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
177 32984 : DO iatom = 1, natom
178 29530 : ikind = kind_of(iatom)
179 29530 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
180 29530 : CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
181 177180 : huckel(:, iatom) = 0.0_dp
182 89144 : DO i = 1, nshell
183 89144 : huckel(i, iatom) = hena(i)*(1._dp + kcnlk(lval(i), ikind)*cnumbers(iatom))
184 : END DO
185 62514 : IF (calculate_forces) THEN
186 24564 : dhuckel(:, iatom) = 0.0_dp
187 12502 : DO i = 1, nshell
188 12502 : dhuckel(i, iatom) = hena(i)*kcnlk(lval(i), ikind)
189 : END DO
190 : END IF
191 : END DO
192 :
193 3454 : DEALLOCATE (kcnlk)
194 :
195 6908 : END SUBROUTINE gfn1_huckel
196 :
197 : ! **************************************************************************************************
198 : !> \brief ...
199 : !> \param qs_env ...
200 : !> \param kijab ...
201 : ! **************************************************************************************************
202 1702 : SUBROUTINE gfn0_kpair(qs_env, kijab)
203 : TYPE(qs_environment_type), POINTER :: qs_env
204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
205 :
206 : INTEGER :: i, ikind, j, jkind, la, lb, maxs, na, &
207 : natorb_a, natorb_b, nb, nkind, za, zb
208 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
209 : LOGICAL :: defined
210 : REAL(KIND=dp) :: ben, den, etaa, etab, kab, kd, kden, &
211 : kdiff, ken, kia, kjb, km, kp, kpen, &
212 : ks, ksen, ksp, xijab, yijab
213 : REAL(KIND=dp), DIMENSION(0:3) :: ke, kl
214 : REAL(KIND=dp), DIMENSION(5) :: zetaa, zetab
215 : TYPE(dft_control_type), POINTER :: dft_control
216 1702 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
217 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
218 : TYPE(xtb_control_type), POINTER :: xtb_control
219 :
220 : CALL get_qs_env(qs_env=qs_env, &
221 : qs_kind_set=qs_kind_set, &
222 1702 : dft_control=dft_control)
223 1702 : xtb_control => dft_control%qs_control%xtb_control
224 :
225 1702 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
226 1702 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
227 :
228 1702 : ks = xtb_control%ks
229 1702 : kp = xtb_control%kp
230 1702 : kd = xtb_control%kd
231 1702 : ksp = xtb_control%ksp
232 1702 : ksen = xtb_control%ksen
233 1702 : kpen = xtb_control%kpen
234 1702 : kden = xtb_control%kden
235 1702 : ben = xtb_control%ben
236 1702 : kdiff = xtb_control%k2sh
237 :
238 1702 : kl(0) = ks
239 1702 : kl(1) = kp
240 1702 : kl(2) = kd
241 1702 : kl(3) = 0.0_dp
242 :
243 1702 : ke(0) = ksen
244 1702 : ke(1) = kpen
245 1702 : ke(2) = kden
246 1702 : ke(3) = 0.0_dp
247 :
248 : ! Calculate KAB parameters and electronegativity correction
249 10212 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
250 457762 : kijab = 0.0_dp
251 :
252 5954 : DO ikind = 1, nkind
253 4252 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
254 4252 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
255 4252 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
256 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, &
257 4252 : en=etaa, zeta=zetaa)
258 21342 : DO jkind = 1, nkind
259 11136 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
260 11136 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
261 11136 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
262 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, &
263 11136 : en=etab, zeta=zetab)
264 : ! Kab
265 11136 : kab = pp_gfn0(za, zb)
266 58258 : DO j = 1, natorb_b
267 47122 : lb = laob(j)
268 47122 : nb = naob(j)
269 282844 : DO i = 1, natorb_a
270 224586 : la = laoa(i)
271 224586 : na = naoa(i)
272 224586 : kia = kl(la)
273 224586 : kjb = kl(lb)
274 224586 : km = 0.5_dp*(kia + kjb)*kab
275 224586 : IF (za == 1 .AND. na == 2) THEN
276 9506 : IF (zb == 1 .AND. nb == 2) THEN
277 : km = 0._dp
278 : ELSE
279 8542 : km = km*kdiff
280 : END IF
281 215080 : ELSEIF (zb == 1 .AND. nb == 2) THEN
282 8542 : km = km*kdiff
283 : END IF
284 271708 : kijab(i, j, ikind, jkind) = km
285 : END DO
286 : END DO
287 : ! Yab
288 58258 : DO j = 1, natorb_b
289 47122 : nb = naob(j)
290 47122 : kjb = zetab(nb)
291 282844 : DO i = 1, natorb_a
292 224586 : na = naoa(i)
293 224586 : kia = zetaa(na)
294 224586 : yijab = 2.0_dp*SQRT(kia*kjb)/(kia + kjb)
295 271708 : kijab(i, j, ikind, jkind) = kijab(i, j, ikind, jkind)*yijab
296 : END DO
297 : END DO
298 : ! X
299 11136 : den = etaa - etab
300 84782 : DO j = 1, natorb_b
301 47122 : lb = laob(j)
302 47122 : kjb = ke(lb)
303 282844 : DO i = 1, natorb_a
304 224586 : la = laoa(i)
305 224586 : kia = ke(la)
306 224586 : ken = 0.5_dp*(kia + kjb)
307 224586 : xijab = 1.0_dp + ken*den**2 + ken*ben*den**4
308 271708 : kijab(i, j, ikind, jkind) = kijab(i, j, ikind, jkind)*xijab
309 : END DO
310 : END DO
311 : END DO
312 : END DO
313 :
314 3404 : END SUBROUTINE gfn0_kpair
315 :
316 : ! **************************************************************************************************
317 : !> \brief ...
318 : !> \param qs_env ...
319 : !> \param kijab ...
320 : ! **************************************************************************************************
321 3454 : SUBROUTINE gfn1_kpair(qs_env, kijab)
322 : TYPE(qs_environment_type), POINTER :: qs_env
323 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
324 :
325 : INTEGER :: i, ikind, j, jkind, la, lb, maxs, na, &
326 : natorb_a, natorb_b, nb, nkind, za, zb
327 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
328 : LOGICAL :: defined
329 : REAL(KIND=dp) :: ena, enb, fen, k2sh, kab, kd, ken, kia, &
330 : kjb, kp, ks, ksp
331 : REAL(KIND=dp), DIMENSION(0:3) :: kl
332 : TYPE(dft_control_type), POINTER :: dft_control
333 3454 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
334 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
335 : TYPE(xtb_control_type), POINTER :: xtb_control
336 :
337 : CALL get_qs_env(qs_env=qs_env, &
338 : qs_kind_set=qs_kind_set, &
339 3454 : dft_control=dft_control)
340 3454 : xtb_control => dft_control%qs_control%xtb_control
341 :
342 3454 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
343 3454 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
344 :
345 3454 : ks = xtb_control%ks
346 3454 : kp = xtb_control%kp
347 3454 : kd = xtb_control%kd
348 3454 : ksp = xtb_control%ksp
349 3454 : k2sh = xtb_control%k2sh
350 3454 : ken = xtb_control%ken
351 :
352 3454 : kl(0) = ks
353 3454 : kl(1) = kp
354 3454 : kl(2) = kd
355 3454 : kl(3) = 0.0_dp
356 :
357 : ! Calculate KAB parameters and electronegativity correction
358 : ! kijab -> K_l_l'[A,B] * X_l_l'[ENa, ENb] * Y[xia, xib]
359 20724 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
360 646114 : kijab = 0.0_dp
361 12082 : DO ikind = 1, nkind
362 8628 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
363 8628 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
364 8628 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
365 8626 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
366 44812 : DO jkind = 1, nkind
367 24104 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
368 24104 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
369 24104 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
370 24102 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
371 : ! get Fen = (1+ken*deltaEN^2)
372 24102 : fen = 1.0_dp + ken*(ena - enb)**2
373 : ! Kab
374 24102 : kab = xtb_set_kab(za, zb, xtb_control)
375 141490 : DO j = 1, natorb_b
376 84656 : lb = laob(j)
377 84656 : nb = naob(j)
378 426994 : DO i = 1, natorb_a
379 318234 : la = laoa(i)
380 318234 : na = naoa(i)
381 318234 : kia = kl(la)
382 318234 : kjb = kl(lb)
383 318234 : IF (zb == 1 .AND. nb == 2) kjb = k2sh
384 318234 : IF (za == 1 .AND. na == 2) kia = k2sh
385 402890 : IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
386 53378 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
387 : ELSE
388 264856 : IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
389 93156 : kijab(i, j, ikind, jkind) = ksp*kab*fen
390 : ELSE
391 171700 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)*kab*fen
392 : END IF
393 : END IF
394 : END DO
395 : END DO
396 : END DO
397 : END DO
398 :
399 6908 : END SUBROUTINE gfn1_kpair
400 :
401 : END MODULE xtb_hcore
402 :
|
