LCOV - code coverage report
Current view: top level - src - xtb_coulomb.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 488 498 98.0 %
Date: 2024-03-28 07:31:50 Functions: 4 4 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculation of Coulomb contributions in xTB
      10             : !> \author JGH
      11             : ! **************************************************************************************************
      12             : MODULE 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             : 
      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       21138 :    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       21138 :       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       21138 :       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       21138 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: xgamma, zeffk
     128       21138 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gammab, gcij, gmcharge
     129       21138 :       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       21138 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: dsblock, ksblock, pblock, sblock
     134       21138 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: dsint
     135       21138 :       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       21138 :       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       21138 :          DIMENSION(:), POINTER                           :: nl_iterator
     148             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     149       21138 :          POINTER                                         :: n_list
     150       21138 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     151       21138 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     152       21138 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     153       21138 :       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       21138 :       CALL timeset(routineN, handle)
     159             : 
     160       21138 :       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       21138 :                       dft_control=dft_control)
     169             : 
     170       21138 :       xtb_control => dft_control%qs_control%xtb_control
     171             : 
     172       21138 :       use_virial = .FALSE.
     173       21138 :       IF (calculate_forces) THEN
     174         788 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     175             :       END IF
     176             : 
     177       21138 :       do_gamma_stress = .FALSE.
     178       21138 :       IF (.NOT. just_energy .AND. use_virial) THEN
     179          48 :          IF (dft_control%nimages == 1) do_gamma_stress = .TRUE.
     180             :       END IF
     181             : 
     182       21138 :       IF (atprop%energy) THEN
     183         172 :          CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
     184         172 :          natom = SIZE(particle_set)
     185         172 :          CALL atprop_array_init(atprop%atecoul, natom)
     186             :       END IF
     187             : 
     188       21138 :       IF (calculate_forces) THEN
     189             :          nmat = 4
     190             :       ELSE
     191       20720 :          nmat = 1
     192             :       END IF
     193             : 
     194       21138 :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
     195       84552 :       ALLOCATE (gchrg(natom, 5, nmat))
     196     1228700 :       gchrg = 0._dp
     197       84552 :       ALLOCATE (gmcharge(natom, nmat))
     198      258172 :       gmcharge = 0._dp
     199             : 
     200             :       ! short range contribution (gamma)
     201             :       ! loop over all atom pairs (sab_xtbe)
     202       21138 :       kg = xtb_control%kg
     203       21138 :       NULLIFY (n_list)
     204       21138 :       IF (xtb_control%old_coulomb_damping) THEN
     205           0 :          CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
     206             :       ELSE
     207       21138 :          CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
     208             :       END IF
     209       21138 :       CALL neighbor_list_iterator_create(nl_iterator, n_list)
     210     7405617 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     211             :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     212     7384479 :                                 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
     213     7384479 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     214     7384479 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
     215     7384479 :          IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     216     7384472 :          CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     217     7384472 :          CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
     218     7384472 :          IF (.NOT. defined .OR. natorb_b < 1) CYCLE
     219             :          ! atomic parameters
     220     7384465 :          CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
     221     7384465 :          CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
     222             :          ! gamma matrix
     223     7384465 :          ni = lmaxa + 1
     224     7384465 :          nj = lmaxb + 1
     225    29537860 :          ALLOCATE (gammab(ni, nj))
     226     7384465 :          rcut = rcuta + rcutb
     227    29537860 :          dr = SQRT(SUM(rij(:)**2))
     228     7384465 :          CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
     229    90116413 :          gchrg(iatom, 1:ni, 1) = gchrg(iatom, 1:ni, 1) + MATMUL(gammab, charges(jatom, 1:nj))
     230     7384465 :          IF (iatom /= jatom) THEN
     231    90292506 :             gchrg(jatom, 1:nj, 1) = gchrg(jatom, 1:nj, 1) + MATMUL(charges(iatom, 1:ni), gammab)
     232             :          END IF
     233     7384465 :          IF (calculate_forces) THEN
     234      280771 :             IF (dr > 1.e-6_dp) THEN
     235      278831 :                CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
     236     1115324 :                DO i = 1, 3
     237             :                   gchrg(iatom, 1:ni, i + 1) = gchrg(iatom, 1:ni, i + 1) &
     238    10194414 :                                               + MATMUL(gammab, charges(jatom, 1:nj))*rij(i)/dr
     239     1115324 :                   IF (iatom /= jatom) THEN
     240             :                      gchrg(jatom, 1:nj, i + 1) = gchrg(jatom, 1:nj, i + 1) &
     241    10533915 :                                                  - MATMUL(charges(iatom, 1:ni), gammab)*rij(i)/dr
     242             :                   END IF
     243             :                END DO
     244      278831 :                IF (use_virial) THEN
     245     1161150 :                   gcint(1:ni) = MATMUL(gammab, charges(jatom, 1:nj))
     246      516084 :                   DO i = 1, 3
     247     1092960 :                      fij(i) = -SUM(charges(iatom, 1:ni)*gcint(1:ni))*rij(i)/dr
     248             :                   END DO
     249      129021 :                   fi = 1.0_dp
     250      129021 :                   IF (iatom == jatom) fi = 0.5_dp
     251      129021 :                   CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
     252      129021 :                   IF (atprop%stress) THEN
     253       97654 :                      CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
     254       97654 :                      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    29537881 :          DEALLOCATE (gammab)
     260             :       END DO
     261       21138 :       CALL neighbor_list_iterator_release(nl_iterator)
     262             : 
     263             :       ! 1/R contribution
     264             : 
     265       21138 :       IF (xtb_control%coulomb_lr) THEN
     266       21138 :          do_ewald = xtb_control%do_ewald
     267       21138 :          IF (do_ewald) THEN
     268             :             ! Ewald sum
     269        8570 :             NULLIFY (ewald_env, ewald_pw)
     270             :             CALL get_qs_env(qs_env=qs_env, &
     271        8570 :                             ewald_env=ewald_env, ewald_pw=ewald_pw)
     272        8570 :             CALL get_cell(cell=cell, periodic=periodic, deth=deth)
     273        8570 :             CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
     274        8570 :             CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
     275             :             CALL tb_ewald_overlap(gmcharge, mcharge, alpha, n_list, virial, &
     276        8570 :                                   use_virial, atprop=atprop)
     277           0 :             SELECT CASE (ewald_type)
     278             :             CASE DEFAULT
     279           0 :                CPABORT("Invalid Ewald type")
     280             :             CASE (do_ewald_none)
     281           0 :                CPABORT("Not allowed with DFTB")
     282             :             CASE (do_ewald_ewald)
     283           0 :                CPABORT("Standard Ewald not implemented in DFTB")
     284             :             CASE (do_ewald_pme)
     285           0 :                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        8570 :                                      use_virial, atprop=atprop)
     290             :             END SELECT
     291             :          ELSE
     292             :             ! direct sum
     293             :             CALL get_qs_env(qs_env=qs_env, &
     294       12568 :                             local_particles=local_particles)
     295       48342 :             DO ikind = 1, SIZE(local_particles%n_el)
     296      115343 :                DO ia = 1, local_particles%n_el(ikind)
     297       67001 :                   iatom = local_particles%list(ikind)%array(ia)
     298      832721 :                   DO jatom = 1, iatom - 1
     299     2919784 :                      rij = particle_set(iatom)%r - particle_set(jatom)%r
     300     2919784 :                      rij = pbc(rij, cell)
     301     2919784 :                      dr = SQRT(SUM(rij(:)**2))
     302      796947 :                      IF (dr > 1.e-6_dp) THEN
     303      729946 :                         gmcharge(iatom, 1) = gmcharge(iatom, 1) + mcharge(jatom)/dr
     304      729946 :                         gmcharge(jatom, 1) = gmcharge(jatom, 1) + mcharge(iatom)/dr
     305      735106 :                         DO i = 2, nmat
     306        5160 :                            gmcharge(iatom, i) = gmcharge(iatom, i) + rij(i - 1)*mcharge(jatom)/dr**3
     307      735106 :                            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       12568 :             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       21138 :                       force=force, para_env=para_env)
     321       21138 :       CALL para_env%sum(gmcharge(:, 1))
     322       21138 :       CALL para_env%sum(gchrg(:, :, 1))
     323             : 
     324       21138 :       IF (xtb_control%coulomb_lr) THEN
     325       21138 :          IF (do_ewald) THEN
     326             :             ! add self charge interaction and background charge contribution
     327       81446 :             gmcharge(:, 1) = gmcharge(:, 1) - 2._dp*alpha*oorootpi*mcharge(:)
     328        9506 :             IF (ANY(periodic(:) == 1)) THEN
     329       79886 :                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       21138 :                                atom_of_kind=atom_of_kind)
     338       21138 :       ecsr = 0.0_dp
     339      225058 :       DO iatom = 1, natom
     340      203920 :          ikind = kind_of(iatom)
     341      203920 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     342      203920 :          CALL get_xtb_atom_param(xtb_kind, lmax=ni)
     343      203920 :          ni = ni + 1
     344      540804 :          ecsr = ecsr + SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, 1))
     345             :       END DO
     346             : 
     347       21138 :       energy%hartree = energy%hartree + 0.5_dp*ecsr
     348      225058 :       energy%hartree = energy%hartree + 0.5_dp*SUM(mcharge(:)*gmcharge(:, 1))
     349             : 
     350       21138 :       IF (atprop%energy) THEN
     351         172 :          CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
     352         748 :          DO ikind = 1, SIZE(local_particles%n_el)
     353         576 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     354         576 :             CALL get_xtb_atom_param(xtb_kind, lmax=ni, occupation=occ)
     355         576 :             ni = ni + 1
     356        3456 :             zeff = SUM(REAL(occ, KIND=dp))
     357        4360 :             DO ia = 1, local_particles%n_el(ikind)
     358        3036 :                iatom = local_particles%list(ikind)%array(ia)
     359             :                atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
     360        7258 :                                        0.5_dp*SUM(REAL(occ(1:ni), KIND=dp)*gchrg(iatom, 1:ni, 1))
     361             :                atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
     362        3612 :                                        0.5_dp*zeff*gmcharge(iatom, 1)
     363             :             END DO
     364             :          END DO
     365             :       END IF
     366             : 
     367       21138 :       IF (calculate_forces) THEN
     368        3992 :          DO iatom = 1, natom
     369        3574 :             ikind = kind_of(iatom)
     370        3574 :             atom_i = atom_of_kind(iatom)
     371        3574 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     372        3574 :             CALL get_xtb_atom_param(xtb_kind, lmax=ni)
     373             :             ! short range
     374        3574 :             ni = ni + 1
     375       14296 :             DO i = 1, 3
     376       29914 :                fij(i) = SUM(charges(iatom, 1:ni)*gchrg(iatom, 1:ni, i + 1))
     377             :             END DO
     378        3574 :             force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
     379        3574 :             force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
     380        3574 :             force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
     381             :             ! long range
     382       14296 :             DO i = 1, 3
     383       14296 :                fij(i) = gmcharge(iatom, i + 1)*mcharge(iatom)
     384             :             END DO
     385        3574 :             force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
     386        3574 :             force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
     387        7566 :             force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
     388             :          END DO
     389             :       END IF
     390             : 
     391       21138 :       IF (.NOT. just_energy) THEN
     392       21096 :          CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
     393       21096 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
     394             : 
     395       21096 :          nimg = dft_control%nimages
     396       21096 :          NULLIFY (cell_to_index)
     397       21096 :          IF (nimg > 1) THEN
     398        2644 :             NULLIFY (kpoints)
     399        2644 :             CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
     400        2644 :             CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     401             :          END IF
     402             : 
     403       21096 :          IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
     404         432 :             DO img = 1, nimg
     405             :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     406         432 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     407             :             END DO
     408             :          END IF
     409             : 
     410       21096 :          NULLIFY (sap_int)
     411       21096 :          IF (do_gamma_stress) THEN
     412             :             ! derivative overlap integral (non collapsed)
     413          34 :             CALL xtb_dsint_list(qs_env, sap_int)
     414             :          END IF
     415             : 
     416       21096 :          IF (nimg == 1) THEN
     417             :             ! no k-points; all matrices have been transformed to periodic bsf
     418       18452 :             CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
     419     1394940 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
     420     1376488 :                CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
     421     1376488 :                ikind = kind_of(irow)
     422     1376488 :                jkind = kind_of(icol)
     423             : 
     424             :                ! atomic parameters
     425     1376488 :                CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     426     1376488 :                CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     427     1376488 :                CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
     428     1376488 :                CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
     429             : 
     430     1376488 :                ni = SIZE(sblock, 1)
     431     1376488 :                nj = SIZE(sblock, 2)
     432     5505952 :                ALLOCATE (gcij(ni, nj))
     433     5575307 :                DO i = 1, ni
     434    18551113 :                   DO j = 1, nj
     435    12975806 :                      la = laoa(i) + 1
     436    12975806 :                      lb = laob(j) + 1
     437    17174625 :                      gcij(i, j) = 0.5_dp*(gchrg(irow, la, 1) + gchrg(icol, lb, 1))
     438             :                   END DO
     439             :                END DO
     440     1376488 :                gmij = 0.5_dp*(gmcharge(irow, 1) + gmcharge(icol, 1))
     441     2761352 :                DO is = 1, SIZE(ks_matrix, 1)
     442     1384864 :                   NULLIFY (ksblock)
     443             :                   CALL dbcsr_get_block_p(matrix=ks_matrix(is, 1)%matrix, &
     444     1384864 :                                          row=irow, col=icol, block=ksblock, found=found)
     445     1384864 :                   CPASSERT(found)
     446    35938332 :                   ksblock = ksblock - gcij*sblock
     447    38699684 :                   ksblock = ksblock - gmij*sblock
     448             :                END DO
     449     1376488 :                IF (calculate_forces) THEN
     450       45310 :                   atom_i = atom_of_kind(irow)
     451       45310 :                   atom_j = atom_of_kind(icol)
     452       45310 :                   NULLIFY (pblock)
     453             :                   CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
     454       45310 :                                          row=irow, col=icol, block=pblock, found=found)
     455       45310 :                   CPASSERT(found)
     456      181240 :                   DO i = 1, 3
     457      135930 :                      NULLIFY (dsblock)
     458             :                      CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
     459      135930 :                                             row=irow, col=icol, block=dsblock, found=found)
     460      135930 :                      CPASSERT(found)
     461      135930 :                      fij(i) = 0.0_dp
     462             :                      ! short range
     463     1537413 :                      fi = -2.0_dp*SUM(pblock*dsblock*gcij)
     464      135930 :                      force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     465      135930 :                      force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     466      135930 :                      fij(i) = fij(i) + fi
     467             :                      ! long range
     468     1537413 :                      fi = -2.0_dp*gmij*SUM(pblock*dsblock)
     469      135930 :                      force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     470      135930 :                      force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     471      317170 :                      fij(i) = fij(i) + fi
     472             :                   END DO
     473             :                END IF
     474     4147916 :                DEALLOCATE (gcij)
     475             :             END DO
     476       18452 :             CALL dbcsr_iterator_stop(iter)
     477             :             ! stress tensor (needs recalculation of overlap integrals)
     478       18452 :             IF (do_gamma_stress) THEN
     479         118 :                DO ikind = 1, nkind
     480         358 :                   DO jkind = 1, nkind
     481         240 :                      iac = ikind + nkind*(jkind - 1)
     482         240 :                      IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
     483             :                      ! atomic parameters
     484         138 :                      CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     485         138 :                      CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     486         138 :                      CALL get_xtb_atom_param(xtb_atom_a, lao=laoa, natorb=ni)
     487         138 :                      CALL get_xtb_atom_param(xtb_atom_b, lao=laob, natorb=nj)
     488        1196 :                      DO ia = 1, sap_int(iac)%nalist
     489         974 :                         IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
     490         956 :                         iatom = sap_int(iac)%alist(ia)%aatom
     491       70471 :                         DO ic = 1, sap_int(iac)%alist(ia)%nclist
     492       69275 :                            jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
     493      277100 :                            rij = sap_int(iac)%alist(ia)%clist(ic)%rac
     494      277100 :                            dr = SQRT(SUM(rij(:)**2))
     495       70249 :                            IF (dr > 1.e-6_dp) THEN
     496       68803 :                               dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
     497      275212 :                               ALLOCATE (gcij(ni, nj))
     498      283604 :                               DO i = 1, ni
     499     1177695 :                                  DO j = 1, nj
     500      894091 :                                     la = laoa(i) + 1
     501      894091 :                                     lb = laob(j) + 1
     502     1108892 :                                     gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
     503             :                                  END DO
     504             :                               END DO
     505       68803 :                               gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
     506       68803 :                               icol = MAX(iatom, jatom)
     507       68803 :                               irow = MIN(iatom, jatom)
     508       68803 :                               NULLIFY (pblock)
     509             :                               CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
     510       68803 :                                                      row=irow, col=icol, block=pblock, found=found)
     511       68803 :                               CPASSERT(found)
     512       68803 :                               fij = 0.0_dp
     513      275212 :                               DO i = 1, 3
     514             :                                  ! short/long range
     515      206409 :                                  IF (irow == iatom) THEN
     516     1956816 :                                     f1 = -2.0_dp*SUM(pblock*dsint(:, :, i)*gcij)
     517     1956816 :                                     f2 = -2.0_dp*gmij*SUM(pblock*dsint(:, :, i))
     518             :                                  ELSE
     519     1575975 :                                     f1 = -2.0_dp*SUM(TRANSPOSE(pblock)*dsint(:, :, i)*gcij)
     520     1575975 :                                     f2 = -2.0_dp*gmij*SUM(TRANSPOSE(pblock)*dsint(:, :, i))
     521             :                                  END IF
     522      275212 :                                  fij(i) = f1 + f2
     523             :                               END DO
     524       68803 :                               DEALLOCATE (gcij)
     525       68803 :                               fi = 1.0_dp
     526       68803 :                               IF (iatom == jatom) fi = 0.5_dp
     527       68803 :                               CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
     528      137606 :                               IF (atprop%stress) THEN
     529       51271 :                                  CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
     530       51271 :                                  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        2644 :             NULLIFY (n_list)
     540        2644 :             CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
     541        2644 :             CALL neighbor_list_iterator_create(nl_iterator, n_list)
     542     1959732 :             DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     543             :                CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     544     1957088 :                                       iatom=iatom, jatom=jatom, r=rij, cell=cellind)
     545             : 
     546     1957088 :                icol = MAX(iatom, jatom)
     547     1957088 :                irow = MIN(iatom, jatom)
     548             : 
     549     1957088 :                ic = cell_to_index(cellind(1), cellind(2), cellind(3))
     550     1957088 :                CPASSERT(ic > 0)
     551             : 
     552     1957088 :                NULLIFY (sblock)
     553             :                CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
     554     1957088 :                                       row=irow, col=icol, block=sblock, found=found)
     555     1957088 :                CPASSERT(found)
     556             : 
     557             :                ! atomic parameters
     558     1957088 :                CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     559     1957088 :                CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     560     1957088 :                CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
     561     1957088 :                CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
     562             : 
     563     1957088 :                ni = SIZE(sblock, 1)
     564     1957088 :                nj = SIZE(sblock, 2)
     565     7828352 :                ALLOCATE (gcij(ni, nj))
     566    11425467 :                DO i = 1, ni
     567    72875334 :                   DO j = 1, nj
     568    70918246 :                      IF (irow == iatom) THEN
     569    34938761 :                         la = laoa(i) + 1
     570    34938761 :                         lb = laob(j) + 1
     571    34938761 :                         gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
     572             :                      ELSE
     573    26511106 :                         la = laoa(j) + 1
     574    26511106 :                         lb = laob(i) + 1
     575    26511106 :                         gcij(i, j) = 0.5_dp*(gchrg(iatom, la, 1) + gchrg(jatom, lb, 1))
     576             :                      END IF
     577             :                   END DO
     578             :                END DO
     579     1957088 :                gmij = 0.5_dp*(gmcharge(iatom, 1) + gmcharge(jatom, 1))
     580     4098426 :                DO is = 1, SIZE(ks_matrix, 1)
     581     2141338 :                   NULLIFY (ksblock)
     582             :                   CALL dbcsr_get_block_p(matrix=ks_matrix(is, ic)%matrix, &
     583     2141338 :                                          row=irow, col=icol, block=ksblock, found=found)
     584     2141338 :                   CPASSERT(found)
     585   146953790 :                   ksblock = ksblock - gcij*sblock
     586   151052216 :                   ksblock = ksblock - gmij*sblock
     587             :                END DO
     588             : 
     589     1957088 :                IF (calculate_forces) THEN
     590       29400 :                   atom_i = atom_of_kind(iatom)
     591       29400 :                   atom_j = atom_of_kind(jatom)
     592       29400 :                   IF (irow /= iatom) THEN
     593       11957 :                      gmij = -gmij
     594      759507 :                      gcij = -gcij
     595             :                   END IF
     596       29400 :                   NULLIFY (pblock)
     597             :                   CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
     598       29400 :                                          row=irow, col=icol, block=pblock, found=found)
     599       29400 :                   CPASSERT(found)
     600      117600 :                   DO i = 1, 3
     601       88200 :                      NULLIFY (dsblock)
     602             :                      CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
     603       88200 :                                             row=irow, col=icol, block=dsblock, found=found)
     604       88200 :                      CPASSERT(found)
     605       88200 :                      fij(i) = 0.0_dp
     606             :                      ! short range
     607     5846208 :                      fi = -2.0_dp*SUM(pblock*dsblock*gcij)
     608       88200 :                      force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     609       88200 :                      force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     610       88200 :                      fij(i) = fij(i) + fi
     611             :                      ! long range
     612     5846208 :                      fi = -2.0_dp*gmij*SUM(pblock*dsblock)
     613       88200 :                      force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     614       88200 :                      force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     615      205800 :                      fij(i) = fij(i) + fi
     616             :                   END DO
     617       29400 :                   IF (use_virial) THEN
     618       73828 :                      dr = SQRT(SUM(rij(:)**2))
     619       18457 :                      IF (dr > 1.e-6_dp) THEN
     620       18393 :                         fi = 1.0_dp
     621       18393 :                         IF (iatom == jatom) fi = 0.5_dp
     622       18393 :                         CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
     623       18393 :                         IF (atprop%stress) THEN
     624           0 :                            CALL virial_pair_force(atprop%atstress(:, :, iatom), fi*0.5_dp, fij, rij)
     625           0 :                            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     5873908 :                DEALLOCATE (gcij)
     631             : 
     632             :             END DO
     633        2644 :             CALL neighbor_list_iterator_release(nl_iterator)
     634             :          END IF
     635             : 
     636       21096 :          IF (calculate_forces .AND. SIZE(matrix_p, 1) == 2) THEN
     637         432 :             DO img = 1, nimg
     638             :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     639         432 :                               alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
     640             :             END DO
     641             :          END IF
     642             :       END IF
     643             : 
     644       21138 :       IF (xtb_control%tb3_interaction) THEN
     645       21138 :          CALL get_qs_env(qs_env, nkind=nkind)
     646      105690 :          ALLOCATE (zeffk(nkind), xgamma(nkind))
     647       76960 :          DO ikind = 1, nkind
     648       55822 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     649       76960 :             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       21138 :                                    sap_int, calculate_forces, just_energy)
     654       21138 :          DEALLOCATE (zeffk, xgamma)
     655             :       END IF
     656             : 
     657             :       ! QMMM
     658       21138 :       IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
     659             :          CALL build_tb_coulomb_qmqm(qs_env, ks_matrix, rho, mcharge, energy, &
     660         862 :                                     calculate_forces, just_energy)
     661             :       END IF
     662             : 
     663       21138 :       IF (do_gamma_stress) THEN
     664          34 :          CALL release_sap_int(sap_int)
     665             :       END IF
     666             : 
     667       21138 :       CALL timestop(handle)
     668             : 
     669       42276 :    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     7606992 :    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     7606992 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: eta
     707             : 
     708    30427968 :       ALLOCATE (eta(nla, nlb))
     709    36778073 :       eta = 0.0_dp
     710             : 
     711    18801744 :       DO j = 1, nlb
     712    36778073 :          DO i = 1, nla
     713    17976329 :             eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
     714    29171081 :             eta(i, j) = 2._dp/eta(i, j)
     715             :          END DO
     716             :       END DO
     717             : 
     718    36778073 :       gmat = 0.0_dp
     719     7606992 :       IF (rab < 1.e-6_dp) THEN
     720             :          ! on site terms
     721      546996 :          gmat(:, :) = eta(:, :)
     722     7502466 :       ELSEIF (rab > rcut) THEN
     723             :          ! do nothing
     724             :       ELSE
     725     7502466 :          rk = rab**kg
     726    36231077 :          eta = eta**(-kg)
     727     7502466 :          IF (rab < rcut - rsmooth) THEN
     728             :             fcut = 1.0_dp
     729             :          ELSE
     730      832596 :             r = rab - (rcut - rsmooth)
     731      832596 :             x = r/rsmooth
     732      832596 :             fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
     733             :          END IF
     734    36231077 :          gmat(:, :) = fcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg) - fcut/rab
     735             :       END IF
     736             : 
     737     7606992 :       DEALLOCATE (eta)
     738             : 
     739     7606992 :    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      303408 :    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      303408 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: eta
     777             : 
     778     1213632 :       ALLOCATE (eta(nla, nlb))
     779             : 
     780      742804 :       DO j = 1, nlb
     781     1452211 :          DO i = 1, nla
     782      709407 :             eta(i, j) = 1._dp/(etaa*(1._dp + kappaa(i))) + 1._dp/(etab*(1._dp + kappab(j)))
     783     1148803 :             eta(i, j) = 2._dp/eta(i, j)
     784             :          END DO
     785             :       END DO
     786             : 
     787      303408 :       IF (rab < 1.e-6) THEN
     788             :          ! on site terms
     789           0 :          dgmat(:, :) = 0.0_dp
     790      303408 :       ELSEIF (rab > rcut) THEN
     791           0 :          dgmat(:, :) = 0.0_dp
     792             :       ELSE
     793     1452211 :          eta = eta**(-kg)
     794      303408 :          rk = rab**kg
     795      303408 :          IF (rab < rcut - rsmooth) THEN
     796             :             fcut = 1.0_dp
     797             :             dfcut = 0.0_dp
     798             :          ELSE
     799       40640 :             r = rab - (rcut - rsmooth)
     800       40640 :             x = r/rsmooth
     801       40640 :             fcut = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
     802       40640 :             dfcut = -30._dp*x**4 + 60._dp*x**3 - 30._dp*x**2
     803       40640 :             dfcut = dfcut/rsmooth
     804             :          END IF
     805     1452211 :          dgmat(:, :) = dfcut*(1._dp/(rk + eta(:, :)))**(1._dp/kg)
     806     1452211 :          dgmat(:, :) = dgmat(:, :) - dfcut/rab + fcut/rab**2
     807     1452211 :          dgmat(:, :) = dgmat(:, :) - fcut/(rk + eta(:, :))*(1._dp/(rk + eta(:, :)))**(1._dp/kg)*rk/rab
     808             :       END IF
     809             : 
     810      303408 :       DEALLOCATE (eta)
     811             : 
     812      303408 :    END SUBROUTINE dgamma_rab_sr
     813             : 
     814             : ! **************************************************************************************************
     815             : !> \brief ...
     816             : !> \param qs_env ...
     817             : !> \param sap_int ...
     818             : ! **************************************************************************************************
     819          38 :    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          38 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     830          38 :                                                             npgfb, nsgfa, nsgfb
     831          38 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     832             :       LOGICAL                                            :: defined
     833             :       REAL(KIND=dp)                                      :: dr
     834          38 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: owork
     835          38 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: oint, sint
     836             :       REAL(KIND=dp), DIMENSION(3)                        :: rij
     837          38 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     838          38 :       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          38 :       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          38 :          DIMENSION(:), POINTER                           :: nl_iterator
     845             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     846          38 :          POINTER                                         :: sab_orb
     847          38 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     848             :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b
     849             : 
     850          38 :       CALL timeset(routineN, handle)
     851             : 
     852          38 :       CALL get_qs_env(qs_env=qs_env, nkind=nkind)
     853          38 :       CPASSERT(.NOT. ASSOCIATED(sap_int))
     854         390 :       ALLOCATE (sap_int(nkind*nkind))
     855         314 :       DO i = 1, nkind*nkind
     856         276 :          NULLIFY (sap_int(i)%alist, sap_int(i)%asort, sap_int(i)%aindex)
     857         314 :          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          38 :                       sab_orb=sab_orb)
     864             : 
     865             :       ! set up basis set lists
     866         210 :       ALLOCATE (basis_set_list(nkind))
     867          38 :       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          38 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     871       70653 :       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       70615 :                                 inode=jneighbor, cell=cell, r=rij)
     875       70615 :          iac = ikind + nkind*(jkind - 1)
     876             :          !
     877       70615 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     878       70615 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
     879       70615 :          IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     880       70615 :          CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     881       70615 :          CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
     882       70615 :          IF (.NOT. defined .OR. natorb_b < 1) CYCLE
     883             : 
     884      282460 :          dr = SQRT(SUM(rij(:)**2))
     885             : 
     886             :          ! integral list
     887       70615 :          IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) THEN
     888         156 :             sap_int(iac)%a_kind = ikind
     889         156 :             sap_int(iac)%p_kind = jkind
     890         156 :             sap_int(iac)%nalist = nlist
     891        1460 :             ALLOCATE (sap_int(iac)%alist(nlist))
     892        1148 :             DO i = 1, nlist
     893         992 :                NULLIFY (sap_int(iac)%alist(i)%clist)
     894         992 :                sap_int(iac)%alist(i)%aatom = 0
     895        1148 :                sap_int(iac)%alist(i)%nclist = 0
     896             :             END DO
     897             :          END IF
     898       70615 :          IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ilist)%clist)) THEN
     899         974 :             sap_int(iac)%alist(ilist)%aatom = iatom
     900         974 :             sap_int(iac)%alist(ilist)%nclist = nneighbor
     901       79381 :             ALLOCATE (sap_int(iac)%alist(ilist)%clist(nneighbor))
     902       71589 :             DO i = 1, nneighbor
     903       71589 :                sap_int(iac)%alist(ilist)%clist(i)%catom = 0
     904             :             END DO
     905             :          END IF
     906       70615 :          clist => sap_int(iac)%alist(ilist)%clist(jneighbor)
     907       70615 :          clist%catom = jatom
     908      282460 :          clist%cell = cell
     909      282460 :          clist%rac = rij
     910      353075 :          ALLOCATE (clist%acint(natorb_a, natorb_b, 3))
     911       70615 :          NULLIFY (clist%achint)
     912     3673726 :          clist%acint = 0._dp
     913       70615 :          clist%nsgf_cnt = 0
     914       70615 :          NULLIFY (clist%sgf_list)
     915             : 
     916             :          ! overlap
     917       70615 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     918       70615 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     919       70615 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     920       70615 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     921             :          ! basis ikind
     922       70615 :          first_sgfa => basis_set_a%first_sgf
     923       70615 :          la_max => basis_set_a%lmax
     924       70615 :          la_min => basis_set_a%lmin
     925       70615 :          npgfa => basis_set_a%npgf
     926       70615 :          nseta = basis_set_a%nset
     927       70615 :          nsgfa => basis_set_a%nsgf_set
     928       70615 :          rpgfa => basis_set_a%pgf_radius
     929       70615 :          set_radius_a => basis_set_a%set_radius
     930       70615 :          scon_a => basis_set_a%scon
     931       70615 :          zeta => basis_set_a%zet
     932             :          ! basis jkind
     933       70615 :          first_sgfb => basis_set_b%first_sgf
     934       70615 :          lb_max => basis_set_b%lmax
     935       70615 :          lb_min => basis_set_b%lmin
     936       70615 :          npgfb => basis_set_b%npgf
     937       70615 :          nsetb = basis_set_b%nset
     938       70615 :          nsgfb => basis_set_b%nsgf_set
     939       70615 :          rpgfb => basis_set_b%pgf_radius
     940       70615 :          set_radius_b => basis_set_b%set_radius
     941       70615 :          scon_b => basis_set_b%scon
     942       70615 :          zetb => basis_set_b%zet
     943             : 
     944       70615 :          ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
     945      564920 :          ALLOCATE (oint(ldsab, ldsab, 4), owork(ldsab, ldsab))
     946      353075 :          ALLOCATE (sint(natorb_a, natorb_b, 4))
     947     4874763 :          sint = 0.0_dp
     948             : 
     949      218232 :          DO iset = 1, nseta
     950      147617 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     951      147617 :             n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     952      147617 :             sgfa = first_sgfa(1, iset)
     953      531797 :             DO jset = 1, nsetb
     954      313565 :                IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
     955      257686 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     956      257686 :                n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     957      257686 :                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      257686 :                                rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
     961             :                ! Contraction
     962     1436047 :                DO i = 1, 4
     963             :                   CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
     964     1030744 :                                    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     1344309 :                                  sgfa, sgfb, trans=.FALSE.)
     967             :                END DO
     968             :             END DO
     969             :          END DO
     970             :          ! update dS/dR matrix
     971     3673726 :          clist%acint(1:natorb_a, 1:natorb_b, 1:3) = sint(1:natorb_a, 1:natorb_b, 2:4)
     972             : 
     973      211883 :          DEALLOCATE (oint, owork, sint)
     974             : 
     975             :       END DO
     976          38 :       CALL neighbor_list_iterator_release(nl_iterator)
     977             : 
     978          38 :       DEALLOCATE (basis_set_list)
     979             : 
     980          38 :       CALL timestop(handle)
     981             : 
     982          76 :    END SUBROUTINE xtb_dsint_list
     983             : 
     984    15350076 : END MODULE xtb_coulomb
     985             : 

Generated by: LCOV version 1.15