(git:374b731)
Loading...
Searching...
No Matches
xtb_coulomb.F
Go to the documentation of this file.
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 xTB
10!> \author JGH
11! **************************************************************************************************
12MODULE xtb_coulomb
13 USE ai_contraction, ONLY: block_add,&
14 contraction
15 USE ai_overlap, ONLY: overlap_ab
16 USE atomic_kind_types, ONLY: atomic_kind_type,&
17 get_atomic_kind_set
18 USE atprop_types, ONLY: atprop_array_init,&
19 atprop_type
20 USE basis_set_types, ONLY: gto_basis_set_p_type,&
21 gto_basis_set_type
22 USE cell_types, ONLY: cell_type,&
23 get_cell,&
24 pbc
25 USE cp_control_types, ONLY: dft_control_type,&
26 xtb_control_type
27 USE dbcsr_api, ONLY: dbcsr_add,&
28 dbcsr_get_block_p,&
29 dbcsr_iterator_blocks_left,&
30 dbcsr_iterator_next_block,&
31 dbcsr_iterator_start,&
32 dbcsr_iterator_stop,&
33 dbcsr_iterator_type,&
34 dbcsr_p_type
35 USE distribution_1d_types, ONLY: distribution_1d_type
36 USE ewald_environment_types, ONLY: ewald_env_get,&
37 ewald_environment_type
38 USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
39 tb_spme_evaluate
40 USE ewald_pw_types, ONLY: ewald_pw_type
41 USE kinds, ONLY: dp
42 USE kpoint_types, ONLY: get_kpoint_info,&
43 kpoint_type
44 USE mathconstants, ONLY: oorootpi,&
45 pi
46 USE message_passing, ONLY: mp_para_env_type
47 USE orbital_pointers, ONLY: ncoset
48 USE particle_types, ONLY: particle_type
49 USE pw_poisson_types, ONLY: do_ewald_ewald,&
50 do_ewald_none,&
51 do_ewald_pme,&
52 do_ewald_spme
53 USE qmmm_tb_coulomb, ONLY: build_tb_coulomb_qmqm
54 USE qs_dftb3_methods, ONLY: build_dftb3_diagonal
55 USE qs_energy_types, ONLY: qs_energy_type
56 USE qs_environment_types, ONLY: get_qs_env,&
57 qs_environment_type
58 USE qs_force_types, ONLY: qs_force_type
59 USE qs_integral_utils, ONLY: basis_set_list_setup,&
60 get_memory_usage
61 USE qs_kind_types, ONLY: get_qs_kind,&
62 qs_kind_type
63 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
64 neighbor_list_iterate,&
65 neighbor_list_iterator_create,&
66 neighbor_list_iterator_p_type,&
67 neighbor_list_iterator_release,&
68 neighbor_list_set_p_type
69 USE qs_rho_types, ONLY: qs_rho_get,&
70 qs_rho_type
71 USE sap_kind_types, ONLY: clist_type,&
72 release_sap_int,&
73 sap_int_type
74 USE virial_methods, ONLY: virial_pair_force
75 USE virial_types, ONLY: virial_type
76 USE xtb_types, ONLY: get_xtb_atom_param,&
77 xtb_atom_type
78#include "./base/base_uses.f90"
79
80 IMPLICIT NONE
81
82 PRIVATE
83
84 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_coulomb'
85
86 PUBLIC :: build_xtb_coulomb, gamma_rab_sr, dgamma_rab_sr, xtb_dsint_list
87
88CONTAINS
89
90! **************************************************************************************************
91!> \brief ...
92!> \param qs_env ...
93!> \param ks_matrix ...
94!> \param rho ...
95!> \param charges ...
96!> \param mcharge ...
97!> \param energy ...
98!> \param calculate_forces ...
99!> \param just_energy ...
100! **************************************************************************************************
101 SUBROUTINE build_xtb_coulomb(qs_env, ks_matrix, rho, charges, 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(:, :), INTENT(in) :: charges
108 REAL(dp), DIMENSION(:), INTENT(in) :: mcharge
109 TYPE(qs_energy_type), POINTER :: energy
110 LOGICAL, INTENT(in) :: calculate_forces, just_energy
111
112 CHARACTER(len=*), PARAMETER :: routinen = 'build_xtb_coulomb'
113
114 INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iac, iatom, ic, icol, ikind, img, &
115 irow, is, j, jatom, jkind, la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nimg, nj, &
116 nkind, nmat, za, zb
117 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
118 INTEGER, DIMENSION(25) :: laoa, laob
119 INTEGER, DIMENSION(3) :: cellind, periodic
120 INTEGER, DIMENSION(5) :: occ
121 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
122 LOGICAL :: defined, do_ewald, do_gamma_stress, &
123 found, use_virial
124 REAL(kind=dp) :: alpha, deth, dr, ecsr, etaa, etab, f1, &
125 f2, fi, gmij, kg, rcut, rcuta, rcutb, &
126 zeff
127 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: xgamma, zeffk
128 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gammab, gcij, gmcharge
129 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gchrg
130 REAL(kind=dp), DIMENSION(25) :: gcint
131 REAL(kind=dp), DIMENSION(3) :: fij, rij
132 REAL(kind=dp), DIMENSION(5) :: kappaa, kappab
133 REAL(kind=dp), DIMENSION(:, :), POINTER :: dsblock, ksblock, pblock, sblock
134 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dsint
135 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
136 TYPE(atprop_type), POINTER :: atprop
137 TYPE(cell_type), POINTER :: cell
138 TYPE(dbcsr_iterator_type) :: iter
139 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
140 TYPE(dft_control_type), POINTER :: dft_control
141 TYPE(distribution_1d_type), POINTER :: local_particles
142 TYPE(ewald_environment_type), POINTER :: ewald_env
143 TYPE(ewald_pw_type), POINTER :: ewald_pw
144 TYPE(kpoint_type), POINTER :: kpoints
145 TYPE(mp_para_env_type), POINTER :: para_env
146 TYPE(neighbor_list_iterator_p_type), &
147 DIMENSION(:), POINTER :: nl_iterator
148 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149 POINTER :: n_list
150 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
152 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
153 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
154 TYPE(virial_type), POINTER :: virial
155 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b, xtb_kind
156 TYPE(xtb_control_type), POINTER :: xtb_control
157
158 CALL timeset(routinen, handle)
159
160 NULLIFY (matrix_p, matrix_s, virial, atprop, dft_control)
161
162 CALL get_qs_env(qs_env, &
163 qs_kind_set=qs_kind_set, &
164 particle_set=particle_set, &
165 cell=cell, &
166 virial=virial, &
167 atprop=atprop, &
168 dft_control=dft_control)
169
170 xtb_control => dft_control%qs_control%xtb_control
171
172 use_virial = .false.
173 IF (calculate_forces) THEN
174 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
175 END IF
176
177 do_gamma_stress = .false.
178 IF (.NOT. just_energy .AND. use_virial) THEN
179 IF (dft_control%nimages == 1) do_gamma_stress = .true.
180 END IF
181
182 IF (atprop%energy) THEN
183 CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
184 natom = SIZE(particle_set)
185 CALL atprop_array_init(atprop%atecoul, natom)
186 END IF
187
188 IF (calculate_forces) THEN
189 nmat = 4
190 ELSE
191 nmat = 1
192 END IF
193
194 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
195 ALLOCATE (gchrg(natom, 5, nmat))
196 gchrg = 0._dp
197 ALLOCATE (gmcharge(natom, nmat))
198 gmcharge = 0._dp
199
200 ! short range contribution (gamma)
201 ! loop over all atom pairs (sab_xtbe)
202 kg = xtb_control%kg
203 NULLIFY (n_list)
204 IF (xtb_control%old_coulomb_damping) THEN
205 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
206 ELSE
207 CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
208 END IF
209 CALL neighbor_list_iterator_create(nl_iterator, n_list)
210 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
211 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
212 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
213 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
214 CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
215 IF (.NOT. defined .OR. natorb_a < 1) cycle
216 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
217 CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
218 IF (.NOT. defined .OR. natorb_b < 1) cycle
219 ! atomic parameters
220 CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
221 CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
222 ! gamma matrix
223 ni = lmaxa + 1
224 nj = lmaxb + 1
225 ALLOCATE (gammab(ni, nj))
226 rcut = rcuta + rcutb
227 dr = sqrt(sum(rij(:)**2))
228 CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
229 gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + matmul(gammab, charges(jatom, 1:nj))
230 IF (iatom /= jatom) THEN
231 gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + matmul(charges(iatom, 1:ni), gammab)
232 END IF
233 IF (calculate_forces) THEN
234 IF (dr > 1.e-6_dp) THEN
235 CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
236 DO i = 1, 3
237 gchrg(iatom, 1:ni, i + 1) = gchrg(iatom, 1:ni, i + 1) &
238 + matmul(gammab, charges(jatom, 1:nj))*rij(i)/dr
239 IF (iatom /= jatom) THEN
240 gchrg(jatom, 1:nj, i + 1) = gchrg(jatom, 1:nj, i + 1) &
241 - matmul(charges(iatom, 1:ni), gammab)*rij(i)/dr
242 END IF
243 END DO
244 IF (use_virial) THEN
245 gcint(1:ni) = matmul(gammab, charges(jatom, 1:nj))
246 DO i = 1, 3
247 fij(i) = -sum(charges(iatom, 1:ni)*gcint(1:ni))*rij(i)/dr
248 END DO
249 fi = 1.0_dp
250 IF (iatom == jatom) fi = 0.5_dp
251 CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
252 IF (atprop%stress) THEN
253 CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
254 CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
255 END IF
256 END IF
257 END IF
258 END IF
259 DEALLOCATE (gammab)
260 END DO
261 CALL neighbor_list_iterator_release(nl_iterator)
262
263 ! 1/R contribution
264
265 IF (xtb_control%coulomb_lr) THEN
266 do_ewald = xtb_control%do_ewald
267 IF (do_ewald) THEN
268 ! Ewald sum
269 NULLIFY (ewald_env, ewald_pw)
270 CALL get_qs_env(qs_env=qs_env, &
271 ewald_env=ewald_env, ewald_pw=ewald_pw)
272 CALL get_cell(cell=cell, periodic=periodic, deth=deth)
273 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
274 CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
275 CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, &
276 use_virial, atprop=atprop)
277 SELECT CASE (ewald_type)
278 CASE DEFAULT
279 cpabort("Invalid Ewald type")
280 CASE (do_ewald_none)
281 cpabort("Not allowed with DFTB")
282 CASE (do_ewald_ewald)
283 cpabort("Standard Ewald not implemented in DFTB")
284 CASE (do_ewald_pme)
285 cpabort("PME not implemented in DFTB")
286 CASE (do_ewald_spme)
287 CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
288 gmcharge, mcharge, calculate_forces, virial, &
289 use_virial, atprop=atprop)
290 END SELECT
291 ELSE
292 ! direct sum
293 CALL get_qs_env(qs_env=qs_env, &
294 local_particles=local_particles)
295 DO ikind = 1, SIZE(local_particles%n_el)
296 DO ia = 1, local_particles%n_el(ikind)
297 iatom = local_particles%list(ikind)%array(ia)
298 DO jatom = 1, iatom - 1
299 rij = particle_set(iatom)%r - particle_set(jatom)%r
300 rij = pbc(rij, cell)
301 dr = sqrt(sum(rij(:)**2))
302 IF (dr > 1.e-6_dp) THEN
303 gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
304 gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
305 DO i = 2, nmat
306 gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
307 gmcharge(jatom, i) = gmcharge(jatom, i) - rij(i - 1)*mcharge(iatom)/dr**3
308 END DO
309 END IF
310 END DO
311 END DO
312 END DO
313 cpassert(.NOT. use_virial)
314 END IF
315 END IF
316
317 ! global sum of gamma*p arrays
318 CALL get_qs_env(qs_env=qs_env, &
319 atomic_kind_set=atomic_kind_set, &
320 force=force, para_env=para_env)
321 CALL para_env%sum(gmcharge(:, 1))
322 CALL para_env%sum(gchrg(:, :, 1))
323
324 IF (xtb_control%coulomb_lr) THEN
325 IF (do_ewald) THEN
326 ! add self charge interaction and background charge contribution
327 gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
328 IF (any(periodic(:) == 1)) THEN
329 gmcharge(:, 1) = gmcharge(:, 1) - pi/alpha**2/deth
330 END IF
331 END IF
332 END IF
333
334 ! energy
335 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
336 kind_of=kind_of, &
337 atom_of_kind=atom_of_kind)
338 ecsr = 0.0_dp
339 DO iatom = 1, natom
340 ikind = kind_of(iatom)
341 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
342 CALL get_xtb_atom_param(xtb_kind, lmax=ni)
343 ni = ni + 1
344 ecsr = ecsr + sum(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
345 END DO
346
347 energy%hartree = energy%hartree + 0.5_dp*ecsr
348 energy%hartree = energy%hartree + 0.5_dp*sum(mcharge(:)*gmcharge(:, 1))
349
350 IF (atprop%energy) THEN
351 CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
352 DO ikind = 1, SIZE(local_particles%n_el)
353 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
354 CALL get_xtb_atom_param(xtb_kind, lmax=ni, occupation=occ)
355 ni = ni + 1
356 zeff = sum(real(occ, kind=dp))
357 DO ia = 1, local_particles%n_el(ikind)
358 iatom = local_particles%list(ikind)%array(ia)
359 atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
360 0.5_dp*sum(real(occ(1:ni), kind=dp)*gchrg(iatom, 1:ni, 1))
361 atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
362 0.5_dp*zeff*gmcharge(iatom, 1)
363 END DO
364 END DO
365 END IF
366
367 IF (calculate_forces) THEN
368 DO iatom = 1, natom
369 ikind = kind_of(iatom)
370 atom_i = atom_of_kind(iatom)
371 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
372 CALL get_xtb_atom_param(xtb_kind, lmax=ni)
373 ! short range
374 ni = ni + 1
375 DO i = 1, 3
376 fij(i) = sum(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i + 1))
377 END DO
378 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
379 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
380 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
381 ! long range
382 DO i = 1, 3
383 fij(i) = gmcharge(iatom, i + 1)*mcharge(iatom)
384 END DO
385 force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
386 force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
387 force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
388 END DO
389 END IF
390
391 IF (.NOT. just_energy) THEN
392 CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
393 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
394
395 nimg = dft_control%nimages
396 NULLIFY (cell_to_index)
397 IF (nimg > 1) THEN
398 NULLIFY (kpoints)
399 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
400 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
401 END IF
402
403 IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
404 DO img = 1, nimg
405 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
406 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
407 END DO
408 END IF
409
410 NULLIFY (sap_int)
411 IF (do_gamma_stress) THEN
412 ! derivative overlap integral (non collapsed)
413 CALL xtb_dsint_list(qs_env, sap_int)
414 END IF
415
416 IF (nimg == 1) THEN
417 ! no k-points; all matrices have been transformed to periodic bsf
418 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
419 DO WHILE (dbcsr_iterator_blocks_left(iter))
420 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
421 ikind = kind_of(irow)
422 jkind = kind_of(icol)
423
424 ! atomic parameters
425 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
426 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
427 CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
428 CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
429
430 ni = SIZE(sblock, 1)
431 nj = SIZE(sblock, 2)
432 ALLOCATE (gcij(ni, nj))
433 DO i = 1, ni
434 DO j = 1, nj
435 la = laoa(i) + 1
436 lb = laob(j) + 1
437 gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
438 END DO
439 END DO
440 gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
441 DO is = 1, SIZE(ks_matrix, 1)
442 NULLIFY (ksblock)
443 CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
444 row=irow, col=icol, block=ksblock, found=found)
445 cpassert(found)
446 ksblock = ksblock - gcij*sblock
447 ksblock = ksblock - gmij*sblock
448 END DO
449 IF (calculate_forces) THEN
450 atom_i = atom_of_kind(irow)
451 atom_j = atom_of_kind(icol)
452 NULLIFY (pblock)
453 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
454 row=irow, col=icol, block=pblock, found=found)
455 cpassert(found)
456 DO i = 1, 3
457 NULLIFY (dsblock)
458 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
459 row=irow, col=icol, block=dsblock, found=found)
460 cpassert(found)
461 fij(i) = 0.0_dp
462 ! short range
463 fi = -2.0_dp*sum(pblock*dsblock*gcij)
464 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
465 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
466 fij(i) = fij(i) + fi
467 ! long range
468 fi = -2.0_dp*gmij*sum(pblock*dsblock)
469 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
470 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
471 fij(i) = fij(i) + fi
472 END DO
473 END IF
474 DEALLOCATE (gcij)
475 END DO
476 CALL dbcsr_iterator_stop(iter)
477 ! stress tensor (needs recalculation of overlap integrals)
478 IF (do_gamma_stress) THEN
479 DO ikind = 1, nkind
480 DO jkind = 1, nkind
481 iac = ikind + nkind*(jkind - 1)
482 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
483 ! atomic parameters
484 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
485 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
486 CALL get_xtb_atom_param(xtb_atom_a, lao=laoa, natorb=ni)
487 CALL get_xtb_atom_param(xtb_atom_b, lao=laob, natorb=nj)
488 DO ia = 1, sap_int(iac)%nalist
489 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
490 iatom = sap_int(iac)%alist(ia)%aatom
491 DO ic = 1, sap_int(iac)%alist(ia)%nclist
492 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
493 rij = sap_int(iac)%alist(ia)%clist(ic)%rac
494 dr = sqrt(sum(rij(:)**2))
495 IF (dr > 1.e-6_dp) THEN
496 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
497 ALLOCATE (gcij(ni, nj))
498 DO i = 1, ni
499 DO j = 1, nj
500 la = laoa(i) + 1
501 lb = laob(j) + 1
502 gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
503 END DO
504 END DO
505 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
506 icol = max(iatom, jatom)
507 irow = min(iatom, jatom)
508 NULLIFY (pblock)
509 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
510 row=irow, col=icol, block=pblock, found=found)
511 cpassert(found)
512 fij = 0.0_dp
513 DO i = 1, 3
514 ! short/long range
515 IF (irow == iatom) THEN
516 f1 = -2.0_dp*sum(pblock*dsint(:, :, i)*gcij)
517 f2 = -2.0_dp*gmij*sum(pblock*dsint(:, :, i))
518 ELSE
519 f1 = -2.0_dp*sum(transpose(pblock)*dsint(:, :, i)*gcij)
520 f2 = -2.0_dp*gmij*sum(transpose(pblock)*dsint(:, :, i))
521 END IF
522 fij(i) = f1 + f2
523 END DO
524 DEALLOCATE (gcij)
525 fi = 1.0_dp
526 IF (iatom == jatom) fi = 0.5_dp
527 CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
528 IF (atprop%stress) THEN
529 CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
530 CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
531 END IF
532 END IF
533 END DO
534 END DO
535 END DO
536 END DO
537 END IF
538 ELSE
539 NULLIFY (n_list)
540 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
541 CALL neighbor_list_iterator_create(nl_iterator, n_list)
542 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
543 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
544 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
545
546 icol = max(iatom, jatom)
547 irow = min(iatom, jatom)
548
549 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
550 cpassert(ic > 0)
551
552 NULLIFY (sblock)
553 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
554 row=irow, col=icol, block=sblock, found=found)
555 cpassert(found)
556
557 ! atomic parameters
558 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
559 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
560 CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
561 CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
562
563 ni = SIZE(sblock, 1)
564 nj = SIZE(sblock, 2)
565 ALLOCATE (gcij(ni, nj))
566 DO i = 1, ni
567 DO j = 1, nj
568 IF (irow == iatom) THEN
569 la = laoa(i) + 1
570 lb = laob(j) + 1
571 gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
572 ELSE
573 la = laoa(j) + 1
574 lb = laob(i) + 1
575 gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
576 END IF
577 END DO
578 END DO
579 gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
580 DO is = 1, SIZE(ks_matrix, 1)
581 NULLIFY (ksblock)
582 CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
583 row=irow, col=icol, block=ksblock, found=found)
584 cpassert(found)
585 ksblock = ksblock - gcij*sblock
586 ksblock = ksblock - gmij*sblock
587 END DO
588
589 IF (calculate_forces) THEN
590 atom_i = atom_of_kind(iatom)
591 atom_j = atom_of_kind(jatom)
592 IF (irow /= iatom) THEN
593 gmij = -gmij
594 gcij = -gcij
595 END IF
596 NULLIFY (pblock)
597 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
598 row=irow, col=icol, block=pblock, found=found)
599 cpassert(found)
600 DO i = 1, 3
601 NULLIFY (dsblock)
602 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
603 row=irow, col=icol, block=dsblock, found=found)
604 cpassert(found)
605 fij(i) = 0.0_dp
606 ! short range
607 fi = -2.0_dp*sum(pblock*dsblock*gcij)
608 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
609 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
610 fij(i) = fij(i) + fi
611 ! long range
612 fi = -2.0_dp*gmij*sum(pblock*dsblock)
613 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
614 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
615 fij(i) = fij(i) + fi
616 END DO
617 IF (use_virial) THEN
618 dr = sqrt(sum(rij(:)**2))
619 IF (dr > 1.e-6_dp) THEN
620 fi = 1.0_dp
621 IF (iatom == jatom) fi = 0.5_dp
622 CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
623 IF (atprop%stress) THEN
624 CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
625 CALL virial_pair_force(atprop%atstress(:, :, jatom), fi*0.5_dp, fij, rij)
626 END IF
627 END IF
628 END IF
629 END IF
630 DEALLOCATE (gcij)
631
632 END DO
633 CALL neighbor_list_iterator_release(nl_iterator)
634 END IF
635
636 IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
637 DO img = 1, nimg
638 CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
639 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
640 END DO
641 END IF
642 END IF
643
644 IF (xtb_control%tb3_interaction) THEN
645 CALL get_qs_env(qs_env, nkind=nkind)
646 ALLOCATE (zeffk(nkind), xgamma(nkind))
647 DO ikind = 1, nkind
648 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
649 CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind), zeff=zeffk(ikind))
650 END DO
651 ! Diagonal 3rd order correction (DFTB3)
652 CALL build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeffk, &
653 sap_int, calculate_forces, just_energy)
654 DEALLOCATE (zeffk, xgamma)
655 END IF
656
657 ! QMMM
658 IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
659 CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
660 calculate_forces, just_energy)
661 END IF
662
663 IF (do_gamma_stress) THEN
664 CALL release_sap_int(sap_int)
665 END IF
666
667 CALL timestop(handle)
668
669 END SUBROUTINE build_xtb_coulomb
670
671! **************************************************************************************************
672!> \brief Computes the short-range gamma parameter from
673!> Nataga-Mishimoto-Ohno-Klopman formula for xTB
674!> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
675!> behaviour. We use a cutoff function to smoothly remove this part.
676!> However, this will change energies and effect final results.
677!>
678!> \param gmat ...
679!> \param rab ...
680!> \param nla ...
681!> \param kappaa ...
682!> \param etaa ...
683!> \param nlb ...
684!> \param kappab ...
685!> \param etab ...
686!> \param kg ...
687!> \param rcut ...
688!> \par History
689!> 10.2018 JGH
690!> \version 1.1
691! **************************************************************************************************
692 SUBROUTINE gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
693 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: gmat
694 REAL(dp), INTENT(IN) :: rab
695 INTEGER, INTENT(IN) :: nla
696 REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
697 REAL(dp), INTENT(IN) :: etaa
698 INTEGER, INTENT(IN) :: nlb
699 REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
700 REAL(dp), INTENT(IN) :: etab, kg, rcut
701
702 REAL(kind=dp), PARAMETER :: rsmooth = 1.0_dp
703
704 INTEGER :: i, j
705 REAL(kind=dp) :: fcut, r, rk, x
706 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
707
708 ALLOCATE (eta(nla, nlb))
709 eta = 0.0_dp
710
711 DO j = 1, nlb
712 DO i = 1, nla
713 eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
714 eta(i, j) = 2._dp/eta(i, j)
715 END DO
716 END DO
717
718 gmat = 0.0_dp
719 IF (rab < 1.e-6_dp) THEN
720 ! on site terms
721 gmat(:, :) = eta(:, :)
722 ELSEIF (rab > rcut) THEN
723 ! do nothing
724 ELSE
725 rk = rab**kg
726 eta = eta**(-kg)
727 IF (rab < rcut - rsmooth) THEN
728 fcut = 1.0_dp
729 ELSE
730 r = rab - (rcut - rsmooth)
731 x = r/rsmooth
732 fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
733 END IF
734 gmat(:, :) = fcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg) - fcut/rab
735 END IF
736
737 DEALLOCATE (eta)
738
739 END SUBROUTINE gamma_rab_sr
740
741! **************************************************************************************************
742!> \brief Computes the derivative of the short-range gamma parameter from
743!> Nataga-Mishimoto-Ohno-Klopman formula for xTB
744!> WARNING: The xTB function (gamma - 1/r) has still an l-dependent longrange (1/r^3)
745!> behaviour. We use a cutoff function to smoothly remove this part.
746!> However, this will change energies and effect final results.
747!>
748!> \param dgmat ...
749!> \param rab ...
750!> \param nla ...
751!> \param kappaa ...
752!> \param etaa ...
753!> \param nlb ...
754!> \param kappab ...
755!> \param etab ...
756!> \param kg ...
757!> \param rcut ...
758!> \par History
759!> 10.2018 JGH
760!> \version 1.1
761! **************************************************************************************************
762 SUBROUTINE dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
763 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: dgmat
764 REAL(dp), INTENT(IN) :: rab
765 INTEGER, INTENT(IN) :: nla
766 REAL(dp), DIMENSION(:), INTENT(IN) :: kappaa
767 REAL(dp), INTENT(IN) :: etaa
768 INTEGER, INTENT(IN) :: nlb
769 REAL(dp), DIMENSION(:), INTENT(IN) :: kappab
770 REAL(dp), INTENT(IN) :: etab, kg, rcut
771
772 REAL(kind=dp), PARAMETER :: rsmooth = 1.0_dp
773
774 INTEGER :: i, j
775 REAL(kind=dp) :: dfcut, fcut, r, rk, x
776 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: eta
777
778 ALLOCATE (eta(nla, nlb))
779
780 DO j = 1, nlb
781 DO i = 1, nla
782 eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
783 eta(i, j) = 2._dp/eta(i, j)
784 END DO
785 END DO
786
787 IF (rab < 1.e-6) THEN
788 ! on site terms
789 dgmat(:, :) = 0.0_dp
790 ELSEIF (rab > rcut) THEN
791 dgmat(:, :) = 0.0_dp
792 ELSE
793 eta = eta**(-kg)
794 rk = rab**kg
795 IF (rab < rcut - rsmooth) THEN
796 fcut = 1.0_dp
797 dfcut = 0.0_dp
798 ELSE
799 r = rab - (rcut - rsmooth)
800 x = r/rsmooth
801 fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
802 dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
803 dfcut = dfcut/rsmooth
804 END IF
805 dgmat(:, :) = dfcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg)
806 dgmat(:, :) = dgmat(:, :) - dfcut/rab + fcut/rab**2
807 dgmat(:, :) = dgmat(:, :) - fcut/(rk + eta(:, :))*(1._dp/(rk + eta(:, :)))**(1._dp/kg)*rk/rab
808 END IF
809
810 DEALLOCATE (eta)
811
812 END SUBROUTINE dgamma_rab_sr
813
814! **************************************************************************************************
815!> \brief ...
816!> \param qs_env ...
817!> \param sap_int ...
818! **************************************************************************************************
819 SUBROUTINE xtb_dsint_list(qs_env, sap_int)
820
821 TYPE(qs_environment_type), POINTER :: qs_env
822 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
823
824 CHARACTER(LEN=*), PARAMETER :: routinen = 'xtb_dsint_list'
825
826 INTEGER :: handle, i, iac, iatom, ikind, ilist, iset, jatom, jkind, jneighbor, jset, ldsab, &
827 n1, n2, natorb_a, natorb_b, ncoa, ncob, nkind, nlist, nneighbor, nseta, nsetb, sgfa, sgfb
828 INTEGER, DIMENSION(3) :: cell
829 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
830 npgfb, nsgfa, nsgfb
831 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
832 LOGICAL :: defined
833 REAL(kind=dp) :: dr
834 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: owork
835 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: oint, sint
836 REAL(kind=dp), DIMENSION(3) :: rij
837 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
838 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
839 TYPE(clist_type), POINTER :: clist
840 TYPE(dft_control_type), POINTER :: dft_control
841 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
842 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
843 TYPE(neighbor_list_iterator_p_type), &
844 DIMENSION(:), POINTER :: nl_iterator
845 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
846 POINTER :: sab_orb
847 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
848 TYPE(xtb_atom_type), POINTER :: xtb_atom_a, xtb_atom_b
849
850 CALL timeset(routinen, handle)
851
852 CALL get_qs_env(qs_env=qs_env, nkind=nkind)
853 cpassert(.NOT. ASSOCIATED(sap_int))
854 ALLOCATE (sap_int(nkind*nkind))
855 DO i = 1, nkind*nkind
856 NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
857 sap_int(i)%nalist = 0
858 END DO
859
860 CALL get_qs_env(qs_env=qs_env, &
861 qs_kind_set=qs_kind_set, &
862 dft_control=dft_control, &
863 sab_orb=sab_orb)
864
865 ! set up basis set lists
866 ALLOCATE (basis_set_list(nkind))
867 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
868
869 ! loop over all atom pairs with a non-zero overlap (sab_orb)
870 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
871 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
872 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, &
873 jatom=jatom, nlist=nlist, ilist=ilist, nnode=nneighbor, &
874 inode=jneighbor, cell=cell, r=rij)
875 iac = ikind + nkind*(jkind - 1)
876 !
877 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
878 CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
879 IF (.NOT. defined .OR. natorb_a < 1) cycle
880 CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
881 CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
882 IF (.NOT. defined .OR. natorb_b < 1) cycle
883
884 dr = sqrt(sum(rij(:)**2))
885
886 ! integral list
887 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
888 sap_int(iac)%a_kind = ikind
889 sap_int(iac)%p_kind = jkind
890 sap_int(iac)%nalist = nlist
891 ALLOCATE (sap_int(iac)%alist(nlist))
892 DO i = 1, nlist
893 NULLIFY (sap_int(iac)%alist(i)%clist)
894 sap_int(iac)%alist(i)%aatom = 0
895 sap_int(iac)%alist(i)%nclist = 0
896 END DO
897 END IF
898 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
899 sap_int(iac)%alist(ilist)%aatom = iatom
900 sap_int(iac)%alist(ilist)%nclist = nneighbor
901 ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
902 DO i = 1, nneighbor
903 sap_int(iac)%alist(ilist)%clist(i)%catom = 0
904 END DO
905 END IF
906 clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
907 clist%catom = jatom
908 clist%cell = cell
909 clist%rac = rij
910 ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
911 NULLIFY (clist%achint)
912 clist%acint = 0._dp
913 clist%nsgf_cnt = 0
914 NULLIFY (clist%sgf_list)
915
916 ! overlap
917 basis_set_a => basis_set_list(ikind)%gto_basis_set
918 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
919 basis_set_b => basis_set_list(jkind)%gto_basis_set
920 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
921 ! basis ikind
922 first_sgfa => basis_set_a%first_sgf
923 la_max => basis_set_a%lmax
924 la_min => basis_set_a%lmin
925 npgfa => basis_set_a%npgf
926 nseta = basis_set_a%nset
927 nsgfa => basis_set_a%nsgf_set
928 rpgfa => basis_set_a%pgf_radius
929 set_radius_a => basis_set_a%set_radius
930 scon_a => basis_set_a%scon
931 zeta => basis_set_a%zet
932 ! basis jkind
933 first_sgfb => basis_set_b%first_sgf
934 lb_max => basis_set_b%lmax
935 lb_min => basis_set_b%lmin
936 npgfb => basis_set_b%npgf
937 nsetb = basis_set_b%nset
938 nsgfb => basis_set_b%nsgf_set
939 rpgfb => basis_set_b%pgf_radius
940 set_radius_b => basis_set_b%set_radius
941 scon_b => basis_set_b%scon
942 zetb => basis_set_b%zet
943
944 ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
945 ALLOCATE (oint(ldsab, ldsab, 4), owork(ldsab, ldsab))
946 ALLOCATE (sint(natorb_a, natorb_b, 4))
947 sint = 0.0_dp
948
949 DO iset = 1, nseta
950 ncoa = npgfa(iset)*ncoset(la_max(iset))
951 n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
952 sgfa = first_sgfa(1, iset)
953 DO jset = 1, nsetb
954 IF (set_radius_a(iset) + set_radius_b(jset) < dr) cycle
955 ncob = npgfb(jset)*ncoset(lb_max(jset))
956 n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
957 sgfb = first_sgfb(1, jset)
958 CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
959 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
960 rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
961 ! Contraction
962 DO i = 1, 4
963 CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
964 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.false.)
965 CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), &
966 sgfa, sgfb, trans=.false.)
967 END DO
968 END DO
969 END DO
970 ! update dS/dR matrix
971 clist%acint(1:natorb_a, 1:natorb_b, 1:3) = sint(1:natorb_a, 1:natorb_b, 2:4)
972
973 DEALLOCATE (oint, owork, sint)
974
975 END DO
976 CALL neighbor_list_iterator_release(nl_iterator)
977
978 DEALLOCATE (basis_set_list)
979
980 CALL timestop(handle)
981
982 END SUBROUTINE xtb_dsint_list
983
984END MODULE xtb_coulomb
985
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
cell_types::pbc
Definition cell_types.F:92
qs_integral_utils::get_memory_usage
Definition qs_integral_utils.F:34
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
ai_overlap
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition ai_overlap.F:18
ai_overlap::overlap_ab
subroutine, public overlap_ab(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, sab, dab, ddab)
Calculation of the two-center overlap integrals [a|b] over Cartesian Gaussian-type functions....
Definition ai_overlap.F:680
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
atprop_types::atprop_array_init
subroutine, public atprop_array_init(atarray, natom)
...
Definition atprop_types.F:98
basis_set_types
Definition basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:123
ewald_methods_tb
Calculation of Ewald contributions in DFTB.
Definition ewald_methods_tb.F:12
ewald_methods_tb::tb_ewald_overlap
subroutine, public tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial, atprop)
...
Definition ewald_methods_tb.F:364
ewald_methods_tb::tb_spme_evaluate
subroutine, public tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
...
Definition ewald_methods_tb.F:90
ewald_pw_types
pw_types
Definition ewald_pw_types.F:12
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:333
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_poisson_types::do_ewald_pme
integer, parameter, public do_ewald_pme
Definition pw_poisson_types.F:75
pw_poisson_types::do_ewald_ewald
integer, parameter, public do_ewald_ewald
Definition pw_poisson_types.F:75
pw_poisson_types::do_ewald_none
integer, parameter, public do_ewald_none
Definition pw_poisson_types.F:75
pw_poisson_types::do_ewald_spme
integer, parameter, public do_ewald_spme
Definition pw_poisson_types.F:75
qmmm_tb_coulomb
Calculation of QMMM Coulomb contributions in TB.
Definition qmmm_tb_coulomb.F:12
qmmm_tb_coulomb::build_tb_coulomb_qmqm
subroutine, public build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Definition qmmm_tb_coulomb.F:73
qs_dftb3_methods
Calculation of DFTB3 Terms.
Definition qs_dftb3_methods.F:12
qs_dftb3_methods::build_dftb3_diagonal
subroutine, public build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, sap_int, calculate_forces, just_energy)
...
Definition qs_dftb3_methods.F:74
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_force_types
Definition qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
sap_kind_types
General overlap type integrals containers.
Definition sap_kind_types.F:12
sap_kind_types::release_sap_int
subroutine, public release_sap_int(sap_int)
...
Definition sap_kind_types.F:86
virial_methods
Definition virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition virial_methods.F:139
virial_types
Definition virial_types.F:13
xtb_coulomb
Calculation of Coulomb contributions in xTB.
Definition xtb_coulomb.F:12
xtb_coulomb::gamma_rab_sr
subroutine, public gamma_rab_sr(gmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
Computes the short-range gamma parameter from Nataga-Mishimoto-Ohno-Klopman formula for xTB WARNING: ...
Definition xtb_coulomb.F:693
xtb_coulomb::dgamma_rab_sr
subroutine, public dgamma_rab_sr(dgmat, rab, nla, kappaa, etaa, nlb, kappab, etab, kg, rcut)
Computes the derivative of the short-range gamma parameter from Nataga-Mishimoto-Ohno-Klopman formula...
Definition xtb_coulomb.F:763
xtb_coulomb::build_xtb_coulomb
subroutine, public build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, calculate_forces, just_energy)
...
Definition xtb_coulomb.F:103
xtb_coulomb::xtb_dsint_list
subroutine, public xtb_dsint_list(qs_env, sap_int)
...
Definition xtb_coulomb.F:820
xtb_types
Definition of the xTB parameter types.
Definition xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...
Definition xtb_types.F:175
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_control_types::xtb_control_type
Definition cp_control_types.F:127
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:55
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:64
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:138
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
sap_kind_types::clist_type
Definition sap_kind_types.F:43
sap_kind_types::sap_int_type
Definition sap_kind_types.F:64
virial_types::virial_type
Definition virial_types.F:26
xtb_types::xtb_atom_type
Definition xtb_types.F:42