Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 2104 : 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 2104 : 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 2104 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
68 : TYPE(dft_control_type), POINTER :: dft_control
69 2104 : 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 2104 : dft_control=dft_control)
77 2104 : xtb_control => dft_control%qs_control%xtb_control
78 :
79 2104 : CALL get_qs_env(qs_env=qs_env, natom=natom)
80 :
81 6312 : ALLOCATE (huckel(5, natom))
82 2104 : IF (calculate_forces) THEN
83 378 : ALLOCATE (dhuckel(5, natom), dqhuckel(5, natom))
84 : END IF
85 2104 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
86 12048 : DO iatom = 1, natom
87 9944 : ikind = kind_of(iatom)
88 9944 : 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 9944 : kcn=kcn, kq=kq, kqat2=kqat2, hen=hena)
91 59664 : kcn = kcn/evolt
92 59664 : kq = kq/evolt
93 9944 : kqat2 = kqat2/evolt
94 59664 : huckel(:, iatom) = 0.0_dp
95 33494 : DO i = 1, nshell
96 23550 : l = lval(i) + 1
97 : huckel(i, iatom) = hena(i) - kcn(l)*cnumbers(iatom) &
98 33494 : - kq(l)*charges(iatom) - kqat2*charges(iatom)**2
99 : END DO
100 21992 : IF (calculate_forces) THEN
101 3480 : dhuckel(:, iatom) = 0.0_dp
102 3480 : dqhuckel(:, iatom) = 0.0_dp
103 1896 : DO i = 1, nshell
104 1316 : l = lval(i) + 1
105 1316 : dhuckel(i, iatom) = -kcn(l)
106 1896 : dqhuckel(i, iatom) = -kq(l) - 2.0_dp*kqat2*charges(iatom)
107 : END DO
108 : END IF
109 : END DO
110 :
111 4208 : 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 3986 : 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 3986 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
129 : INTEGER, DIMENSION(25) :: lval
130 : REAL(KIND=dp) :: kcnd, kcnp, kcns
131 3986 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: kcnlk
132 : REAL(KIND=dp), DIMENSION(5) :: hena
133 3986 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
134 : TYPE(dft_control_type), POINTER :: dft_control
135 3986 : 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 3986 : dft_control=dft_control)
143 3986 : xtb_control => dft_control%qs_control%xtb_control
144 :
145 3986 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom)
146 :
147 3986 : kcns = xtb_control%kcns
148 3986 : kcnp = xtb_control%kcnp
149 3986 : kcnd = xtb_control%kcnd
150 :
151 : ! Calculate Huckel parameters
152 : ! Eq 12
153 : ! huckel(nshell,natom)
154 11958 : ALLOCATE (kcnlk(0:3, nkind))
155 13634 : DO ikind = 1, nkind
156 9648 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
157 13634 : IF (metal(za)) THEN
158 560 : kcnlk(0:3, ikind) = 0.0_dp
159 9536 : 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 9536 : kcnlk(0, ikind) = kcns
166 9536 : kcnlk(1, ikind) = kcnp
167 9536 : kcnlk(2, ikind) = kcnd
168 9536 : kcnlk(3, ikind) = 0.0_dp
169 : END IF
170 : END DO
171 :
172 11958 : ALLOCATE (huckel(5, natom))
173 3986 : IF (calculate_forces) THEN
174 1084 : ALLOCATE (dhuckel(5, natom))
175 : END IF
176 3986 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
177 36000 : DO iatom = 1, natom
178 32014 : ikind = kind_of(iatom)
179 32014 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
180 32014 : CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
181 192084 : huckel(:, iatom) = 0.0_dp
182 97366 : DO i = 1, nshell
183 97366 : huckel(i, iatom) = hena(i)*(1._dp + kcnlk(lval(i), ikind)*cnumbers(iatom))
184 : END DO
185 68014 : IF (calculate_forces) THEN
186 27744 : dhuckel(:, iatom) = 0.0_dp
187 14348 : DO i = 1, nshell
188 14348 : dhuckel(i, iatom) = hena(i)*kcnlk(lval(i), ikind)
189 : END DO
190 : END IF
191 : END DO
192 :
193 3986 : DEALLOCATE (kcnlk)
194 :
195 7972 : END SUBROUTINE gfn1_huckel
196 :
197 : ! **************************************************************************************************
198 : !> \brief ...
199 : !> \param qs_env ...
200 : !> \param kijab ...
201 : ! **************************************************************************************************
202 2104 : 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 2104 : 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 2104 : dft_control=dft_control)
223 2104 : xtb_control => dft_control%qs_control%xtb_control
224 :
225 2104 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
226 2104 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
227 :
228 2104 : ks = xtb_control%ks
229 2104 : kp = xtb_control%kp
230 2104 : kd = xtb_control%kd
231 2104 : ksp = xtb_control%ksp
232 2104 : ksen = xtb_control%ksen
233 2104 : kpen = xtb_control%kpen
234 2104 : kden = xtb_control%kden
235 2104 : ben = xtb_control%ben
236 2104 : kdiff = xtb_control%k2sh
237 :
238 2104 : kl(0) = ks
239 2104 : kl(1) = kp
240 2104 : kl(2) = kd
241 2104 : kl(3) = 0.0_dp
242 :
243 2104 : ke(0) = ksen
244 2104 : ke(1) = kpen
245 2104 : ke(2) = kden
246 2104 : ke(3) = 0.0_dp
247 :
248 : ! Calculate KAB parameters and electronegativity correction
249 12624 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
250 2104 : kijab = 0.0_dp
251 :
252 7012 : DO ikind = 1, nkind
253 4908 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
254 4908 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
255 4908 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
256 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, &
257 4908 : en=etaa, zeta=zetaa)
258 24228 : DO jkind = 1, nkind
259 12308 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
260 12308 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
261 12308 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
262 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, &
263 12308 : en=etab, zeta=zetab)
264 : ! Kab
265 12308 : kab = pp_gfn0(za, zb)
266 63870 : DO j = 1, natorb_b
267 51562 : lb = laob(j)
268 51562 : nb = naob(j)
269 310024 : DO i = 1, natorb_a
270 246154 : la = laoa(i)
271 246154 : na = naoa(i)
272 246154 : kia = kl(la)
273 246154 : kjb = kl(lb)
274 246154 : km = 0.5_dp*(kia + kjb)*kab
275 246154 : IF (za == 1 .AND. na == 2) THEN
276 11022 : IF (zb == 1 .AND. nb == 2) THEN
277 : km = 0._dp
278 : ELSE
279 9808 : km = km*kdiff
280 : END IF
281 235132 : ELSEIF (zb == 1 .AND. nb == 2) THEN
282 9808 : km = km*kdiff
283 : END IF
284 297716 : kijab(i, j, ikind, jkind) = km
285 : END DO
286 : END DO
287 : ! Yab
288 63870 : DO j = 1, natorb_b
289 51562 : nb = naob(j)
290 51562 : kjb = zetab(nb)
291 310024 : DO i = 1, natorb_a
292 246154 : na = naoa(i)
293 246154 : kia = zetaa(na)
294 246154 : yijab = 2.0_dp*SQRT(kia*kjb)/(kia + kjb)
295 297716 : kijab(i, j, ikind, jkind) = kijab(i, j, ikind, jkind)*yijab
296 : END DO
297 : END DO
298 : ! X
299 12308 : den = etaa - etab
300 93394 : DO j = 1, natorb_b
301 51562 : lb = laob(j)
302 51562 : kjb = ke(lb)
303 310024 : DO i = 1, natorb_a
304 246154 : la = laoa(i)
305 246154 : kia = ke(la)
306 246154 : ken = 0.5_dp*(kia + kjb)
307 246154 : xijab = 1.0_dp + ken*den**2 + ken*ben*den**4
308 297716 : 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 4208 : END SUBROUTINE gfn0_kpair
315 :
316 : ! **************************************************************************************************
317 : !> \brief ...
318 : !> \param qs_env ...
319 : !> \param kijab ...
320 : ! **************************************************************************************************
321 3986 : 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 3986 : 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 3986 : dft_control=dft_control)
340 3986 : xtb_control => dft_control%qs_control%xtb_control
341 :
342 3986 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
343 3986 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
344 :
345 3986 : ks = xtb_control%ks
346 3986 : kp = xtb_control%kp
347 3986 : kd = xtb_control%kd
348 3986 : ksp = xtb_control%ksp
349 3986 : k2sh = xtb_control%k2sh
350 3986 : ken = xtb_control%ken
351 :
352 3986 : kl(0) = ks
353 3986 : kl(1) = kp
354 3986 : kl(2) = kd
355 3986 : 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 23916 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
360 3986 : kijab = 0.0_dp
361 13634 : DO ikind = 1, nkind
362 9648 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
363 9648 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
364 9648 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
365 9646 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
366 49464 : DO jkind = 1, nkind
367 26184 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
368 26184 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
369 26184 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
370 26182 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
371 : ! get Fen = (1+ken*deltaEN^2)
372 26182 : fen = 1.0_dp + ken*(ena - enb)**2
373 : ! Kab
374 26182 : kab = xtb_set_kab(za, zb, xtb_control)
375 155092 : DO j = 1, natorb_b
376 93078 : lb = laob(j)
377 93078 : nb = naob(j)
378 477836 : DO i = 1, natorb_a
379 358574 : la = laoa(i)
380 358574 : na = naoa(i)
381 358574 : kia = kl(la)
382 358574 : kjb = kl(lb)
383 358574 : IF (zb == 1 .AND. nb == 2) kjb = k2sh
384 358574 : IF (za == 1 .AND. na == 2) kia = k2sh
385 451652 : IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
386 57410 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
387 : ELSE
388 301164 : IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
389 101352 : kijab(i, j, ikind, jkind) = ksp*kab*fen
390 : ELSE
391 199812 : 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 7972 : END SUBROUTINE gfn1_kpair
400 :
401 : END MODULE xtb_hcore
402 :
|
