LCOV - code coverage report
Current view: top level - src - qs_dispersion_d3.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.2 % 639 634
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 5 5

            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 D3 dispersion
      10              : !> \author JGH
      11              : ! **************************************************************************************************
      12              : MODULE qs_dispersion_d3
      13              : 
      14              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15              :                                               get_atomic_kind,&
      16              :                                               get_atomic_kind_set
      17              :    USE atprop_types,                    ONLY: atprop_array_init,&
      18              :                                               atprop_type
      19              :    USE cell_types,                      ONLY: cell_type,&
      20              :                                               get_cell,&
      21              :                                               pbc,&
      22              :                                               plane_distance
      23              :    USE cp_files,                        ONLY: close_file,&
      24              :                                               open_file
      25              :    USE input_constants,                 ONLY: vdw_pairpot_dftd3,&
      26              :                                               vdw_pairpot_dftd3bj
      27              :    USE kinds,                           ONLY: dp
      28              :    USE mathconstants,                   ONLY: twopi
      29              :    USE message_passing,                 ONLY: mp_para_env_type
      30              :    USE molecule_types,                  ONLY: molecule_type
      31              :    USE particle_types,                  ONLY: particle_type
      32              :    USE physcon,                         ONLY: kcalmol
      33              :    USE qs_dispersion_cnum,              ONLY: d3_cnumber,&
      34              :                                               dcnum_distribute,&
      35              :                                               dcnum_type,&
      36              :                                               exclude_d3_kind_pair
      37              :    USE qs_dispersion_types,             ONLY: dftd3_pp,&
      38              :                                               qs_atom_dispersion_type,&
      39              :                                               qs_dispersion_type
      40              :    USE qs_dispersion_utils,             ONLY: cellhash
      41              :    USE qs_environment_types,            ONLY: get_qs_env,&
      42              :                                               qs_environment_type
      43              :    USE qs_force_types,                  ONLY: qs_force_type
      44              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      45              :                                               qs_kind_type
      46              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      47              :                                               neighbor_list_iterate,&
      48              :                                               neighbor_list_iterator_create,&
      49              :                                               neighbor_list_iterator_p_type,&
      50              :                                               neighbor_list_iterator_release,&
      51              :                                               neighbor_list_set_p_type
      52              :    USE virial_methods,                  ONLY: virial_pair_force
      53              :    USE virial_types,                    ONLY: virial_type
      54              : #include "./base/base_uses.f90"
      55              : 
      56              :    IMPLICIT NONE
      57              : 
      58              :    PRIVATE
      59              : 
      60              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d3'
      61              : 
      62              :    PUBLIC :: calculate_dispersion_d3_pairpot, dftd3_c6_param
      63              : 
      64              : ! **************************************************************************************************
      65              : 
      66              : CONTAINS
      67              : 
      68              : ! **************************************************************************************************
      69              : !> \brief ...
      70              : !> \param qs_env ...
      71              : !> \param dispersion_env ...
      72              : !> \param evdw ...
      73              : !> \param calculate_forces ...
      74              : !> \param unit_nr ...
      75              : !> \param atevdw ...
      76              : ! **************************************************************************************************
      77         4724 :    SUBROUTINE calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
      78         4724 :                                               unit_nr, atevdw)
      79              : 
      80              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      81              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
      82              :       REAL(KIND=dp), INTENT(OUT)                         :: evdw
      83              :       LOGICAL, INTENT(IN)                                :: calculate_forces
      84              :       INTEGER, INTENT(IN)                                :: unit_nr
      85              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atevdw
      86              : 
      87              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d3_pairpot'
      88              : 
      89              :       INTEGER :: atom_a, atom_b, atom_c, atom_d, handle, hashb, hashc, i, ia, iat, iatom, icx, &
      90              :          icy, icz, idmp, ikind, ilist, imol, jatom, jkind, katom, kkind, kstart, latom, lkind, &
      91              :          max_elem, maxc, mepos, na, natom, nb, nc, nkind, num_pe, za, zb, zc
      92         4724 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, atomnumber, kind_of
      93              :       INTEGER, DIMENSION(3)                              :: cell_b, cell_c, ncell, periodic
      94         4724 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
      95              :       LOGICAL :: atenergy, atex, debugall, domol, exclude_pair, floating_a, floating_b, &
      96              :          floating_c, ghost_a, ghost_b, ghost_c, is000, use_virial
      97         4724 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: dodisp, exclude
      98              :       REAL(KIND=dp) :: a1, a2, alp6, alp8, ang, c6, c8, c9, cc6ab, cc6bc, cc6ca, cnum, dc6a, dc6b, &
      99              :          dc8a, dc8b, dcc6aba, dcc6abb, dcc6bcb, dcc6bcc, dcc6caa, dcc6cac, de6, de8, de91, de921, &
     100              :          de922, dea, dfdab6, dfdab8, dfdabc, dr, drk, e6, e6tot, e8, e8tot, e9, e9tot, elrc6, &
     101              :          elrc8, elrc9, eps_cn, esrb, f0ab, fac, fac0, fdab6, fdab8, fdabc, gsrb, kgc8, nab, nabc, &
     102              :          r0, r0ab, r2ab, r2abc, r2bc, r2ca, r6, r8, rabc, rc2, rcc, rcut, rs6, rs8, s1, s2, s3, &
     103              :          s6, s8, s8i, s9, srbe, ssrb, t1srb, t2srb
     104         4724 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: atom2mol, c6d2, cnkind, cnumbers, &
     105         4724 :                                                             cnumfix, radd2
     106         4724 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: rcpbc
     107              :       REAL(KIND=dp), DIMENSION(3)                        :: fdij, fdik, ra, rab, rb, rb0, rbc, rc, &
     108              :                                                             rc0, rca, rij, rik, sab_max
     109              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_virial_thread
     110         4724 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: atener
     111         4724 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     112              :       TYPE(atprop_type), POINTER                         :: atprop
     113              :       TYPE(cell_type), POINTER                           :: cell
     114         4724 :       TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:)        :: dcnum
     115         4724 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     116              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     117              :       TYPE(neighbor_list_iterator_p_type), &
     118         4724 :          DIMENSION(:), POINTER                           :: nl_iterator
     119              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     120         4724 :          POINTER                                         :: sab_vdw
     121         4724 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     122              :       TYPE(qs_atom_dispersion_type), POINTER             :: disp_a, disp_b, disp_c
     123         4724 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     124         4724 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     125              :       TYPE(virial_type), POINTER                         :: virial
     126              : 
     127         4724 :       CALL timeset(routineN, handle)
     128              : 
     129         4724 :       evdw = 0._dp
     130              : 
     131         4724 :       NULLIFY (atomic_kind_set, qs_kind_set, sab_vdw)
     132              : 
     133              :       CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
     134         4724 :                       qs_kind_set=qs_kind_set, cell=cell, virial=virial, para_env=para_env, atprop=atprop)
     135              : 
     136         4724 :       debugall = dispersion_env%verbose
     137              : 
     138              :       ! atomic energy and stress arrays
     139         4724 :       atenergy = atprop%energy
     140         4724 :       IF (atenergy) THEN
     141           88 :          CALL atprop_array_init(atprop%atevdw, natom)
     142           88 :          atener => atprop%atevdw
     143              :       END IF
     144              :       ! external atomic energy
     145         4724 :       atex = .FALSE.
     146         4724 :       IF (PRESENT(atevdw)) THEN
     147            2 :          atex = .TRUE.
     148              :       END IF
     149              : 
     150         4724 :       NULLIFY (particle_set)
     151         4724 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
     152         4724 :       natom = SIZE(particle_set)
     153              : 
     154         4724 :       NULLIFY (force)
     155         4724 :       CALL get_qs_env(qs_env=qs_env, force=force)
     156         4724 :       CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     157         4724 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     158         4724 :       pv_virial_thread(:, :) = 0._dp
     159              : 
     160        37792 :       ALLOCATE (dodisp(nkind), exclude(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
     161        16306 :       DO ikind = 1, nkind
     162        11582 :          CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
     163        11582 :          CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
     164        11582 :          dodisp(ikind) = disp_a%defined
     165        11582 :          exclude(ikind) = ghost_a .OR. floating_a
     166        11582 :          atomnumber(ikind) = za
     167        11582 :          c6d2(ikind) = disp_a%c6
     168        27888 :          radd2(ikind) = disp_a%vdw_radii
     169              :       END DO
     170              : 
     171        14172 :       ALLOCATE (rcpbc(3, natom))
     172        49464 :       DO iatom = 1, natom
     173        49464 :          rcpbc(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
     174              :       END DO
     175              : 
     176         4724 :       rcut = 2._dp*dispersion_env%rc_disp
     177         4724 :       rc2 = rcut*rcut
     178              : 
     179         4724 :       maxc = SIZE(dispersion_env%c6ab, 3)
     180         4724 :       max_elem = SIZE(dispersion_env%c6ab, 1)
     181         4724 :       alp6 = dispersion_env%alp
     182         4724 :       alp8 = alp6 + 2._dp
     183         4724 :       s6 = dispersion_env%s6
     184         4724 :       s8 = dispersion_env%s8
     185         4724 :       s9 = dispersion_env%s6
     186         4724 :       rs6 = dispersion_env%sr6
     187         4724 :       rs8 = 1._dp
     188         4724 :       a1 = dispersion_env%a1
     189         4724 :       a2 = dispersion_env%a2
     190         4724 :       eps_cn = dispersion_env%eps_cn
     191         4724 :       e6tot = 0._dp
     192         4724 :       e8tot = 0._dp
     193         4724 :       e9tot = 0._dp
     194         4724 :       esrb = 0._dp
     195         4724 :       domol = dispersion_env%domol
     196              :       ! molecule correction
     197         4724 :       kgc8 = dispersion_env%kgc8
     198         4724 :       IF (domol) THEN
     199            2 :          CALL get_qs_env(qs_env=qs_env, molecule_set=molecule_set)
     200            6 :          ALLOCATE (atom2mol(natom))
     201            6 :          DO imol = 1, SIZE(molecule_set)
     202           16 :             DO iat = molecule_set(imol)%first_atom, molecule_set(imol)%last_atom
     203           14 :                atom2mol(iat) = imol
     204              :             END DO
     205              :          END DO
     206              :       END IF
     207              :       ! damping type
     208         4724 :       idmp = 0
     209         4724 :       IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
     210              :          idmp = 1
     211         4464 :       ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
     212         4464 :          idmp = 2
     213              :       END IF
     214              :       ! SRB parameters
     215         4724 :       ssrb = dispersion_env%srb_params(1)
     216         4724 :       gsrb = dispersion_env%srb_params(2)
     217         4724 :       t1srb = dispersion_env%srb_params(3)
     218         4724 :       t2srb = dispersion_env%srb_params(4)
     219              : 
     220         4724 :       IF (unit_nr > 0) THEN
     221            8 :          WRITE (unit_nr, *) " Scaling parameter (s6) ", s6
     222            8 :          WRITE (unit_nr, *) " Scaling parameter (s8) ", s8
     223            8 :          IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
     224            4 :             WRITE (unit_nr, *) " Zero Damping parameter (sr6)", rs6
     225            4 :             WRITE (unit_nr, *) " Zero Damping parameter (sr8)", rs8
     226            4 :          ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
     227            4 :             WRITE (unit_nr, *) " BJ Damping parameter (a1) ", a1
     228            4 :             WRITE (unit_nr, *) " BJ Damping parameter (a2) ", a2
     229              :          END IF
     230            8 :          WRITE (unit_nr, *) " Cutoff coordination numbers", eps_cn
     231            8 :          IF (dispersion_env%lrc) THEN
     232            1 :             WRITE (unit_nr, *) " Apply a long range correction"
     233              :          END IF
     234            8 :          IF (dispersion_env%srb) THEN
     235            2 :             WRITE (unit_nr, *) " Apply a short range bond correction"
     236            2 :             WRITE (unit_nr, *) " SRB parameters (s,g,t1,t2) ", ssrb, gsrb, t1srb, t2srb
     237              :          END IF
     238            8 :          IF (domol) THEN
     239            1 :             WRITE (unit_nr, *) " Inter-molecule scaling parameter (s8) ", kgc8
     240              :          END IF
     241              :       END IF
     242              :       ! Calculate coordination numbers
     243         4724 :       NULLIFY (particle_set)
     244         4724 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
     245         4724 :       natom = SIZE(particle_set)
     246        14172 :       ALLOCATE (cnumbers(natom))
     247        49464 :       cnumbers = 0._dp
     248              : 
     249         4724 :       IF (calculate_forces .OR. debugall) THEN
     250        11918 :          ALLOCATE (dcnum(natom))
     251        10394 :          dcnum(:)%neighbors = 0
     252        10394 :          DO iatom = 1, natom
     253        10394 :             ALLOCATE (dcnum(iatom)%nlist(10), dcnum(iatom)%dvals(10), dcnum(iatom)%rik(3, 10))
     254              :          END DO
     255              :       ELSE
     256         7924 :          ALLOCATE (dcnum(1))
     257              :       END IF
     258              : 
     259              :       CALL d3_cnumber(qs_env, dispersion_env, cnumbers, dcnum, exclude, atomnumber, &
     260         4724 :                       calculate_forces, 1)
     261              : 
     262         4724 :       CALL para_env%sum(cnumbers)
     263              :       ! for parallel runs we have to update dcnum on all processors
     264         4724 :       IF (calculate_forces .OR. debugall) THEN
     265          762 :          CALL dcnum_distribute(dcnum, para_env)
     266          762 :          IF (unit_nr > 0 .AND. SIZE(dcnum) > 0) THEN
     267            3 :             WRITE (unit_nr, *)
     268            3 :             WRITE (unit_nr, *) "  ATOM       Coordination   Neighbors"
     269           15 :             DO i = 1, natom
     270           15 :                WRITE (unit_nr, "(I6,F20.10,I12)") i, cnumbers(i), dcnum(i)%neighbors
     271              :             END DO
     272            3 :             WRITE (unit_nr, *)
     273              :          END IF
     274              :       END IF
     275              : 
     276         4724 :       nab = 0._dp
     277         4724 :       nabc = 0._dp
     278         4724 :       IF (dispersion_env%doabc) THEN
     279          104 :          rcc = 2._dp*dispersion_env%rc_disp
     280          104 :          CALL get_cell(cell=cell, periodic=periodic)
     281          104 :          sab_max(1) = rcc/plane_distance(1, 0, 0, cell)
     282          104 :          sab_max(2) = rcc/plane_distance(0, 1, 0, cell)
     283          104 :          sab_max(3) = rcc/plane_distance(0, 0, 1, cell)
     284          416 :          ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
     285          104 :          IF (unit_nr > 0) THEN
     286            3 :             WRITE (unit_nr, *) " Calculate C9 Terms"
     287            3 :             WRITE (unit_nr, "(A,T20,I3,A,I3)") "  Search in cells ", -ncell(1), ":", ncell(1)
     288            3 :             WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(2), ":", ncell(2)
     289            3 :             WRITE (unit_nr, "(T20,I3,A,I3)") - ncell(3), ":", ncell(3)
     290            3 :             WRITE (unit_nr, *)
     291              :          END IF
     292          104 :          IF (dispersion_env%c9cnst) THEN
     293           62 :             IF (unit_nr > 0) WRITE (unit_nr, *) " Use reference coordination numbers for C9 term"
     294          186 :             ALLOCATE (cnumfix(natom))
     295          306 :             cnumfix = 0._dp
     296              :             ! first use the default values
     297          306 :             DO iatom = 1, natom
     298          244 :                ikind = kind_of(iatom)
     299          244 :                CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
     300          306 :                cnumfix(iatom) = dispersion_env%cn(za)
     301              :             END DO
     302              :             ! now check for changes from default
     303           62 :             IF (ASSOCIATED(dispersion_env%cnkind)) THEN
     304           12 :                DO i = 1, SIZE(dispersion_env%cnkind)
     305            6 :                   ikind = dispersion_env%cnkind(i)%kind
     306            6 :                   cnum = dispersion_env%cnkind(i)%cnum
     307            6 :                   CPASSERT(ikind <= nkind)
     308            6 :                   CPASSERT(ikind > 0)
     309            6 :                   CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, atom_list=atom_list)
     310           32 :                   DO ia = 1, na
     311           20 :                      iatom = atom_list(ia)
     312           26 :                      cnumfix(iatom) = cnum
     313              :                   END DO
     314              :                END DO
     315              :             END IF
     316           62 :             IF (ASSOCIATED(dispersion_env%cnlist)) THEN
     317            6 :                DO i = 1, SIZE(dispersion_env%cnlist)
     318           14 :                   DO ilist = 1, dispersion_env%cnlist(i)%natom
     319            8 :                      iatom = dispersion_env%cnlist(i)%atom(ilist)
     320           12 :                      cnumfix(iatom) = dispersion_env%cnlist(i)%cnum
     321              :                   END DO
     322              :                END DO
     323              :             END IF
     324           62 :             IF (unit_nr > 0) THEN
     325            5 :                DO i = 1, natom
     326            5 :                   IF (ABS(cnumbers(i) - cnumfix(i)) > 0.5_dp) THEN
     327            0 :                      WRITE (unit_nr, "(A,T20,A,I6,T41,2F10.3)") "  Difference in CN ", "Atom:", &
     328            0 :                         i, cnumbers(i), cnumfix(i)
     329              :                   END IF
     330              :                END DO
     331            1 :                WRITE (unit_nr, *)
     332              :             END IF
     333              :          END IF
     334              :       END IF
     335              : 
     336         4724 :       sab_vdw => dispersion_env%sab_vdw
     337              : 
     338         4724 :       num_pe = 1
     339         4724 :       CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
     340              : 
     341         4724 :       mepos = 0
     342      8693554 :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     343      8688830 :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
     344              : 
     345      8688830 :          IF (exclude(ikind) .OR. exclude(jkind)) CYCLE
     346              : 
     347      8688773 :          IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
     348              : 
     349      8688539 :          za = atomnumber(ikind)
     350      8688539 :          zb = atomnumber(jkind)
     351              :          ! vdW potential
     352     34754156 :          dr = SQRT(SUM(rij(:)**2))
     353      8693263 :          IF (dr <= rcut) THEN
     354      8688539 :             nab = nab + 1._dp
     355      8688539 :             fac = 1._dp
     356      8688539 :             IF (iatom == jatom) fac = 0.5_dp
     357      8688539 :             IF (disp_a%type == dftd3_pp .AND. dr > 0.001_dp) THEN
     358      8665866 :                IF (dispersion_env%nd3_exclude_pair > 0) THEN
     359              :                   CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
     360          320 :                                             exclude=exclude_pair)
     361          320 :                   IF (exclude_pair) CYCLE
     362              :                END IF
     363              :                ! C6 coefficient
     364              :                CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
     365      8665634 :                           cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, c6, dc6a, dc6b)
     366      8665634 :                c8 = 3.0d0*c6*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
     367      8665634 :                dc8a = 3.0d0*dc6a*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
     368      8665634 :                dc8b = 3.0d0*dc6b*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
     369      8665634 :                r6 = dr**6
     370      8665634 :                r8 = r6*dr*dr
     371      8665634 :                s8i = s8
     372      8665634 :                IF (domol) THEN
     373         3568 :                   IF (atom2mol(iatom) /= atom2mol(jatom)) THEN
     374         1452 :                      s8i = kgc8
     375              :                   END IF
     376              :                END IF
     377              :                ! damping
     378      8665634 :                IF (idmp == 1) THEN
     379              :                   ! zero
     380      3878299 :                   CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs6, alp6, rcut, fdab6, dfdab6)
     381      3878299 :                   CALL damping_d3(dr, dispersion_env%r0ab(za, zb), rs8, alp8, rcut, fdab8, dfdab8)
     382      3878299 :                   e6 = s6*fac*c6*fdab6/r6
     383      3878299 :                   e8 = s8i*fac*c8*fdab8/r8
     384      4787335 :                ELSE IF (idmp == 2) THEN
     385              :                   ! BJ
     386      4787335 :                   r0ab = SQRT(3.0d0*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb))
     387      4787335 :                   f0ab = a1*r0ab + a2
     388      4787335 :                   fdab6 = 1.0_dp/(r6 + f0ab**6)
     389      4787335 :                   fdab8 = 1.0_dp/(r8 + f0ab**8)
     390      4787335 :                   e6 = s6*fac*c6*fdab6
     391      4787335 :                   e8 = s8i*fac*c8*fdab8
     392              :                ELSE
     393            0 :                   CPABORT("Unknown DFT-D3 damping function:")
     394              :                END IF
     395      8665634 :                IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
     396           69 :                   srbe = ssrb*(REAL((za*zb), KIND=dp))**t1srb*EXP(-gsrb*dr*dispersion_env%r0ab(za, zb)**t2srb)
     397           69 :                   esrb = esrb + srbe
     398           69 :                   evdw = evdw - srbe
     399              :                ELSE
     400              :                   srbe = 0.0_dp
     401              :                END IF
     402      8665634 :                evdw = evdw - e6 - e8
     403      8665634 :                e6tot = e6tot - e6
     404      8665634 :                e8tot = e8tot - e8
     405      8665634 :                IF (atenergy) THEN
     406      2765504 :                   atener(iatom) = atener(iatom) - 0.5_dp*(e6 + e8 + srbe)
     407      2765504 :                   atener(jatom) = atener(jatom) - 0.5_dp*(e6 + e8 + srbe)
     408              :                END IF
     409      8665634 :                IF (atex) THEN
     410          850 :                   atevdw(iatom) = atevdw(iatom) - 0.5_dp*(e6 + e8 + srbe)
     411          850 :                   atevdw(jatom) = atevdw(jatom) - 0.5_dp*(e6 + e8 + srbe)
     412              :                END IF
     413      8665634 :                IF (calculate_forces) THEN
     414              :                   ! damping
     415      3159756 :                   IF (idmp == 1) THEN
     416              :                      ! zero
     417      2002625 :                      de6 = -s6*c6/r6*(dfdab6 - 6._dp*fdab6/dr)
     418      2002625 :                      de8 = -s8i*c8/r8*(dfdab8 - 8._dp*fdab8/dr)
     419      1157131 :                   ELSE IF (idmp == 2) THEN
     420              :                      ! BJ
     421      1157131 :                      de6 = s6*c6*fdab6*fdab6*6.0_dp*dr**5
     422      1157131 :                      de8 = s8i*c8*fdab8*fdab8*8.0_dp*dr**7
     423              :                   ELSE
     424            0 :                      CPABORT("Unknown DFT-D3 damping function:")
     425              :                   END IF
     426     12639024 :                   fdij(:) = (de6 + de8)*rij(:)/dr*fac
     427      3159756 :                   IF (dispersion_env%srb .AND. dr .LT. 30.0d0) THEN
     428           92 :                      fdij(:) = fdij(:) + srbe*gsrb*dispersion_env%r0ab(za, zb)**t2srb*rij(:)/dr
     429              :                   END IF
     430      3159756 :                   atom_a = atom_of_kind(iatom)
     431      3159756 :                   atom_b = atom_of_kind(jatom)
     432     12639024 :                   force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
     433     12639024 :                   force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
     434      3159756 :                   IF (use_virial) THEN
     435      1999276 :                      CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
     436              :                   END IF
     437              :                   ! forces from the r-dependence of the coordination numbers
     438      3159756 :                   IF (idmp == 1) THEN
     439              :                      ! zero
     440      2002625 :                      dc6a = -s6*fac*dc6a*fdab6/r6
     441      2002625 :                      dc6b = -s6*fac*dc6b*fdab6/r6
     442      2002625 :                      dc8a = -s8i*fac*dc8a*fdab8/r8
     443      2002625 :                      dc8b = -s8i*fac*dc8b*fdab8/r8
     444      1157131 :                   ELSE IF (idmp == 2) THEN
     445              :                      ! BJ
     446      1157131 :                      dc6a = -s6*fac*dc6a*fdab6
     447      1157131 :                      dc6b = -s6*fac*dc6b*fdab6
     448      1157131 :                      dc8a = -s8i*fac*dc8a*fdab8
     449      1157131 :                      dc8b = -s8i*fac*dc8b*fdab8
     450              :                   ELSE
     451            0 :                      CPABORT("Unknown DFT-D3 damping function:")
     452              :                   END IF
     453     40796377 :                   DO i = 1, dcnum(iatom)%neighbors
     454     37636621 :                      katom = dcnum(iatom)%nlist(i)
     455     37636621 :                      kkind = kind_of(katom)
     456    150546484 :                      rik = dcnum(iatom)%rik(:, i)
     457    150546484 :                      drk = SQRT(SUM(rik(:)**2))
     458    150546484 :                      fdik(:) = (dc6a + dc8a)*dcnum(iatom)%dvals(i)*rik(:)/drk
     459     37636621 :                      atom_c = atom_of_kind(katom)
     460    150546484 :                      force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
     461    150546484 :                      force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
     462     40796377 :                      IF (use_virial) THEN
     463     26619815 :                         CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     464              :                      END IF
     465              :                   END DO
     466     40794825 :                   DO i = 1, dcnum(jatom)%neighbors
     467     37635069 :                      katom = dcnum(jatom)%nlist(i)
     468     37635069 :                      kkind = kind_of(katom)
     469    150540276 :                      rik = dcnum(jatom)%rik(:, i)
     470    150540276 :                      drk = SQRT(SUM(rik(:)**2))
     471    150540276 :                      fdik(:) = (dc6b + dc8b)*dcnum(jatom)%dvals(i)*rik(:)/drk
     472     37635069 :                      atom_c = atom_of_kind(katom)
     473    150540276 :                      force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
     474    150540276 :                      force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdik(:)
     475     40794825 :                      IF (use_virial) THEN
     476     26617067 :                         CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     477              :                      END IF
     478              :                   END DO
     479              :                END IF
     480      8665634 :                IF (dispersion_env%doabc) THEN
     481        16052 :                   CALL get_iterator_info(nl_iterator, cell=cell_b)
     482        16052 :                   hashb = cellhash(cell_b, ncell)
     483        24294 :                   is000 = (ALL(cell_b == 0))
     484       256832 :                   rb0(:) = MATMUL(cell%hmat, cell_b)
     485        16052 :                   ra(:) = pbc(particle_set(iatom)%r(:), cell)
     486        80260 :                   rb(:) = pbc(particle_set(jatom)%r(:), cell) + rb0
     487        64208 :                   r2ab = SUM((ra - rb)**2)
     488       103248 :                   DO icx = -ncell(1), ncell(1)
     489       626724 :                      DO icy = -ncell(2), ncell(2)
     490      4092908 :                         DO icz = -ncell(3), ncell(3)
     491      3482236 :                            cell_c(1) = icx
     492      3482236 :                            cell_c(2) = icy
     493      3482236 :                            cell_c(3) = icz
     494      3482236 :                            hashc = cellhash(cell_c, ncell)
     495      4108960 :                            IF (is000 .AND. (ALL(cell_c == 0))) THEN
     496              :                               ! CASE 1: all atoms in (000), use only ordered triples
     497          874 :                               kstart = MAX(jatom + 1, iatom + 1)
     498          874 :                               fac0 = 1.0_dp
     499      3481362 :                            ELSE IF (is000) THEN
     500              :                               ! CASE 2: AB in (000), C in other cell
     501              :                               !         This case covers also all instances with BC in same
     502              :                               !         cell not (000)
     503              :                               kstart = 1
     504              :                               fac0 = 1.0_dp
     505              :                            ELSE
     506              :                               ! These are case 2 again, cycle
     507      3446242 :                               IF (hashc == hashb) CYCLE
     508      4042074 :                               IF (ALL(cell_c == 0)) CYCLE
     509              :                               ! CASE 3: A in (000) and B and C in different cells
     510              :                               kstart = 1
     511              :                               fac0 = 1.0_dp/3.0_dp
     512              :                            END IF
     513     55230080 :                            rc0 = MATMUL(cell%hmat, cell_c)
     514     14809501 :                            DO katom = kstart, natom
     515     10834145 :                               kkind = kind_of(katom)
     516     10834145 :                               IF (exclude(kkind) .OR. .NOT. dodisp(kkind)) CYCLE
     517     43249972 :                               rc(:) = rcpbc(:, katom) + rc0(:)
     518     43249972 :                               r2bc = SUM((rb - rc)**2)
     519     10812493 :                               IF (r2bc >= rc2) CYCLE
     520      9114324 :                               r2ca = SUM((rc - ra)**2)
     521      2278581 :                               IF (r2ca >= rc2) CYCLE
     522      1168819 :                               r2abc = r2ab*r2bc*r2ca
     523      1168819 :                               IF (r2abc <= 0.001_dp) CYCLE
     524      1168819 :                               IF (dispersion_env%nd3_exclude_pair > 0) THEN
     525              :                                  CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
     526         5066 :                                                            kkind, exclude_pair)
     527         5066 :                                  IF (exclude_pair) CYCLE
     528              :                               END IF
     529              :                               ! this is an empirical scaling
     530      4617895 :                               IF (r2abc <= 0.01*rc2*rc2*rc2) THEN
     531        47775 :                                  kkind = kind_of(katom)
     532        47775 :                                  atom_c = atom_of_kind(katom)
     533        47775 :                                  zc = atomnumber(kkind)
     534              :                                  ! avoid double counting!
     535        47775 :                                  fac = 1._dp
     536        47775 :                                  IF (iatom == jatom .OR. iatom == katom .OR. jatom == katom) fac = 0.5_dp
     537        47775 :                                  IF (iatom == jatom .AND. iatom == katom) fac = 1._dp/3._dp
     538        47775 :                                  fac = fac*fac0
     539        47775 :                                  nabc = nabc + 1._dp
     540        47775 :                                  IF (dispersion_env%c9cnst) THEN
     541              :                                     CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
     542        32520 :                                                cnumfix(iatom), cnumfix(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
     543              :                                     CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
     544        32520 :                                                cnumfix(jatom), cnumfix(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
     545              :                                     CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
     546        32520 :                                                cnumfix(katom), cnumfix(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
     547              :                                  ELSE
     548              :                                     CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
     549        15255 :                                                cnumbers(iatom), cnumbers(jatom), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
     550              :                                     CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
     551        15255 :                                                cnumbers(jatom), cnumbers(katom), dispersion_env%k3, cc6bc, dcc6bcb, dcc6bcc)
     552              :                                     CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
     553        15255 :                                                cnumbers(katom), cnumbers(iatom), dispersion_env%k3, cc6ca, dcc6cac, dcc6caa)
     554              :                                  END IF
     555        47775 :                                  c9 = -SQRT(cc6ab*cc6bc*cc6ca)
     556        47775 :                                  rabc = r2abc**(1._dp/6._dp)
     557              :                                  r0 = (dispersion_env%r0ab(za, zb)*dispersion_env%r0ab(zb, zc)* &
     558        47775 :                                        dispersion_env%r0ab(zc, za))**(1._dp/3._dp)
     559              :                                  ! bug fixed 3.10.2017
     560              :                                  ! correct value from alp6=14 to 16 as used in original paper
     561        47775 :                                  CALL damping_d3(rabc, r0, 4._dp/3._dp, 16.0_dp, rcut, fdabc, dfdabc)
     562        47775 :                                  s1 = r2ab + r2bc - r2ca
     563        47775 :                                  s2 = r2ab + r2ca - r2bc
     564        47775 :                                  s3 = r2ca + r2bc - r2ab
     565        47775 :                                  ang = 0.375_dp*s1*s2*s3/r2abc + 1.0_dp
     566              : 
     567        47775 :                                  e9 = s9*fac*fdabc*c9*ang/r2abc**1.50d0
     568        47775 :                                  evdw = evdw - e9
     569        47775 :                                  e9tot = e9tot - e9
     570        47775 :                                  IF (atenergy) THEN
     571        20568 :                                     atener(iatom) = atener(iatom) - e9/3._dp
     572        20568 :                                     atener(jatom) = atener(jatom) - e9/3._dp
     573        20568 :                                     atener(katom) = atener(katom) - e9/3._dp
     574              :                                  END IF
     575        47775 :                                  IF (atex) THEN
     576        10284 :                                     atevdw(iatom) = atevdw(iatom) - e9/3._dp
     577        10284 :                                     atevdw(jatom) = atevdw(jatom) - e9/3._dp
     578        10284 :                                     atevdw(katom) = atevdw(katom) - e9/3._dp
     579              :                                  END IF
     580              : 
     581        47775 :                                  IF (calculate_forces) THEN
     582       230390 :                                     rab = rb - ra; rbc = rc - rb; rca = ra - rc
     583        23039 :                                     de91 = s9*fac*dfdabc*c9*ang/r2abc**1.50d0
     584        23039 :                                     de91 = -de91/3._dp*rabc + 3._dp*e9
     585        23039 :                                     dea = s9*fac*fdabc*c9/r2abc**2.50d0*0.75_dp
     586        92156 :                                     fdij(:) = de91*rab(:)/r2ab
     587        92156 :                                     fdij(:) = fdij(:) + dea*s1*s2*s3*rab(:)/r2ab
     588        92156 :                                     fdij(:) = fdij(:) - dea*(s2*s3 + s1*s3 - s1*s2)*rab(:)
     589        92156 :                                     force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
     590        92156 :                                     force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
     591        23039 :                                     IF (use_virial) THEN
     592        15262 :                                        CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rab)
     593              :                                     END IF
     594        92156 :                                     fdij(:) = de91*rbc(:)/r2bc
     595        92156 :                                     fdij(:) = fdij(:) + dea*s1*s2*s3*rbc(:)/r2bc
     596        92156 :                                     fdij(:) = fdij(:) - dea*(s2*s3 - s1*s3 + s1*s2)*rbc(:)
     597        92156 :                                     force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdij(:)
     598        92156 :                                     force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) + fdij(:)
     599        23039 :                                     IF (use_virial) THEN
     600        15262 :                                        CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rbc)
     601              :                                     END IF
     602        92156 :                                     fdij(:) = de91*rca(:)/r2ca
     603        92156 :                                     fdij(:) = fdij(:) + dea*s1*s2*s3*rca(:)/r2ca
     604        92156 :                                     fdij(:) = fdij(:) - dea*(-s2*s3 + s1*s3 + s1*s2)*rca(:)
     605        92156 :                                     force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdij(:)
     606        92156 :                                     force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) + fdij(:)
     607        23039 :                                     IF (use_virial) THEN
     608        15262 :                                        CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rca)
     609              :                                     END IF
     610              : 
     611        23039 :                                     IF (.NOT. dispersion_env%c9cnst) THEN
     612              :                                        ! forces from the r-dependence of the coordination numbers
     613              :                                        ! atomic stress not implemented
     614              : 
     615        11087 :                                        de91 = 0.5_dp*e9/cc6ab
     616        11087 :                                        de921 = de91*dcc6aba
     617        11087 :                                        de922 = de91*dcc6abb
     618        38781 :                                        DO i = 1, dcnum(iatom)%neighbors
     619        27694 :                                           latom = dcnum(iatom)%nlist(i)
     620        27694 :                                           lkind = kind_of(latom)
     621        27694 :                                           rik(1) = dcnum(iatom)%rik(1, i)
     622        27694 :                                           rik(2) = dcnum(iatom)%rik(2, i)
     623        27694 :                                           rik(3) = dcnum(iatom)%rik(3, i)
     624        27694 :                                           drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
     625       110776 :                                           fdik(:) = -de921*dcnum(iatom)%dvals(i)*rik(:)/drk
     626        27694 :                                           atom_d = atom_of_kind(latom)
     627       110776 :                                           force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
     628       110776 :                                           force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
     629        38781 :                                           IF (use_virial) THEN
     630        27612 :                                              CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     631              :                                           END IF
     632              :                                        END DO
     633        38781 :                                        DO i = 1, dcnum(jatom)%neighbors
     634        27694 :                                           latom = dcnum(jatom)%nlist(i)
     635        27694 :                                           lkind = kind_of(latom)
     636        27694 :                                           rik(1) = dcnum(jatom)%rik(1, i)
     637        27694 :                                           rik(2) = dcnum(jatom)%rik(2, i)
     638        27694 :                                           rik(3) = dcnum(jatom)%rik(3, i)
     639        27694 :                                           drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
     640       110776 :                                           fdik(:) = -de922*dcnum(jatom)%dvals(i)*rik(:)/drk
     641        27694 :                                           atom_d = atom_of_kind(latom)
     642       110776 :                                           force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
     643       110776 :                                           force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
     644        38781 :                                           IF (use_virial) THEN
     645        27612 :                                              CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     646              :                                           END IF
     647              :                                        END DO
     648              : 
     649        11087 :                                        de91 = 0.5_dp*e9/cc6bc
     650        11087 :                                        de921 = de91*dcc6bcb
     651        11087 :                                        de922 = de91*dcc6bcc
     652        38781 :                                        DO i = 1, dcnum(jatom)%neighbors
     653        27694 :                                           latom = dcnum(jatom)%nlist(i)
     654        27694 :                                           lkind = kind_of(latom)
     655        27694 :                                           rik(1) = dcnum(jatom)%rik(1, i)
     656        27694 :                                           rik(2) = dcnum(jatom)%rik(2, i)
     657        27694 :                                           rik(3) = dcnum(jatom)%rik(3, i)
     658        27694 :                                           drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
     659       110776 :                                           fdik(:) = -de921*dcnum(jatom)%dvals(i)*rik(:)/drk
     660        27694 :                                           atom_d = atom_of_kind(latom)
     661       110776 :                                           force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) - fdik(:)
     662       110776 :                                           force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
     663        38781 :                                           IF (use_virial) THEN
     664        27612 :                                              CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     665              :                                           END IF
     666              :                                        END DO
     667        38781 :                                        DO i = 1, dcnum(katom)%neighbors
     668        27694 :                                           latom = dcnum(katom)%nlist(i)
     669        27694 :                                           lkind = kind_of(latom)
     670        27694 :                                           rik(1) = dcnum(katom)%rik(1, i)
     671        27694 :                                           rik(2) = dcnum(katom)%rik(2, i)
     672        27694 :                                           rik(3) = dcnum(katom)%rik(3, i)
     673        27694 :                                           drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
     674       110776 :                                           fdik(:) = -de922*dcnum(katom)%dvals(i)*rik(:)/drk
     675        27694 :                                           atom_d = atom_of_kind(latom)
     676       110776 :                                           force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
     677       110776 :                                           force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
     678        38781 :                                           IF (use_virial) THEN
     679        27612 :                                              CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     680              :                                           END IF
     681              :                                        END DO
     682              : 
     683        11087 :                                        de91 = 0.5_dp*e9/cc6ca
     684        11087 :                                        de921 = de91*dcc6cac
     685        11087 :                                        de922 = de91*dcc6caa
     686        38781 :                                        DO i = 1, dcnum(katom)%neighbors
     687        27694 :                                           latom = dcnum(katom)%nlist(i)
     688        27694 :                                           lkind = kind_of(latom)
     689        27694 :                                           rik(1) = dcnum(katom)%rik(1, i)
     690        27694 :                                           rik(2) = dcnum(katom)%rik(2, i)
     691        27694 :                                           rik(3) = dcnum(katom)%rik(3, i)
     692        27694 :                                           drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
     693       110776 :                                           fdik(:) = -de921*dcnum(katom)%dvals(i)*rik(:)/drk
     694        27694 :                                           atom_d = atom_of_kind(latom)
     695       110776 :                                           force(kkind)%dispersion(:, atom_c) = force(kkind)%dispersion(:, atom_c) - fdik(:)
     696       110776 :                                           force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
     697        38781 :                                           IF (use_virial) THEN
     698        27612 :                                              CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     699              :                                           END IF
     700              :                                        END DO
     701        38781 :                                        DO i = 1, dcnum(iatom)%neighbors
     702        27694 :                                           latom = dcnum(iatom)%nlist(i)
     703        27694 :                                           lkind = kind_of(latom)
     704        27694 :                                           rik(1) = dcnum(iatom)%rik(1, i)
     705        27694 :                                           rik(2) = dcnum(iatom)%rik(2, i)
     706        27694 :                                           rik(3) = dcnum(iatom)%rik(3, i)
     707        27694 :                                           drk = SQRT(rik(1)*rik(1) + rik(2)*rik(2) + rik(3)*rik(3))
     708       110776 :                                           fdik(:) = -de922*dcnum(iatom)%dvals(i)*rik(:)/drk
     709        27694 :                                           atom_d = atom_of_kind(latom)
     710       110776 :                                           force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdik(:)
     711       110776 :                                           force(lkind)%dispersion(:, atom_d) = force(lkind)%dispersion(:, atom_d) + fdik(:)
     712        38781 :                                           IF (use_virial) THEN
     713        27612 :                                              CALL virial_pair_force(pv_virial_thread, -1._dp, fdik, rik)
     714              :                                           END IF
     715              :                                        END DO
     716              :                                     END IF
     717              : 
     718              :                                  END IF
     719              : 
     720              :                               END IF
     721              :                            END DO
     722              :                         END DO
     723              :                      END DO
     724              :                   END DO
     725              : 
     726              :                END IF
     727              :             END IF
     728              :          END IF
     729              :       END DO
     730              : 
     731        61412 :       virial%pv_virial = virial%pv_virial + pv_virial_thread
     732              : 
     733         4724 :       CALL neighbor_list_iterator_release(nl_iterator)
     734              : 
     735         4724 :       elrc6 = 0._dp
     736         4724 :       elrc8 = 0._dp
     737         4724 :       elrc9 = 0._dp
     738              :       ! Long range correction (atomic contributions not implemented)
     739         4724 :       IF (dispersion_env%lrc) THEN
     740          144 :          ALLOCATE (cnkind(nkind))
     741          106 :          cnkind = 0._dp
     742              :          ! first use the default values
     743          106 :          DO ikind = 1, nkind
     744           58 :             CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
     745          106 :             cnkind(ikind) = dispersion_env%cn(za)
     746              :          END DO
     747              :          ! now check for changes from default
     748           48 :          IF (ASSOCIATED(dispersion_env%cnkind)) THEN
     749           12 :             DO i = 1, SIZE(dispersion_env%cnkind)
     750            6 :                ikind = dispersion_env%cnkind(i)%kind
     751           12 :                cnkind(ikind) = dispersion_env%cnkind(i)%cnum
     752              :             END DO
     753              :          END IF
     754          106 :          DO ikind = 1, nkind
     755           58 :             CALL get_atomic_kind(atomic_kind_set(ikind), natom=na, z=za)
     756           58 :             CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
     757           58 :             IF (.NOT. disp_a%defined .OR. ghost_a .OR. floating_a) CYCLE
     758          236 :             DO jkind = 1, nkind
     759           74 :                CALL get_atomic_kind(atomic_kind_set(jkind), natom=nb, z=zb)
     760           74 :                CALL get_qs_kind(qs_kind_set(jkind), dispersion=disp_b, ghost=ghost_b, floating=floating_b)
     761           74 :                IF (.NOT. disp_b%defined .OR. ghost_b .OR. floating_b) CYCLE
     762           72 :                IF (dispersion_env%nd3_exclude_pair > 0) THEN
     763              :                   CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, &
     764            8 :                                             exclude=exclude_pair)
     765            8 :                   IF (exclude_pair) CYCLE
     766              :                END IF
     767              :                CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
     768           66 :                           cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
     769           66 :                elrc6 = elrc6 - s6*twopi*REAL(na*nb, KIND=dp)*cc6ab/(3._dp*rcut**3*cell%deth)
     770           66 :                c8 = 3.0d0*cc6ab*dispersion_env%r2r4(za)*dispersion_env%r2r4(zb)
     771           66 :                elrc8 = elrc8 - s8*twopi*REAL(na*nb, KIND=dp)*c8/(5._dp*rcut**5*cell%deth)
     772          198 :                IF (dispersion_env%doabc) THEN
     773          160 :                   DO kkind = 1, nkind
     774           94 :                      CALL get_atomic_kind(atomic_kind_set(kkind), natom=nc, z=zc)
     775           94 :                      CALL get_qs_kind(qs_kind_set(kkind), dispersion=disp_c, ghost=ghost_c, floating=floating_c)
     776           94 :                      IF (.NOT. disp_c%defined .OR. ghost_c .OR. floating_c) CYCLE
     777           92 :                      IF (dispersion_env%nd3_exclude_pair > 0) THEN
     778              :                         CALL exclude_d3_kind_pair(dispersion_env%d3_exclude_pair, ikind, jkind, kkind, &
     779            4 :                                                   exclude_pair)
     780            4 :                         IF (exclude_pair) CYCLE
     781              :                      END IF
     782              :                      CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, za, zb, &
     783           90 :                                 cnkind(ikind), cnkind(jkind), dispersion_env%k3, cc6ab, dcc6aba, dcc6abb)
     784              :                      CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zc, za, &
     785           90 :                                 cnkind(kkind), cnkind(ikind), dispersion_env%k3, cc6ca, dcc6aba, dcc6abb)
     786              :                      CALL getc6(maxc, max_elem, dispersion_env%c6ab, dispersion_env%maxci, zb, zc, &
     787           90 :                                 cnkind(jkind), cnkind(kkind), dispersion_env%k3, cc6bc, dcc6aba, dcc6abb)
     788           90 :                      c9 = -SQRT(cc6ab*cc6bc*cc6ca)
     789          250 :                      elrc9 = elrc9 - s9*64._dp*twopi*REAL(na*nb*nc, KIND=dp)*c9/(6._dp*rcut**3*cell%deth**2)
     790              :                   END DO
     791              :                END IF
     792              :             END DO
     793              :          END DO
     794           48 :          IF (use_virial) THEN
     795           34 :             IF (para_env%is_source()) THEN
     796           68 :                DO i = 1, 3
     797           68 :                   virial%pv_virial(i, i) = virial%pv_virial(i, i) + (elrc6 + elrc8 + 2._dp*elrc9)
     798              :                END DO
     799              :             END IF
     800              :          END IF
     801           48 :          DEALLOCATE (cnkind)
     802              :       END IF
     803              : 
     804         4724 :       DEALLOCATE (cnumbers)
     805         4724 :       IF (dispersion_env%doabc .AND. dispersion_env%c9cnst) THEN
     806           62 :          DEALLOCATE (cnumfix)
     807              :       END IF
     808         4724 :       IF (calculate_forces .OR. debugall) THEN
     809        10394 :          DO iatom = 1, natom
     810        10394 :             DEALLOCATE (dcnum(iatom)%nlist, dcnum(iatom)%dvals, dcnum(iatom)%rik)
     811              :          END DO
     812          762 :          DEALLOCATE (dcnum)
     813              :       ELSE
     814         3962 :          DEALLOCATE (dcnum)
     815              :       END IF
     816              : 
     817              :       ! set dispersion energy
     818              : 
     819         4724 :       CALL para_env%sum(e6tot)
     820         4724 :       CALL para_env%sum(e8tot)
     821         4724 :       CALL para_env%sum(e9tot)
     822         4724 :       CALL para_env%sum(evdw)
     823         4724 :       CALL para_env%sum(nab)
     824         4724 :       CALL para_env%sum(nabc)
     825         4724 :       e6tot = e6tot + elrc6
     826         4724 :       e8tot = e8tot + elrc8
     827         4724 :       e9tot = e9tot + elrc9
     828              :       ! For printing, we need all contributions
     829         4724 :       evdw = evdw + (elrc6 + elrc8 + elrc9)
     830         4724 :       IF (unit_nr > 0) THEN
     831            8 :          WRITE (unit_nr, "(A,F20.0)") "  E6 vdW terms              :", nab
     832            8 :          WRITE (unit_nr, *) " E6 vdW energy [au/kcal]   :", e6tot, e6tot*kcalmol
     833            8 :          WRITE (unit_nr, *) " E8 vdW energy [au/kcal]   :", e8tot, e8tot*kcalmol
     834            8 :          WRITE (unit_nr, *) " %E8 on total vdW energy   :", e8tot/evdw*100._dp
     835            8 :          WRITE (unit_nr, "(A,F20.0)") "  E9 vdW terms              :", nabc
     836            8 :          WRITE (unit_nr, *) " E9 vdW energy [au/kcal]   :", e9tot, e9tot*kcalmol
     837            8 :          WRITE (unit_nr, *) " %E9 on total vdW energy   :", e9tot/evdw*100._dp
     838            8 :          IF (dispersion_env%lrc) THEN
     839            1 :             WRITE (unit_nr, *) " E LRC C6 [au/kcal]        :", elrc6, elrc6*kcalmol
     840            1 :             WRITE (unit_nr, *) " E LRC C8 [au/kcal]        :", elrc8, elrc8*kcalmol
     841            1 :             WRITE (unit_nr, *) " E LRC C9 [au/kcal]        :", elrc9, elrc9*kcalmol
     842              :          END IF
     843              :       END IF
     844              : 
     845         4724 :       evdw = evdw/para_env%num_pe
     846              : 
     847         4724 :       DEALLOCATE (dodisp, exclude, atomnumber, rcpbc, radd2, c6d2)
     848              : 
     849         4724 :       IF (domol) THEN
     850            2 :          DEALLOCATE (atom2mol)
     851              :       END IF
     852              : 
     853         4724 :       CALL timestop(handle)
     854              : 
     855         9448 :    END SUBROUTINE calculate_dispersion_d3_pairpot
     856              : 
     857              : ! **************************************************************************************************
     858              : !> \brief ...
     859              : !> \param c6ab ...
     860              : !> \param maxci ...
     861              : !> \param filename ...
     862              : !> \param para_env ...
     863              : ! **************************************************************************************************
     864          392 :    SUBROUTINE dftd3_c6_param(c6ab, maxci, filename, para_env)
     865              : 
     866              :       REAL(KIND=dp), DIMENSION(:, :, :, :, :)            :: c6ab
     867              :       INTEGER, DIMENSION(:)                              :: maxci
     868              :       CHARACTER(LEN=*)                                   :: filename
     869              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     870              : 
     871              :       INTEGER                                            :: funit, iadr, iat, jadr, jat, kk, nl, &
     872              :                                                             nlines, nn
     873          392 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: pars
     874              : 
     875          392 :       IF (para_env%is_source()) THEN
     876              :          ! Read the DFT-D3 C6AB parameters from file "filename"
     877          201 :          CALL open_file(file_name=filename, unit_number=funit, file_form="FORMATTED")
     878          201 :          READ (funit, *) nl, nlines
     879              :       END IF
     880          392 :       CALL para_env%bcast(nl)
     881          392 :       CALL para_env%bcast(nlines)
     882         1176 :       ALLOCATE (pars(nl))
     883          392 :       IF (para_env%is_source()) THEN
     884          201 :          READ (funit, *) pars(1:nl)
     885          201 :          CALL close_file(unit_number=funit)
     886              :       END IF
     887          392 :       CALL para_env%bcast(pars)
     888              : 
     889              :       ! Store C6AB coefficients in an array
     890    262578848 :       c6ab = -1._dp
     891        37240 :       maxci = 0
     892          392 :       kk = 1
     893     12695312 :       DO nn = 1, nlines
     894     12694920 :          iat = NINT(pars(kk + 1))
     895     12694920 :          jat = NINT(pars(kk + 2))
     896     12694920 :          CALL limit(iat, jat, iadr, jadr)
     897     12694920 :          maxci(iat) = MAX(maxci(iat), iadr)
     898     12694920 :          maxci(jat) = MAX(maxci(jat), jadr)
     899     12694920 :          c6ab(iat, jat, iadr, jadr, 1) = pars(kk)
     900     12694920 :          c6ab(iat, jat, iadr, jadr, 2) = pars(kk + 3)
     901     12694920 :          c6ab(iat, jat, iadr, jadr, 3) = pars(kk + 4)
     902              : 
     903     12694920 :          c6ab(jat, iat, jadr, iadr, 1) = pars(kk)
     904     12694920 :          c6ab(jat, iat, jadr, iadr, 2) = pars(kk + 4)
     905     12694920 :          c6ab(jat, iat, jadr, iadr, 3) = pars(kk + 3)
     906     12695312 :          kk = (nn*5) + 1
     907              :       END DO
     908              : 
     909          392 :       DEALLOCATE (pars)
     910              : 
     911          392 :    END SUBROUTINE dftd3_c6_param
     912              : 
     913              : ! **************************************************************************************************
     914              : !> \brief ...
     915              : !> \param iat ...
     916              : !> \param jat ...
     917              : !> \param iadr ...
     918              : !> \param jadr ...
     919              : ! **************************************************************************************************
     920     12694920 :    SUBROUTINE limit(iat, jat, iadr, jadr)
     921              :       INTEGER                                            :: iat, jat, iadr, jadr
     922              : 
     923              :       INTEGER                                            :: i
     924              : 
     925     12694920 :       iadr = 1
     926     12694920 :       jadr = 1
     927     12694920 :       i = 100
     928     32660264 :       DO WHILE (iat .GT. 100)
     929     19965344 :          iat = iat - 100
     930     32660264 :          iadr = iadr + 1
     931              :       END DO
     932              : 
     933     18919096 :       i = 100
     934     18919096 :       DO WHILE (jat .GT. 100)
     935      6224176 :          jat = jat - 100
     936      6224176 :          jadr = jadr + 1
     937              :       END DO
     938     12694920 :    END SUBROUTINE limit
     939              : 
     940              : ! **************************************************************************************************
     941              : !> \brief ...
     942              : !> \param rab ...
     943              : !> \param rcutab ...
     944              : !> \param srn ...
     945              : !> \param alpn ...
     946              : !> \param rcut ...
     947              : !> \param fdab ...
     948              : !> \param dfdab ...
     949              : ! **************************************************************************************************
     950      7804373 :    SUBROUTINE damping_d3(rab, rcutab, srn, alpn, rcut, fdab, dfdab)
     951              : 
     952              :       REAL(KIND=dp), INTENT(IN)                          :: rab, rcutab, srn, alpn, rcut
     953              :       REAL(KIND=dp), INTENT(OUT)                         :: fdab, dfdab
     954              : 
     955              :       REAL(KIND=dp)                                      :: a, b, c, d, dd, dfab, dfcc, dz, fab, &
     956              :                                                             fcc, rl, rr, ru, z, zz
     957              : 
     958      7804373 :       rl = rcut - 1._dp
     959      7804373 :       ru = rcut
     960      7804373 :       IF (rab >= ru) THEN
     961              :          fcc = 0._dp
     962              :          dfcc = 0._dp
     963      7804373 :       ELSEIF (rab <= rl) THEN
     964              :          fcc = 1._dp
     965              :          dfcc = 0._dp
     966              :       ELSE
     967       712324 :          z = rab*rab - rl*rl
     968       712324 :          dz = 2._dp*rab
     969       712324 :          zz = z*z*z
     970       712324 :          d = ru*ru - rl*rl
     971       712324 :          dd = d*d*d
     972       712324 :          a = -10._dp/dd
     973       712324 :          b = 15._dp/(dd*d)
     974       712324 :          c = -6._dp/(dd*d*d)
     975       712324 :          fcc = 1._dp + zz*(a + b*z + c*z*z)
     976       712324 :          dfcc = zz*dz/z*(3._dp*a + 4._dp*b*z + 5._dp*c*z*z)
     977              :       END IF
     978              : 
     979      7804373 :       rr = 6._dp*(rab/(srn*rcutab))**(-alpn)
     980      7804373 :       fab = 1._dp/(1._dp + rr)
     981      7804373 :       dfab = fab*fab*rr*alpn/rab
     982      7804373 :       fdab = fab*fcc
     983      7804373 :       dfdab = dfab*fcc + fab*dfcc
     984              : 
     985      7804373 :    END SUBROUTINE damping_d3
     986              : 
     987              : ! **************************************************************************************************
     988              : !> \brief ...
     989              : !> \param maxc ...
     990              : !> \param max_elem ...
     991              : !> \param c6ab ...
     992              : !> \param mxc ...
     993              : !> \param iat ...
     994              : !> \param jat ...
     995              : !> \param nci ...
     996              : !> \param ncj ...
     997              : !> \param k3 ...
     998              : !> \param c6 ...
     999              : !> \param dc6a ...
    1000              : !> \param dc6b ...
    1001              : ! **************************************************************************************************
    1002      8809295 :    SUBROUTINE getc6(maxc, max_elem, c6ab, mxc, iat, jat, nci, ncj, k3, c6, dc6a, dc6b)
    1003              : 
    1004              :       INTEGER, INTENT(IN)                                :: maxc, max_elem
    1005              :       REAL(KIND=dp), INTENT(IN) :: c6ab(max_elem, max_elem, maxc, maxc, 3)
    1006              :       INTEGER, INTENT(IN)                                :: mxc(max_elem), iat, jat
    1007              :       REAL(KIND=dp), INTENT(IN)                          :: nci, ncj, k3
    1008              :       REAL(KIND=dp), INTENT(OUT)                         :: c6, dc6a, dc6b
    1009              : 
    1010              :       INTEGER                                            :: i, j
    1011              :       REAL(KIND=dp)                                      :: c6mem, cn1, cn2, csum, dtmpa, dtmpb, &
    1012              :                                                             dwa, dwb, dza, dzb, r, rsave, rsum, &
    1013              :                                                             tmp1
    1014              : 
    1015              : ! the exponential is sensitive to numerics
    1016              : ! when nci or ncj is much larger than cn1/cn2
    1017              : 
    1018      8809295 :       c6mem = -1.0e+99_dp
    1019      8809295 :       rsave = 1.0e+99_dp
    1020      8809295 :       rsum = 0.0_dp
    1021      8809295 :       csum = 0.0_dp
    1022      8809295 :       dza = 0.0_dp
    1023      8809295 :       dzb = 0.0_dp
    1024      8809295 :       dwa = 0.0_dp
    1025      8809295 :       dwb = 0.0_dp
    1026      8809295 :       c6 = 0.0_dp
    1027     32102402 :       DO i = 1, mxc(iat)
    1028     98870466 :          DO j = 1, mxc(jat)
    1029     66768064 :             c6 = c6ab(iat, jat, i, j, 1)
    1030     90061171 :             IF (c6 > 0.0_dp) THEN
    1031     66768064 :                cn1 = c6ab(iat, jat, i, j, 2)
    1032     66768064 :                cn2 = c6ab(iat, jat, i, j, 3)
    1033              :                ! distance
    1034     66768064 :                r = (cn1 - nci)**2 + (cn2 - ncj)**2
    1035     66768064 :                IF (r < rsave) THEN
    1036     33154538 :                   rsave = r
    1037     33154538 :                   c6mem = c6
    1038              :                END IF
    1039     66768064 :                tmp1 = EXP(k3*r)
    1040     66768064 :                dtmpa = -2.0_dp*k3*(cn1 - nci)*tmp1
    1041     66768064 :                dtmpb = -2.0_dp*k3*(cn2 - ncj)*tmp1
    1042     66768064 :                rsum = rsum + tmp1
    1043     66768064 :                csum = csum + tmp1*c6
    1044     66768064 :                dza = dza + dtmpa*c6
    1045     66768064 :                dwa = dwa + dtmpa
    1046     66768064 :                dzb = dzb + dtmpb*c6
    1047     66768064 :                dwb = dwb + dtmpb
    1048              :             END IF
    1049              :          END DO
    1050              :       END DO
    1051              : 
    1052      8809295 :       IF (c6 == 0.0_dp) c6mem = 0.0_dp
    1053              : 
    1054      8809295 :       IF (rsum > 1.0e-66_dp) THEN
    1055      8794927 :          c6 = csum/rsum
    1056      8794927 :          dc6a = (dza - c6*dwa)/rsum
    1057      8794927 :          dc6b = (dzb - c6*dwb)/rsum
    1058              :       ELSE
    1059        14368 :          c6 = c6mem
    1060        14368 :          dc6a = 0._dp
    1061        14368 :          dc6b = 0._dp
    1062              :       END IF
    1063              : 
    1064      8809295 :    END SUBROUTINE getc6
    1065              : 
    1066              : END MODULE qs_dispersion_d3
        

Generated by: LCOV version 2.0-1