qs_dispersion_d4.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Calculation of dispersion using pair potentials
!> \author Johann Pototschnig
! **************************************************************************************************
MODULE qs_dispersion_d4
   USE atomic_kind_types, ONLY: atomic_kind_type, &
                                get_atomic_kind, &
                                get_atomic_kind_set
   USE distribution_1d_types, ONLY: distribution_1d_type
   USE eeq_method, ONLY: eeq_charges, eeq_forces
   USE machine, ONLY: m_flush, &
                      m_walltime
   USE cell_types, ONLY: cell_type, &
                         plane_distance, &
                         pbc, &
                         get_cell
   USE qs_environment_types, ONLY: get_qs_env, &
                                   qs_environment_type
   USE qs_force_types, ONLY: qs_force_type
   USE qs_kind_types, ONLY: get_qs_kind, &
                            qs_kind_type, &
                            set_qs_kind
   USE qs_neighbor_list_types, ONLY: get_iterator_info, &
                                     neighbor_list_iterate, &
                                     neighbor_list_iterator_create, &
                                     neighbor_list_iterator_p_type, &
                                     neighbor_list_iterator_release, &
                                     neighbor_list_set_p_type
   USE virial_methods, ONLY: virial_pair_force
   USE virial_types, ONLY: virial_type
   USE kinds, ONLY: dp
   USE particle_types, ONLY: particle_type
   USE qs_dispersion_types, ONLY: qs_dispersion_type
   USE qs_dispersion_utils, ONLY: cellhash
   USE qs_dispersion_cnum, ONLY: cnumber_init, dcnum_type, cnumber_release
   USE message_passing, ONLY: mp_para_env_type
 
#if defined(__DFTD4)
!&<
   USE dftd4,                           ONLY: d4_model, &
                                              damping_param, &
                                              get_dispersion, &
                                              get_rational_damping, &
                                              new, &
                                              new_d4_model, &
                                              realspace_cutoff, &
                                              structure_type, &
                                              rational_damping_param, &
                                              get_coordination_number, &
                                              get_lattice_points
   USE multicharge,                     ONLY: get_charges
   USE mctc_env,                        ONLY: error_type
!&>
#endif
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d4'
 
   PUBLIC :: calculate_dispersion_d4_pairpot
 
! **************************************************************************************************
 
CONTAINS
 
#if defined(__DFTD4)
! **************************************************************************************************
!> \brief Read smooth-cutoff width override from environment.
!> \param primary primary environment variable name
!> \param fallback fallback environment variable name
!> \param cutoff real-space cutoff used to validate the width
!> \param width smooth-cutoff width to update
! **************************************************************************************************
   SUBROUTINE d4_smooth_width_from_env(primary, fallback, cutoff, width)
      CHARACTER(LEN=*), INTENT(IN)                       :: primary, fallback
      REAL(KIND=dp), INTENT(IN)                          :: cutoff
      REAL(KIND=dp), INTENT(INOUT)                       :: width
 
      CHARACTER(LEN=64)                                  :: env
      INTEGER                                            :: io, stat
      REAL(KIND=dp)                                      :: env_width
 
      CALL get_environment_variable(primary, env, status=stat)
      IF (stat /= 0 .OR. len_trim(env) == 0) THEN
         CALL get_environment_variable(fallback, env, status=stat)
      END IF
      IF (stat == 0 .AND. len_trim(env) > 0) THEN
         READ (env, *, iostat=io) env_width
         IF (io == 0) THEN
            IF (env_width > 0.0_dp .AND. env_width < cutoff) THEN
               width = env_width
            ELSE
               width = 0.0_dp
            END IF
         END IF
      END IF
 
   END SUBROUTINE d4_smooth_width_from_env
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param dispersion_env ...
!> \param evdw ...
!> \param calculate_forces ...
!> \param iw ...
!> \param atomic_energy ...
! **************************************************************************************************
   SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw, &
                                              atomic_energy)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(qs_dispersion_type), INTENT(IN), POINTER      :: dispersion_env
      REAL(KIND=dp), INTENT(INOUT)                       :: evdw
      LOGICAL, INTENT(IN)                                :: calculate_forces
      INTEGER, INTENT(IN)                                :: iw
      REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atomic_energy
 
      CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d4_pairpot'
 
      INTEGER                                            :: atoma, cnfun, enshift, handle, i, iatom, &
                                                            ifull, ikind, mref, natom, natom_full, &
                                                            ncoup, nghost
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_map, atom_map_back, atom_of_kind, &
                                                            atomtype, kind_of, species, t_atomtype
      INTEGER, DIMENSION(3)                              :: periodic
      LOGICAL                                            :: debug, grad, ifloating, ighost, &
                                                            use_virial
      LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: a_ghost, exclude_ghost
      LOGICAL, DIMENSION(3)                              :: lperiod
      REAL(KIND=dp)                                      :: ed2, ed3, ev1, ev2, ev3, ev4, pd2, pd3, &
                                                            ta, tb, tc, td, te, ts
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cn, cn_red, cnd, dedcn, dedq, edcn, edq, &
                                                            enerd2, enerd3, energies, energies3, &
                                                            q_red, qd
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ga, gradient, t_xyz, tvec, xyz
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gdeb, gwdcn, gwdq, gwvec
      REAL(KIND=dp), DIMENSION(3, 3)                     :: sigma, stress
      REAL(KIND=dp), DIMENSION(3, 3, 4)                  :: sdeb
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:)        :: dcnum
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(virial_type), POINTER                         :: virial
 
      CLASS(damping_param), ALLOCATABLE                  :: param
      TYPE(d4_model)                                     :: disp
      TYPE(structure_type)                               :: mol
      TYPE(realspace_cutoff)                             :: cutoff
 
      TYPE(error_type), ALLOCATABLE                      :: error
 
      CALL timeset(routinen, handle)
 
      debug = dispersion_env%d4_debug
 
      CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
                      cell=cell, force=force, virial=virial, para_env=para_env)
      CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
 
      !get information about particles
      natom_full = SIZE(particle_set)
      nghost = 0
      ALLOCATE (t_xyz(3, natom_full), t_atomtype(natom_full), a_ghost(natom_full))
      CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
      DO iatom = 1, natom_full
         t_xyz(:, iatom) = particle_set(iatom)%r(:)
         ikind = kind_of(iatom)
         CALL get_qs_kind(qs_kind_set(ikind), zatom=t_atomtype(iatom), ghost=ighost, floating=ifloating)
         a_ghost(iatom) = ighost .OR. ifloating
         IF (a_ghost(iatom)) nghost = nghost + 1
      END DO
 
      natom = natom_full - nghost
      ! Build atom mapping: full index -> reduced index (0 for ghost)
      ALLOCATE (atom_map(natom_full), atom_map_back(natom))
      atom_map = 0
      iatom = 0
      ALLOCATE (xyz(3, natom), atomtype(natom))
      DO i = 1, natom_full
         IF (.NOT. a_ghost(i)) THEN
            iatom = iatom + 1
            atom_map(i) = iatom
            atom_map_back(iatom) = i
            xyz(:, iatom) = t_xyz(:, i)
            atomtype(iatom) = t_atomtype(i)
         END IF
      END DO
      DEALLOCATE (a_ghost, t_xyz, t_atomtype)
 
      !get information about cell / lattice
      CALL get_cell(cell=cell, periodic=periodic)
      lperiod(1) = periodic(1) == 1
      lperiod(2) = periodic(2) == 1
      lperiod(3) = periodic(3) == 1
      ! enforce en shift method 1 (original/molecular)
      ! method 2 from paper on PBC seems not to work
      enshift = 1
      !IF (ALL(periodic == 0)) enshift = 1
 
      !prepare for the call to the dispersion function
      CALL new(mol, atomtype, xyz, lattice=cell%hmat, periodic=lperiod)
      CALL new_d4_model(error, disp, mol)
      IF (ALLOCATED(error)) THEN
         cpabort(error%message)
      END IF
 
      ! Number of coupling
      ncoup = disp%ncoup
 
      ! Build species mapping: full atom index -> D4 species ID (0 for ghost)
      ALLOCATE (species(natom_full))
      species = 0
      DO i = 1, natom
         species(atom_map_back(i)) = mol%id(i)
      END DO
 
      ! Build per-kind exclusion mask for EEQ (ghost/floating kinds excluded)
      ALLOCATE (exclude_ghost(SIZE(qs_kind_set)))
      exclude_ghost = .false.
      DO i = 1, SIZE(qs_kind_set)
         CALL get_qs_kind(qs_kind_set(i), ghost=ighost, floating=ifloating)
         exclude_ghost(i) = ighost .OR. ifloating
      END DO
 
      IF (dispersion_env%ref_functional == "none") THEN
         CALL get_rational_damping("pbe", param, s9=0.0_dp)
         IF (.NOT. ALLOCATED(param)) THEN
            cpabort("D4: Failed to get rational damping parameters for default functional")
         END IF
         SELECT TYPE (param)
         TYPE is (rational_damping_param)
            param%s6 = dispersion_env%s6
            param%s8 = dispersion_env%s8
            param%a1 = dispersion_env%a1
            param%a2 = dispersion_env%a2
            param%alp = dispersion_env%alp
         END SELECT
      ELSE
         CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
         IF (.NOT. ALLOCATED(param)) THEN
            cpabort("D4: Unknown reference functional '"//trim(dispersion_env%ref_functional)//"'")
         END IF
         SELECT TYPE (param)
         TYPE is (rational_damping_param)
            dispersion_env%s6 = param%s6
            dispersion_env%s8 = param%s8
            dispersion_env%a1 = param%a1
            dispersion_env%a2 = param%a2
            dispersion_env%alp = param%alp
         END SELECT
      END IF
 
      ! Coordination number cutoff
      cutoff%cn = dispersion_env%rc_cn
      ! Two-body interaction cutoff
      cutoff%disp2 = dispersion_env%rc_d4*2._dp
      ! Three-body interaction cutoff
      cutoff%disp3 = dispersion_env%rc_disp*2._dp
      IF (cutoff%disp3 > cutoff%disp2) THEN
         cpabort("D4: Three-body cutoff should be smaller than two-body cutoff")
      END IF
      ! DFTD4 4.2 introduced explicit smooth-cutoff widths.
      ! Keep the same narrow two-body smoothing used by the toolchain patch.
      cutoff%width2 = 0.05_dp
      cutoff%width3 = 0.0_dp
      CALL d4_smooth_width_from_env("DFTD4_DISP2_SMOOTH_WIDTH", "TBLITE_D4_DISP2_SMOOTH_WIDTH", &
                                    cutoff%disp2, cutoff%width2)
      CALL d4_smooth_width_from_env("DFTD4_DISP3_SMOOTH_WIDTH", "TBLITE_D4_DISP3_SMOOTH_WIDTH", &
                                    cutoff%disp3, cutoff%width3)
 
      IF (calculate_forces) THEN
         grad = .true.
         use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
      ELSE
         grad = .false.
         use_virial = .false.
      END IF
 
      IF (dispersion_env%d4_reference_code) THEN
 
         !> Wrapper to handle the evaluation of dispersion energy and derivatives
         IF (.NOT. dispersion_env%doabc) THEN
            cpwarn("Using D4_REFERENCE_CODE enforces calculation of C9 term.")
         END IF
         IF (grad) THEN
            ALLOCATE (gradient(3, natom))
            CALL get_dispersion(mol, disp, param, cutoff, evdw, gradient, stress)
            IF (calculate_forces) THEN
               IF (use_virial) THEN
                  virial%pv_virial = virial%pv_virial - stress/para_env%num_pe
               END IF
               DO iatom = 1, natom
                  ifull = atom_map_back(iatom)
                  ikind = kind_of(ifull)
                  atoma = atom_of_kind(ifull)
                  force(ikind)%dispersion(:, atoma) = &
                     force(ikind)%dispersion(:, atoma) + gradient(:, iatom)/para_env%num_pe
               END DO
            END IF
            DEALLOCATE (gradient)
         ELSE
            CALL get_dispersion(mol, disp, param, cutoff, evdw)
         END IF
         !dispersion energy is computed by every MPI process
         evdw = evdw/para_env%num_pe
         IF (dispersion_env%ext_charges) dispersion_env%dcharges = 0.0_dp
         IF (PRESENT(atomic_energy)) THEN
            cpwarn("Atomic energies not available for D4 reference code")
            atomic_energy = 0.0_dp
         END IF
 
      ELSE
 
         IF (iw > 0) THEN
            WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
            WRITE (iw, fmt="(T32,A)") "DEBUG D4 DISPERSION"
            WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
            WRITE (iw, '(A,T71,A10)') " DEBUG D4| Reference functional   ", trim(dispersion_env%ref_functional)
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s6) ", dispersion_env%s6
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s8) ", dispersion_env%s8
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a1) ", dispersion_env%a1
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a2) ", dispersion_env%a2
            WRITE (iw, '(A,T71,E10.4)') " DEBUG D4| Cutoff value coordination numbers ", dispersion_env%eps_cn
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius coordination numbers ", dispersion_env%rc_cn
            WRITE (iw, '(A,T71,I10)') " DEBUG D4| Coordination number function type ", dispersion_env%cnfun
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 2-body terms [bohr]", 2._dp*dispersion_env%rc_d4
            WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 3-body terms [bohr]", 2._dp*dispersion_env%rc_disp
         END IF
 
         td = 0.0_dp
         IF (debug .AND. iw > 0) THEN
            ts = m_walltime()
            CALL refd4_debug(param, disp, mol, cutoff, grad, dispersion_env%doabc, &
                             enerd2, enerd3, cnd, qd, edcn, edq, gdeb, sdeb)
            te = m_walltime()
            td = te - ts
         END IF
 
         tc = 0.0_dp
         ts = m_walltime()
 
         mref = maxval(disp%ref)
         ! Coordination numbers (full-size from qs_env; ghosts excluded internally)
         cnfun = dispersion_env%cnfun
         CALL cnumber_init(qs_env, cn, dcnum, cnfun, grad)
         ! cn has size natom_full; ghost entries are 0
 
         ! Filter CN to reduced space for D4 model
         ALLOCATE (cn_red(natom))
         DO i = 1, natom
            cn_red(i) = cn(atom_map_back(i))
         END DO
         IF (debug .AND. iw > 0) THEN
            WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (max)", maxval(abs(cn_red - cnd))
            WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (ave)", sum(abs(cn_red - cnd))/natom
         END IF
 
         ! EEQ charges
         ! Use CP2K's MPI-parallel EEQ solver with ghost exclusion.
         ! Ghost kinds get huge hardness + zero coupling -> q_ghost = 0 exactly,
         ! without influencing real atom charges. Preserves full MPI parallelism.
         IF (dispersion_env%ext_charges) THEN
            ALLOCATE (q_red(natom))
            q_red(1:natom) = dispersion_env%charges(1:natom)
         ELSE
            ! eeq_charges writes natom_full entries (from qs_env), so allocate full-size
            ALLOCATE (q_red(natom_full))
            CALL eeq_charges(qs_env, q_red, dispersion_env%eeq_sparam, 2, enshift, &
                             exclude=exclude_ghost, cn_max=8.0_dp)
            ! Filter to reduced size for D4 model
            block
               REAL(KIND=dp), ALLOCATABLE :: q_tmp(:)
               ALLOCATE (q_tmp(natom))
               DO i = 1, natom
                  q_tmp(i) = q_red(atom_map_back(i))
               END DO
               DEALLOCATE (q_red)
               CALL move_alloc(q_tmp, q_red)
            END block
         END IF
         IF (debug .AND. iw > 0) THEN
            WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (max)", maxval(abs(q_red - qd))
            WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (ave)", sum(abs(q_red - qd))/natom
         END IF
         ! Weights for C6 calculation (reduced space)
         ALLOCATE (gwvec(mref, natom, ncoup))
         IF (grad) ALLOCATE (gwdcn(mref, natom, ncoup), gwdq(mref, natom, ncoup))
         CALL disp%weight_references(mol, cn_red, q_red, gwvec, gwdcn, gwdq)
 
         ! Energies and derivatives (full-size for CP2K infrastructure compatibility)
         ALLOCATE (energies(natom_full))
         energies(:) = 0.0_dp
         IF (grad) THEN
            ALLOCATE (gradient(3, natom_full), ga(3, natom_full))
            ALLOCATE (dedcn(natom_full), dedq(natom_full))
            dedcn(:) = 0.0_dp; dedq(:) = 0.0_dp
            ga(:, :) = 0.0_dp
            sigma(:, :) = 0.0_dp
         END IF
         CALL dispersion_2b(dispersion_env, cutoff%disp2, disp%r4r2, &
                            gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
                            energies, dedcn, dedq, grad, ga, sigma, &
                            atom_map, species)
         IF (grad) THEN
            gradient(1:3, 1:natom_full) = ga(1:3, 1:natom_full)
            stress = sigma
            IF (debug) THEN
               CALL para_env%sum(ga)
               CALL para_env%sum(sigma)
               IF (iw > 0) THEN
                  CALL gerror(ga, gdeb(:, :, 1), ev1, ev2, ev3, ev4)
                  WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [2B]", ev1, ev2, " %"
                  WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [2B]", ev3, ev4, " %"
                  IF (use_virial) THEN
                     CALL serror(sigma, sdeb(:, :, 1), ev1, ev2)
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [2B]", ev1, ev2, " %"
                  END IF
               END IF
            END IF
         END IF
         ! no contribution from dispersion_3b as q=0 (but q is changed!)
         ! so we calculate this here
         IF (grad) THEN
            IF (dispersion_env%ext_charges) THEN
               dispersion_env%dcharges = dedq
            ELSE
               CALL para_env%sum(dedq)
               ! Reconstruct full-sized charges for eeq_forces
               block
                  REAL(KIND=dp), ALLOCATABLE :: q_full(:)
                  ALLOCATE (q_full(natom_full))
                  q_full = 0.0_dp
                  DO i = 1, natom
                     q_full(atom_map_back(i)) = q_red(i)
                  END DO
                  ga(:, :) = 0.0_dp
                  sigma = 0.0_dp
                  CALL eeq_forces(qs_env, q_full, dedq, ga, sigma, dispersion_env%eeq_sparam, &
                                  2, enshift, response_only=.true., exclude=exclude_ghost, &
                                  cn_max=8.0_dp)
                  DEALLOCATE (q_full)
               END block
               gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
               stress = stress + sigma
               IF (debug) THEN
                  CALL para_env%sum(ga)
                  CALL para_env%sum(sigma)
                  IF (iw > 0) THEN
                     CALL verror(dedq, edq, ev1, ev2)
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdq", ev1, ev2, " %"
                     CALL gerror(ga, gdeb(:, :, 2), ev1, ev2, ev3, ev4)
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdq]", ev1, ev2, " %"
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdq]", ev3, ev4, " %"
                     IF (use_virial) THEN
                        CALL serror(sigma, sdeb(:, :, 2), ev1, ev2)
                        WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdq]", ev1, ev2, " %"
                     END IF
                  END IF
               END IF
            END IF
         END IF
 
         IF (dispersion_env%doabc) THEN
            ALLOCATE (energies3(natom_full))
            energies3(:) = 0.0_dp
            q_red(:) = 0.0_dp
            ! i.e. dc6dq = dEdq = 0
            CALL disp%weight_references(mol, cn_red, q_red, gwvec, gwdcn, gwdq)
            !
            IF (grad) THEN
               gwdq = 0.0_dp
               ga(:, :) = 0.0_dp
               sigma = 0.0_dp
            END IF
            CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, tvec)
            CALL dispersion_3b(qs_env, dispersion_env, tvec, cutoff%disp3, disp%r4r2, &
                               gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
                               energies3, dedcn, dedq, grad, ga, sigma, &
                               atom_map, species)
            IF (grad) THEN
               gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
               stress = stress + sigma
               IF (debug) THEN
                  CALL para_env%sum(ga)
                  CALL para_env%sum(sigma)
                  IF (iw > 0) THEN
                     CALL gerror(ga, gdeb(:, :, 3), ev1, ev2, ev3, ev4)
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [3B]", ev1, ev2, " %"
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [3B]", ev3, ev4, " %"
                     IF (use_virial) THEN
                        CALL serror(sigma, sdeb(:, :, 3), ev1, ev2)
                        WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [3B]", ev1, ev2, " %"
                     END IF
                  END IF
               END IF
            END IF
         END IF
 
         IF (grad) THEN
            CALL para_env%sum(dedcn)
            ga(:, :) = 0.0_dp
            sigma = 0.0_dp
            CALL dedcn_force(qs_env, dedcn, dcnum, ga, sigma)
            gradient(1:3, 1:natom_full) = gradient(1:3, 1:natom_full) + ga(1:3, 1:natom_full)
            stress = stress + sigma
            IF (debug) THEN
               CALL para_env%sum(ga)
               CALL para_env%sum(sigma)
               IF (iw > 0) THEN
                  CALL verror(dedcn, edcn, ev1, ev2)
                  WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdcn", ev1, ev2, " %"
                  CALL gerror(ga, gdeb(:, :, 4), ev1, ev2, ev3, ev4)
                  WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdcn]", ev1, ev2, " %"
                  WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdcn]", ev3, ev4, " %"
                  IF (use_virial) THEN
                     CALL serror(sigma, sdeb(:, :, 4), ev1, ev2)
                     WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdcn]", ev1, ev2, " %"
                  END IF
               END IF
            END IF
         END IF
         DEALLOCATE (q_red, cn_red)
         CALL cnumber_release(cn, dcnum, grad)
         te = m_walltime()
         tc = tc + te - ts
 
         IF (debug) THEN
            ta = sum(energies)
            CALL para_env%sum(ta)
            IF (iw > 0) THEN
               tb = sum(enerd2)
               ed2 = ta - tb
               pd2 = abs(ed2)/abs(tb)*100.
               WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 2-body", ed2, pd2, " %"
            END IF
            IF (dispersion_env%doabc) THEN
               ta = sum(energies3)
               CALL para_env%sum(ta)
               IF (iw > 0) THEN
                  tb = sum(enerd3)
                  ed3 = ta - tb
                  pd3 = abs(ed3)/abs(tb)*100.
                  WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 3-body", ed3, pd3, " %"
               END IF
            END IF
            IF (iw > 0) THEN
               WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for reference code [s]", td
               WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for production code [s]", tc
            END IF
         END IF
 
         IF (dispersion_env%doabc) THEN
            energies(:) = energies(:) + energies3(:)
         END IF
         evdw = sum(energies)
         IF (PRESENT(atomic_energy)) THEN
            atomic_energy(1:natom_full) = energies(1:natom_full)
         END IF
 
         IF (use_virial .AND. calculate_forces) THEN
            virial%pv_virial = virial%pv_virial - stress
         END IF
         IF (calculate_forces) THEN
            DO iatom = 1, natom_full
               IF (atom_map(iatom) == 0) cycle
               ikind = kind_of(iatom)
               atoma = atom_of_kind(iatom)
               force(ikind)%dispersion(:, atoma) = &
                  force(ikind)%dispersion(:, atoma) + gradient(:, iatom)
            END DO
         END IF
 
         DEALLOCATE (energies)
         IF (dispersion_env%doabc) DEALLOCATE (energies3)
         IF (grad) THEN
            DEALLOCATE (gradient, ga)
         END IF
 
      END IF
 
      DEALLOCATE (xyz, atomtype, atom_map, atom_map_back, species, exclude_ghost)
 
      CALL timestop(handle)
 
   END SUBROUTINE calculate_dispersion_d4_pairpot
 
! **************************************************************************************************
!> \brief ...
!> \param param ...
!> \param disp ...
!> \param mol ...
!> \param cutoff ...
!> \param grad ...
!> \param doabc ...
!> \param enerd2 ...
!> \param enerd3 ...
!> \param cnd ...
!> \param qd ...
!> \param dEdcn ...
!> \param dEdq ...
!> \param gradient ...
!> \param stress ...
! **************************************************************************************************
   SUBROUTINE refd4_debug(param, disp, mol, cutoff, grad, doabc, &
                          enerd2, enerd3, cnd, qd, dEdcn, dEdq, gradient, stress)
      CLASS(damping_param)                               :: param
      TYPE(d4_model)                                     :: disp
      TYPE(structure_type)                               :: mol
      TYPE(realspace_cutoff)                             :: cutoff
      TYPE(error_type), ALLOCATABLE                      :: error
      LOGICAL, INTENT(IN)                                :: grad, doabc
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: enerd2, enerd3, cnd, qd, dedcn, dedq
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gradient
      REAL(KIND=dp), DIMENSION(3, 3, 4)                  :: stress
 
      INTEGER                                            :: mref, natom, i, ncoup
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: q, qq
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: lattr, c6, dc6dcn, dc6dq
      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: cndr, cndl, qdr, qdl, gwdcn, gwdq, gwvec
 
      mref = maxval(disp%ref)
      natom = mol%nat
      ncoup = disp%ncoup
 
      ! Coordination numbers
      ALLOCATE (cnd(natom))
      IF (grad) ALLOCATE (cndr(3, natom, natom), cndl(3, 3, natom))
      CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%cn, lattr)
      CALL get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, disp%en, &
                                   cnd, cndr, cndl)
      ! EEQ charges
      ALLOCATE (qd(natom))
      IF (grad) ALLOCATE (qdr(3, natom, natom), qdl(3, 3, natom))
      CALL get_charges(disp%mchrg, mol, error, qd, qdr, qdl)
      IF (ALLOCATED(error)) THEN
         cpabort(error%message)
      END IF
      ! C6 interpolation
      ALLOCATE (gwvec(mref, natom, ncoup))
      IF (grad) ALLOCATE (gwdcn(mref, natom, ncoup), gwdq(mref, natom, ncoup))
      CALL disp%weight_references(mol, cnd, qd, gwvec, gwdcn, gwdq)
      ALLOCATE (c6(natom, natom))
      IF (grad) ALLOCATE (dc6dcn(natom, natom), dc6dq(natom, natom))
      CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
      CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp2, lattr)
      !
      IF (grad) THEN
         ALLOCATE (gradient(3, natom, 4))
         gradient = 0.0_dp
         stress = 0.0_dp
      END IF
      !
      ALLOCATE (enerd2(natom))
      enerd2(:) = 0.0_dp
      IF (grad) THEN
         ALLOCATE (dedcn(natom), dedq(natom))
         dedcn(:) = 0.0_dp; dedq(:) = 0.0_dp
      END IF
      CALL param%get_dispersion2(mol, lattr, cutoff%disp2, cutoff%width2, disp%r4r2, c6, &
                                 dc6dcn, dc6dq, enerd2, dedcn, dedq, gradient(:, :, 1), &
                                 stress(:, :, 1))
      !
      IF (grad) THEN
         DO i = 1, 3
            gradient(i, :, 2) = matmul(qdr(i, :, :), dedq(:))
            stress(i, :, 2) = matmul(qdl(i, :, :), dedq(:))
         END DO
      END IF
      !
      IF (doabc) THEN
         ALLOCATE (q(natom), qq(natom))
         q(:) = 0.0_dp; qq(:) = 0.0_dp
         ALLOCATE (enerd3(natom))
         enerd3(:) = 0.0_dp
         CALL disp%weight_references(mol, cnd, q, gwvec, gwdcn, gwdq)
         CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
         CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, lattr)
         CALL param%get_dispersion3(mol, lattr, cutoff%disp3, cutoff%width3, disp%r4r2, c6, &
                                    dc6dcn, dc6dq, enerd3, dedcn, qq, gradient(:, :, 3), &
                                    stress(:, :, 3))
      END IF
      IF (grad) THEN
         DO i = 1, 3
            gradient(i, :, 4) = matmul(cndr(i, :, :), dedcn(:))
            stress(i, :, 4) = matmul(cndl(i, :, :), dedcn(:))
         END DO
      END IF
 
   END SUBROUTINE refd4_debug
 
#else
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param dispersion_env ...
!> \param evdw ...
!> \param calculate_forces ...
!> \param iw ...
!> \param atomic_energy ...
! **************************************************************************************************
   SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
                                              iw, atomic_energy)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(qs_dispersion_type), INTENT(IN), POINTER      :: dispersion_env
      REAL(kind=dp), INTENT(INOUT)                       :: evdw
      LOGICAL, INTENT(IN)                                :: calculate_forces
      INTEGER, INTENT(IN)                                :: iw
      REAL(kind=dp), DIMENSION(:), OPTIONAL              :: atomic_energy
 
      mark_used(qs_env)
      mark_used(dispersion_env)
      mark_used(evdw)
      mark_used(calculate_forces)
      mark_used(iw)
      mark_used(atomic_energy)
 
      cpabort("CP2K build without DFTD4")
 
   END SUBROUTINE calculate_dispersion_d4_pairpot
 
#endif
 
! **************************************************************************************************
!> \brief ...
!> \param dispersion_env ...
!> \param cutoff ...
!> \param r4r2 ...
!> \param gwvec ...
!> \param gwdcn ...
!> \param gwdq ...
!> \param c6ref ...
!> \param mrefs ...
!> \param energies ...
!> \param dEdcn ...
!> \param dEdq ...
!> \param calculate_forces ...
!> \param gradient ...
!> \param stress ...
!> \param atom_map ...
!> \param species ...
! **************************************************************************************************
   SUBROUTINE dispersion_2b(dispersion_env, cutoff, r4r2, &
                            gwvec, gwdcn, gwdq, c6ref, mrefs, &
                            energies, dEdcn, dEdq, &
                            calculate_forces, gradient, stress, &
                            atom_map, species)
      TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: r4r2
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: gwvec, gwdcn, gwdq
      REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: energies, dedcn, dedq
      LOGICAL, INTENT(IN)                                :: calculate_forces
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: gradient, stress
      INTEGER, DIMENSION(:), INTENT(IN)                  :: atom_map, species
 
      INTEGER                                            :: ia, iatom, ik, ikind, ja, jatom, jk, &
                                                            jkind, mepos, num_pe
      REAL(kind=dp)                                      :: a1, a2, c6ij, cutoff2, d6, d8, de, dr2, &
                                                            edisp, fac, gdisp, r0ij, rrij, s6, s8, &
                                                            t6, t8
      REAL(kind=dp), DIMENSION(2)                        :: dcdcn, dcdq
      REAL(kind=dp), DIMENSION(3)                        :: dg, rij
      REAL(kind=dp), DIMENSION(3, 3)                     :: ds
      TYPE(neighbor_list_iterator_p_type), &
         DIMENSION(:), POINTER                           :: nl_iterator
      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
         POINTER                                         :: sab_vdw
 
      a1 = dispersion_env%a1
      a2 = dispersion_env%a2
      s6 = dispersion_env%s6
      s8 = dispersion_env%s8
      cutoff2 = cutoff*cutoff
 
      sab_vdw => dispersion_env%sab_vdw
 
      num_pe = 1
      CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
 
      mepos = 0
      DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
         CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
                                iatom=iatom, jatom=jatom, r=rij)
         ! Skip ghost/floating atoms
         ia = atom_map(iatom)
         ja = atom_map(jatom)
         IF (ia == 0 .OR. ja == 0) cycle
         ! D4 species indices
         ik = species(iatom)
         jk = species(jatom)
         ! vdW potential
         dr2 = sum(rij(:)**2)
         IF (dr2 <= cutoff2 .AND. dr2 > 0.0000001_dp) THEN
            rrij = 3._dp*r4r2(ik)*r4r2(jk)
            r0ij = a1*sqrt(rrij) + a2
            IF (calculate_forces) THEN
               CALL get_c6derivs(c6ij, dcdcn, dcdq, ia, ja, ik, jk, &
                                 gwvec, gwdcn, gwdq, c6ref, mrefs)
            ELSE
               CALL get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
            END IF
            fac = 1._dp
            IF (iatom == jatom) fac = 0.5_dp
            t6 = 1.0_dp/(dr2**3 + r0ij**6)
            t8 = 1.0_dp/(dr2**4 + r0ij**8)
 
            edisp = (s6*t6 + s8*rrij*t8)*fac
            de = -c6ij*edisp
            energies(iatom) = energies(iatom) + de*0.5_dp
            energies(jatom) = energies(jatom) + de*0.5_dp
 
            IF (calculate_forces) THEN
               d6 = -6.0_dp*dr2**2*t6**2
               d8 = -8.0_dp*dr2**3*t8**2
               gdisp = (s6*d6 + s8*rrij*d8)*fac
               dg(:) = -c6ij*gdisp*rij(:)
               gradient(:, iatom) = gradient(:, iatom) - dg
               gradient(:, jatom) = gradient(:, jatom) + dg
               ds(:, :) = spread(dg, 1, 3)*spread(rij, 2, 3)
               stress(:, :) = stress(:, :) + ds(:, :)
               dedcn(iatom) = dedcn(iatom) - dcdcn(1)*edisp
               dedq(iatom) = dedq(iatom) - dcdq(1)*edisp
               dedcn(jatom) = dedcn(jatom) - dcdcn(2)*edisp
               dedq(jatom) = dedq(jatom) - dcdq(2)*edisp
            END IF
         END IF
      END DO
 
      CALL neighbor_list_iterator_release(nl_iterator)
 
   END SUBROUTINE dispersion_2b
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param dispersion_env ...
!> \param tvec ...
!> \param cutoff ...
!> \param r4r2 ...
!> \param gwvec ...
!> \param gwdcn ...
!> \param gwdq ...
!> \param c6ref ...
!> \param mrefs ...
!> \param energies ...
!> \param dEdcn ...
!> \param dEdq ...
!> \param calculate_forces ...
!> \param gradient ...
!> \param stress ...
!> \param atom_map ...
!> \param species ...
! **************************************************************************************************
   SUBROUTINE dispersion_3b(qs_env, dispersion_env, tvec, cutoff, r4r2, &
                            gwvec, gwdcn, gwdq, c6ref, mrefs, &
                            energies, dEdcn, dEdq, &
                            calculate_forces, gradient, stress, &
                            atom_map, species)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: tvec
      REAL(kind=dp), INTENT(IN)                          :: cutoff
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: r4r2
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: gwvec, gwdcn, gwdq
      REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
      REAL(kind=dp), DIMENSION(:), INTENT(INOUT)         :: energies, dedcn, dedq
      LOGICAL, INTENT(IN)                                :: calculate_forces
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: gradient, stress
      INTEGER, DIMENSION(:), INTENT(IN)                  :: atom_map, species
 
      INTEGER                                            :: ia, iatom, ik, ikind, ja, jatom, jk, &
                                                            jkind, ka, katom, kk, ktr, mepos, &
                                                            natom, num_pe
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
      INTEGER, DIMENSION(3)                              :: cell_b
      REAL(kind=dp)                                      :: a1, a2, alp, ang, c6ij, c6ik, c6jk, c9, &
                                                            cutoff2, dang, de, dfdmp, fac, fdmp, &
                                                            r0, r0ij, r0ik, r0jk, r1, r2, r2ij, &
                                                            r2ik, r2jk, r3, r5, rr, s6, s8, s9
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: rcpbc
      REAL(kind=dp), DIMENSION(2)                        :: dc6dcnij, dc6dcnik, dc6dcnjk, dc6dqij, &
                                                            dc6dqik, dc6dqjk
      REAL(kind=dp), DIMENSION(3)                        :: dgij, dgik, dgjk, ra, rb, rb0, rij, vij, &
                                                            vik, vjk
      REAL(kind=dp), DIMENSION(3, 3)                     :: ds
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(neighbor_list_iterator_p_type), &
         DIMENSION(:), POINTER                           :: nl_iterator
      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
         POINTER                                         :: sab_vdw
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      CALL get_qs_env(qs_env=qs_env, natom=natom, cell=cell, &
                      atomic_kind_set=atomic_kind_set, particle_set=particle_set)
 
      ALLOCATE (rcpbc(3, natom))
      DO iatom = 1, natom
         rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
      END DO
      CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
 
      a1 = dispersion_env%a1
      a2 = dispersion_env%a2
      s6 = dispersion_env%s6
      s8 = dispersion_env%s8
      s9 = dispersion_env%s9
      alp = dispersion_env%alp
 
      cutoff2 = cutoff**2
 
      sab_vdw => dispersion_env%sab_vdw
 
      num_pe = 1
      CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
 
      mepos = 0
      DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
         CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
 
         ! Skip ghost/floating atoms
         ia = atom_map(iatom)
         ja = atom_map(jatom)
         IF (ia == 0 .OR. ja == 0) cycle
         ik = species(iatom)
         jk = species(jatom)
 
         r2ij = sum(rij(:)**2)
         IF (calculate_forces) THEN
            CALL get_c6derivs(c6ij, dc6dcnij, dc6dqij, ia, ja, ik, jk, &
                              gwvec, gwdcn, gwdq, c6ref, mrefs)
         ELSE
            CALL get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
         END IF
         r0ij = a1*sqrt(3._dp*r4r2(jk)*r4r2(ik)) + a2
         IF (r2ij <= cutoff2 .AND. r2ij > epsilon(1._dp)) THEN
            CALL get_iterator_info(nl_iterator, cell=cell_b)
            rb0(:) = matmul(cell%hmat, cell_b)
            ra(:) = rcpbc(:, iatom)
            rb(:) = rcpbc(:, jatom) + rb0
            vij(:) = rb(:) - ra(:)
 
            DO katom = 1, min(iatom, jatom)
               ka = atom_map(katom)
               IF (ka == 0) cycle
               kk = species(katom)
               IF (calculate_forces) THEN
                  CALL get_c6derivs(c6ik, dc6dcnik, dc6dqik, ka, ia, kk, ik, &
                                    gwvec, gwdcn, gwdq, c6ref, mrefs)
                  CALL get_c6derivs(c6jk, dc6dcnjk, dc6dqjk, ka, ja, kk, jk, &
                                    gwvec, gwdcn, gwdq, c6ref, mrefs)
               ELSE
                  CALL get_c6value(c6ik, ka, ia, kk, ik, gwvec, c6ref, mrefs)
                  CALL get_c6value(c6jk, ka, ja, kk, jk, gwvec, c6ref, mrefs)
               END IF
               c9 = -s9*sqrt(abs(c6ij*c6ik*c6jk))
               r0ik = a1*sqrt(3._dp*r4r2(kk)*r4r2(ik)) + a2
               r0jk = a1*sqrt(3._dp*r4r2(kk)*r4r2(jk)) + a2
               r0 = r0ij*r0ik*r0jk
               fac = triple_scale(iatom, jatom, katom)
               DO ktr = 1, SIZE(tvec, 2)
                  vik(:) = rcpbc(:, katom) + tvec(:, ktr) - rcpbc(:, iatom)
                  r2ik = vik(1)*vik(1) + vik(2)*vik(2) + vik(3)*vik(3)
                  IF (r2ik > cutoff2 .OR. r2ik < epsilon(1.0_dp)) cycle
                  vjk(:) = rcpbc(:, katom) + tvec(:, ktr) - rb(:)
                  r2jk = vjk(1)*vjk(1) + vjk(2)*vjk(2) + vjk(3)*vjk(3)
                  IF (r2jk > cutoff2 .OR. r2jk < epsilon(1.0_dp)) cycle
                  r2 = r2ij*r2ik*r2jk
                  r1 = sqrt(r2)
                  r3 = r2*r1
                  r5 = r3*r2
 
                  fdmp = 1.0_dp/(1.0_dp + 6.0_dp*(r0/r1)**(alp/3.0_dp))
                  ang = 0.375_dp*(r2ij + r2jk - r2ik)*(r2ij - r2jk + r2ik)* &
                        (-r2ij + r2jk + r2ik)/r5 + 1.0_dp/r3
 
                  rr = ang*fdmp
                  de = rr*c9*fac
                  energies(iatom) = energies(iatom) - de/3._dp
                  energies(jatom) = energies(jatom) - de/3._dp
                  energies(katom) = energies(katom) - de/3._dp
 
                  IF (calculate_forces) THEN
 
                     dfdmp = -2.0_dp*alp*(r0/r1)**(alp/3.0_dp)*fdmp**2
 
                     ! d/drij
                     dang = -0.375_dp*(r2ij**3 + r2ij**2*(r2jk + r2ik) &
                                       + r2ij*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ik &
                                               + 3.0_dp*r2ik**2) &
                                       - 5.0_dp*(r2jk - r2ik)**2*(r2jk + r2ik))/r5
                     dgij(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ij*vij
 
                     ! d/drik
                     dang = -0.375_dp*(r2ik**3 + r2ik**2*(r2jk + r2ij) &
                                       + r2ik*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ij &
                                               + 3.0_dp*r2ij**2) &
                                       - 5.0_dp*(r2jk - r2ij)**2*(r2jk + r2ij))/r5
                     dgik(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ik*vik
 
                     ! d/drjk
                     dang = -0.375_dp*(r2jk**3 + r2jk**2*(r2ik + r2ij) &
                                       + r2jk*(3.0_dp*r2ik**2 + 2.0_dp*r2ik*r2ij &
                                               + 3.0_dp*r2ij**2) &
                                       - 5.0_dp*(r2ik - r2ij)**2*(r2ik + r2ij))/r5
                     dgjk(:) = c9*(-dang*fdmp + ang*dfdmp)/r2jk*vjk
 
                     gradient(:, iatom) = gradient(:, iatom) - dgij - dgik
                     gradient(:, jatom) = gradient(:, jatom) + dgij - dgjk
                     gradient(:, katom) = gradient(:, katom) + dgik + dgjk
 
                     ds(:, :) = spread(dgij, 1, 3)*spread(vij, 2, 3) &
                                + spread(dgik, 1, 3)*spread(vik, 2, 3) &
                                + spread(dgjk, 1, 3)*spread(vjk, 2, 3)
 
                     stress(:, :) = stress + ds*fac
 
                     dedcn(iatom) = dedcn(iatom) - de*0.5_dp &
                                    *(dc6dcnij(1)/c6ij + dc6dcnik(2)/c6ik)
                     dedcn(jatom) = dedcn(jatom) - de*0.5_dp &
                                    *(dc6dcnij(2)/c6ij + dc6dcnjk(2)/c6jk)
                     dedcn(katom) = dedcn(katom) - de*0.5_dp &
                                    *(dc6dcnik(1)/c6ik + dc6dcnjk(1)/c6jk)
 
                     dedq(iatom) = dedq(iatom) - de*0.5_dp &
                                   *(dc6dqij(1)/c6ij + dc6dqik(2)/c6ik)
                     dedq(jatom) = dedq(jatom) - de*0.5_dp &
                                   *(dc6dqij(2)/c6ij + dc6dqjk(2)/c6jk)
                     dedq(katom) = dedq(katom) - de*0.5_dp &
                                   *(dc6dqik(1)/c6ik + dc6dqjk(1)/c6jk)
 
                  END IF
 
               END DO
            END DO
         END IF
      END DO
 
      CALL neighbor_list_iterator_release(nl_iterator)
 
      DEALLOCATE (rcpbc)
 
   END SUBROUTINE dispersion_3b
 
! **************************************************************************************************
!> \brief ...
!> \param ii ...
!> \param jj ...
!> \param kk ...
!> \return ...
! **************************************************************************************************
   FUNCTION triple_scale(ii, jj, kk) RESULT(triple)
      INTEGER, INTENT(IN)                                :: ii, jj, kk
      REAL(kind=dp)                                      :: triple
 
      IF (ii == jj) THEN
         IF (ii == kk) THEN
            ! ii'i" -> 1/6
            triple = 1.0_dp/6.0_dp
         ELSE
            ! ii'j -> 1/2
            triple = 0.5_dp
         END IF
      ELSE
         IF (ii /= kk .AND. jj /= kk) THEN
            ! ijk -> 1 (full)
            triple = 1.0_dp
         ELSE
            ! ijj' and iji' -> 1/2
            triple = 0.5_dp
         END IF
      END IF
 
   END FUNCTION triple_scale
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param dEdcn ...
!> \param dcnum ...
!> \param gradient ...
!> \param stress ...
! **************************************************************************************************
   SUBROUTINE dedcn_force(qs_env, dEdcn, dcnum, gradient, stress)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: dedcn
      TYPE(dcnum_type), DIMENSION(:), INTENT(IN)         :: dcnum
      REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT)      :: gradient
      REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT)      :: stress
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'dEdcn_force'
 
      INTEGER                                            :: handle, i, ia, iatom, ikind, katom, &
                                                            natom, nkind
      LOGICAL                                            :: use_virial
      REAL(kind=dp)                                      :: drk
      REAL(kind=dp), DIMENSION(3)                        :: fdik, rik
      TYPE(distribution_1d_type), POINTER                :: local_particles
      TYPE(virial_type), POINTER                         :: virial
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, nkind=nkind, natom=natom, &
                      local_particles=local_particles, &
                      virial=virial)
      use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
 
      DO ikind = 1, nkind
         DO ia = 1, local_particles%n_el(ikind)
            iatom = local_particles%list(ikind)%array(ia)
            DO i = 1, dcnum(iatom)%neighbors
               katom = dcnum(iatom)%nlist(i)
               rik = dcnum(iatom)%rik(:, i)
               drk = sqrt(sum(rik(:)**2))
               fdik(:) = -(dedcn(iatom) + dedcn(katom))*dcnum(iatom)%dvals(i)*rik(:)/drk
               gradient(:, iatom) = gradient(:, iatom) + fdik(:)
               IF (use_virial) THEN
                  CALL virial_pair_force(stress, -0.5_dp, fdik, rik)
               END IF
            END DO
         END DO
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE dedcn_force
 
! **************************************************************************************************
!> \brief ...
!> \param c6ij ...
!> \param ia ...
!> \param ja ...
!> \param ik ...
!> \param jk ...
!> \param gwvec ...
!> \param c6ref ...
!> \param mrefs ...
! **************************************************************************************************
   SUBROUTINE get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
      REAL(kind=dp), INTENT(OUT)                         :: c6ij
      INTEGER, INTENT(IN)                                :: ia, ja, ik, jk
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: gwvec
      REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
 
      INTEGER                                            :: iref, jref
      REAL(kind=dp)                                      :: refc6
 
      c6ij = 0.0_dp
      DO jref = 1, mrefs(jk)
         DO iref = 1, mrefs(ik)
            refc6 = c6ref(iref, jref, ik, jk)
            c6ij = c6ij + gwvec(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
         END DO
      END DO
 
   END SUBROUTINE get_c6value
 
! **************************************************************************************************
!> \brief ...
!> \param c6ij ...
!> \param dc6dcn ...
!> \param dc6dq ...
!> \param ia ...
!> \param ja ...
!> \param ik ...
!> \param jk ...
!> \param gwvec ...
!> \param gwdcn ...
!> \param gwdq ...
!> \param c6ref ...
!> \param mrefs ...
! **************************************************************************************************
   SUBROUTINE get_c6derivs(c6ij, dc6dcn, dc6dq, ia, ja, ik, jk, &
                           gwvec, gwdcn, gwdq, c6ref, mrefs)
      REAL(kind=dp), INTENT(OUT)                         :: c6ij
      REAL(kind=dp), DIMENSION(2), INTENT(OUT)           :: dc6dcn, dc6dq
      INTEGER, INTENT(IN)                                :: ia, ja, ik, jk
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: gwvec, gwdcn, gwdq
      REAL(kind=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
      INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
 
      INTEGER                                            :: iref, jref
      REAL(kind=dp)                                      :: refc6
 
      c6ij = 0.0_dp
      dc6dcn = 0.0_dp
      dc6dq = 0.0_dp
      DO jref = 1, mrefs(jk)
         DO iref = 1, mrefs(ik)
            refc6 = c6ref(iref, jref, ik, jk)
            c6ij = c6ij + gwvec(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
            dc6dcn(1) = dc6dcn(1) + gwdcn(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
            dc6dcn(2) = dc6dcn(2) + gwvec(iref, ia, 1)*gwdcn(jref, ja, 1)*refc6
            dc6dq(1) = dc6dq(1) + gwdq(iref, ia, 1)*gwvec(jref, ja, 1)*refc6
            dc6dq(2) = dc6dq(2) + gwvec(iref, ia, 1)*gwdq(jref, ja, 1)*refc6
         END DO
      END DO
 
   END SUBROUTINE get_c6derivs
 
! **************************************************************************************************
!> \brief ...
!> \param ga ...
!> \param gd ...
!> \param ev1 ...
!> \param ev2 ...
!> \param ev3 ...
!> \param ev4 ...
! **************************************************************************************************
   SUBROUTINE gerror(ga, gd, ev1, ev2, ev3, ev4)
      REAL(kind=dp), DIMENSION(:, :)                     :: ga, gd
      REAL(kind=dp), INTENT(OUT)                         :: ev1, ev2, ev3, ev4
 
      INTEGER                                            :: na, np(2)
 
      na = SIZE(ga, 2)
      ev1 = sqrt(sum((gd - ga)**2)/na)
      ev2 = ev1/sqrt(sum(gd**2)/na)*100._dp
      np = maxloc(abs(gd - ga))
      ev3 = abs(gd(np(1), np(2)) - ga(np(1), np(2)))
      ev4 = abs(gd(np(1), np(2)))
      IF (ev4 > 1.e-6) THEN
         ev4 = ev3/ev4*100._dp
      ELSE
         ev4 = 100.0_dp
      END IF
 
   END SUBROUTINE gerror
 
! **************************************************************************************************
!> \brief ...
!> \param sa ...
!> \param sd ...
!> \param ev1 ...
!> \param ev2 ...
! **************************************************************************************************
   SUBROUTINE serror(sa, sd, ev1, ev2)
      REAL(kind=dp), DIMENSION(3, 3)                     :: sa, sd
      REAL(kind=dp), INTENT(OUT)                         :: ev1, ev2
 
      INTEGER                                            :: i, j
      REAL(kind=dp)                                      :: rel
 
      ev1 = maxval(abs(sd - sa))
      ev2 = 0.0_dp
      DO i = 1, 3
         DO j = 1, 3
            IF (abs(sd(i, j)) > 1.e-6_dp) THEN
               rel = abs(sd(i, j) - sa(i, j))/abs(sd(i, j))*100._dp
               ev2 = max(ev2, rel)
            END IF
         END DO
      END DO
 
   END SUBROUTINE serror
 
! **************************************************************************************************
!> \brief ...
!> \param va ...
!> \param vd ...
!> \param ev1 ...
!> \param ev2 ...
! **************************************************************************************************
   SUBROUTINE verror(va, vd, ev1, ev2)
      REAL(kind=dp), DIMENSION(:)                        :: va, vd
      REAL(kind=dp), INTENT(OUT)                         :: ev1, ev2
 
      INTEGER                                            :: i, na
      REAL(kind=dp)                                      :: rel
 
      na = SIZE(va)
      ev1 = maxval(abs(vd - va))
      ev2 = 0.0_dp
      DO i = 1, na
         IF (abs(vd(i)) > 1.e-8_dp) THEN
            rel = abs(vd(i) - va(i))/abs(vd(i))*100._dp
            ev2 = max(ev2, rel)
         END IF
      END DO
 
   END SUBROUTINE verror
 
END MODULE qs_dispersion_d4