Line data Source code
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 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 1436 : 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 1436 : 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 1436 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
68 : TYPE(dft_control_type), POINTER :: dft_control
69 1436 : 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 1436 : dft_control=dft_control)
77 1436 : xtb_control => dft_control%qs_control%xtb_control
78 :
79 1436 : CALL get_qs_env(qs_env=qs_env, natom=natom)
80 :
81 4308 : ALLOCATE (huckel(5, natom))
82 1436 : IF (calculate_forces) THEN
83 140 : ALLOCATE (dhuckel(5, natom), dqhuckel(5, natom))
84 : END IF
85 1436 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
86 8346 : DO iatom = 1, natom
87 6910 : ikind = kind_of(iatom)
88 6910 : 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 6910 : kcn=kcn, kq=kq, kqat2=kqat2, hen=hena)
91 41460 : kcn = kcn/evolt
92 41460 : kq = kq/evolt
93 6910 : kqat2 = kqat2/evolt
94 41460 : huckel(:, iatom) = 0.0_dp
95 23176 : DO i = 1, nshell
96 16266 : l = lval(i) + 1
97 : huckel(i, iatom) = hena(i) - kcn(l)*cnumbers(iatom) &
98 23176 : - kq(l)*charges(iatom) - kqat2*charges(iatom)**2
99 : END DO
100 15256 : IF (calculate_forces) THEN
101 864 : dhuckel(:, iatom) = 0.0_dp
102 864 : dqhuckel(:, iatom) = 0.0_dp
103 476 : DO i = 1, nshell
104 332 : l = lval(i) + 1
105 332 : dhuckel(i, iatom) = -kcn(l)
106 476 : dqhuckel(i, iatom) = -kq(l) - 2.0_dp*kqat2*charges(iatom)
107 : END DO
108 : END IF
109 : END DO
110 :
111 2872 : 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 3202 : 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 3202 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
129 : INTEGER, DIMENSION(25) :: lval
130 : REAL(KIND=dp) :: kcnd, kcnp, kcns
131 3202 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: kcnlk
132 : REAL(KIND=dp), DIMENSION(5) :: hena
133 3202 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
134 : TYPE(dft_control_type), POINTER :: dft_control
135 3202 : 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 3202 : dft_control=dft_control)
143 3202 : xtb_control => dft_control%qs_control%xtb_control
144 :
145 3202 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom)
146 :
147 3202 : kcns = xtb_control%kcns
148 3202 : kcnp = xtb_control%kcnp
149 3202 : kcnd = xtb_control%kcnd
150 :
151 : ! Calculate Huckel parameters
152 : ! Eq 12
153 : ! huckel(nshell,natom)
154 9606 : ALLOCATE (kcnlk(0:3, nkind))
155 11074 : DO ikind = 1, nkind
156 7872 : CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
157 11074 : IF (metal(za)) THEN
158 250 : kcnlk(0:3, ikind) = 0.0_dp
159 7822 : 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 7822 : kcnlk(0, ikind) = kcns
166 7822 : kcnlk(1, ikind) = kcnp
167 7822 : kcnlk(2, ikind) = kcnd
168 7822 : kcnlk(3, ikind) = 0.0_dp
169 : END IF
170 : END DO
171 :
172 9606 : ALLOCATE (huckel(5, natom))
173 3202 : IF (calculate_forces) THEN
174 888 : ALLOCATE (dhuckel(5, natom))
175 : END IF
176 3202 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
177 31724 : DO iatom = 1, natom
178 28522 : ikind = kind_of(iatom)
179 28522 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
180 28522 : CALL get_xtb_atom_param(xtb_atom_a, nshell=nshell, lval=lval, hen=hena)
181 171132 : huckel(:, iatom) = 0.0_dp
182 86120 : DO i = 1, nshell
183 86120 : huckel(i, iatom) = hena(i)*(1._dp + kcnlk(lval(i), ikind)*cnumbers(iatom))
184 : END DO
185 60246 : IF (calculate_forces) THEN
186 24180 : dhuckel(:, iatom) = 0.0_dp
187 12310 : DO i = 1, nshell
188 12310 : dhuckel(i, iatom) = hena(i)*kcnlk(lval(i), ikind)
189 : END DO
190 : END IF
191 : END DO
192 :
193 3202 : DEALLOCATE (kcnlk)
194 :
195 6404 : END SUBROUTINE gfn1_huckel
196 :
197 : ! **************************************************************************************************
198 : !> \brief ...
199 : !> \param qs_env ...
200 : !> \param kijab ...
201 : ! **************************************************************************************************
202 1436 : 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 1436 : 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 1436 : dft_control=dft_control)
223 1436 : xtb_control => dft_control%qs_control%xtb_control
224 :
225 1436 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
226 1436 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
227 :
228 1436 : ks = xtb_control%ks
229 1436 : kp = xtb_control%kp
230 1436 : kd = xtb_control%kd
231 1436 : ksp = xtb_control%ksp
232 1436 : ksen = xtb_control%ksen
233 1436 : kpen = xtb_control%kpen
234 1436 : kden = xtb_control%kden
235 1436 : ben = xtb_control%ben
236 1436 : kdiff = xtb_control%k2sh
237 :
238 1436 : kl(0) = ks
239 1436 : kl(1) = kp
240 1436 : kl(2) = kd
241 1436 : kl(3) = 0.0_dp
242 :
243 1436 : ke(0) = ksen
244 1436 : ke(1) = kpen
245 1436 : ke(2) = kden
246 1436 : ke(3) = 0.0_dp
247 :
248 : ! Calculate KAB parameters and electronegativity correction
249 8616 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
250 406424 : kijab = 0.0_dp
251 :
252 4890 : DO ikind = 1, nkind
253 3454 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
254 3454 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
255 3454 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
256 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, &
257 3454 : en=etaa, zeta=zetaa)
258 17086 : DO jkind = 1, nkind
259 8742 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
260 8742 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
261 8742 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
262 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, &
263 8742 : en=etab, zeta=zetab)
264 : ! Kab
265 8742 : kab = pp_gfn0(za, zb)
266 47884 : DO j = 1, natorb_b
267 39142 : lb = laob(j)
268 39142 : nb = naob(j)
269 245870 : DO i = 1, natorb_a
270 197986 : la = laoa(i)
271 197986 : na = naoa(i)
272 197986 : kia = kl(la)
273 197986 : kjb = kl(lb)
274 197986 : km = 0.5_dp*(kia + kjb)*kab
275 197986 : IF (za == 1 .AND. na == 2) THEN
276 6846 : IF (zb == 1 .AND. nb == 2) THEN
277 : km = 0._dp
278 : ELSE
279 6148 : km = km*kdiff
280 : END IF
281 191140 : ELSEIF (zb == 1 .AND. nb == 2) THEN
282 6148 : km = km*kdiff
283 : END IF
284 237128 : kijab(i, j, ikind, jkind) = km
285 : END DO
286 : END DO
287 : ! Yab
288 47884 : DO j = 1, natorb_b
289 39142 : nb = naob(j)
290 39142 : kjb = zetab(nb)
291 245870 : DO i = 1, natorb_a
292 197986 : na = naoa(i)
293 197986 : kia = zetaa(na)
294 197986 : yijab = 2.0_dp*SQRT(kia*kjb)/(kia + kjb)
295 237128 : kijab(i, j, ikind, jkind) = kijab(i, j, ikind, jkind)*yijab
296 : END DO
297 : END DO
298 : ! X
299 8742 : den = etaa - etab
300 68822 : DO j = 1, natorb_b
301 39142 : lb = laob(j)
302 39142 : kjb = ke(lb)
303 245870 : DO i = 1, natorb_a
304 197986 : la = laoa(i)
305 197986 : kia = ke(la)
306 197986 : ken = 0.5_dp*(kia + kjb)
307 197986 : xijab = 1.0_dp + ken*den**2 + ken*ben*den**4
308 237128 : 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 2872 : END SUBROUTINE gfn0_kpair
315 :
316 : ! **************************************************************************************************
317 : !> \brief ...
318 : !> \param qs_env ...
319 : !> \param kijab ...
320 : ! **************************************************************************************************
321 3202 : 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 3202 : 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 3202 : dft_control=dft_control)
340 3202 : xtb_control => dft_control%qs_control%xtb_control
341 :
342 3202 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
343 3202 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxsgf=maxs, basis_type="ORB")
344 :
345 3202 : ks = xtb_control%ks
346 3202 : kp = xtb_control%kp
347 3202 : kd = xtb_control%kd
348 3202 : ksp = xtb_control%ksp
349 3202 : k2sh = xtb_control%k2sh
350 3202 : ken = xtb_control%ken
351 :
352 3202 : kl(0) = ks
353 3202 : kl(1) = kp
354 3202 : kl(2) = kd
355 3202 : 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 19212 : ALLOCATE (kijab(maxs, maxs, nkind, nkind))
360 597478 : kijab = 0.0_dp
361 11074 : DO ikind = 1, nkind
362 7872 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
363 7872 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
364 7872 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
365 7870 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, electronegativity=ena)
366 40780 : DO jkind = 1, nkind
367 21836 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
368 21836 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
369 21836 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
370 21834 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, electronegativity=enb)
371 : ! get Fen = (1+ken*deltaEN^2)
372 21834 : fen = 1.0_dp + ken*(ena - enb)**2
373 : ! Kab
374 21834 : kab = xtb_set_kab(za, zb, xtb_control)
375 128638 : DO j = 1, natorb_b
376 77096 : lb = laob(j)
377 77096 : nb = naob(j)
378 391966 : DO i = 1, natorb_a
379 293034 : la = laoa(i)
380 293034 : na = naoa(i)
381 293034 : kia = kl(la)
382 293034 : kjb = kl(lb)
383 293034 : IF (zb == 1 .AND. nb == 2) kjb = k2sh
384 293034 : IF (za == 1 .AND. na == 2) kia = k2sh
385 370130 : IF ((zb == 1 .AND. nb == 2) .OR. (za == 1 .AND. na == 2)) THEN
386 48590 : kijab(i, j, ikind, jkind) = 0.5_dp*(kia + kjb)
387 : ELSE
388 244444 : IF ((la == 0 .AND. lb == 1) .OR. (la == 1 .AND. lb == 0)) THEN
389 84084 : kijab(i, j, ikind, jkind) = ksp*kab*fen
390 : ELSE
391 160360 : 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 6404 : END SUBROUTINE gfn1_kpair
400 :
401 : END MODULE xtb_hcore
402 :
|
