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 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 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 :
86 : PUBLIC :: build_dftb_coulomb, gamma_rab_sr
87 :
88 : CONTAINS
89 :
90 : ! **************************************************************************************************
91 : !> \brief ...
92 : !> \param qs_env ...
93 : !> \param ks_matrix ...
94 : !> \param rho ...
95 : !> \param mcharge ...
96 : !> \param energy ...
97 : !> \param calculate_forces ...
98 : !> \param just_energy ...
99 : ! **************************************************************************************************
100 11918 : SUBROUTINE build_dftb_coulomb(qs_env, ks_matrix, rho, mcharge, energy, &
101 : calculate_forces, just_energy)
102 :
103 : TYPE(qs_environment_type), POINTER :: qs_env
104 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
105 : TYPE(qs_rho_type), POINTER :: rho
106 : REAL(dp), DIMENSION(:) :: mcharge
107 : TYPE(qs_energy_type), POINTER :: energy
108 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
109 :
110 : CHARACTER(len=*), PARAMETER :: routineN = 'build_dftb_coulomb'
111 :
112 : INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, &
113 : irow, is, jatom, jkind, natom, natorb_a, natorb_b, nimg, nkind, nmat
114 11918 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
115 : INTEGER, DIMENSION(3) :: cellind, periodic
116 11918 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
117 : LOGICAL :: defined, do_ewald, do_gamma_stress, &
118 : found, hb_sr_damp, use_virial
119 : REAL(KIND=dp) :: alpha, ddr, deth, dgam, dr, drm, drp, &
120 : fi, ga, gb, gmat, gmij, hb_para, zeff
121 11918 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
122 : REAL(KIND=dp), DIMENSION(0:3) :: eta_a, eta_b
123 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
124 11918 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, gmcharge, ksblock, pblock, &
125 11918 : sblock
126 11918 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
127 11918 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
128 : TYPE(atprop_type), POINTER :: atprop
129 : TYPE(cell_type), POINTER :: cell
130 : TYPE(dbcsr_iterator_type) :: iter
131 11918 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
132 : TYPE(dft_control_type), POINTER :: dft_control
133 : TYPE(distribution_1d_type), POINTER :: local_particles
134 : TYPE(ewald_environment_type), POINTER :: ewald_env
135 : TYPE(ewald_pw_type), POINTER :: ewald_pw
136 : TYPE(kpoint_type), POINTER :: kpoints
137 : TYPE(mp_para_env_type), POINTER :: para_env
138 : TYPE(neighbor_list_iterator_p_type), &
139 11918 : DIMENSION(:), POINTER :: nl_iterator
140 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
141 11918 : POINTER :: n_list
142 11918 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
143 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind, dftb_kind_a, dftb_kind_b
144 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
145 11918 : POINTER :: dftb_potential
146 11918 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
147 11918 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
148 11918 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
149 : TYPE(virial_type), POINTER :: virial
150 :
151 11918 : CALL timeset(routineN, handle)
152 :
153 11918 : NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
154 :
155 11918 : use_virial = .FALSE.
156 :
157 11918 : IF (calculate_forces) THEN
158 : nmat = 4
159 : ELSE
160 11302 : nmat = 1
161 : END IF
162 :
163 11918 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
164 47672 : ALLOCATE (gmcharge(natom, nmat))
165 223608 : gmcharge = 0._dp
166 :
167 : CALL get_qs_env(qs_env, &
168 : particle_set=particle_set, &
169 : cell=cell, &
170 : virial=virial, &
171 : atprop=atprop, &
172 : dft_control=dft_control, &
173 : atomic_kind_set=atomic_kind_set, &
174 : qs_kind_set=qs_kind_set, &
175 11918 : force=force, para_env=para_env)
176 :
177 11918 : IF (calculate_forces) THEN
178 1126 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
179 : END IF
180 :
181 11918 : NULLIFY (dftb_potential)
182 11918 : CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
183 11918 : NULLIFY (n_list)
184 11918 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
185 11918 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
186 2275516 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
187 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
188 2263598 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
189 :
190 2263598 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
191 : CALL get_dftb_atom_param(dftb_kind_a, &
192 2263598 : defined=defined, eta=eta_a, natorb=natorb_a)
193 2263598 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
194 2263598 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
195 : CALL get_dftb_atom_param(dftb_kind_b, &
196 2263598 : defined=defined, eta=eta_b, natorb=natorb_b)
197 2263598 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
198 :
199 : ! gamma matrix
200 2263598 : hb_sr_damp = dft_control%qs_control%dftb_control%hb_sr_damp
201 2263598 : IF (hb_sr_damp) THEN
202 : ! short range correction enabled only when iatom XOR jatom are hydrogens
203 : hb_sr_damp = is_hydrogen(particle_set(iatom)%atomic_kind) .NEQV. &
204 949716 : is_hydrogen(particle_set(jatom)%atomic_kind)
205 : END IF
206 2263598 : IF (hb_sr_damp) THEN
207 380582 : hb_para = dft_control%qs_control%dftb_control%hb_sr_para
208 : ELSE
209 1883016 : hb_para = 0.0_dp
210 : END IF
211 2263598 : ga = eta_a(0)
212 2263598 : gb = eta_b(0)
213 9054392 : dr = SQRT(SUM(rij(:)**2))
214 2263598 : gmat = gamma_rab_sr(dr, ga, gb, hb_para)
215 2263598 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + gmat*mcharge(iatom)
216 2263598 : IF (iatom /= jatom) THEN
217 2123947 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + gmat*mcharge(jatom)
218 : END IF
219 2275516 : IF (calculate_forces .AND. (iatom /= jatom .OR. dr > 0.001_dp)) THEN
220 253110 : ddr = 0.1_dp*dftb_potential(ikind, jkind)%dgrd
221 253110 : drp = dr + ddr
222 253110 : drm = dr - ddr
223 253110 : dgam = 0.5_dp*(gamma_rab_sr(drp, ga, gb, hb_para) - gamma_rab_sr(drm, ga, gb, hb_para))/ddr
224 1012440 : DO i = 1, 3
225 759330 : gmcharge(jatom, i + 1) = gmcharge(jatom, i + 1) - dgam*mcharge(iatom)*rij(i)/dr
226 759330 : IF (dr > 0.001_dp) THEN
227 759330 : gmcharge(iatom, i + 1) = gmcharge(iatom, i + 1) + dgam*mcharge(jatom)*rij(i)/dr
228 : END IF
229 1012440 : IF (use_virial) THEN
230 507639 : fij(i) = -mcharge(iatom)*mcharge(jatom)*dgam*rij(i)/dr
231 : END IF
232 : END DO
233 253110 : IF (use_virial) THEN
234 169213 : fi = 1.0_dp
235 169213 : IF (iatom == jatom) fi = 0.5_dp
236 169213 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
237 169213 : IF (atprop%stress) THEN
238 70678 : CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
239 70678 : CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
240 : END IF
241 : END IF
242 : END IF
243 : END DO
244 11918 : CALL neighbor_list_iterator_release(nl_iterator)
245 :
246 11918 : IF (atprop%energy) THEN
247 466 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
248 466 : natom = SIZE(particle_set)
249 466 : CALL atprop_array_init(atprop%atecoul, natom)
250 : END IF
251 :
252 : ! 1/R contribution
253 11918 : do_ewald = dft_control%qs_control%dftb_control%do_ewald
254 11918 : IF (do_ewald) THEN
255 : ! Ewald sum
256 6420 : NULLIFY (ewald_env, ewald_pw)
257 : CALL get_qs_env(qs_env=qs_env, &
258 6420 : ewald_env=ewald_env, ewald_pw=ewald_pw)
259 6420 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
260 6420 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
261 6420 : CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
262 : CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, &
263 6420 : virial, use_virial, atprop=atprop)
264 0 : SELECT CASE (ewald_type)
265 : CASE DEFAULT
266 0 : CPABORT("Invalid Ewald type")
267 : CASE (do_ewald_none)
268 0 : CPABORT("Not allowed with DFTB")
269 : CASE (do_ewald_ewald)
270 0 : CPABORT("Standard Ewald not implemented in DFTB")
271 : CASE (do_ewald_pme)
272 0 : CPABORT("PME not implemented in DFTB")
273 : CASE (do_ewald_spme)
274 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
275 : gmcharge, mcharge, calculate_forces, virial, &
276 6420 : use_virial, atprop=atprop)
277 : END SELECT
278 : ELSE
279 : ! direct sum
280 : CALL get_qs_env(qs_env=qs_env, &
281 5498 : local_particles=local_particles)
282 17846 : DO ikind = 1, SIZE(local_particles%n_el)
283 26764 : DO ia = 1, local_particles%n_el(ikind)
284 8918 : iatom = local_particles%list(ikind)%array(ia)
285 31531 : DO jatom = 1, iatom - 1
286 41060 : rij = particle_set(iatom)%r - particle_set(jatom)%r
287 41060 : rij = pbc(rij, cell)
288 41060 : dr = SQRT(SUM(rij(:)**2))
289 10265 : gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
290 10265 : gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
291 20077 : DO i = 2, nmat
292 894 : gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
293 11159 : gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
294 : END DO
295 : END DO
296 : END DO
297 : END DO
298 5498 : CPASSERT(.NOT. use_virial)
299 : END IF
300 :
301 322482 : CALL para_env%sum(gmcharge(:, 1))
302 :
303 11918 : IF (do_ewald) THEN
304 : ! add self charge interaction and background charge contribution
305 143866 : gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
306 7242 : IF (ANY(periodic(:) == 1)) THEN
307 142502 : gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
308 : END IF
309 : END IF
310 :
311 167200 : energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1))
312 11918 : IF (atprop%energy) THEN
313 : CALL get_qs_env(qs_env=qs_env, &
314 466 : local_particles=local_particles)
315 1526 : DO ikind = 1, SIZE(local_particles%n_el)
316 1060 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
317 1060 : CALL get_dftb_atom_param(dftb_kind, zeff=zeff)
318 18006 : DO ia = 1, local_particles%n_el(ikind)
319 16480 : iatom = local_particles%list(ikind)%array(ia)
320 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
321 17540 : 0.5_dp*zeff*gmcharge(iatom, 1)
322 : END DO
323 : END DO
324 : END IF
325 :
326 11918 : IF (calculate_forces) THEN
327 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
328 : kind_of=kind_of, &
329 616 : atom_of_kind=atom_of_kind)
330 :
331 14830 : gmcharge(:, 2) = gmcharge(:, 2)*mcharge(:)
332 14830 : gmcharge(:, 3) = gmcharge(:, 3)*mcharge(:)
333 14830 : gmcharge(:, 4) = gmcharge(:, 4)*mcharge(:)
334 14830 : DO iatom = 1, natom
335 14214 : ikind = kind_of(iatom)
336 14214 : atom_i = atom_of_kind(iatom)
337 14214 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - gmcharge(iatom, 2)
338 14214 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - gmcharge(iatom, 3)
339 14830 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - gmcharge(iatom, 4)
340 : END DO
341 : END IF
342 :
343 11918 : do_gamma_stress = .FALSE.
344 11918 : IF (.NOT. just_energy .AND. use_virial) THEN
345 106 : IF (dft_control%nimages == 1) do_gamma_stress = .TRUE.
346 : END IF
347 :
348 11918 : IF (.NOT. just_energy) THEN
349 11714 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
350 11714 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
351 :
352 11714 : nimg = dft_control%nimages
353 11714 : NULLIFY (cell_to_index)
354 11714 : IF (nimg > 1) THEN
355 1988 : NULLIFY (kpoints)
356 1988 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
357 1988 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
358 : END IF
359 :
360 11714 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
361 328 : DO img = 1, nimg
362 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
363 328 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
364 : END DO
365 : END IF
366 :
367 11714 : NULLIFY (sap_int)
368 11714 : IF (do_gamma_stress) THEN
369 : ! derivative overlap integral (non collapsed)
370 88 : CALL dftb_dsint_list(qs_env, sap_int)
371 : END IF
372 :
373 11714 : IF (nimg == 1) THEN
374 : ! no k-points; all matrices have been transformed to periodic bsf
375 9726 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
376 1674188 : DO WHILE (dbcsr_iterator_blocks_left(iter))
377 1664462 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
378 1664462 : gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
379 3330194 : DO is = 1, SIZE(ks_matrix, 1)
380 1665732 : NULLIFY (ksblock)
381 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
382 1665732 : row=irow, col=icol, block=ksblock, found=found)
383 1665732 : CPASSERT(found)
384 25083404 : ksblock = ksblock - gmij*sblock
385 : END DO
386 1674188 : IF (calculate_forces) THEN
387 230484 : ikind = kind_of(irow)
388 230484 : atom_i = atom_of_kind(irow)
389 230484 : jkind = kind_of(icol)
390 230484 : atom_j = atom_of_kind(icol)
391 230484 : NULLIFY (pblock)
392 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
393 230484 : row=irow, col=icol, block=pblock, found=found)
394 230484 : CPASSERT(found)
395 921936 : DO i = 1, 3
396 691452 : NULLIFY (dsblock)
397 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
398 691452 : row=irow, col=icol, block=dsblock, found=found)
399 691452 : CPASSERT(found)
400 4854195 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
401 691452 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
402 691452 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
403 1613388 : fij(i) = fi
404 : END DO
405 : END IF
406 : END DO
407 9726 : CALL dbcsr_iterator_stop(iter)
408 : !
409 9726 : IF (do_gamma_stress) THEN
410 264 : DO ikind = 1, nkind
411 616 : DO jkind = 1, nkind
412 352 : iac = ikind + nkind*(jkind - 1)
413 352 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
414 8596 : DO ia = 1, sap_int(iac)%nalist
415 8076 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
416 8074 : iatom = sap_int(iac)%alist(ia)%aatom
417 177568 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
418 169142 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
419 676568 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
420 676568 : dr = SQRT(SUM(rij(:)**2))
421 177218 : IF (dr > 1.e-6_dp) THEN
422 165104 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
423 165104 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
424 165104 : icol = MAX(iatom, jatom)
425 165104 : irow = MIN(iatom, jatom)
426 165104 : NULLIFY (pblock)
427 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
428 165104 : row=irow, col=icol, block=pblock, found=found)
429 165104 : CPASSERT(found)
430 660416 : DO i = 1, 3
431 660416 : IF (irow == iatom) THEN
432 1721034 : fij(i) = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i))
433 : ELSE
434 1730961 : fij(i) = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
435 : END IF
436 : END DO
437 165104 : fi = 1.0_dp
438 165104 : IF (iatom == jatom) fi = 0.5_dp
439 165104 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
440 165104 : IF (atprop%stress) THEN
441 70678 : CALL virial_pair_force(atprop%atstress(:, :, iatom), 0.5_dp*fi, fij, rij)
442 70678 : CALL virial_pair_force(atprop%atstress(:, :, jatom), 0.5_dp*fi, fij, rij)
443 : END IF
444 : END IF
445 : END DO
446 : END DO
447 : END DO
448 : END DO
449 : END IF
450 : ELSE
451 1988 : NULLIFY (n_list)
452 1988 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
453 1988 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
454 245594 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
455 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
456 243606 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
457 :
458 243606 : icol = MAX(iatom, jatom)
459 243606 : irow = MIN(iatom, jatom)
460 243606 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
461 243606 : CPASSERT(ic > 0)
462 :
463 243606 : gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
464 :
465 243606 : NULLIFY (sblock)
466 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
467 243606 : row=irow, col=icol, block=sblock, found=found)
468 243606 : CPASSERT(found)
469 487212 : DO is = 1, SIZE(ks_matrix, 1)
470 243606 : NULLIFY (ksblock)
471 : CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
472 243606 : row=irow, col=icol, block=ksblock, found=found)
473 243606 : CPASSERT(found)
474 4849788 : ksblock = ksblock - gmij*sblock
475 : END DO
476 :
477 245594 : IF (calculate_forces) THEN
478 4182 : ikind = kind_of(iatom)
479 4182 : atom_i = atom_of_kind(iatom)
480 4182 : jkind = kind_of(jatom)
481 4182 : atom_j = atom_of_kind(jatom)
482 4182 : IF (irow == jatom) gmij = -gmij
483 4182 : NULLIFY (pblock)
484 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
485 4182 : row=irow, col=icol, block=pblock, found=found)
486 4182 : CPASSERT(found)
487 16728 : DO i = 1, 3
488 12546 : NULLIFY (dsblock)
489 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
490 12546 : row=irow, col=icol, block=dsblock, found=found)
491 12546 : CPASSERT(found)
492 227889 : fi = -2.0_dp*gmij*SUM(pblock*dsblock)
493 12546 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
494 12546 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
495 29274 : fij(i) = fi
496 : END DO
497 4182 : IF (use_virial) THEN
498 4182 : fi = 1.0_dp
499 4182 : IF (iatom == jatom) fi = 0.5_dp
500 4182 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
501 4182 : IF (atprop%stress) THEN
502 0 : CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
503 0 : CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
504 : END IF
505 : END IF
506 : END IF
507 : END DO
508 1988 : CALL neighbor_list_iterator_release(nl_iterator)
509 : END IF
510 :
511 11714 : IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
512 328 : DO img = 1, nimg
513 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
514 328 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
515 : END DO
516 : END IF
517 : END IF
518 :
519 11918 : IF (dft_control%qs_control%dftb_control%dftb3_diagonal) THEN
520 2818 : CALL get_qs_env(qs_env, nkind=nkind)
521 14090 : ALLOCATE (zeffk(nkind), xgamma(nkind))
522 8354 : DO ikind = 1, nkind
523 5536 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind)
524 8354 : CALL get_dftb_atom_param(dftb_kind, dudq=xgamma(ikind), zeff=zeffk(ikind))
525 : END DO
526 : ! Diagonal 3rd order correction (DFTB3)
527 : CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
528 2818 : sap_int, calculate_forces, just_energy)
529 2818 : DEALLOCATE (zeffk, xgamma)
530 : END IF
531 :
532 : ! QMMM
533 11918 : IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
534 : CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
535 1626 : calculate_forces, just_energy)
536 : END IF
537 :
538 11918 : IF (do_gamma_stress) THEN
539 88 : CALL release_sap_int(sap_int)
540 : END IF
541 :
542 11918 : DEALLOCATE (gmcharge)
543 :
544 11918 : CALL timestop(handle)
545 :
546 23836 : END SUBROUTINE build_dftb_coulomb
547 :
548 : ! **************************************************************************************************
549 : !> \brief Computes the short-range gamma parameter from exact Coulomb
550 : !> interaction of normalized exp(-a*r) charge distribution - 1/r
551 : !> \param r ...
552 : !> \param ga ...
553 : !> \param gb ...
554 : !> \param hb_para ...
555 : !> \return ...
556 : !> \par Literature
557 : !> Elstner et al, PRB 58 (1998) 7260
558 : !> \par History
559 : !> 10.2008 Axel Kohlmeyer - adding sr_damp
560 : !> 08.2014 JGH - adding flexible exponent for damping
561 : !> \version 1.1
562 : ! **************************************************************************************************
563 2799914 : FUNCTION gamma_rab_sr(r, ga, gb, hb_para) RESULT(gamma)
564 : REAL(dp), INTENT(in) :: r, ga, gb, hb_para
565 : REAL(dp) :: gamma
566 :
567 : REAL(dp) :: a, b, fac, g_sum
568 :
569 2799914 : gamma = 0.0_dp
570 2799914 : a = 3.2_dp*ga ! 3.2 = 16/5 in Eq. 18 and ff.
571 2799914 : b = 3.2_dp*gb
572 2799914 : g_sum = a + b
573 2799914 : IF (g_sum < tol_gamma) RETURN ! hardness screening
574 2799914 : IF (r < rtiny) THEN ! This is for short distances but non-onsite terms
575 : ! This gives also correct diagonal elements (a=b, r=0)
576 77641 : gamma = 0.5_dp*(a*b/g_sum + (a*b)**2/g_sum**3)
577 77641 : RETURN
578 : END IF
579 : !
580 : ! distinguish two cases: Gamma's are very close, e.g. for the same atom type,
581 : ! and Gamma's are different
582 : !
583 2722273 : IF (ABS(a - b) < rtiny) THEN
584 1536622 : fac = 1.6_dp*r*a*b/g_sum*(1.0_dp + a*b/g_sum**2)
585 1536622 : gamma = -(48.0_dp + 33._dp*fac + (9.0_dp + fac)*fac**2)*EXP(-fac)/(48._dp*r)
586 : ELSE
587 : gamma = -EXP(-a*r)*(0.5_dp*a*b**4/(a**2 - b**2)**2 - &
588 : (b**6 - 3._dp*a**2*b**4)/(r*(a**2 - b**2)**3)) - & ! a-> b
589 : EXP(-b*r)*(0.5_dp*b*a**4/(b**2 - a**2)**2 - &
590 1185651 : (a**6 - 3._dp*b**2*a**4)/(r*(b**2 - a**2)**3)) ! b-> a
591 : END IF
592 : !
593 : ! damping function for better short range hydrogen bonds.
594 : ! functional form from Hu H. et al., J. Phys. Chem. A 2007, 111, 5685-5691
595 : ! according to Elstner M, Theor. Chem. Acc. 2006, 116, 316-325,
596 : ! this should only be applied to a-b pairs involving hydrogen.
597 2722273 : IF (hb_para > 0.0_dp) THEN
598 438246 : gamma = gamma*EXP(-(0.5_dp*(ga + gb))**hb_para*r*r)
599 : END IF
600 : END FUNCTION gamma_rab_sr
601 :
602 : ! **************************************************************************************************
603 : !> \brief ...
604 : !> \param qs_env ...
605 : !> \param sap_int ...
606 : ! **************************************************************************************************
607 88 : SUBROUTINE dftb_dsint_list(qs_env, sap_int)
608 :
609 : TYPE(qs_environment_type), POINTER :: qs_env
610 : TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
611 :
612 : CHARACTER(LEN=*), PARAMETER :: routineN = 'dftb_dsint_list'
613 :
614 : INTEGER :: handle, i, iac, iatom, ikind, ilist, jatom, jkind, jneighbor, llm, lmaxi, lmaxj, &
615 : n1, n2, natorb_a, natorb_b, ngrd, ngrdcut, nkind, nlist, nneighbor
616 : INTEGER, DIMENSION(3) :: cell
617 : LOGICAL :: defined
618 : REAL(KIND=dp) :: ddr, dgrd, dr
619 : REAL(KIND=dp), DIMENSION(3) :: drij, rij
620 88 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, dsblockm, smatij, smatji
621 : TYPE(clist_type), POINTER :: clist
622 : TYPE(dft_control_type), POINTER :: dft_control
623 : TYPE(dftb_control_type), POINTER :: dftb_control
624 : TYPE(neighbor_list_iterator_p_type), &
625 88 : DIMENSION(:), POINTER :: nl_iterator
626 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
627 88 : POINTER :: sab_orb
628 : TYPE(qs_dftb_atom_type), POINTER :: dftb_kind_a, dftb_kind_b
629 : TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
630 88 : POINTER :: dftb_potential
631 : TYPE(qs_dftb_pairpot_type), POINTER :: dftb_param_ij, dftb_param_ji
632 88 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
633 :
634 88 : CALL timeset(routineN, handle)
635 :
636 88 : CALL get_qs_env(qs_env=qs_env, nkind=nkind)
637 88 : CPASSERT(.NOT. ASSOCIATED(sap_int))
638 616 : ALLOCATE (sap_int(nkind*nkind))
639 440 : DO i = 1, nkind*nkind
640 352 : NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
641 440 : sap_int(i)%nalist = 0
642 : END DO
643 :
644 : CALL get_qs_env(qs_env=qs_env, &
645 : qs_kind_set=qs_kind_set, &
646 : dft_control=dft_control, &
647 88 : sab_orb=sab_orb)
648 :
649 88 : dftb_control => dft_control%qs_control%dftb_control
650 :
651 88 : NULLIFY (dftb_potential)
652 : CALL get_qs_env(qs_env=qs_env, &
653 88 : dftb_potential=dftb_potential)
654 :
655 : ! loop over all atom pairs with a non-zero overlap (sab_orb)
656 88 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
657 169230 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
658 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, &
659 : jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
660 169142 : inode=jneighbor, cell=cell, r=rij)
661 169142 : iac = ikind + nkind*(jkind - 1)
662 : !
663 169142 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
664 : CALL get_dftb_atom_param(dftb_kind_a, &
665 169142 : defined=defined, lmax=lmaxi, natorb=natorb_a)
666 169142 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
667 169142 : CALL get_qs_kind(qs_kind_set(jkind), dftb_parameter=dftb_kind_b)
668 : CALL get_dftb_atom_param(dftb_kind_b, &
669 169142 : defined=defined, lmax=lmaxj, natorb=natorb_b)
670 169142 : IF (.NOT. defined .OR. natorb_b < 1) CYCLE
671 :
672 676568 : dr = SQRT(SUM(rij(:)**2))
673 :
674 : ! integral list
675 169142 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
676 344 : sap_int(iac)%a_kind = ikind
677 344 : sap_int(iac)%p_kind = jkind
678 344 : sap_int(iac)%nalist = nlist
679 9108 : ALLOCATE (sap_int(iac)%alist(nlist))
680 8420 : DO i = 1, nlist
681 8076 : NULLIFY (sap_int(iac)%alist(i)%clist)
682 8076 : sap_int(iac)%alist(i)%aatom = 0
683 8420 : sap_int(iac)%alist(i)%nclist = 0
684 : END DO
685 : END IF
686 169142 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
687 8074 : sap_int(iac)%alist(ilist)%aatom = iatom
688 8074 : sap_int(iac)%alist(ilist)%nclist = nneighbor
689 241808 : ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
690 177216 : DO i = 1, nneighbor
691 177216 : sap_int(iac)%alist(ilist)%clist(i)%catom = 0
692 : END DO
693 : END IF
694 169142 : clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
695 169142 : clist%catom = jatom
696 676568 : clist%cell = cell
697 676568 : clist%rac = rij
698 845710 : ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
699 169142 : NULLIFY (clist%achint)
700 3730163 : clist%acint = 0._dp
701 169142 : clist%nsgf_cnt = 0
702 169142 : NULLIFY (clist%sgf_list)
703 :
704 : ! retrieve information on S matrix
705 169142 : dftb_param_ij => dftb_potential(ikind, jkind)
706 169142 : dftb_param_ji => dftb_potential(jkind, ikind)
707 : ! assume table size and type is symmetric
708 169142 : ngrd = dftb_param_ij%ngrd
709 169142 : ngrdcut = dftb_param_ij%ngrdcut
710 169142 : dgrd = dftb_param_ij%dgrd
711 169142 : ddr = dgrd*0.1_dp
712 169142 : CPASSERT(dftb_param_ij%llm == dftb_param_ji%llm)
713 169142 : llm = dftb_param_ij%llm
714 169142 : smatij => dftb_param_ij%smat
715 169142 : smatji => dftb_param_ji%smat
716 :
717 507514 : IF (NINT(dr/dgrd) <= ngrdcut) THEN
718 168768 : IF (iatom == jatom .AND. dr < 0.001_dp) THEN
719 : ! diagonal block
720 : ELSE
721 : ! off-diagonal block
722 164730 : n1 = natorb_a
723 164730 : n2 = natorb_b
724 1153110 : ALLOCATE (dsblock(n1, n2), dsblockm(n1, n2))
725 658920 : DO i = 1, 3
726 3443355 : dsblock = 0._dp
727 3443355 : dsblockm = 0.0_dp
728 494190 : drij = rij
729 :
730 494190 : drij(i) = rij(i) - ddr
731 : CALL compute_block_sk(dsblockm, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
732 494190 : llm, lmaxi, lmaxj, iatom, iatom)
733 :
734 494190 : drij(i) = rij(i) + ddr
735 : CALL compute_block_sk(dsblock, smatij, smatji, drij, ngrd, ngrdcut, dgrd, &
736 494190 : llm, lmaxi, lmaxj, iatom, iatom)
737 :
738 6392520 : dsblock = dsblock - dsblockm
739 3443355 : dsblock = dsblock/(2.0_dp*ddr)
740 :
741 6557250 : clist%acint(1:n1, 1:n2, i) = -dsblock(1:n1, 1:n2)
742 : END DO
743 164730 : DEALLOCATE (dsblock, dsblockm)
744 : END IF
745 : END IF
746 :
747 : END DO
748 88 : CALL neighbor_list_iterator_release(nl_iterator)
749 :
750 88 : CALL timestop(handle)
751 :
752 88 : END SUBROUTINE dftb_dsint_list
753 :
754 : END MODULE qs_dftb_coulomb
755 :
|
