(git:f7c1873)
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 ! **************************************************************************************************
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,&
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
37  ewald_environment_type
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,&
51  do_ewald_pme,&
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
60  get_memory_usage
61  USE qs_kind_types, ONLY: get_qs_kind,&
62  qs_kind_type
66  neighbor_list_iterator_p_type,&
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,&
73  sap_int_type
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 
87 
88 CONTAINS
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 
984 END MODULE xtb_coulomb
985 
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Calculation of the overlap integrals over Cartesian Gaussian-type functions.
Definition: ai_overlap.F:18
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
Define the atomic kind types and their sub types.
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.
Holds information on atomic properties.
Definition: atprop_types.F:14
subroutine, public atprop_array_init(atarray, natom)
...
Definition: atprop_types.F:98
Handles all functions related to the CELL.
Definition: cell_types.F:15
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
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
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.
Calculation of Ewald contributions in DFTB.
subroutine, public tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, use_virial, atprop)
...
subroutine, public tb_spme_evaluate(ewald_env, ewald_pw, particle_set, box, gmcharge, mcharge, calculate_forces, virial, use_virial, atprop)
...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
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
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public oorootpi
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
Define the data structure for the particle information.
functions related to the poisson solver on regular grids
integer, parameter, public do_ewald_pme
integer, parameter, public do_ewald_ewald
integer, parameter, public do_ewald_none
integer, parameter, public do_ewald_spme
Calculation of QMMM Coulomb contributions in TB.
subroutine, public build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, calculate_forces, just_energy)
...
Calculation of DFTB3 Terms.
subroutine, public build_dftb3_diagonal(qs_env, ks_matrix, rho, mcharge, energy, xgamma, zeff, sap_int, calculate_forces, just_energy)
...
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.
Some utility functions for the calculation of integrals.
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.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
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.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
General overlap type integrals containers.
subroutine, public release_sap_int(sap_int)
...
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Calculation of Coulomb contributions in xTB.
Definition: xtb_coulomb.F:12
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
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
subroutine, public build_xtb_coulomb(qs_env, ks_matrix, rho, charges, mcharge, energy, calculate_forces, just_energy)
...
Definition: xtb_coulomb.F:103
subroutine, public xtb_dsint_list(qs_env, sap_int)
...
Definition: xtb_coulomb.F:820
Definition of the xTB parameter types.
Definition: xtb_types.F:20
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