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 xTB Hamiltonian derivative
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_hab_force
16 : USE ai_contraction, ONLY: block_add,&
17 : contraction
18 : USE ai_overlap, ONLY: overlap_ab
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind_set
21 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
22 : gto_basis_set_type
23 : USE block_p_types, ONLY: block_p_type
24 : USE cp_control_types, ONLY: dft_control_type,&
25 : xtb_control_type
26 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
27 : dbcsr_finalize,&
28 : dbcsr_get_block_p,&
29 : dbcsr_p_type,&
30 : dbcsr_type
31 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
32 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
33 : dbcsr_deallocate_matrix_set
34 : USE cp_log_handling, ONLY: cp_get_default_logger,&
35 : cp_logger_type
36 : USE kinds, ONLY: dp
37 : USE message_passing, ONLY: mp_para_env_type
38 : USE orbital_pointers, ONLY: ncoset
39 : USE particle_types, ONLY: particle_type
40 : USE qs_dispersion_cnum, ONLY: cnumber_init,&
41 : cnumber_release,&
42 : dcnum_type
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_force_types, ONLY: qs_force_type
46 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
47 : get_memory_usage
48 : USE qs_kind_types, ONLY: get_qs_kind,&
49 : qs_kind_type
50 : USE qs_ks_types, ONLY: qs_ks_env_type
51 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
52 : neighbor_list_iterate,&
53 : neighbor_list_iterator_create,&
54 : neighbor_list_iterator_p_type,&
55 : neighbor_list_iterator_release,&
56 : neighbor_list_set_p_type
57 : USE qs_overlap, ONLY: create_sab_matrix
58 : USE xtb_hcore, ONLY: gfn1_huckel,&
59 : gfn1_kpair
60 : USE xtb_types, ONLY: get_xtb_atom_param,&
61 : xtb_atom_type
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_hab_force'
69 :
70 : PUBLIC :: build_xtb_hab_force
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param qs_env ...
77 : !> \param p_matrix ...
78 : ! **************************************************************************************************
79 16 : SUBROUTINE build_xtb_hab_force(qs_env, p_matrix)
80 :
81 : TYPE(qs_environment_type), POINTER :: qs_env
82 : TYPE(dbcsr_type), POINTER :: p_matrix
83 :
84 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_xtb_hab_force'
85 :
86 : INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
87 : j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, maxder, n1, n2, na, natom, natorb_a, &
88 : natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, nsetb, sgfa, sgfb, &
89 : za, zb
90 32 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
91 : INTEGER, DIMENSION(25) :: laoa, laob, naoa, naob
92 : INTEGER, DIMENSION(3) :: cell
93 16 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
94 16 : npgfb, nsgfa, nsgfb
95 16 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
96 : LOGICAL :: defined, diagblock, found, use_virial
97 : REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, f0, fhua, fhub, &
98 : fhud, foab, hij, rcova, rcovab, rcovb, &
99 : rrab
100 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cnumbers
101 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dfblock, dhuckel, huckel, owork
102 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
103 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: kijab
104 : REAL(KIND=dp), DIMENSION(3) :: fdik, fdika, fdikb, force_ab, rij, rik
105 : REAL(KIND=dp), DIMENSION(5) :: dpia, dpib, kpolya, kpolyb, pia, pib
106 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
107 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: fblock, pblock, rpgfa, rpgfb, sblock, &
108 16 : scon_a, scon_b, zeta, zetb
109 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
110 64 : TYPE(block_p_type), DIMENSION(2:4) :: dsblocks
111 : TYPE(cp_logger_type), POINTER :: logger
112 16 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_s
113 16 : TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:) :: dcnum
114 : TYPE(dft_control_type), POINTER :: dft_control
115 16 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
116 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
117 : TYPE(mp_para_env_type), POINTER :: para_env
118 : TYPE(neighbor_list_iterator_p_type), &
119 16 : DIMENSION(:), POINTER :: nl_iterator
120 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
121 16 : POINTER :: sab_orb
122 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 16 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
124 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
125 : TYPE(qs_ks_env_type), POINTER :: ks_env
126 : TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
127 : TYPE(xtb_control_type), POINTER :: xtb_control
128 :
129 16 : CALL timeset(routineN, handle)
130 :
131 16 : NULLIFY (logger)
132 16 : logger => cp_get_default_logger()
133 :
134 16 : NULLIFY (matrix_h, matrix_s, atomic_kind_set, qs_kind_set, sab_orb)
135 :
136 : CALL get_qs_env(qs_env=qs_env, &
137 : atomic_kind_set=atomic_kind_set, &
138 : qs_kind_set=qs_kind_set, &
139 : dft_control=dft_control, &
140 : para_env=para_env, &
141 16 : sab_orb=sab_orb)
142 :
143 16 : CPASSERT(dft_control%qs_control%xtb_control%gfn_type == 1)
144 :
145 16 : nkind = SIZE(atomic_kind_set)
146 16 : xtb_control => dft_control%qs_control%xtb_control
147 16 : nimg = dft_control%nimages
148 16 : nderivatives = 1
149 16 : maxder = ncoset(nderivatives)
150 :
151 16 : NULLIFY (particle_set)
152 16 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
153 16 : natom = SIZE(particle_set)
154 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
155 16 : atom_of_kind=atom_of_kind, kind_of=kind_of)
156 :
157 16 : NULLIFY (force)
158 16 : CALL get_qs_env(qs_env=qs_env, force=force)
159 16 : use_virial = .FALSE.
160 16 : CPASSERT(nimg == 1)
161 :
162 : ! set up basis set lists
163 86 : ALLOCATE (basis_set_list(nkind))
164 16 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
165 :
166 : ! allocate overlap matrix
167 16 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
168 16 : CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
169 : CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
170 16 : sab_orb, .TRUE.)
171 : ! initialize H matrix
172 16 : CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
173 32 : DO img = 1, nimg
174 16 : ALLOCATE (matrix_h(1, img)%matrix)
175 : CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
176 16 : name="HAMILTONIAN MATRIX")
177 32 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
178 : END DO
179 :
180 : ! Calculate coordination numbers
181 : ! needed for effective atomic energy levels (Eq. 12)
182 : ! code taken from D3 dispersion energy
183 16 : CALL cnumber_init(qs_env, cnumbers, dcnum, 1, .TRUE.)
184 :
185 : ! Calculate Huckel parameters
186 16 : CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, .TRUE.)
187 :
188 : ! Calculate KAB parameters and electronegativity correction
189 16 : CALL gfn1_kpair(qs_env, kijab)
190 :
191 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
192 16 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
193 13003 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
194 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
195 12987 : iatom=iatom, jatom=jatom, r=rij, cell=cell)
196 12987 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
197 12987 : CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
198 12987 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
199 12987 : CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
200 12987 : CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
201 12987 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
202 :
203 51948 : dr = SQRT(SUM(rij(:)**2))
204 :
205 : ! atomic parameters
206 : CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, &
207 12987 : nshell=nsa, kpoly=kpolya)
208 : CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, &
209 12987 : nshell=nsb, kpoly=kpolyb)
210 :
211 12987 : ic = 1
212 12987 : icol = MAX(iatom, jatom)
213 12987 : irow = MIN(iatom, jatom)
214 12987 : NULLIFY (sblock, fblock)
215 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
216 12987 : row=irow, col=icol, BLOCK=sblock, found=found)
217 12987 : CPASSERT(found)
218 : CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
219 12987 : row=irow, col=icol, BLOCK=fblock, found=found)
220 12987 : CPASSERT(found)
221 :
222 12987 : NULLIFY (pblock)
223 : CALL dbcsr_get_block_p(matrix=p_matrix, &
224 12987 : row=irow, col=icol, block=pblock, found=found)
225 12987 : CPASSERT(ASSOCIATED(pblock))
226 51948 : DO i = 2, 4
227 38961 : NULLIFY (dsblocks(i)%block)
228 : CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
229 38961 : row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
230 51948 : CPASSERT(found)
231 : END DO
232 :
233 : ! overlap
234 12987 : basis_set_a => basis_set_list(ikind)%gto_basis_set
235 12987 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
236 12987 : basis_set_b => basis_set_list(jkind)%gto_basis_set
237 12987 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
238 12987 : atom_a = atom_of_kind(iatom)
239 12987 : atom_b = atom_of_kind(jatom)
240 : ! basis ikind
241 12987 : first_sgfa => basis_set_a%first_sgf
242 12987 : la_max => basis_set_a%lmax
243 12987 : la_min => basis_set_a%lmin
244 12987 : npgfa => basis_set_a%npgf
245 12987 : nseta = basis_set_a%nset
246 12987 : nsgfa => basis_set_a%nsgf_set
247 12987 : rpgfa => basis_set_a%pgf_radius
248 12987 : set_radius_a => basis_set_a%set_radius
249 12987 : scon_a => basis_set_a%scon
250 12987 : zeta => basis_set_a%zet
251 : ! basis jkind
252 12987 : first_sgfb => basis_set_b%first_sgf
253 12987 : lb_max => basis_set_b%lmax
254 12987 : lb_min => basis_set_b%lmin
255 12987 : npgfb => basis_set_b%npgf
256 12987 : nsetb = basis_set_b%nset
257 12987 : nsgfb => basis_set_b%nsgf_set
258 12987 : rpgfb => basis_set_b%pgf_radius
259 12987 : set_radius_b => basis_set_b%set_radius
260 12987 : scon_b => basis_set_b%scon
261 12987 : zetb => basis_set_b%zet
262 :
263 12987 : ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
264 103896 : ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
265 64935 : ALLOCATE (sint(natorb_a, natorb_b, maxder))
266 515591 : sint = 0.0_dp
267 :
268 38961 : DO iset = 1, nseta
269 25974 : ncoa = npgfa(iset)*ncoset(la_max(iset))
270 25974 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
271 25974 : sgfa = first_sgfa(1, iset)
272 90909 : DO jset = 1, nsetb
273 51948 : IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
274 45815 : ncob = npgfb(jset)*ncoset(lb_max(jset))
275 45815 : n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
276 45815 : sgfb = first_sgfb(1, jset)
277 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
278 : lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
279 45815 : rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
280 : ! Contraction
281 255049 : DO i = 1, 4
282 : CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
283 183260 : cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
284 235208 : CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
285 : END DO
286 : END DO
287 : END DO
288 : ! update S matrix
289 12987 : IF (iatom <= jatom) THEN
290 61008 : sblock(:, :) = sblock(:, :) + sint(:, :, 1)
291 : ELSE
292 59865 : sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
293 : END IF
294 51948 : DO i = 2, 4
295 51948 : IF (iatom <= jatom) THEN
296 183024 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
297 : ELSE
298 179595 : dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
299 : END IF
300 : END DO
301 :
302 : ! Calculate Pi = Pia * Pib (Eq. 11)
303 12987 : rcovab = rcova + rcovb
304 12987 : rrab = SQRT(dr/rcovab)
305 38961 : pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
306 38961 : pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
307 12987 : IF (dr > 1.e-6_dp) THEN
308 12867 : drx = 0.5_dp/rrab/rcovab
309 : ELSE
310 : drx = 0.0_dp
311 : END IF
312 38961 : dpia(1:nsa) = drx*kpolya(1:nsa)
313 38961 : dpib(1:nsb) = drx*kpolyb(1:nsb)
314 :
315 : ! diagonal block
316 12987 : diagblock = .FALSE.
317 12987 : IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
318 : !
319 : ! Eq. 10
320 : !
321 : IF (diagblock) THEN
322 444 : DO i = 1, natorb_a
323 324 : na = naoa(i)
324 444 : fblock(i, i) = fblock(i, i) + huckel(na, iatom)
325 : END DO
326 : ELSE
327 44819 : DO j = 1, natorb_b
328 31952 : nb = naob(j)
329 124223 : DO i = 1, natorb_a
330 79404 : na = naoa(i)
331 79404 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
332 111356 : IF (iatom <= jatom) THEN
333 39648 : fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
334 : ELSE
335 39756 : fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
336 : END IF
337 : END DO
338 : END DO
339 : END IF
340 :
341 12987 : f0 = 1.0_dp
342 12987 : IF (irow == iatom) f0 = -1.0_dp
343 : ! Derivative wrt coordination number
344 12987 : fhua = 0.0_dp
345 12987 : fhub = 0.0_dp
346 12987 : fhud = 0.0_dp
347 12987 : IF (diagblock) THEN
348 444 : DO i = 1, natorb_a
349 324 : la = laoa(i)
350 324 : na = naoa(i)
351 444 : fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
352 : END DO
353 : ELSE
354 44819 : DO j = 1, natorb_b
355 31952 : lb = laob(j)
356 31952 : nb = naob(j)
357 124223 : DO i = 1, natorb_a
358 79404 : la = laoa(i)
359 79404 : na = naoa(i)
360 79404 : hij = 0.5_dp*pia(na)*pib(nb)
361 111356 : IF (iatom <= jatom) THEN
362 39648 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
363 39648 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
364 : ELSE
365 39756 : fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
366 39756 : fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
367 : END IF
368 : END DO
369 : END DO
370 12867 : IF (iatom /= jatom) THEN
371 12825 : fhua = 2.0_dp*fhua
372 12825 : fhub = 2.0_dp*fhub
373 : END IF
374 : END IF
375 : ! iatom
376 12987 : atom_a = atom_of_kind(iatom)
377 184926 : DO i = 1, dcnum(iatom)%neighbors
378 171939 : katom = dcnum(iatom)%nlist(i)
379 171939 : kkind = kind_of(katom)
380 171939 : atom_c = atom_of_kind(katom)
381 687756 : rik = dcnum(iatom)%rik(:, i)
382 687756 : drk = SQRT(SUM(rik(:)**2))
383 184926 : IF (drk > 1.e-3_dp) THEN
384 687756 : fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
385 687756 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
386 687756 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
387 687756 : fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
388 687756 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
389 687756 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
390 : END IF
391 : END DO
392 : ! jatom
393 12987 : atom_b = atom_of_kind(jatom)
394 184342 : DO i = 1, dcnum(jatom)%neighbors
395 171355 : katom = dcnum(jatom)%nlist(i)
396 171355 : kkind = kind_of(katom)
397 171355 : atom_c = atom_of_kind(katom)
398 685420 : rik = dcnum(jatom)%rik(:, i)
399 685420 : drk = SQRT(SUM(rik(:)**2))
400 184342 : IF (drk > 1.e-3_dp) THEN
401 685420 : fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
402 685420 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
403 685420 : force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
404 : END IF
405 : END DO
406 12987 : IF (diagblock) THEN
407 120 : force_ab = 0._dp
408 : ELSE
409 : ! force from R dendent Huckel element
410 12867 : n1 = SIZE(fblock, 1)
411 12867 : n2 = SIZE(fblock, 2)
412 51468 : ALLOCATE (dfblock(n1, n2))
413 119445 : dfblock = 0.0_dp
414 44819 : DO j = 1, natorb_b
415 31952 : lb = laob(j)
416 31952 : nb = naob(j)
417 124223 : DO i = 1, natorb_a
418 79404 : la = laoa(i)
419 79404 : na = naoa(i)
420 79404 : dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
421 111356 : IF (iatom <= jatom) THEN
422 39648 : dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
423 : ELSE
424 39756 : dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
425 : END IF
426 : END DO
427 : END DO
428 119445 : dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
429 51468 : DO ir = 1, 3
430 38601 : foab = 2.0_dp*dfp*rij(ir)/dr
431 : ! force from overlap matrix contribution to H
432 134457 : DO j = 1, natorb_b
433 95856 : lb = laob(j)
434 95856 : nb = naob(j)
435 372669 : DO i = 1, natorb_a
436 238212 : la = laoa(i)
437 238212 : na = naoa(i)
438 238212 : hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
439 334068 : IF (iatom <= jatom) THEN
440 118944 : foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
441 : ELSE
442 119268 : foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
443 : END IF
444 : END DO
445 : END DO
446 51468 : force_ab(ir) = foab
447 : END DO
448 12867 : DEALLOCATE (dfblock)
449 : END IF
450 :
451 12987 : atom_a = atom_of_kind(iatom)
452 12987 : atom_b = atom_of_kind(jatom)
453 32565 : IF (irow == iatom) force_ab = -force_ab
454 51948 : force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
455 51948 : force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
456 :
457 77938 : DEALLOCATE (oint, owork, sint)
458 :
459 : END DO
460 16 : CALL neighbor_list_iterator_release(nl_iterator)
461 :
462 32 : DO i = 1, SIZE(matrix_h, 1)
463 48 : DO img = 1, nimg
464 16 : CALL dbcsr_finalize(matrix_h(i, img)%matrix)
465 32 : CALL dbcsr_finalize(matrix_s(i, img)%matrix)
466 : END DO
467 : END DO
468 16 : CALL dbcsr_deallocate_matrix_set(matrix_s)
469 16 : CALL dbcsr_deallocate_matrix_set(matrix_h)
470 :
471 : ! deallocate coordination numbers
472 16 : CALL cnumber_release(cnumbers, dcnum, .TRUE.)
473 :
474 : ! deallocate Huckel parameters
475 16 : DEALLOCATE (huckel, dhuckel)
476 : ! deallocate KAB parameters
477 16 : DEALLOCATE (kijab)
478 :
479 16 : DEALLOCATE (basis_set_list)
480 :
481 16 : CALL timestop(handle)
482 :
483 48 : END SUBROUTINE build_xtb_hab_force
484 :
485 : END MODULE xtb_hab_force
486 :
|
