Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculation of D2 dispersion
      10             : !> \author JGH
      11             : ! **************************************************************************************************
      12             : MODULE qs_dispersion_d2
      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             :    USE kinds,                           ONLY: dp
      21             :    USE physcon,                         ONLY: bohr,&
      22             :                                               kjmol
      23             :    USE qs_dispersion_types,             ONLY: qs_atom_dispersion_type,&
      24             :                                               qs_dispersion_type
      25             :    USE qs_environment_types,            ONLY: get_qs_env,&
      26             :                                               qs_environment_type
      27             :    USE qs_force_types,                  ONLY: qs_force_type
      28             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      29             :                                               qs_kind_type
      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             : #include "./base/base_uses.f90"
      39             : 
      40             :    IMPLICIT NONE
      41             : 
      42             :    PRIVATE
      43             : 
      44             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d2'
      45             : 
      46             :    PUBLIC :: calculate_dispersion_d2_pairpot, dftd2_param
      47             : 
      48             : ! **************************************************************************************************
      49             : 
      50             : CONTAINS
      51             : 
      52             : ! **************************************************************************************************
      53             : !> \brief ...
      54             : !> \param qs_env ...
      55             : !> \param dispersion_env ...
      56             : !> \param evdw ...
      57             : !> \param calculate_forces ...
      58             : !> \param atevdw ...
      59             : ! **************************************************************************************************
      60          58 :    SUBROUTINE calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
      61             : 
      62             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      63             :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
      64             :       REAL(KIND=dp), INTENT(OUT)                         :: evdw
      65             :       LOGICAL, INTENT(IN)                                :: calculate_forces
      66             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atevdw
      67             : 
      68             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d2_pairpot'
      69             : 
      70             :       INTEGER                                            :: atom_a, atom_b, handle, iatom, ikind, &
      71             :                                                             jatom, jkind, mepos, natom, nkind, &
      72             :                                                             num_pe, za, zb
      73          58 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, atomnumber, kind_of
      74             :       LOGICAL                                            :: atenergy, atex, floating_a, ghost_a, &
      75             :                                                             use_virial
      76          58 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: dodisp, floating, ghost
      77             :       REAL(KIND=dp)                                      :: c6, dd, devdw, dfdmp, dr, er, fac, fdmp, &
      78             :                                                             rcc, rcut, s6, xp
      79          58 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: c6d2, radd2
      80             :       REAL(KIND=dp), DIMENSION(3)                        :: fdij, rij
      81             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pv_virial_thread
      82          58 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: atener
      83          58 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      84             :       TYPE(atprop_type), POINTER                         :: atprop
      85             :       TYPE(cell_type), POINTER                           :: cell
      86             :       TYPE(neighbor_list_iterator_p_type), &
      87          58 :          DIMENSION(:), POINTER                           :: nl_iterator
      88             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
      89          58 :          POINTER                                         :: sab_vdw
      90             :       TYPE(qs_atom_dispersion_type), POINTER             :: disp_a
      91          58 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      92          58 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      93             :       TYPE(virial_type), POINTER                         :: virial
      94             : 
      95          58 :       CALL timeset(routineN, handle)
      96             : 
      97          58 :       evdw = 0._dp
      98             : 
      99          58 :       NULLIFY (atomic_kind_set, qs_kind_set, sab_vdw)
     100             : 
     101             :       CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
     102          58 :                       qs_kind_set=qs_kind_set, cell=cell, virial=virial, atprop=atprop)
     103             : 
     104             :       ! atomic energy and stress arrays
     105          58 :       atenergy = atprop%energy
     106          58 :       IF (atenergy) THEN
     107           0 :          CALL atprop_array_init(atprop%atevdw, natom)
     108           0 :          atener => atprop%atevdw
     109             :       END IF
     110             :       ! external atomic energy
     111          58 :       atex = .FALSE.
     112          58 :       IF (PRESENT(atevdw)) THEN
     113           0 :          atex = .TRUE.
     114             :       END IF
     115             : 
     116          58 :       NULLIFY (force)
     117          58 :       CALL get_qs_env(qs_env=qs_env, force=force)
     118          58 :       CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     119          58 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     120          58 :       pv_virial_thread(:, :) = 0._dp
     121             : 
     122         522 :       ALLOCATE (dodisp(nkind), ghost(nkind), floating(nkind), atomnumber(nkind), c6d2(nkind), radd2(nkind))
     123         188 :       DO ikind = 1, nkind
     124         130 :          CALL get_atomic_kind(atomic_kind_set(ikind), z=za)
     125         130 :          CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp_a, ghost=ghost_a, floating=floating_a)
     126         130 :          dodisp(ikind) = disp_a%defined
     127         130 :          ghost(ikind) = ghost_a
     128         130 :          floating(ikind) = floating_a
     129         130 :          atomnumber(ikind) = za
     130         130 :          c6d2(ikind) = disp_a%c6
     131         318 :          radd2(ikind) = disp_a%vdw_radii
     132             :       END DO
     133             : 
     134          58 :       rcut = 2._dp*dispersion_env%rc_disp
     135          58 :       s6 = dispersion_env%scaling
     136          58 :       dd = dispersion_env%exp_pre
     137             : 
     138          58 :       sab_vdw => dispersion_env%sab_vdw
     139          58 :       num_pe = 1
     140          58 :       CALL neighbor_list_iterator_create(nl_iterator, sab_vdw, nthread=num_pe)
     141             : 
     142          58 :       mepos = 0
     143       32703 :       DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
     144       32645 :          CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
     145             : 
     146       32645 :          IF (ghost(ikind) .OR. ghost(jkind) .OR. floating(ikind) .OR. floating(jkind)) CYCLE
     147             : 
     148       32147 :          IF (.NOT. (dodisp(ikind) .AND. dodisp(jkind))) CYCLE
     149             : 
     150       32147 :          za = atomnumber(ikind)
     151       32147 :          zb = atomnumber(jkind)
     152             :          ! vdW potential
     153      128588 :          dr = SQRT(SUM(rij(:)**2))
     154       32205 :          IF (dr <= rcut) THEN
     155       32147 :             fac = 1._dp
     156       32147 :             IF (iatom == jatom) fac = 0.5_dp
     157       32147 :             IF (dr > 0.001_dp) THEN
     158       31997 :                c6 = SQRT(c6d2(ikind)*c6d2(jkind))
     159       31997 :                rcc = radd2(ikind) + radd2(jkind)
     160       31997 :                er = EXP(-dd*(dr/rcc - 1._dp))
     161       31997 :                fdmp = 1._dp/(1._dp + er)
     162       31997 :                xp = s6*c6/dr**6
     163       31997 :                evdw = evdw - xp*fdmp*fac
     164       31997 :                IF (calculate_forces) THEN
     165       31432 :                   dfdmp = dd/rcc*er*fdmp*fdmp
     166       31432 :                   devdw = -xp*(-6._dp*fdmp/dr + dfdmp)
     167      125728 :                   fdij(:) = devdw*rij(:)/dr*fac
     168       31432 :                   atom_a = atom_of_kind(iatom)
     169       31432 :                   atom_b = atom_of_kind(jatom)
     170      125728 :                   force(ikind)%dispersion(:, atom_a) = force(ikind)%dispersion(:, atom_a) - fdij(:)
     171      125728 :                   force(jkind)%dispersion(:, atom_b) = force(jkind)%dispersion(:, atom_b) + fdij(:)
     172       31432 :                   IF (use_virial) THEN
     173        9842 :                      CALL virial_pair_force(pv_virial_thread, -1._dp, fdij, rij)
     174             :                   END IF
     175             :                END IF
     176       31997 :                IF (atenergy) THEN
     177           0 :                   atener(iatom) = atener(iatom) - 0.5_dp*xp*fdmp*fac
     178           0 :                   atener(jatom) = atener(jatom) - 0.5_dp*xp*fdmp*fac
     179             :                END IF
     180       31997 :                IF (atex) THEN
     181           0 :                   atevdw(iatom) = atevdw(iatom) - 0.5_dp*xp*fdmp*fac
     182           0 :                   atevdw(jatom) = atevdw(jatom) - 0.5_dp*xp*fdmp*fac
     183             :                END IF
     184             :             END IF
     185             :          END IF
     186             : 
     187             :       END DO
     188             : 
     189         754 :       virial%pv_virial = virial%pv_virial + pv_virial_thread
     190             : 
     191          58 :       CALL neighbor_list_iterator_release(nl_iterator)
     192             : 
     193          58 :       DEALLOCATE (dodisp, ghost, floating, atomnumber, radd2, c6d2)
     194             : 
     195          58 :       CALL timestop(handle)
     196             : 
     197         116 :    END SUBROUTINE calculate_dispersion_d2_pairpot
     198             : 
     199             : ! **************************************************************************************************
     200             : !> \brief ...
     201             : !> \param z ...
     202             : !> \param c6 ...
     203             : !> \param r ...
     204             : !> \param found ...
     205             : ! **************************************************************************************************
     206          38 :    SUBROUTINE dftd2_param(z, c6, r, found)
     207             : 
     208             :       INTEGER, INTENT(in)                                :: z
     209             :       REAL(KIND=dp), INTENT(inout)                       :: c6, r
     210             :       LOGICAL, INTENT(inout)                             :: found
     211             : 
     212             :       REAL(KIND=dp), DIMENSION(54), PARAMETER :: c6val = (/0.14_dp, 0.08_dp, 1.61_dp, 1.61_dp, &
     213             :          3.13_dp, 1.75_dp, 1.23_dp, 0.70_dp, 0.75_dp, 0.63_dp, 5.71_dp, 5.71_dp, 10.79_dp, 9.23_dp,&
     214             :          7.84_dp, 5.57_dp, 5.07_dp, 4.61_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, &
     215             :          10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 10.80_dp, 16.99_dp, 17.10_dp, &
     216             :          16.37_dp, 12.64_dp, 12.47_dp, 12.01_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, &
     217             :          24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 24.67_dp, 37.32_dp, 38.71_dp, &
     218             :          38.44_dp, 31.74_dp, 31.50_dp, 29.99_dp/)
     219             :       REAL(KIND=dp), DIMENSION(54), PARAMETER :: rval = (/1.001_dp, 1.012_dp, 0.825_dp, 1.408_dp, &
     220             :          1.485_dp, 1.452_dp, 1.397_dp, 1.342_dp, 1.287_dp, 1.243_dp, 1.144_dp, 1.364_dp, 1.639_dp, &
     221             :          1.716_dp, 1.705_dp, 1.683_dp, 1.639_dp, 1.595_dp, 1.485_dp, 1.474_dp, 1.562_dp, 1.562_dp, &
     222             :          1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.562_dp, 1.650_dp, &
     223             :          1.727_dp, 1.760_dp, 1.771_dp, 1.749_dp, 1.727_dp, 1.628_dp, 1.606_dp, 1.639_dp, 1.639_dp, &
     224             :          1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.639_dp, 1.672_dp, &
     225             :          1.804_dp, 1.881_dp, 1.892_dp, 1.892_dp, 1.881_dp/)
     226             : 
     227             : !
     228             : ! GRIMME DISPERSION PARAMETERS
     229             : ! Stefan Grimme, Semiempirical GGA-Type Density Functional Constructed
     230             : !                with a Long-Range Dispersion Correction, J. Comp. Chem. 27: 1787-1799 (2006)
     231             : ! doi:10.1002/jcc.20495
     232             : !
     233             : ! Conversion factor [Jnm^6mol^-1] -> [a.u.] : 17.34527758021901
     234             : ! Conversion factor [A] -> [a.u.] : 1.889726132885643
     235             : !
     236             : ! C6 values in [Jnm^6/mol]
     237             : ! vdW radii [A]
     238             : 
     239          38 :       IF (z > 0 .AND. z <= 54) THEN
     240          38 :          found = .TRUE.
     241          38 :          c6 = c6val(z)*1000._dp*bohr**6/kjmol
     242          38 :          r = rval(z)*bohr
     243             :       ELSE
     244           0 :          found = .FALSE.
     245             :       END IF
     246             : 
     247          38 :    END SUBROUTINE dftd2_param
     248             : 
     249             : ! **************************************************************************************************
     250             : 
     251             : END MODULE qs_dispersion_d2