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 Coulomb contributions in DFTB
10 : !> \author JGH
11 : ! **************************************************************************************************
12 : MODULE qs_dftb_coulomb
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set,&
15 : is_hydrogen
16 : USE atprop_types, ONLY: atprop_array_init,&
17 : atprop_type
18 : USE cell_types, ONLY: cell_type,&
19 : get_cell,&
20 : pbc
21 : USE cp_control_types, ONLY: dft_control_type,&
22 : dftb_control_type
23 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
24 : dbcsr_get_block_p,&
25 : dbcsr_iterator_blocks_left,&
26 : dbcsr_iterator_next_block,&
27 : dbcsr_iterator_start,&
28 : dbcsr_iterator_stop,&
29 : dbcsr_iterator_type,&
30 : dbcsr_p_type
31 : USE distribution_1d_types, ONLY: distribution_1d_type
32 : USE ewald_environment_types, ONLY: ewald_env_get,&
33 : ewald_environment_type
34 : USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
35 : tb_spme_evaluate
36 : USE ewald_pw_types, ONLY: ewald_pw_type
37 : USE kinds, ONLY: dp
38 : USE kpoint_types, ONLY: get_kpoint_info,&
39 : kpoint_type
40 : USE mathconstants, ONLY: oorootpi,&
41 : pi
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE particle_types, ONLY: particle_type
44 : USE pw_poisson_types, ONLY: do_ewald_ewald,&
45 : do_ewald_none,&
46 : do_ewald_pme,&
47 : do_ewald_spme
48 : USE qmmm_tb_coulomb, ONLY: build_tb_coulomb_qmqm
49 : USE qs_dftb3_methods, ONLY: build_dftb3_diagonal
50 : USE qs_dftb_types, ONLY: qs_dftb_atom_type,&
51 : qs_dftb_pairpot_type
52 : USE qs_dftb_utils, ONLY: compute_block_sk,&
53 : get_dftb_atom_param
54 : USE qs_energy_types, ONLY: qs_energy_type
55 : USE qs_environment_types, ONLY: get_qs_env,&
56 : qs_environment_type
57 : USE qs_force_types, ONLY: qs_force_type
58 : USE qs_kind_types, ONLY: get_qs_kind,&
59 : qs_kind_type
60 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
61 : neighbor_list_iterate,&
62 : neighbor_list_iterator_create,&
63 : neighbor_list_iterator_p_type,&
64 : neighbor_list_iterator_release,&
65 : neighbor_list_set_p_type
66 : USE qs_rho_types, ONLY: qs_rho_get,&
67 : qs_rho_type
68 : USE sap_kind_types, ONLY: clist_type,&
69 : release_sap_int,&
70 : sap_int_type
71 : USE virial_methods, ONLY: virial_pair_force
72 : USE virial_types, ONLY: virial_type
73 : #include "./base/base_uses.f90"
74 :
75 : IMPLICIT NONE
76 :
77 : PRIVATE
78 :
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_coulomb'
80 :
81 : ! screening for gamma function
82 : REAL(dp), PARAMETER :: tol_gamma = 1.e-4_dp
83 : ! small real number
84 : REAL(dp), PARAMETER :: rtiny = 1.e-10_dp
85 : REAL(KIND=dp), PARAMETER, PRIVATE :: dftb_fd_deriv_step = 1.0e-3_dp
86 :
87 : PUBLIC :: build_dftb_coulomb, gamma_rab_sr
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief ...
93 : !> \param qs_env ...
94 : !> \param ks_matrix ...
95 : !> \param rho ...
96 : !> \param mcharge ...
97 : !> \param energy ...
98 : !> \param calculate_forces ...
99 : !> \param just_energy ...
100 : ! **************************************************************************************************
101 28054 : SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
102 : calculate_forces, just_energy)
103 :
104 : TYPE(qs_environment_type), POINTER :: qs_env
105 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
106 : TYPE(qs_rho_type), POINTER :: rho
107 : REAL(dp), DIMENSION(:) :: mcharge
108 : TYPE(qs_energy_type), POINTER :: energy
109 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
110 :
111 : CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb_coulomb'
112 :
113 : INTEGER :: atom_i, atom_j, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, &
114 : irow, is, jatom, jkind, llm, lmaxi, lmaxj, n1, n2, natom, natorb_a, natorb_b, ngrd, &
115 : ngrdcut, nimg, nkind, nmat
116 28054 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
117 : INTEGER, DIMENSION(3) :: cellind, periodic
118 28054 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
119 : LOGICAL :: defined, defined_a, defined_b, do_ewald, &
120 : do_gamma_stress, found, hb_sr_damp, &
121 : use_virial
122 : REAL(KIND=dp) :: alpha, background, ddr, deth, dgam, &
123 : dgrd, dr, drm, drp, fi, ga, gb, gmat, &
124 : gmij, gmij_virial, hb_para, zeff
125 28054 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
126 : REAL(KIND=dp), DIMENSION(0:3) :: eta_a, eta_b
127 : REAL(KIND=dp), DIMENSION(3) :: drij, fij, rij
128 28054 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblock_vir, dsblockm_vir, &
129 28054 : gmcharge, ksblock, pblock, sblock, &
130 28054 : smatij, smatji
131 28054 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
132 28054 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
133 : TYPE(atprop_type), POINTER :: atprop
134 : TYPE(cell_type), POINTER :: cell
135 : TYPE(dbcsr_iterator_type) :: iter
136 28054 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
137 : TYPE(dft_control_type), POINTER :: dft_control
138 : TYPE(distribution_1d_type), POINTER :: local_particles
139 : TYPE(ewald_environment_type), POINTER :: ewald_env
140 : TYPE(ewald_pw_type), POINTER :: ewald_pw
141 : TYPE(kpoint_type), POINTER :: kpoints
142 : TYPE(mp_para_env_type), POINTER :: para_env
143 : TYPE(neighbor_list_iterator_p_type), &
144 28054 : DIMENSION(:), POINTER :: nl_iterator
145 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
146 28054 : POINTER :: n_list
147 28054 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
148 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind, dftb_kind_a, dftb_kind_b
149 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
150 28054 : POINTER :: dftb_potential
151 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
152 28054 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
153 28054 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 28054 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
155 : TYPE(virial_type), POINTER :: virial
156 :
157 28054 : CALL timeset(routineN, handle)
158 :
159 28054 : NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
160 :
161 28054 : use_virial = .FALSE.
162 :
163 28054 : IF (calculate_forces) THEN
164 : nmat = 4
165 : ELSE
166 27390 : nmat = 1
167 : END IF
168 :
169 28054 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
170 112216 : ALLOCATE (gmcharge(natom, nmat))
171 320004 : gmcharge = 0._dp
172 :
173 : CALL get_qs_env(qs_env, &
174 : particle_set=particle_set, &
175 : cell=cell, &
176 : virial=virial, &
177 : atprop=atprop, &
178 : dft_control=dft_control, &
179 : atomic_kind_set=atomic_kind_set, &
180 : qs_kind_set=qs_kind_set, &
181 28054 : force=force, para_env=para_env)
182 :
183 28054 : IF (calculate_forces) THEN
184 1198 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
185 : END IF
186 :
187 28054 : NULLIFY (dftb_potential)
188 28054 : CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
189 28054 : NULLIFY (n_list)
190 28054 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
191 28054 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
192 2985263 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
193 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
194 2957209 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
195 :
196 2957209 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
197 : CALL get_dftb_atom_param(dftb_kind_a, &
198 2957209 : defined=defined, eta=eta_a, natorb=natorb_a)
199 2957209 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
200 2957209 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
201 : CALL get_dftb_atom_param(dftb_kind_b, &
202 2957209 : defined=defined, eta=eta_b, natorb=natorb_b)
203 2957209 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
204 :
205 : ! gamma matrix
206 2957209 : hb_sr_damp = dft_control%qs_control%dftb_control%hb_sr_damp
207 2957209 : IF (hb_sr_damp) THEN
208 : ! short range correction enabled only when iatom XOR jatom are hydrogens
209 : hb_sr_damp = is_hydrogen(particle_set(iatom)%atomic_kind) .NEQV. &
210 1183097 : is_hydrogen(particle_set(jatom)%atomic_kind)
211 : END IF
212 1183097 : IF (hb_sr_damp) THEN
213 464808 : hb_para = dft_control%qs_control%dftb_control%hb_sr_para
214 : ELSE
215 2492401 : hb_para = 0.0_dp
216 : END IF
217 2957209 : ga = eta_a(0)
218 2957209 : gb = eta_b(0)
219 11828836 : dr = SQRT(SUM(rij(:)**2))
220 2957209 : gmat = gamma_rab_sr(dr, ga, gb, hb_para)
221 2957209 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + gmat*mcharge(iatom)
222 2957209 : IF (iatom /= jatom) THEN
223 2622548 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + gmat*mcharge(jatom)
224 : END IF
225 2985263 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
226 255211 : ddr = dftb_fd_deriv_step*dftb_potential(ikind, jkind)%dgrd
227 255211 : drp = dr + ddr
228 255211 : drm = dr - ddr
229 255211 : dgam = 0.5_dp*(gamma_rab_sr(drp, ga, gb, hb_para) - gamma_rab_sr(drm, ga, gb, hb_para))/ddr
230 1020844 : DO i = 1, 3
231 765633 : gmcharge(jatom, i + 1) = gmcharge(jatom, i + 1) - dgam*mcharge(iatom)*rij(i)/dr
232 765633 : IF (dr > 0.001_dp) THEN
233 765633 : gmcharge(iatom, i + 1) = gmcharge(iatom, i + 1) + dgam*mcharge(jatom)*rij(i)/dr
234 : END IF
235 1020844 : IF (use_virial) THEN
236 513246 : fij(i) = -mcharge(iatom)*mcharge(jatom)*dgam*rij(i)/dr
237 : END IF
238 : END DO
239 255211 : IF (use_virial) THEN
240 171082 : fi = 1.0_dp
241 171082 : IF (iatom == jatom) fi = 0.5_dp
242 171082 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
243 : END IF
244 : END IF
245 : END DO
246 28054 : CALL neighbor_list_iterator_release(nl_iterator)
247 :
248 28054 : IF (atprop%energy) THEN
249 466 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
250 466 : natom = SIZE(particle_set)
251 466 : CALL atprop_array_init(atprop%atecoul, natom)
252 : END IF
253 :
254 : ! 1/R contribution
255 28054 : do_ewald = dft_control%qs_control%dftb_control%do_ewald
256 28054 : IF (do_ewald) THEN
257 : ! Ewald sum
258 11276 : NULLIFY (ewald_env, ewald_pw)
259 : CALL get_qs_env(qs_env=qs_env, &
260 11276 : ewald_env=ewald_env, ewald_pw=ewald_pw)
261 11276 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
262 11276 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
263 11276 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
264 11276 : CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial)
265 0 : SELECT CASE (ewald_type)
266 : CASE DEFAULT
267 0 : CPABORT("Invalid Ewald type")
268 : CASE (do_ewald_none)
269 0 : CPABORT("Not allowed with DFTB")
270 : CASE (do_ewald_ewald)
271 0 : CPABORT("Standard Ewald not implemented in DFTB")
272 : CASE (do_ewald_pme)
273 0 : CPABORT("PME not implemented in DFTB")
274 : CASE (do_ewald_spme)
275 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
276 11276 : gmcharge, mcharge, calculate_forces, virial, use_virial)
277 : END SELECT
278 : ELSE
279 : ! direct sum
280 : CALL get_qs_env(qs_env=qs_env, &
281 16778 : local_particles=local_particles)
282 54034 : DO ikind = 1, SIZE(local_particles%n_el)
283 81046 : DO ia = 1, local_particles%n_el(ikind)
284 27012 : iatom = local_particles%list(ikind)%array(ia)
285 94978 : DO jatom = 1, iatom - 1
286 122840 : rij = particle_set(iatom)%r - particle_set(jatom)%r
287 122840 : rij = pbc(rij, cell)
288 122840 : dr = SQRT(SUM(rij(:)**2))
289 30710 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
290 30710 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
291 31784 : DO i = 2, nmat
292 1074 : gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
293 31784 : gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
294 : END DO
295 57722 : IF (use_virial .AND. dr > 1.e-6_dp) THEN
296 72 : DO i = 1, 3
297 72 : fij(i) = mcharge(iatom)*mcharge(jatom)*rij(i)/dr**3
298 : END DO
299 18 : CALL virial_pair_force(virial%pv_virial, 1.0_dp, fij, rij)
300 : END IF
301 : END DO
302 : END DO
303 : END DO
304 : END IF
305 :
306 465498 : CALL para_env%sum(gmcharge(:, 1))
307 :
308 28054 : background = 0.0_dp
309 28054 : IF (do_ewald) THEN
310 : ! add self charge interaction and background charge contribution
311 175974 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
312 12098 : IF (ANY(periodic(:) == 1)) THEN
313 11002 : background = pi/alpha**2/deth
314 : END IF
315 : END IF
316 :
317 246776 : energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*(gmcharge(:, 1) - background))
318 28054 : IF (atprop%energy) THEN
319 : CALL get_qs_env(qs_env=qs_env, &
320 466 : local_particles=local_particles)
321 1526 : DO ikind = 1, SIZE(local_particles%n_el)
322 1060 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
323 1060 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
324 18006 : DO ia = 1, local_particles%n_el(ikind)
325 16480 : iatom = local_particles%list(ikind)%array(ia)
326 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
327 17540 : 0.5_dp*zeff*gmcharge(iatom, 1)
328 : END DO
329 : END DO
330 : END IF
331 :
332 28054 : IF (calculate_forces) THEN
333 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
334 : kind_of=kind_of, &
335 664 : atom_of_kind=atom_of_kind)
336 :
337 15058 : gmcharge(:, 2) = gmcharge(:, 2)*mcharge(:)
338 15058 : gmcharge(:, 3) = gmcharge(:, 3)*mcharge(:)
339 15058 : gmcharge(:, 4) = gmcharge(:, 4)*mcharge(:)
340 15058 : DO iatom = 1, natom
341 14394 : ikind = kind_of(iatom)
342 14394 : atom_i = atom_of_kind(iatom)
343 14394 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - gmcharge(iatom, 2)
344 14394 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - gmcharge(iatom, 3)
345 15058 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - gmcharge(iatom, 4)
346 : END DO
347 : END IF
348 :
349 28054 : do_gamma_stress = .FALSE.
350 28054 : IF (.NOT. just_energy .AND. use_virial) THEN
351 130 : IF (dft_control%nimages == 1) do_gamma_stress = .TRUE.
352 : END IF
353 :
354 28054 : IF (.NOT. just_energy) THEN
355 27782 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
356 27782 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
357 :
358 27782 : nimg = dft_control%nimages
359 27782 : NULLIFY (cell_to_index)
360 27782 : IF (nimg > 1) THEN
361 12120 : NULLIFY (kpoints)
362 12120 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
363 12120 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
364 : END IF
365 :
366 27782 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
367 336 : DO img = 1, nimg
368 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
369 336 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
370 : END DO
371 : END IF
372 :
373 27782 : NULLIFY (sap_int)
374 27782 : IF (do_gamma_stress) THEN
375 : ! derivative overlap integral (non collapsed)
376 94 : CALL dftb_dsint_list(qs_env, sap_int)
377 : END IF
378 :
379 27782 : IF (nimg == 1) THEN
380 : ! no k-points; all matrices have been transformed to periodic bsf
381 15662 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
382 1764163 : DO WHILE (dbcsr_iterator_blocks_left(iter))
383 1748501 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
384 1748501 : gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
385 3505699 : DO is = 1, SIZE(ks_matrix, 1)
386 1757198 : NULLIFY (ksblock)
387 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
388 1757198 : row=irow, col=icol, block=ksblock, found=found)
389 1757198 : CPASSERT(found)
390 26544713 : ksblock = ksblock - gmij*sblock
391 : END DO
392 1764163 : IF (calculate_forces) THEN
393 230549 : ikind = kind_of(irow)
394 230549 : atom_i = atom_of_kind(irow)
395 230549 : jkind = kind_of(icol)
396 230549 : atom_j = atom_of_kind(icol)
397 230549 : NULLIFY (pblock)
398 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
399 230549 : row=irow, col=icol, block=pblock, found=found)
400 230549 : CPASSERT(found)
401 922196 : DO i = 1, 3
402 691647 : NULLIFY (dsblock)
403 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
404 691647 : row=irow, col=icol, block=dsblock, found=found)
405 691647 : CPASSERT(found)
406 4855806 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
407 691647 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
408 691647 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
409 1613843 : fij(i) = fi
410 : END DO
411 : END IF
412 : END DO
413 15662 : CALL dbcsr_iterator_stop(iter)
414 : !
415 15662 : IF (do_gamma_stress) THEN
416 282 : DO ikind = 1, nkind
417 658 : DO jkind = 1, nkind
418 376 : iac = ikind + nkind*(jkind - 1)
419 376 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
420 8638 : DO ia = 1, sap_int(iac)%nalist
421 8094 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
422 8089 : iatom = sap_int(iac)%alist(ia)%aatom
423 177713 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
424 169248 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
425 676992 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
426 676992 : dr = SQRT(SUM(rij(:)**2))
427 177342 : IF (dr > 1.e-6_dp) THEN
428 165201 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
429 165201 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
430 165201 : icol = MAX(iatom, jatom)
431 165201 : irow = MIN(iatom, jatom)
432 165201 : NULLIFY (pblock)
433 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
434 165201 : row=irow, col=icol, block=pblock, found=found)
435 165201 : CPASSERT(found)
436 660804 : DO i = 1, 3
437 660804 : IF (irow == iatom) THEN
438 1722762 : fij(i) = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i))
439 : ELSE
440 1731483 : fij(i) = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
441 : END IF
442 : END DO
443 165201 : fi = 1.0_dp
444 165201 : IF (iatom == jatom) fi = 0.5_dp
445 165201 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
446 : END IF
447 : END DO
448 : END DO
449 : END DO
450 : END DO
451 : END IF
452 : ELSE
453 12120 : NULLIFY (n_list)
454 12120 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
455 12120 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
456 721307 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
457 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
458 709187 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
459 :
460 709187 : icol = MAX(iatom, jatom)
461 709187 : irow = MIN(iatom, jatom)
462 709187 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
463 709187 : CPASSERT(ic > 0)
464 :
465 709187 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
466 :
467 709187 : NULLIFY (sblock)
468 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
469 709187 : row=irow, col=icol, block=sblock, found=found)
470 709187 : CPASSERT(found)
471 1418374 : DO is = 1, SIZE(ks_matrix, 1)
472 709187 : NULLIFY (ksblock)
473 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
474 709187 : row=irow, col=icol, block=ksblock, found=found)
475 709187 : CPASSERT(found)
476 46814543 : ksblock = ksblock - gmij*sblock
477 : END DO
478 :
479 721307 : IF (calculate_forces) THEN
480 6218 : ikind = kind_of(iatom)
481 6218 : atom_i = atom_of_kind(iatom)
482 6218 : jkind = kind_of(jatom)
483 6218 : atom_j = atom_of_kind(jatom)
484 6218 : IF (irow == jatom) gmij = -gmij
485 6218 : NULLIFY (pblock)
486 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
487 6218 : row=irow, col=icol, block=pblock, found=found)
488 6218 : CPASSERT(found)
489 24872 : DO i = 1, 3
490 18654 : NULLIFY (dsblock)
491 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
492 18654 : row=irow, col=icol, block=dsblock, found=found)
493 18654 : CPASSERT(found)
494 393381 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
495 18654 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
496 18654 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
497 43526 : fij(i) = fi
498 : END DO
499 6218 : IF (use_virial) THEN
500 : ! The image-resolved overlap derivative blocks are matrix oriented. Rebuild the
501 : ! atom-oriented derivative for the virial contraction, matching the Gamma path.
502 5996 : fij = 0.0_dp
503 5996 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
504 : CALL get_dftb_atom_param(dftb_kind_a, &
505 5996 : defined=defined_a, lmax=lmaxi, natorb=natorb_a)
506 5996 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
507 : CALL get_dftb_atom_param(dftb_kind_b, &
508 5996 : defined=defined_b, lmax=lmaxj, natorb=natorb_b)
509 5996 : IF (defined_a .AND. defined_b .AND. natorb_a > 0 .AND. natorb_b > 0) THEN
510 5996 : dftb_param_ij => dftb_potential(ikind, jkind)
511 5996 : dftb_param_ji => dftb_potential(jkind, ikind)
512 5996 : ngrd = dftb_param_ij%ngrd
513 5996 : ngrdcut = dftb_param_ij%ngrdcut
514 5996 : dgrd = dftb_param_ij%dgrd
515 5996 : ddr = dgrd*dftb_fd_deriv_step
516 5996 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
517 5996 : llm = dftb_param_ij%llm
518 5996 : smatij => dftb_param_ij%smat
519 5996 : smatji => dftb_param_ji%smat
520 23984 : dr = SQRT(SUM(rij(:)**2))
521 5996 : IF (NINT(dr/dgrd) <= ngrdcut .AND. dr > 1.e-6_dp) THEN
522 5870 : n1 = natorb_a
523 5870 : n2 = natorb_b
524 35220 : ALLOCATE (dsblock_vir(n1, n2), dsblockm_vir(n1, n2))
525 5870 : gmij_virial = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
526 23480 : DO i = 1, 3
527 381909 : dsblock_vir = 0._dp
528 381909 : dsblockm_vir = 0._dp
529 17610 : drij = rij
530 17610 : drij(i) = rij(i) - ddr
531 : CALL compute_block_sk(dsblockm_vir, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
532 17610 : llm, lmaxi, lmaxj, iatom, iatom)
533 17610 : drij(i) = rij(i) + ddr
534 : CALL compute_block_sk(dsblock_vir, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
535 17610 : llm, lmaxi, lmaxj, iatom, iatom)
536 746208 : dsblock_vir = (dsblock_vir - dsblockm_vir)/(2.0_dp*ddr)
537 23480 : IF (irow == iatom) THEN
538 209892 : fij(i) = 2.0_dp*gmij_virial*SUM(pblock*dsblock_vir)
539 : ELSE
540 172017 : fij(i) = 2.0_dp*gmij_virial*SUM(TRANSPOSE(pblock)*dsblock_vir)
541 : END IF
542 : END DO
543 5870 : DEALLOCATE (dsblock_vir, dsblockm_vir)
544 5870 : fi = 1.0_dp
545 5870 : IF (iatom == jatom) fi = 0.5_dp
546 5870 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
547 : END IF
548 : END IF
549 : END IF
550 : END IF
551 : END DO
552 12120 : CALL neighbor_list_iterator_release(nl_iterator)
553 : END IF
554 :
555 27782 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
556 336 : DO img = 1, nimg
557 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
558 336 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
559 : END DO
560 : END IF
561 : END IF
562 :
563 28054 : IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN
564 10932 : CALL get_qs_env(qs_env, nkind=nkind)
565 43728 : ALLOCATE (zeffk(nkind), xgamma(nkind))
566 32690 : DO ikind = 1, nkind
567 21758 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
568 32690 : CALL get_dftb_atom_param(dftb_kind, dudq=xgamma(ikind), zeff=zeffk(ikind))
569 : END DO
570 : ! Diagonal 3rd order correction (DFTB3)
571 : CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
572 10932 : sap_int, calculate_forces, just_energy)
573 10932 : DEALLOCATE (zeffk, xgamma)
574 : END IF
575 :
576 : ! QMMM
577 28054 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
578 : CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
579 1626 : calculate_forces, just_energy)
580 : END IF
581 :
582 28054 : IF (do_gamma_stress) THEN
583 94 : CALL release_sap_int(sap_int)
584 : END IF
585 :
586 28054 : DEALLOCATE (gmcharge)
587 :
588 28054 : CALL timestop(handle)
589 :
590 56108 : END SUBROUTINE build_dftb_coulomb
591 :
592 : ! **************************************************************************************************
593 : !> \brief Computes the short-range gamma parameter from exact Coulomb
594 : !> interaction of normalized exp(-a*r) charge distribution - 1/r
595 : !> \param r ...
596 : !> \param ga ...
597 : !> \param gb ...
598 : !> \param hb_para ...
599 : !> \return ...
600 : !> \par Literature
601 : !> Elstner et al, PRB 58 (1998) 7260
602 : !> \par History
603 : !> 10.2008 Axel Kohlmeyer - adding sr_damp
604 : !> 08.2014 JGH - adding flexible exponent for damping
605 : !> \version 1.1
606 : ! **************************************************************************************************
607 3497727 : FUNCTION gamma_rab_sr(r, ga, gb, hb_para) RESULT(gamma)
608 : REAL(dp), INTENT(in) :: r, ga, gb, hb_para
609 : REAL(dp) :: gamma
610 :
611 : REAL(dp) :: a, b, fac, g_sum
612 :
613 3497727 : gamma = 0.0_dp
614 3497727 : a = 3.2_dp*ga ! 3.2 = 16/5 in Eq. 18 and ff.
615 3497727 : b = 3.2_dp*gb
616 3497727 : g_sum = a + b
617 3497727 : IF (g_sum < tol_gamma) RETURN ! hardness screening
618 3497727 : IF (r < rtiny) THEN ! This is for short distances but non-onsite terms
619 : ! This gives also correct diagonal elements (a=b, r=0)
620 109361 : gamma = 0.5_dp*(a*b/g_sum + (a*b)**2/g_sum**3)
621 109361 : RETURN
622 : END IF
623 : !
624 : ! distinguish two cases: Gamma's are very close, e.g. for the same atom type,
625 : ! and Gamma's are different
626 : !
627 3388366 : IF (ABS(a - b) < rtiny) THEN
628 2018628 : fac = 1.6_dp*r*a*b/g_sum*(1.0_dp + a*b/g_sum**2)
629 2018628 : gamma = -(48.0_dp + 33._dp*fac + (9.0_dp + fac)*fac**2)*EXP(-fac)/(48._dp*r)
630 : ELSE
631 : gamma = -EXP(-a*r)*(0.5_dp*a*b**4/(a**2 - b**2)**2 - &
632 : (b**6 - 3._dp*a**2*b**4)/(r*(a**2 - b**2)**3)) - & ! a-> b
633 : EXP(-b*r)*(0.5_dp*b*a**4/(b**2 - a**2)**2 - &
634 1369738 : (a**6 - 3._dp*b**2*a**4)/(r*(b**2 - a**2)**3)) ! b-> a
635 : END IF
636 : !
637 : ! damping function for better short range hydrogen bonds.
638 : ! functional form from Hu H. et al., J. Phys. Chem. A 2007, 111, 5685-5691
639 : ! according to Elstner M, Theor. Chem. Acc. 2006, 116, 316-325,
640 : ! this should only be applied to a-b pairs involving hydrogen.
641 3388366 : IF (hb_para > 0.0_dp) THEN
642 522820 : gamma = gamma*EXP(-(0.5_dp*(ga + gb))**hb_para*r*r)
643 : END IF
644 : END FUNCTION gamma_rab_sr
645 :
646 : ! **************************************************************************************************
647 : !> \brief ...
648 : !> \param qs_env ...
649 : !> \param sap_int ...
650 : ! **************************************************************************************************
651 94 : SUBROUTINE dftb_dsint_list(qs_env, sap_int)
652 :
653 : TYPE(qs_environment_type), POINTER :: qs_env
654 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
655 :
656 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dftb_dsint_list'
657 :
658 : INTEGER :: handle, i, iac, iatom, ikind, ilist, jatom, jkind, jneighbor, llm, lmaxi, lmaxj, &
659 : n1, n2, natorb_a, natorb_b, ngrd, ngrdcut, nkind, nlist, nneighbor
660 : INTEGER, DIMENSION(3) :: cell
661 : LOGICAL :: defined
662 : REAL(KIND=dp) :: ddr, dgrd, dr
663 : REAL(KIND=dp), DIMENSION(3) :: drij, rij
664 94 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, smatij, smatji
665 : TYPE(clist_type), POINTER :: clist
666 : TYPE(dft_control_type), POINTER :: dft_control
667 : TYPE(dftb_control_type), POINTER :: dftb_control
668 : TYPE(neighbor_list_iterator_p_type), &
669 94 : DIMENSION(:), POINTER :: nl_iterator
670 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
671 94 : POINTER :: sab_orb
672 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b
673 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
674 94 : POINTER :: dftb_potential
675 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
676 94 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
677 :
678 94 : CALL timeset(routineN, handle)
679 :
680 94 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
681 94 : CPASSERT(.NOT. ASSOCIATED(sap_int))
682 658 : ALLOCATE (sap_int(nkind*nkind))
683 470 : DO i = 1, nkind*nkind
684 376 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
685 470 : sap_int(i)%nalist = 0
686 : END DO
687 :
688 : CALL get_qs_env(qs_env=qs_env, &
689 : qs_kind_set=qs_kind_set, &
690 : dft_control=dft_control, &
691 94 : sab_orb=sab_orb)
692 :
693 94 : dftb_control => dft_control%qs_control%dftb_control
694 :
695 94 : NULLIFY (dftb_potential)
696 : CALL get_qs_env(qs_env=qs_env, &
697 94 : dftb_potential=dftb_potential)
698 :
699 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
700 94 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
701 169342 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
702 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, &
703 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
704 169248 : inode=jneighbor, cell=cell, r=rij)
705 169248 : iac = ikind + nkind*(jkind - 1)
706 : !
707 169248 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
708 : CALL get_dftb_atom_param(dftb_kind_a, &
709 169248 : defined=defined, lmax=lmaxi, natorb=natorb_a)
710 169248 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
711 169248 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
712 : CALL get_dftb_atom_param(dftb_kind_b, &
713 169248 : defined=defined, lmax=lmaxj, natorb=natorb_b)
714 169248 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
715 :
716 676992 : dr = SQRT(SUM(rij(:)**2))
717 :
718 : ! integral list
719 169248 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
720 356 : sap_int(iac)%a_kind = ikind
721 356 : sap_int(iac)%p_kind = jkind
722 356 : sap_int(iac)%nalist = nlist
723 9162 : ALLOCATE (sap_int(iac)%alist(nlist))
724 8450 : DO i = 1, nlist
725 8094 : NULLIFY (sap_int(iac)%alist(i)%clist)
726 8094 : sap_int(iac)%alist(i)%aatom = 0
727 8450 : sap_int(iac)%alist(i)%nclist = 0
728 : END DO
729 : END IF
730 169248 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
731 8089 : sap_int(iac)%alist(ilist)%aatom = iatom
732 8089 : sap_int(iac)%alist(ilist)%nclist = nneighbor
733 242049 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
734 177337 : DO i = 1, nneighbor
735 177337 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
736 : END DO
737 : END IF
738 169248 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
739 169248 : clist%catom = jatom
740 676992 : clist%cell = cell
741 676992 : clist%rac = rij
742 846240 : ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
743 169248 : NULLIFY (clist%achint)
744 3732762 : clist%acint = 0._dp
745 169248 : clist%nsgf_cnt = 0
746 169248 : NULLIFY (clist%sgf_list)
747 :
748 : ! retrieve information on S matrix
749 169248 : dftb_param_ij => dftb_potential(ikind, jkind)
750 169248 : dftb_param_ji => dftb_potential(jkind, ikind)
751 : ! assume table size and type is symmetric
752 169248 : ngrd = dftb_param_ij%ngrd
753 169248 : ngrdcut = dftb_param_ij%ngrdcut
754 169248 : dgrd = dftb_param_ij%dgrd
755 169248 : ddr = dgrd*dftb_fd_deriv_step
756 169248 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
757 169248 : llm = dftb_param_ij%llm
758 169248 : smatij => dftb_param_ij%smat
759 169248 : smatji => dftb_param_ji%smat
760 :
761 507838 : IF (NINT(dr/dgrd) <= ngrdcut) THEN
762 168874 : IF (iatom == jatom .AND. dr < 0.001_dp) THEN
763 : ! diagonal block
764 : ELSE
765 : ! off-diagonal block
766 164827 : n1 = natorb_a
767 164827 : n2 = natorb_b
768 988962 : ALLOCATE (dsblock(n1, n2), dsblockm(n1, n2))
769 659308 : DO i = 1, 3
770 3445605 : dsblock = 0._dp
771 3445605 : dsblockm = 0.0_dp
772 494481 : drij = rij
773 :
774 494481 : drij(i) = rij(i) - ddr
775 : CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
776 494481 : llm, lmaxi, lmaxj, iatom, iatom)
777 :
778 494481 : drij(i) = rij(i) + ddr
779 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
780 494481 : llm, lmaxi, lmaxj, iatom, iatom)
781 :
782 6396729 : dsblock = dsblock - dsblockm
783 3445605 : dsblock = dsblock/(2.0_dp*ddr)
784 :
785 6561556 : clist%acint(1:n1, 1:n2, i) = -dsblock(1:n1, 1:n2)
786 : END DO
787 164827 : DEALLOCATE (dsblock, dsblockm)
788 : END IF
789 : END IF
790 :
791 : END DO
792 94 : CALL neighbor_list_iterator_release(nl_iterator)
793 :
794 94 : CALL timestop(handle)
795 :
796 94 : END SUBROUTINE dftb_dsint_list
797 :
798 : END MODULE qs_dftb_coulomb
|
