LCOV - code coverage report
Current view: top level - src - qs_dispersion_d4.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 68.8 % 516 355
Test Date: 2025-07-25 12:55:17 Functions: 63.6 % 11 7

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculation of dispersion using pair potentials
      10              : !> \author Johann Pototschnig
      11              : ! **************************************************************************************************
      12              : MODULE qs_dispersion_d4
      13              :    USE atomic_kind_types, ONLY: atomic_kind_type, &
      14              :                                 get_atomic_kind, &
      15              :                                 get_atomic_kind_set
      16              :    USE distribution_1d_types, ONLY: distribution_1d_type
      17              :    USE eeq_method, ONLY: eeq_charges, eeq_forces
      18              :    USE machine, ONLY: m_flush, &
      19              :                       m_walltime
      20              :    USE cell_types, ONLY: cell_type, &
      21              :                          plane_distance, &
      22              :                          pbc, &
      23              :                          get_cell
      24              :    USE qs_environment_types, ONLY: get_qs_env, &
      25              :                                    qs_environment_type
      26              :    USE qs_force_types, ONLY: qs_force_type
      27              :    USE qs_kind_types, ONLY: get_qs_kind, &
      28              :                             qs_kind_type, &
      29              :                             set_qs_kind
      30              :    USE qs_neighbor_list_types, ONLY: get_iterator_info, &
      31              :                                      neighbor_list_iterate, &
      32              :                                      neighbor_list_iterator_create, &
      33              :                                      neighbor_list_iterator_p_type, &
      34              :                                      neighbor_list_iterator_release, &
      35              :                                      neighbor_list_set_p_type
      36              :    USE virial_methods, ONLY: virial_pair_force
      37              :    USE virial_types, ONLY: virial_type
      38              :    USE kinds, ONLY: dp
      39              :    USE particle_types, ONLY: particle_type
      40              :    USE qs_dispersion_types, ONLY: qs_dispersion_type
      41              :    USE qs_dispersion_utils, ONLY: cellhash
      42              :    USE qs_dispersion_cnum, ONLY: cnumber_init, dcnum_type, cnumber_release
      43              :    USE message_passing, ONLY: mp_para_env_type
      44              : 
      45              : #if defined(__DFTD4)
      46              : !&<
      47              :    USE dftd4,                           ONLY: d4_model, &
      48              :                                               damping_param, &
      49              :                                               get_dispersion, &
      50              :                                               get_rational_damping, &
      51              :                                               new, &
      52              :                                               new_d4_model, &
      53              :                                               realspace_cutoff, &
      54              :                                               structure_type, &
      55              :                                               rational_damping_param, &
      56              :                                               get_coordination_number, &
      57              :                                               get_lattice_points
      58              :    USE dftd4_charge,                    ONLY: get_charges
      59              : !&>
      60              : #endif
      61              : #include "./base/base_uses.f90"
      62              : 
      63              :    IMPLICIT NONE
      64              : 
      65              :    PRIVATE
      66              : 
      67              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d4'
      68              : 
      69              :    PUBLIC :: calculate_dispersion_d4_pairpot
      70              : 
      71              : ! **************************************************************************************************
      72              : 
      73              : CONTAINS
      74              : 
      75              : #if defined(__DFTD4)
      76              : ! **************************************************************************************************
      77              : !> \brief ...
      78              : !> \param qs_env ...
      79              : !> \param dispersion_env ...
      80              : !> \param evdw ...
      81              : !> \param calculate_forces ...
      82              : !> \param iw ...
      83              : !> \param atomic_energy ...
      84              : ! **************************************************************************************************
      85          894 :    SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw, &
      86          894 :                                               atomic_energy)
      87              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      88              :       TYPE(qs_dispersion_type), INTENT(IN), POINTER      :: dispersion_env
      89              :       REAL(KIND=dp), INTENT(INOUT)                       :: evdw
      90              :       LOGICAL, INTENT(IN)                                :: calculate_forces
      91              :       INTEGER, INTENT(IN)                                :: iw
      92              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atomic_energy
      93              : 
      94              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d4_pairpot'
      95              : 
      96              :       INTEGER                                            :: atoma, cnfun, enshift, handle, i, iatom, &
      97              :                                                             ikind, mref, natom, nghost
      98          894 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, atomtype, kind_of, &
      99          894 :                                                             t_atomtype
     100              :       INTEGER, DIMENSION(3)                              :: periodic
     101              :       LOGICAL                                            :: debug, grad, ifloating, ighost, &
     102              :                                                             use_virial
     103          894 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: a_ghost
     104              :       LOGICAL, DIMENSION(3)                              :: lperiod
     105              :       REAL(KIND=dp)                                      :: ed2, ed3, ev1, ev2, ev3, ev4, pd2, pd3, &
     106              :                                                             ta, tb, tc, td, te, ts
     107          894 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cn, cnd, dEdcn, dEdq, edcn, edq, enerd2, &
     108          894 :                                                             enerd3, energies, energies3, q, qd
     109          894 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ga, gradient, gwdcn, gwdq, gwvec, t_xyz, &
     110          894 :                                                             tvec, xyz
     111          894 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gdeb
     112              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: sigma, stress
     113              :       REAL(KIND=dp), DIMENSION(3, 3, 4)                  :: sdeb
     114          894 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     115              :       TYPE(cell_type), POINTER                           :: cell
     116          894 :       TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:)        :: dcnum
     117              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     118          894 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     119          894 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     120          894 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     121              :       TYPE(virial_type), POINTER                         :: virial
     122              : 
     123         1788 :       CLASS(damping_param), ALLOCATABLE                  :: param
     124          894 :       TYPE(d4_model)                                     :: disp
     125          894 :       TYPE(structure_type)                               :: mol
     126              :       TYPE(realspace_cutoff)                             :: cutoff
     127              : 
     128          894 :       CALL timeset(routineN, handle)
     129              : 
     130          894 :       debug = dispersion_env%d4_debug
     131              : 
     132              :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
     133          894 :                       cell=cell, force=force, virial=virial, para_env=para_env)
     134          894 :       CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     135              : 
     136              :       !get information about particles
     137          894 :       natom = SIZE(particle_set)
     138          894 :       nghost = 0
     139         5364 :       ALLOCATE (t_xyz(3, natom), t_atomtype(natom), a_ghost(natom))
     140          894 :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
     141         4652 :       DO iatom = 1, natom
     142        15032 :          t_xyz(:, iatom) = particle_set(iatom)%r(:)
     143         3758 :          ikind = kind_of(iatom)
     144         3758 :          CALL get_qs_kind(qs_kind_set(ikind), zatom=t_atomtype(iatom), ghost=ighost, floating=ifloating)
     145         3758 :          a_ghost(iatom) = ighost .OR. ifloating
     146         4652 :          IF (a_ghost(iatom)) nghost = nghost + 1
     147              :       END DO
     148              : 
     149          894 :       natom = natom - nghost
     150          894 :       iatom = 0
     151         4470 :       ALLOCATE (xyz(3, natom), atomtype(natom))
     152         4652 :       DO i = 1, natom + nghost
     153         4652 :          IF (.NOT. a_ghost(i)) THEN
     154         3746 :             iatom = iatom + 1
     155        14984 :             xyz(:, iatom) = t_xyz(:, i)
     156         3746 :             atomtype(iatom) = t_atomtype(i)
     157              :          END IF
     158              :       END DO
     159          894 :       DEALLOCATE (a_ghost, t_xyz, t_atomtype)
     160              : 
     161              :       !get information about cell / lattice
     162          894 :       CALL get_cell(cell=cell, periodic=periodic)
     163          894 :       lperiod(1) = periodic(1) == 1
     164          894 :       lperiod(2) = periodic(2) == 1
     165          894 :       lperiod(3) = periodic(3) == 1
     166              :       ! enforce en shift method 1 (original/molecular)
     167              :       ! method 2 from paper on PBC seems not to work
     168          894 :       enshift = 1
     169              :       !IF (ALL(periodic == 0)) enshift = 1
     170              : 
     171              :       !prepare for the call to the dispersion function
     172          894 :       CALL new(mol, atomtype, xyz, lattice=cell%hmat, periodic=lperiod)
     173          894 :       CALL new_d4_model(disp, mol)
     174              : 
     175          894 :       IF (dispersion_env%ref_functional == "none") THEN
     176          858 :          CALL get_rational_damping("pbe", param, s9=0.0_dp)
     177              :          SELECT TYPE (param)
     178              :          TYPE is (rational_damping_param)
     179          858 :             param%s6 = dispersion_env%s6
     180          858 :             param%s8 = dispersion_env%s8
     181          858 :             param%a1 = dispersion_env%a1
     182          858 :             param%a2 = dispersion_env%a2
     183          858 :             param%alp = dispersion_env%alp
     184              :          END SELECT
     185              :       ELSE
     186              :          CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
     187           36 :          SELECT TYPE (param)
     188              :          TYPE is (rational_damping_param)
     189           36 :             dispersion_env%s6 = param%s6
     190           36 :             dispersion_env%s8 = param%s8
     191           36 :             dispersion_env%a1 = param%a1
     192           36 :             dispersion_env%a2 = param%a2
     193           36 :             dispersion_env%alp = param%alp
     194              :          END SELECT
     195              :       END IF
     196              : 
     197              :       ! Coordination number cutoff
     198          894 :       cutoff%cn = dispersion_env%rc_cn
     199              :       ! Two-body interaction cutoff
     200          894 :       cutoff%disp2 = dispersion_env%rc_d4*2._dp
     201              :       ! Three-body interaction cutoff
     202          894 :       cutoff%disp3 = dispersion_env%rc_disp*2._dp
     203          894 :       IF (cutoff%disp3 > cutoff%disp2) THEN
     204            0 :          CPABORT("D4: Three-body cutoff should be smaller than two-body cutoff")
     205              :       END IF
     206              : 
     207          894 :       IF (calculate_forces) THEN
     208           24 :          grad = .TRUE.
     209           24 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     210              :       ELSE
     211          870 :          grad = .FALSE.
     212          870 :          use_virial = .FALSE.
     213              :       END IF
     214              : 
     215          894 :       IF (dispersion_env%d4_reference_code) THEN
     216              : 
     217              :          !> Wrapper to handle the evaluation of dispersion energy and derivatives
     218           22 :          IF (.NOT. dispersion_env%doabc) THEN
     219            0 :             CPWARN("Using D4_REFERENCE_CODE enforces calculation of C9 term.")
     220              :          END IF
     221           22 :          IF (grad) THEN
     222           12 :             ALLOCATE (gradient(3, natom))
     223            6 :             CALL get_dispersion(mol, disp, param, cutoff, evdw, gradient, stress)
     224            6 :             IF (calculate_forces) THEN
     225            6 :                IF (use_virial) THEN
     226           26 :                   virial%pv_virial = virial%pv_virial - stress/para_env%num_pe
     227              :                END IF
     228           54 :                DO iatom = 1, natom
     229           48 :                   ikind = kind_of(iatom)
     230           48 :                   atoma = atom_of_kind(iatom)
     231              :                   force(ikind)%dispersion(:, atoma) = &
     232          198 :                      force(ikind)%dispersion(:, atoma) + gradient(:, iatom)/para_env%num_pe
     233              :                END DO
     234              :             END IF
     235            6 :             DEALLOCATE (gradient)
     236              :          ELSE
     237           16 :             CALL get_dispersion(mol, disp, param, cutoff, evdw)
     238              :          END IF
     239              :          !dispersion energy is computed by every MPI process
     240           22 :          evdw = evdw/para_env%num_pe
     241           22 :          IF (dispersion_env%ext_charges) dispersion_env%dcharges = 0.0_dp
     242           22 :          IF (PRESENT(atomic_energy)) THEN
     243            0 :             CPWARN("Atomic energies not available for D4 reference code")
     244            0 :             atomic_energy = 0.0_dp
     245              :          END IF
     246              : 
     247              :       ELSE
     248              : 
     249          872 :          IF (iw > 0) THEN
     250            0 :             WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
     251            0 :             WRITE (iw, FMT="(T32,A)") "DEBUG D4 DISPERSION"
     252            0 :             WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
     253            0 :             WRITE (iw, '(A,T71,A10)') " DEBUG D4| Reference functional   ", TRIM(dispersion_env%ref_functional)
     254            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s6) ", dispersion_env%s6
     255            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Scaling parameter (s8) ", dispersion_env%s8
     256            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a1) ", dispersion_env%a1
     257            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| BJ Damping parameter (a2) ", dispersion_env%a2
     258            0 :             WRITE (iw, '(A,T71,E10.4)') " DEBUG D4| Cutoff value coordination numbers ", dispersion_env%eps_cn
     259            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius coordination numbers ", dispersion_env%rc_cn
     260            0 :             WRITE (iw, '(A,T71,I10)') " DEBUG D4| Coordination number function type ", dispersion_env%cnfun
     261            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 2-body terms [bohr]", 2._dp*dispersion_env%rc_d4
     262            0 :             WRITE (iw, '(A,T71,F10.4)') " DEBUG D4| Cutoff radius 3-body terms [bohr]", 2._dp*dispersion_env%rc_disp
     263              :          END IF
     264              : 
     265          872 :          td = 0.0_dp
     266          872 :          IF (debug .AND. iw > 0) THEN
     267            0 :             ts = m_walltime()
     268              :             CALL refd4_debug(param, disp, mol, cutoff, grad, dispersion_env%doabc, &
     269            0 :                              enerd2, enerd3, cnd, qd, Edcn, Edq, gdeb, sdeb)
     270            0 :             te = m_walltime()
     271            0 :             td = te - ts
     272              :          END IF
     273              : 
     274          872 :          tc = 0.0_dp
     275          872 :          ts = m_walltime()
     276              : 
     277         2904 :          mref = MAXVAL(disp%ref)
     278              :          ! Coordination numbers
     279          872 :          cnfun = dispersion_env%cnfun
     280          872 :          CALL cnumber_init(qs_env, cn, dcnum, cnfun, grad)
     281          872 :          IF (debug .AND. iw > 0) THEN
     282            0 :             WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (max)", MAXVAL(ABS(cn - cnd))
     283            0 :             WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| CN differences (ave)", SUM(ABS(cn - cnd))/natom
     284              :          END IF
     285              : 
     286              :          ! EEQ charges
     287         2616 :          ALLOCATE (q(natom))
     288          872 :          IF (dispersion_env%ext_charges) THEN
     289         4320 :             q(1:natom) = dispersion_env%charges(1:natom)
     290              :          ELSE
     291           14 :             CALL eeq_charges(qs_env, q, dispersion_env%eeq_sparam, 2, enshift)
     292              :          END IF
     293          872 :          IF (debug .AND. iw > 0) THEN
     294            0 :             WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (max)", MAXVAL(ABS(q - qd))
     295            0 :             WRITE (iw, '(A,T71,F10.6)') " DEBUG D4| Charge differences (ave)", SUM(ABS(q - qd))/natom
     296              :          END IF
     297              :          ! Weights for C6 calculation
     298         3488 :          ALLOCATE (gwvec(mref, natom))
     299          944 :          IF (grad) ALLOCATE (gwdcn(mref, natom), gwdq(mref, natom))
     300          872 :          CALL disp%weight_references(mol, cn, q, gwvec, gwdcn, gwdq)
     301              : 
     302         1744 :          ALLOCATE (energies(natom))
     303         4474 :          energies(:) = 0.0_dp
     304          872 :          IF (grad) THEN
     305           54 :             ALLOCATE (gradient(3, natom), ga(3, natom))
     306           54 :             ALLOCATE (dEdcn(natom), dEdq(natom))
     307          234 :             dEdcn(:) = 0.0_dp; dEdq(:) = 0.0_dp
     308          450 :             ga(:, :) = 0.0_dp
     309           18 :             sigma(:, :) = 0.0_dp
     310              :          END IF
     311              :          CALL dispersion_2b(dispersion_env, cutoff%disp2, disp%r4r2, &
     312              :                             gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
     313          872 :                             energies, dEdcn, dEdq, grad, ga, sigma)
     314          872 :          IF (grad) THEN
     315          450 :             gradient(1:3, 1:natom) = ga(1:3, 1:natom)
     316           18 :             stress = sigma
     317           18 :             IF (debug) THEN
     318            0 :                CALL para_env%sum(ga)
     319            0 :                CALL para_env%sum(sigma)
     320            0 :                IF (iw > 0) THEN
     321            0 :                   CALL gerror(ga, gdeb(:, :, 1), ev1, ev2, ev3, ev4)
     322            0 :                   WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [2B]", ev1, ev2, " %"
     323            0 :                   WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [2B]", ev3, ev4, " %"
     324            0 :                   IF (use_virial) THEN
     325            0 :                      CALL serror(sigma, sdeb(:, :, 1), ev1, ev2)
     326            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [2B]", ev1, ev2, " %"
     327              :                   END IF
     328              :                END IF
     329              :             END IF
     330              :          END IF
     331              :          ! no contribution from dispersion_3b as q=0 (but q is changed!)
     332              :          ! so we callculate this here
     333          872 :          IF (grad) THEN
     334           18 :             IF (dispersion_env%ext_charges) THEN
     335           60 :                dispersion_env%dcharges = dEdq
     336              :             ELSE
     337            6 :                CALL para_env%sum(dEdq)
     338          246 :                ga(:, :) = 0.0_dp
     339            6 :                sigma = 0.0_dp
     340              :                CALL eeq_forces(qs_env, q, dEdq, ga, sigma, dispersion_env%eeq_sparam, &
     341            6 :                                2, enshift, response_only=.TRUE.)
     342          246 :                gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + ga(1:3, 1:natom)
     343           78 :                stress = stress + sigma
     344            6 :                IF (debug) THEN
     345            0 :                   CALL para_env%sum(ga)
     346            0 :                   CALL para_env%sum(sigma)
     347            0 :                   IF (iw > 0) THEN
     348            0 :                      CALL verror(dEdq, Edq, ev1, ev2)
     349            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdq", ev1, ev2, " %"
     350            0 :                      CALL gerror(ga, gdeb(:, :, 2), ev1, ev2, ev3, ev4)
     351            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdq]", ev1, ev2, " %"
     352            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdq]", ev3, ev4, " %"
     353            0 :                      IF (use_virial) THEN
     354            0 :                         CALL serror(sigma, sdeb(:, :, 2), ev1, ev2)
     355            0 :                         WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdq]", ev1, ev2, " %"
     356              :                      END IF
     357              :                   END IF
     358              :                END IF
     359              :             END IF
     360              :          END IF
     361              : 
     362          872 :          IF (dispersion_env%doabc) THEN
     363           28 :             ALLOCATE (energies3(natom))
     364          154 :             energies3(:) = 0.0_dp
     365          154 :             q(:) = 0.0_dp
     366              :             ! i.e. dc6dq = dEdq = 0
     367           14 :             CALL disp%weight_references(mol, cn, q, gwvec, gwdcn, gwdq)
     368              :             !
     369           14 :             IF (grad) THEN
     370          486 :                gwdq = 0.0_dp
     371          246 :                ga(:, :) = 0.0_dp
     372            6 :                sigma = 0.0_dp
     373              :             END IF
     374              :             CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, tvec)
     375              :             CALL dispersion_3b(qs_env, dispersion_env, tvec, cutoff%disp3, disp%r4r2, &
     376              :                                gwvec, gwdcn, gwdq, disp%c6, disp%ref, &
     377           14 :                                energies3, dEdcn, dEdq, grad, ga, sigma)
     378           28 :             IF (grad) THEN
     379          246 :                gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + ga(1:3, 1:natom)
     380           78 :                stress = stress + sigma
     381            6 :                IF (debug) THEN
     382            0 :                   CALL para_env%sum(ga)
     383            0 :                   CALL para_env%sum(sigma)
     384            0 :                   IF (iw > 0) THEN
     385            0 :                      CALL gerror(ga, gdeb(:, :, 3), ev1, ev2, ev3, ev4)
     386            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [3B]", ev1, ev2, " %"
     387            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [3B]", ev3, ev4, " %"
     388            0 :                      IF (use_virial) THEN
     389            0 :                         CALL serror(sigma, sdeb(:, :, 3), ev1, ev2)
     390            0 :                         WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [3B]", ev1, ev2, " %"
     391              :                      END IF
     392              :                   END IF
     393              :                END IF
     394              :             END IF
     395              :          END IF
     396              : 
     397          872 :          IF (grad) THEN
     398           18 :             CALL para_env%sum(dEdcn)
     399          450 :             ga(:, :) = 0.0_dp
     400           18 :             sigma = 0.0_dp
     401           18 :             CALL dEdcn_force(qs_env, dEdcn, dcnum, ga, sigma)
     402          450 :             gradient(1:3, 1:natom) = gradient(1:3, 1:natom) + ga(1:3, 1:natom)
     403          234 :             stress = stress + sigma
     404           18 :             IF (debug) THEN
     405            0 :                CALL para_env%sum(ga)
     406            0 :                CALL para_env%sum(sigma)
     407            0 :                IF (iw > 0) THEN
     408            0 :                   CALL verror(dEdcn, Edcn, ev1, ev2)
     409            0 :                   WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Derivative dEdcn", ev1, ev2, " %"
     410            0 :                   CALL gerror(ga, gdeb(:, :, 4), ev1, ev2, ev3, ev4)
     411            0 :                   WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| RMS error Gradient [dEdcn]", ev1, ev2, " %"
     412            0 :                   WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Gradient [dEdcn]", ev3, ev4, " %"
     413            0 :                   IF (use_virial) THEN
     414            0 :                      CALL serror(sigma, sdeb(:, :, 4), ev1, ev2)
     415            0 :                      WRITE (iw, '(A,T51,F14.10,T69,F10.4,A)') " DEBUG D4| MAV error Stress [dEdcn]", ev1, ev2, " %"
     416              :                   END IF
     417              :                END IF
     418              :             END IF
     419              :          END IF
     420          872 :          DEALLOCATE (q)
     421          872 :          CALL cnumber_release(cn, dcnum, grad)
     422          872 :          te = m_walltime()
     423          872 :          tc = tc + te - ts
     424              : 
     425          872 :          IF (debug) THEN
     426            0 :             ta = SUM(energies)
     427            0 :             CALL para_env%sum(ta)
     428            0 :             IF (iw > 0) THEN
     429            0 :                tb = SUM(enerd2)
     430            0 :                ed2 = ta - tb
     431            0 :                pd2 = ABS(ed2)/ABS(tb)*100.
     432            0 :                WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 2-body", ed2, pd2, " %"
     433              :             END IF
     434            0 :             IF (dispersion_env%doabc) THEN
     435            0 :                ta = SUM(energies3)
     436            0 :                CALL para_env%sum(ta)
     437            0 :                IF (iw > 0) THEN
     438            0 :                   tb = SUM(enerd3)
     439            0 :                   ed3 = ta - tb
     440            0 :                   pd3 = ABS(ed3)/ABS(tb)*100.
     441            0 :                   WRITE (iw, '(A,T51,F14.8,T69,F10.4,A)') " DEBUG D4| Energy error 3-body", ed3, pd3, " %"
     442              :                END IF
     443              :             END IF
     444            0 :             IF (iw > 0) THEN
     445            0 :                WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for reference code [s]", td
     446            0 :                WRITE (iw, '(A,T67,F14.4)') " DEBUG D4| Time for production code [s]", tc
     447              :             END IF
     448              :          END IF
     449              : 
     450          872 :          IF (dispersion_env%doabc) THEN
     451          154 :             energies(:) = energies(:) + energies3(:)
     452              :          END IF
     453         4474 :          evdw = SUM(energies)
     454          872 :          IF (PRESENT(atomic_energy)) THEN
     455            0 :             atomic_energy(1:natom) = energies(1:natom)
     456              :          END IF
     457              : 
     458          872 :          IF (use_virial .AND. calculate_forces) THEN
     459           52 :             virial%pv_virial = virial%pv_virial - stress
     460              :          END IF
     461          872 :          IF (calculate_forces) THEN
     462          126 :             DO iatom = 1, natom
     463          108 :                ikind = kind_of(iatom)
     464          108 :                atoma = atom_of_kind(iatom)
     465              :                force(ikind)%dispersion(:, atoma) = &
     466          450 :                   force(ikind)%dispersion(:, atoma) + gradient(:, iatom)
     467              :             END DO
     468              :          END IF
     469              : 
     470          872 :          DEALLOCATE (energies)
     471          872 :          IF (dispersion_env%doabc) DEALLOCATE (energies3)
     472          872 :          IF (grad) THEN
     473           18 :             DEALLOCATE (gradient, ga)
     474              :          END IF
     475              : 
     476              :       END IF
     477              : 
     478          894 :       DEALLOCATE (xyz, atomtype)
     479              : 
     480          894 :       CALL timestop(handle)
     481              : 
     482         1788 :    END SUBROUTINE calculate_dispersion_d4_pairpot
     483              : 
     484              : ! **************************************************************************************************
     485              : !> \brief ...
     486              : !> \param param ...
     487              : !> \param disp ...
     488              : !> \param mol ...
     489              : !> \param cutoff ...
     490              : !> \param grad ...
     491              : !> \param doabc ...
     492              : !> \param enerd2 ...
     493              : !> \param enerd3 ...
     494              : !> \param cnd ...
     495              : !> \param qd ...
     496              : !> \param dEdcn ...
     497              : !> \param dEdq ...
     498              : !> \param gradient ...
     499              : !> \param stress ...
     500              : ! **************************************************************************************************
     501            0 :    SUBROUTINE refd4_debug(param, disp, mol, cutoff, grad, doabc, &
     502              :                           enerd2, enerd3, cnd, qd, dEdcn, dEdq, gradient, stress)
     503              :       CLASS(damping_param)                               :: param
     504              :       TYPE(d4_model)                                     :: disp
     505              :       TYPE(structure_type)                               :: mol
     506              :       TYPE(realspace_cutoff)                             :: cutoff
     507              :       LOGICAL, INTENT(IN)                                :: grad, doabc
     508              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: enerd2, enerd3, cnd, qd, dEdcn, dEdq
     509              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gradient
     510              :       REAL(KIND=dp), DIMENSION(3, 3, 4)                  :: stress
     511              : 
     512              :       INTEGER                                            :: mref, natom, i
     513            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: q, qq
     514            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: lattr, gwdcn, gwdq, gwvec, &
     515            0 :                                                             c6, dc6dcn, dc6dq
     516            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: cndr, cndL, qdr, qdL
     517              : 
     518            0 :       mref = MAXVAL(disp%ref)
     519            0 :       natom = mol%nat
     520              : 
     521              :       ! Coordination numbers
     522            0 :       ALLOCATE (cnd(natom))
     523            0 :       IF (grad) ALLOCATE (cndr(3, natom, natom), cndL(3, 3, natom))
     524              :       CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%cn, lattr)
     525              :       CALL get_coordination_number(mol, lattr, cutoff%cn, disp%rcov, disp%en, &
     526            0 :                                    cnd, cndr, cndL)
     527              :       ! EEQ charges
     528            0 :       ALLOCATE (qd(natom))
     529            0 :       IF (grad) ALLOCATE (qdr(3, natom, natom), qdL(3, 3, natom))
     530            0 :       CALL get_charges(mol, qd, qdr, qdL)
     531              :       ! C6 interpolation
     532            0 :       ALLOCATE (gwvec(mref, natom))
     533            0 :       IF (grad) ALLOCATE (gwdcn(mref, natom), gwdq(mref, natom))
     534            0 :       CALL disp%weight_references(mol, cnd, qd, gwvec, gwdcn, gwdq)
     535            0 :       ALLOCATE (c6(natom, natom))
     536            0 :       IF (grad) ALLOCATE (dc6dcn(natom, natom), dc6dq(natom, natom))
     537            0 :       CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
     538            0 :       CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp2, lattr)
     539              :       !
     540            0 :       IF (grad) THEN
     541            0 :          ALLOCATE (gradient(3, natom, 4))
     542            0 :          gradient = 0.0_dp
     543            0 :          stress = 0.0_dp
     544              :       END IF
     545              :       !
     546            0 :       ALLOCATE (enerd2(natom))
     547            0 :       enerd2(:) = 0.0_dp
     548            0 :       IF (grad) THEN
     549            0 :          ALLOCATE (dEdcn(natom), dEdq(natom))
     550            0 :          dEdcn(:) = 0.0_dp; dEdq(:) = 0.0_dp
     551              :       END IF
     552              :       CALL param%get_dispersion2(mol, lattr, cutoff%disp2, disp%r4r2, c6, dc6dcn, dc6dq, &
     553            0 :                                  enerd2, dEdcn, dEdq, gradient(:, :, 1), stress(:, :, 1))
     554              :       !
     555            0 :       IF (grad) THEN
     556            0 :          DO i = 1, 3
     557            0 :             gradient(i, :, 2) = MATMUL(qdr(i, :, :), dEdq(:))
     558            0 :             stress(i, :, 2) = MATMUL(qdL(i, :, :), dEdq(:))
     559              :          END DO
     560              :       END IF
     561              :       !
     562            0 :       IF (doabc) THEN
     563            0 :          ALLOCATE (q(natom), qq(natom))
     564            0 :          q(:) = 0.0_dp; qq(:) = 0.0_dp
     565            0 :          ALLOCATE (enerd3(natom))
     566            0 :          enerd3(:) = 0.0_dp
     567            0 :          CALL disp%weight_references(mol, cnd, q, gwvec, gwdcn, gwdq)
     568            0 :          CALL disp%get_atomic_c6(mol, gwvec, gwdcn, gwdq, c6, dc6dcn, dc6dq)
     569            0 :          CALL get_lattice_points(mol%periodic, mol%lattice, cutoff%disp3, lattr)
     570              :          CALL param%get_dispersion3(mol, lattr, cutoff%disp3, disp%r4r2, c6, dc6dcn, dc6dq, &
     571            0 :                                     enerd3, dEdcn, qq, gradient(:, :, 3), stress(:, :, 3))
     572              :       END IF
     573            0 :       IF (grad) THEN
     574            0 :          DO i = 1, 3
     575            0 :             gradient(i, :, 4) = MATMUL(cndr(i, :, :), dEdcn(:))
     576            0 :             stress(i, :, 4) = MATMUL(cndL(i, :, :), dEdcn(:))
     577              :          END DO
     578              :       END IF
     579              : 
     580            0 :    END SUBROUTINE refd4_debug
     581              : 
     582              : #else
     583              : 
     584              : ! **************************************************************************************************
     585              : !> \brief ...
     586              : !> \param qs_env ...
     587              : !> \param dispersion_env ...
     588              : !> \param evdw ...
     589              : !> \param calculate_forces ...
     590              : !> \param iw ...
     591              : !> \param atomic_energy ...
     592              : ! **************************************************************************************************
     593              :    SUBROUTINE calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
     594              :                                               iw, atomic_energy)
     595              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     596              :       TYPE(qs_dispersion_type), INTENT(IN), POINTER      :: dispersion_env
     597              :       REAL(KIND=dp), INTENT(INOUT)                       :: evdw
     598              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     599              :       INTEGER, INTENT(IN)                                :: iw
     600              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atomic_energy
     601              : 
     602              :       MARK_USED(qs_env)
     603              :       MARK_USED(dispersion_env)
     604              :       MARK_USED(evdw)
     605              :       MARK_USED(calculate_forces)
     606              :       MARK_USED(iw)
     607              :       MARK_USED(atomic_energy)
     608              : 
     609              :       CPABORT("CP2K build without DFTD4")
     610              : 
     611              :    END SUBROUTINE calculate_dispersion_d4_pairpot
     612              : 
     613              : #endif
     614              : 
     615              : ! **************************************************************************************************
     616              : !> \brief ...
     617              : !> \param dispersion_env ...
     618              : !> \param cutoff ...
     619              : !> \param r4r2 ...
     620              : !> \param gwvec ...
     621              : !> \param gwdcn ...
     622              : !> \param gwdq ...
     623              : !> \param c6ref ...
     624              : !> \param mrefs ...
     625              : !> \param energies ...
     626              : !> \param dEdcn ...
     627              : !> \param dEdq ...
     628              : !> \param calculate_forces ...
     629              : !> \param gradient ...
     630              : !> \param stress ...
     631              : ! **************************************************************************************************
     632          872 :    SUBROUTINE dispersion_2b(dispersion_env, cutoff, r4r2, &
     633         2580 :                             gwvec, gwdcn, gwdq, c6ref, mrefs, &
     634         1744 :                             energies, dEdcn, dEdq, &
     635         1726 :                             calculate_forces, gradient, stress)
     636              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     637              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     638              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: r4r2
     639              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: gwvec, gwdcn, gwdq
     640              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
     641              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
     642              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: energies, dEdcn, dEdq
     643              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     644              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: gradient, stress
     645              : 
     646              :       INTEGER                                            :: iatom, ikind, jatom, jkind, mepos, num_pe
     647              :       REAL(KINd=dp)                                      :: a1, a2, c6ij, cutoff2, d6, d8, dE, dr2, &
     648              :                                                             edisp, fac, gdisp, r0ij, rrij, s6, s8, &
     649              :                                                             t6, t8
     650              :       REAL(KINd=dp), DIMENSION(2)                        :: dcdcn, dcdq
     651              :       REAL(KINd=dp), DIMENSION(3)                        :: dG, rij
     652              :       REAL(KINd=dp), DIMENSION(3, 3)                     :: dS
     653              :       TYPE(neighbor_list_iterator_p_type), &
     654          872 :          DIMENSION(:), POINTER                           :: nl_iterator
     655              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     656              :          POINTER                                         :: sab_vdw
     657              : 
     658          872 :       a1 = dispersion_env%a1
     659          872 :       a2 = dispersion_env%a2
     660          872 :       s6 = dispersion_env%s6
     661          872 :       s8 = dispersion_env%s8
     662          872 :       cutoff2 = cutoff*cutoff
     663              : 
     664          872 :       sab_vdw => dispersion_env%sab_vdw
     665              : 
     666          872 :       num_pe = 1
     667          872 :       CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
     668              : 
     669          872 :       mepos = 0
     670       192688 :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     671              :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
     672       191816 :                                 iatom=iatom, jatom=jatom, r=rij)
     673              :          ! vdW potential
     674       767264 :          dr2 = SUM(rij(:)**2)
     675       192688 :          IF (dr2 <= cutoff2 .AND. dr2 > 0.0000001_dp) THEN
     676       190015 :             rrij = 3._dp*r4r2(ikind)*r4r2(jkind)
     677       190015 :             r0ij = a1*SQRT(rrij) + a2
     678       190015 :             IF (calculate_forces) THEN
     679              :                CALL get_c6derivs(c6ij, dcdcn, dcdq, iatom, jatom, ikind, jkind, &
     680        91168 :                                  gwvec, gwdcn, gwdq, c6ref, mrefs)
     681              :             ELSE
     682        98847 :                CALL get_c6value(c6ij, iatom, jatom, ikind, jkind, gwvec, c6ref, mrefs)
     683              :             END IF
     684       190015 :             fac = 1._dp
     685       190015 :             IF (iatom == jatom) fac = 0.5_dp
     686       190015 :             t6 = 1.0_dp/(dr2**3 + r0ij**6)
     687       190015 :             t8 = 1.0_dp/(dr2**4 + r0ij**8)
     688              : 
     689       190015 :             edisp = (s6*t6 + s8*rrij*t8)*fac
     690       190015 :             dE = -c6ij*edisp
     691       190015 :             energies(iatom) = energies(iatom) + dE*0.5_dp
     692       190015 :             energies(jatom) = energies(jatom) + dE*0.5_dp
     693              : 
     694       190015 :             IF (calculate_forces) THEN
     695        91168 :                d6 = -6.0_dp*dr2**2*t6**2
     696        91168 :                d8 = -8.0_dp*dr2**3*t8**2
     697        91168 :                gdisp = (s6*d6 + s8*rrij*d8)*fac
     698       364672 :                dG(:) = -c6ij*gdisp*rij(:)
     699       364672 :                gradient(:, iatom) = gradient(:, iatom) - dG
     700       364672 :                gradient(:, jatom) = gradient(:, jatom) + dG
     701      1185184 :                dS(:, :) = SPREAD(dG, 1, 3)*SPREAD(rij, 2, 3)
     702      1185184 :                stress(:, :) = stress(:, :) + dS(:, :)
     703        91168 :                dEdcn(iatom) = dEdcn(iatom) - dcdcn(1)*edisp
     704        91168 :                dEdq(iatom) = dEdq(iatom) - dcdq(1)*edisp
     705        91168 :                dEdcn(jatom) = dEdcn(jatom) - dcdcn(2)*edisp
     706        91168 :                dEdq(jatom) = dEdq(jatom) - dcdq(2)*edisp
     707              :             END IF
     708              :          END IF
     709              :       END DO
     710              : 
     711          872 :       CALL neighbor_list_iterator_release(nl_iterator)
     712              : 
     713          872 :    END SUBROUTINE dispersion_2b
     714              : 
     715              : ! **************************************************************************************************
     716              : !> \brief ...
     717              : !> \param qs_env ...
     718              : !> \param dispersion_env ...
     719              : !> \param tvec ...
     720              : !> \param cutoff ...
     721              : !> \param r4r2 ...
     722              : !> \param gwvec ...
     723              : !> \param gwdcn ...
     724              : !> \param gwdq ...
     725              : !> \param c6ref ...
     726              : !> \param mrefs ...
     727              : !> \param energies ...
     728              : !> \param dEdcn ...
     729              : !> \param dEdq ...
     730              : !> \param calculate_forces ...
     731              : !> \param gradient ...
     732              : !> \param stress ...
     733              : ! **************************************************************************************************
     734           14 :    SUBROUTINE dispersion_3b(qs_env, dispersion_env, tvec, cutoff, r4r2, &
     735           30 :                             gwvec, gwdcn, gwdq, c6ref, mrefs, &
     736           28 :                             energies, dEdcn, dEdq, &
     737           22 :                             calculate_forces, gradient, stress)
     738              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     739              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     740              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: tvec
     741              :       REAL(KIND=dp), INTENT(IN)                          :: cutoff
     742              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: r4r2
     743              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: gwvec, gwdcn, gwdq
     744              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
     745              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
     746              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: energies, dEdcn, dEdq
     747              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     748              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: gradient, stress
     749              : 
     750              :       INTEGER                                            :: iatom, ikind, jatom, jkind, katom, &
     751              :                                                             kkind, ktr, mepos, natom, num_pe
     752           14 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     753              :       INTEGER, DIMENSION(3)                              :: cell_b
     754              :       REAL(KINd=dp)                                      :: a1, a2, alp, ang, c6ij, c6ik, c6jk, c9, &
     755              :                                                             cutoff2, dang, dE, dfdmp, fac, fdmp, &
     756              :                                                             r0, r0ij, r0ik, r0jk, r1, r2, r2ij, &
     757              :                                                             r2ik, r2jk, r3, r5, rr, s6, s8, s9
     758           14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: rcpbc
     759              :       REAL(KINd=dp), DIMENSION(2)                        :: dc6dcnij, dc6dcnik, dc6dcnjk, dc6dqij, &
     760              :                                                             dc6dqik, dc6dqjk
     761              :       REAL(KINd=dp), DIMENSION(3)                        :: dGij, dGik, dGjk, ra, rb, rb0, rij, vij, &
     762              :                                                             vik, vjk
     763              :       REAL(KINd=dp), DIMENSION(3, 3)                     :: dS
     764           14 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     765              :       TYPE(cell_type), POINTER                           :: cell
     766              :       TYPE(neighbor_list_iterator_p_type), &
     767           14 :          DIMENSION(:), POINTER                           :: nl_iterator
     768              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     769           14 :          POINTER                                         :: sab_vdw
     770           14 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     771              : 
     772              :       CALL get_qs_env(qs_env=qs_env, natom=natom, cell=cell, &
     773           14 :                       atomic_kind_set=atomic_kind_set, particle_set=particle_set)
     774              : 
     775           42 :       ALLOCATE (rcpbc(3, natom))
     776          154 :       DO iatom = 1, natom
     777          154 :          rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
     778              :       END DO
     779           14 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     780              : 
     781           14 :       a1 = dispersion_env%a1
     782           14 :       a2 = dispersion_env%a2
     783           14 :       s6 = dispersion_env%s6
     784           14 :       s8 = dispersion_env%s8
     785           14 :       s9 = dispersion_env%s9
     786           14 :       alp = dispersion_env%alp
     787              : 
     788           14 :       cutoff2 = cutoff**2
     789              : 
     790           14 :       sab_vdw => dispersion_env%sab_vdw
     791              : 
     792           14 :       num_pe = 1
     793           14 :       CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
     794              : 
     795           14 :       mepos = 0
     796       129748 :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     797       129734 :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
     798              : 
     799       518936 :          r2ij = SUM(rij(:)**2)
     800       129734 :          IF (calculate_forces) THEN
     801              :             CALL get_c6derivs(c6ij, dc6dcnij, dc6dqij, iatom, jatom, ikind, jkind, &
     802        88868 :                               gwvec, gwdcn, gwdq, c6ref, mrefs)
     803              :          ELSE
     804        40866 :             CALL get_c6value(c6ij, iatom, jatom, ikind, jkind, gwvec, c6ref, mrefs)
     805              :          END IF
     806       129734 :          r0ij = a1*SQRT(3._dp*r4r2(jkind)*r4r2(ikind)) + a2
     807       129748 :          IF (r2ij <= cutoff2 .AND. r2ij > EPSILON(1._dp)) THEN
     808        25564 :             CALL get_iterator_info(nl_iterator, cell=cell_b)
     809       409024 :             rb0(:) = MATMUL(cell%hmat, cell_b)
     810       102256 :             ra(:) = rcpbc(:, iatom)
     811       102256 :             rb(:) = rcpbc(:, jatom) + rb0
     812       102256 :             vij(:) = rb(:) - ra(:)
     813              : 
     814       127741 :             DO katom = 1, MIN(iatom, jatom)
     815       102177 :                kkind = kind_of(katom)
     816       102177 :                IF (calculate_forces) THEN
     817              :                   CALL get_c6derivs(c6ik, dc6dcnik, dc6dqik, katom, iatom, kkind, ikind, &
     818        88383 :                                     gwvec, gwdcn, gwdq, c6ref, mrefs)
     819              :                   CALL get_c6derivs(c6jk, dc6dcnjk, dc6dqjk, katom, jatom, kkind, jkind, &
     820        88383 :                                     gwvec, gwdcn, gwdq, c6ref, mrefs)
     821              :                ELSE
     822        13794 :                   CALL get_c6value(c6ik, katom, iatom, kkind, ikind, gwvec, c6ref, mrefs)
     823        13794 :                   CALL get_c6value(c6jk, katom, jatom, kkind, jkind, gwvec, c6ref, mrefs)
     824              :                END IF
     825       102177 :                c9 = -s9*SQRT(ABS(c6ij*c6ik*c6jk))
     826       102177 :                r0ik = a1*SQRT(3._dp*r4r2(kkind)*r4r2(ikind)) + a2
     827       102177 :                r0jk = a1*SQRT(3._dp*r4r2(kkind)*r4r2(jkind)) + a2
     828       102177 :                r0 = r0ij*r0ik*r0jk
     829       102177 :                fac = triple_scale(iatom, jatom, katom)
     830    111110602 :                DO ktr = 1, SIZE(tvec, 2)
     831    443931444 :                   vik(:) = rcpbc(:, katom) + tvec(:, ktr) - rcpbc(:, iatom)
     832    110982861 :                   r2ik = vik(1)*vik(1) + vik(2)*vik(2) + vik(3)*vik(3)
     833    110982861 :                   IF (r2ik > cutoff2 .OR. r2ik < EPSILON(1.0_dp)) CYCLE
     834    123226284 :                   vjk(:) = rcpbc(:, katom) + tvec(:, ktr) - rb(:)
     835     30806571 :                   r2jk = vjk(1)*vjk(1) + vjk(2)*vjk(2) + vjk(3)*vjk(3)
     836     30806571 :                   IF (r2jk > cutoff2 .OR. r2jk < EPSILON(1.0_dp)) CYCLE
     837     14392504 :                   r2 = r2ij*r2ik*r2jk
     838     14392504 :                   r1 = SQRT(r2)
     839     14392504 :                   r3 = r2*r1
     840     14392504 :                   r5 = r3*r2
     841              : 
     842     14392504 :                   fdmp = 1.0_dp/(1.0_dp + 6.0_dp*(r0/r1)**(alp/3.0_dp))
     843              :                   ang = 0.375_dp*(r2ij + r2jk - r2ik)*(r2ij - r2jk + r2ik)* &
     844     14392504 :                         (-r2ij + r2jk + r2ik)/r5 + 1.0_dp/r3
     845              : 
     846     14392504 :                   rr = ang*fdmp
     847     14392504 :                   dE = rr*c9*fac
     848     14392504 :                   energies(iatom) = energies(iatom) - dE/3._dp
     849     14392504 :                   energies(jatom) = energies(jatom) - dE/3._dp
     850     14392504 :                   energies(katom) = energies(katom) - dE/3._dp
     851              : 
     852     14494681 :                   IF (calculate_forces) THEN
     853              : 
     854     14199360 :                      dfdmp = -2.0_dp*alp*(r0/r1)**(alp/3.0_dp)*fdmp**2
     855              : 
     856              :                      ! d/drij
     857              :                      dang = -0.375_dp*(r2ij**3 + r2ij**2*(r2jk + r2ik) &
     858              :                                        + r2ij*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ik &
     859              :                                                + 3.0_dp*r2ik**2) &
     860     14199360 :                                        - 5.0_dp*(r2jk - r2ik)**2*(r2jk + r2ik))/r5
     861     56797440 :                      dGij(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ij*vij
     862              : 
     863              :                      ! d/drik
     864              :                      dang = -0.375_dp*(r2ik**3 + r2ik**2*(r2jk + r2ij) &
     865              :                                        + r2ik*(3.0_dp*r2jk**2 + 2.0_dp*r2jk*r2ij &
     866              :                                                + 3.0_dp*r2ij**2) &
     867     14199360 :                                        - 5.0_dp*(r2jk - r2ij)**2*(r2jk + r2ij))/r5
     868     56797440 :                      dGik(:) = c9*(-dang*fdmp + ang*dfdmp)/r2ik*vik
     869              : 
     870              :                      ! d/drjk
     871              :                      dang = -0.375_dp*(r2jk**3 + r2jk**2*(r2ik + r2ij) &
     872              :                                        + r2jk*(3.0_dp*r2ik**2 + 2.0_dp*r2ik*r2ij &
     873              :                                                + 3.0_dp*r2ij**2) &
     874     14199360 :                                        - 5.0_dp*(r2ik - r2ij)**2*(r2ik + r2ij))/r5
     875     56797440 :                      dGjk(:) = c9*(-dang*fdmp + ang*dfdmp)/r2jk*vjk
     876              : 
     877     56797440 :                      gradient(:, iatom) = gradient(:, iatom) - dGij - dGik
     878     56797440 :                      gradient(:, jatom) = gradient(:, jatom) + dGij - dGjk
     879     56797440 :                      gradient(:, katom) = gradient(:, katom) + dGik + dGjk
     880              : 
     881              :                      dS(:, :) = SPREAD(dGij, 1, 3)*SPREAD(vij, 2, 3) &
     882              :                                 + SPREAD(dGik, 1, 3)*SPREAD(vik, 2, 3) &
     883    184591680 :                                 + SPREAD(dGjk, 1, 3)*SPREAD(vjk, 2, 3)
     884              : 
     885    184591680 :                      stress(:, :) = stress + dS*fac
     886              : 
     887              :                      dEdcn(iatom) = dEdcn(iatom) - dE*0.5_dp &
     888     14199360 :                                     *(dc6dcnij(1)/c6ij + dc6dcnik(2)/c6ik)
     889              :                      dEdcn(jatom) = dEdcn(jatom) - dE*0.5_dp &
     890     14199360 :                                     *(dc6dcnij(2)/c6ij + dc6dcnjk(2)/c6jk)
     891              :                      dEdcn(katom) = dEdcn(katom) - dE*0.5_dp &
     892     14199360 :                                     *(dc6dcnik(1)/c6ik + dc6dcnjk(1)/c6jk)
     893              : 
     894              :                      dEdq(iatom) = dEdq(iatom) - dE*0.5_dp &
     895     14199360 :                                    *(dc6dqij(1)/c6ij + dc6dqik(2)/c6ik)
     896              :                      dEdq(jatom) = dEdq(jatom) - dE*0.5_dp &
     897     14199360 :                                    *(dc6dqij(2)/c6ij + dc6dqjk(2)/c6jk)
     898              :                      dEdq(katom) = dEdq(katom) - dE*0.5_dp &
     899     14199360 :                                    *(dc6dqik(1)/c6ik + dc6dqjk(1)/c6jk)
     900              : 
     901              :                   END IF
     902              : 
     903              :                END DO
     904              :             END DO
     905              :          END IF
     906              :       END DO
     907              : 
     908           14 :       CALL neighbor_list_iterator_release(nl_iterator)
     909              : 
     910           14 :       DEALLOCATE (rcpbc)
     911              : 
     912           28 :    END SUBROUTINE dispersion_3b
     913              : 
     914              : ! **************************************************************************************************
     915              : !> \brief ...
     916              : !> \param ii ...
     917              : !> \param jj ...
     918              : !> \param kk ...
     919              : !> \return ...
     920              : ! **************************************************************************************************
     921       102177 :    FUNCTION triple_scale(ii, jj, kk) RESULT(triple)
     922              :       INTEGER, INTENT(IN)                                :: ii, jj, kk
     923              :       REAL(KIND=dp)                                      :: triple
     924              : 
     925       102177 :       IF (ii == jj) THEN
     926        25081 :          IF (ii == kk) THEN
     927              :             ! ii'i" -> 1/6
     928              :             triple = 1.0_dp/6.0_dp
     929              :          ELSE
     930              :             ! ii'j -> 1/2
     931        20517 :             triple = 0.5_dp
     932              :          END IF
     933              :       ELSE
     934        77096 :          IF (ii /= kk .AND. jj /= kk) THEN
     935              :             ! ijk -> 1 (full)
     936              :             triple = 1.0_dp
     937              :          ELSE
     938              :             ! ijj' and iji' -> 1/2
     939        21000 :             triple = 0.5_dp
     940              :          END IF
     941              :       END IF
     942              : 
     943       102177 :    END FUNCTION triple_scale
     944              : 
     945              : ! **************************************************************************************************
     946              : !> \brief ...
     947              : !> \param qs_env ...
     948              : !> \param dEdcn ...
     949              : !> \param dcnum ...
     950              : !> \param gradient ...
     951              : !> \param stress ...
     952              : ! **************************************************************************************************
     953           18 :    SUBROUTINE dEdcn_force(qs_env, dEdcn, dcnum, gradient, stress)
     954              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     955              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: dEdcn
     956              :       TYPE(dcnum_type), DIMENSION(:), INTENT(IN)         :: dcnum
     957              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: gradient
     958              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(INOUT)      :: stress
     959              : 
     960              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'dEdcn_force'
     961              : 
     962              :       INTEGER                                            :: handle, i, ia, iatom, ikind, katom, &
     963              :                                                             natom, nkind
     964              :       LOGICAL                                            :: use_virial
     965              :       REAL(KIND=dp)                                      :: drk
     966              :       REAL(KIND=dp), DIMENSION(3)                        :: fdik, rik
     967              :       TYPE(distribution_1d_type), POINTER                :: local_particles
     968              :       TYPE(virial_type), POINTER                         :: virial
     969              : 
     970           18 :       CALL timeset(routineN, handle)
     971              : 
     972              :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom, &
     973              :                       local_particles=local_particles, &
     974           18 :                       virial=virial)
     975           18 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     976              : 
     977           66 :       DO ikind = 1, nkind
     978          120 :          DO ia = 1, local_particles%n_el(ikind)
     979           54 :             iatom = local_particles%list(ikind)%array(ia)
     980          186 :             DO i = 1, dcnum(iatom)%neighbors
     981           84 :                katom = dcnum(iatom)%nlist(i)
     982          336 :                rik = dcnum(iatom)%rik(:, i)
     983          336 :                drk = SQRT(SUM(rik(:)**2))
     984          336 :                fdik(:) = -(dEdcn(iatom) + dEdcn(katom))*dcnum(iatom)%dvals(i)*rik(:)/drk
     985          336 :                gradient(:, iatom) = gradient(:, iatom) + fdik(:)
     986          138 :                IF (use_virial) THEN
     987           22 :                   CALL virial_pair_force(stress, -0.5_dp, fdik, rik)
     988              :                END IF
     989              :             END DO
     990              :          END DO
     991              :       END DO
     992              : 
     993           18 :       CALL timestop(handle)
     994              : 
     995           18 :    END SUBROUTINE dEdcn_force
     996              : 
     997              : ! **************************************************************************************************
     998              : !> \brief ...
     999              : !> \param c6ij ...
    1000              : !> \param ia ...
    1001              : !> \param ja ...
    1002              : !> \param ik ...
    1003              : !> \param jk ...
    1004              : !> \param gwvec ...
    1005              : !> \param c6ref ...
    1006              : !> \param mrefs ...
    1007              : ! **************************************************************************************************
    1008       167301 :    SUBROUTINE get_c6value(c6ij, ia, ja, ik, jk, gwvec, c6ref, mrefs)
    1009              :       REAL(KIND=dp), INTENT(OUT)                         :: c6ij
    1010              :       INTEGER, INTENT(IN)                                :: ia, ja, ik, jk
    1011              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: gwvec
    1012              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
    1013              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
    1014              : 
    1015              :       INTEGER                                            :: iref, jref
    1016              :       REAL(KIND=dp)                                      :: refc6
    1017              : 
    1018       167301 :       c6ij = 0.0_dp
    1019       685103 :       DO jref = 1, mrefs(jk)
    1020      2606798 :          DO iref = 1, mrefs(ik)
    1021      1921695 :             refc6 = c6ref(iref, jref, ik, jk)
    1022      2439497 :             c6ij = c6ij + gwvec(iref, ia)*gwvec(jref, ja)*refc6
    1023              :          END DO
    1024              :       END DO
    1025              : 
    1026       167301 :    END SUBROUTINE get_c6value
    1027              : 
    1028              : ! **************************************************************************************************
    1029              : !> \brief ...
    1030              : !> \param c6ij ...
    1031              : !> \param dc6dcn ...
    1032              : !> \param dc6dq ...
    1033              : !> \param ia ...
    1034              : !> \param ja ...
    1035              : !> \param ik ...
    1036              : !> \param jk ...
    1037              : !> \param gwvec ...
    1038              : !> \param gwdcn ...
    1039              : !> \param gwdq ...
    1040              : !> \param c6ref ...
    1041              : !> \param mrefs ...
    1042              : ! **************************************************************************************************
    1043       356802 :    SUBROUTINE get_c6derivs(c6ij, dc6dcn, dc6dq, ia, ja, ik, jk, &
    1044       356802 :                            gwvec, gwdcn, gwdq, c6ref, mrefs)
    1045              :       REAL(KIND=dp), INTENT(OUT)                         :: c6ij
    1046              :       REAL(KIND=dp), DIMENSION(2), INTENT(OUT)           :: dc6dcn, dc6dq
    1047              :       INTEGER, INTENT(IN)                                :: ia, ja, ik, jk
    1048              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: gwvec, gwdcn, gwdq
    1049              :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: c6ref
    1050              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: mrefs
    1051              : 
    1052              :       INTEGER                                            :: iref, jref
    1053              :       REAL(KIND=dp)                                      :: refc6
    1054              : 
    1055       356802 :       c6ij = 0.0_dp
    1056       356802 :       dc6dcn = 0.0_dp
    1057       356802 :       dc6dq = 0.0_dp
    1058      1315528 :       DO jref = 1, mrefs(jk)
    1059      4969737 :          DO iref = 1, mrefs(ik)
    1060      3654209 :             refc6 = c6ref(iref, jref, ik, jk)
    1061      3654209 :             c6ij = c6ij + gwvec(iref, ia)*gwvec(jref, ja)*refc6
    1062      3654209 :             dc6dcn(1) = dc6dcn(1) + gwdcn(iref, ia)*gwvec(jref, ja)*refc6
    1063      3654209 :             dc6dcn(2) = dc6dcn(2) + gwvec(iref, ia)*gwdcn(jref, ja)*refc6
    1064      3654209 :             dc6dq(1) = dc6dq(1) + gwdq(iref, ia)*gwvec(jref, ja)*refc6
    1065      4612935 :             dc6dq(2) = dc6dq(2) + gwvec(iref, ia)*gwdq(jref, ja)*refc6
    1066              :          END DO
    1067              :       END DO
    1068              : 
    1069       356802 :    END SUBROUTINE get_c6derivs
    1070              : 
    1071              : ! **************************************************************************************************
    1072              : !> \brief ...
    1073              : !> \param ga ...
    1074              : !> \param gd ...
    1075              : !> \param ev1 ...
    1076              : !> \param ev2 ...
    1077              : !> \param ev3 ...
    1078              : !> \param ev4 ...
    1079              : ! **************************************************************************************************
    1080            0 :    SUBROUTINE gerror(ga, gd, ev1, ev2, ev3, ev4)
    1081              :       REAL(KIND=dp), DIMENSION(:, :)                     :: ga, gd
    1082              :       REAL(KIND=dp), INTENT(OUT)                         :: ev1, ev2, ev3, ev4
    1083              : 
    1084              :       INTEGER                                            :: na, np(2)
    1085              : 
    1086            0 :       na = SIZE(ga, 2)
    1087            0 :       ev1 = SQRT(SUM((gd - ga)**2)/na)
    1088            0 :       ev2 = ev1/SQRT(SUM(gd**2)/na)*100._dp
    1089            0 :       np = MAXLOC(ABS(gd - ga))
    1090            0 :       ev3 = ABS(gd(np(1), np(2)) - ga(np(1), np(2)))
    1091            0 :       ev4 = ABS(gd(np(1), np(2)))
    1092            0 :       IF (ev4 > 1.E-6) THEN
    1093            0 :          ev4 = ev3/ev4*100._dp
    1094              :       ELSE
    1095            0 :          ev4 = 100.0_dp
    1096              :       END IF
    1097              : 
    1098            0 :    END SUBROUTINE gerror
    1099              : 
    1100              : ! **************************************************************************************************
    1101              : !> \brief ...
    1102              : !> \param sa ...
    1103              : !> \param sd ...
    1104              : !> \param ev1 ...
    1105              : !> \param ev2 ...
    1106              : ! **************************************************************************************************
    1107            0 :    SUBROUTINE serror(sa, sd, ev1, ev2)
    1108              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: sa, sd
    1109              :       REAL(KIND=dp), INTENT(OUT)                         :: ev1, ev2
    1110              : 
    1111              :       INTEGER                                            :: i, j
    1112              :       REAL(KIND=dp)                                      :: rel
    1113              : 
    1114            0 :       ev1 = MAXVAL(ABS(sd - sa))
    1115            0 :       ev2 = 0.0_dp
    1116            0 :       DO i = 1, 3
    1117            0 :          DO j = 1, 3
    1118            0 :             IF (ABS(sd(i, j)) > 1.E-6_dp) THEN
    1119            0 :                rel = ABS(sd(i, j) - sa(i, j))/ABS(sd(i, j))*100._dp
    1120            0 :                ev2 = MAX(ev2, rel)
    1121              :             END IF
    1122              :          END DO
    1123              :       END DO
    1124              : 
    1125            0 :    END SUBROUTINE serror
    1126              : 
    1127              : ! **************************************************************************************************
    1128              : !> \brief ...
    1129              : !> \param va ...
    1130              : !> \param vd ...
    1131              : !> \param ev1 ...
    1132              : !> \param ev2 ...
    1133              : ! **************************************************************************************************
    1134            0 :    SUBROUTINE verror(va, vd, ev1, ev2)
    1135              :       REAL(KIND=dp), DIMENSION(:)                        :: va, vd
    1136              :       REAL(KIND=dp), INTENT(OUT)                         :: ev1, ev2
    1137              : 
    1138              :       INTEGER                                            :: i, na
    1139              :       REAL(KIND=dp)                                      :: rel
    1140              : 
    1141            0 :       na = SIZE(va)
    1142            0 :       ev1 = MAXVAL(ABS(vd - va))
    1143            0 :       ev2 = 0.0_dp
    1144            0 :       DO i = 1, na
    1145            0 :          IF (ABS(vd(i)) > 1.E-8_dp) THEN
    1146            0 :             rel = ABS(vd(i) - va(i))/ABS(vd(i))*100._dp
    1147            0 :             ev2 = MAX(ev2, rel)
    1148              :          END IF
    1149              :       END DO
    1150              : 
    1151            0 :    END SUBROUTINE verror
    1152              : 
    1153          894 : END MODULE qs_dispersion_d4
        

Generated by: LCOV version 2.0-1