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 dispersion using pair potentials
      10             : !> \author Johann Pototschnig
      11             : ! **************************************************************************************************
      12             : MODULE qs_dispersion_d4
      13             : 
      14             :    USE atomic_kind_types, ONLY: get_atomic_kind
      15             :    USE cell_types, ONLY: cell_type, &
      16             :                          get_cell
      17             : #if defined(__DFTD4)
      18             :    USE dftd4, ONLY: d4_model, &
      19             :                     damping_param, &
      20             :                     get_dispersion, &
      21             :                     get_rational_damping, &
      22             :                     new, &
      23             :                     new_d4_model, &
      24             :                     realspace_cutoff, &
      25             :                     structure_type
      26             : #endif
      27             :    USE kinds, ONLY: dp
      28             :    USE particle_types, ONLY: particle_type
      29             :    USE qs_dispersion_types, ONLY: qs_dispersion_type
      30             :    USE qs_force_types, ONLY: qs_force_type
      31             :    USE message_passing, ONLY: mp_para_env_type
      32             : #include "./base/base_uses.f90"
      33             : 
      34             :    IMPLICIT NONE
      35             : 
      36             :    PRIVATE
      37             : 
      38             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_d4'
      39             : 
      40             :    PUBLIC :: calculate_dispersion_d4_pairpot
      41             : 
      42             : ! **************************************************************************************************
      43             : 
      44             : CONTAINS
      45             : 
      46             : #if defined(__DFTD4)
      47             :    ! **************************************************************************************************
      48             : !> \brief ...
      49             : !> \param dispersion_env ...
      50             : !> \param particle_set ...
      51             : !> \param cell ...
      52             : !> \param para_env ...
      53             : !> \param energy ...
      54             : !> \param force ...
      55             : !> \param atom_of_kind ...
      56             : !> \param virial ...
      57             : ! **************************************************************************************************
      58           6 :    SUBROUTINE calculate_dispersion_d4_pairpot(dispersion_env, particle_set, cell, para_env, energy, force, atom_of_kind, virial)
      59             : 
      60             :       TYPE(qs_dispersion_type), INTENT(IN), POINTER      :: dispersion_env
      61             :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
      62             :          POINTER                                         :: particle_set
      63             :       TYPE(cell_type), INTENT(IN), POINTER               :: cell
      64             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      65             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
      66             :       TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
      67             :          OPTIONAL, POINTER                               :: force
      68             :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: atom_of_kind
      69             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
      70             :          OPTIONAL                                        :: virial
      71             : 
      72             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_d4_pairpot'
      73             : 
      74          12 :       class(damping_param), ALLOCATABLE                  :: param
      75           6 :       TYPE(d4_model)                                     :: disp
      76           6 :       TYPE(structure_type)                               :: mol
      77             :       TYPE(realspace_cutoff)                             :: cutoff
      78             : 
      79             :       INTEGER                                            :: iatom, natom, ind_atom
      80             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: el_num
      81           6 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gradient, xyz
      82             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: stress
      83             :       INTEGER                                            :: handle, ikind
      84             :       !REAL(KIND=dp), DIMENSION(3, 3)                    :: h
      85             :       INTEGER, DIMENSION(3)                              :: periodic
      86             :       LOGICAL, DIMENSION(3)                              :: lperiod
      87             : 
      88           6 :       CALL timeset(routineN, handle)
      89             : 
      90             :       !get information about particles
      91           6 :       natom = SIZE(particle_set)
      92          18 :       ALLOCATE (xyz(3, natom))
      93          18 :       ALLOCATE (el_num(natom))
      94          54 :       DO iatom = 1, natom
      95         192 :          xyz(:, iatom) = particle_set(iatom)%r(:)
      96          48 :          CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
      97          54 :          el_num(iatom) = ikind
      98             :       END DO
      99             : 
     100             :       !get information about cell / lattice
     101           6 :       CALL get_cell(cell=cell, periodic=periodic)
     102           6 :       lperiod(1) = periodic(1) == 1
     103           6 :       lperiod(2) = periodic(2) == 1
     104           6 :       lperiod(3) = periodic(3) == 1
     105             : 
     106             :       !prepare for the call to the dispersion function
     107           6 :       CALL new(mol, el_num, xyz, lattice=cell%hmat, periodic=lperiod)
     108           6 :       CALL new_d4_model(disp, mol)
     109           6 :       CALL get_rational_damping(dispersion_env%ref_functional, param, s9=dispersion_env%s9)
     110             : 
     111             :       ! Coordination number cutoff
     112           6 :       cutoff%cn = dispersion_env%rc_cn
     113             :       ! Two-body interaction cutoff
     114           6 :       cutoff%disp2 = dispersion_env%rc_d4*2._dp
     115             :       ! Three-body interaction cutoff
     116           6 :       cutoff%disp3 = dispersion_env%rc_disp*2._dp
     117             : 
     118             :       !> Wrapper to handle the evaluation of dispersion energy and derivatives
     119           6 :       IF (PRESENT(force)) THEN
     120           4 :          IF (PRESENT(atom_of_kind)) THEN
     121          12 :             ALLOCATE (gradient(3, mol%nat))
     122           4 :             CALL get_dispersion(mol, disp, param, cutoff, energy, gradient, stress)
     123           4 :             IF (PRESENT(virial)) THEN
     124          26 :                virial = -1.00_dp*stress
     125           2 :                IF (para_env%num_pe > 1 .AND. para_env%mepos > 0) virial = 0.00_dp
     126             :             END IF
     127          32 :             DO iatom = 1, natom
     128          28 :                ikind = el_num(iatom)
     129          28 :                ind_atom = atom_of_kind(iatom)
     130         116 :                force(ikind)%dispersion(:, ind_atom) = force(ikind)%dispersion(:, ind_atom) + gradient(:, iatom)
     131             :             END DO
     132           4 :             DEALLOCATE (gradient)
     133             :          ELSE
     134           0 :             CPABORT("missing atom_of_kind")
     135             :          END IF
     136             :       ELSE
     137           2 :          IF (PRESENT(virial)) THEN
     138           0 :             CPABORT("missing force for virial term")
     139             :          END IF
     140           2 :          CALL get_dispersion(mol, disp, param, cutoff, energy)
     141             :       END IF
     142             : 
     143             :       !dispersion energy is computed by every MPI process
     144           6 :       IF (para_env%num_pe > 1 .AND. para_env%mepos > 0) energy = 0.00_dp
     145             : 
     146           6 :       DEALLOCATE (el_num, xyz)
     147             : 
     148           6 :       CALL timestop(handle)
     149             : 
     150          18 :    END SUBROUTINE calculate_dispersion_d4_pairpot
     151             : 
     152             : #else
     153             : 
     154             :    ! **************************************************************************************************
     155             : !> \brief ...
     156             : !> \param dispersion_env ...
     157             : !> \param particle_set ...
     158             : !> \param cell ...
     159             : !> \param para_env ...
     160             : !> \param energy ...
     161             : !> \param force ...
     162             : !> \param atom_of_kind ...
     163             : !> \param virial ...
     164             : ! **************************************************************************************************
     165             :    SUBROUTINE calculate_dispersion_d4_pairpot(dispersion_env, particle_set, cell, para_env, energy, force, atom_of_kind, virial)
     166             : 
     167             :       TYPE(qs_dispersion_type), INTENT(IN), POINTER      :: dispersion_env
     168             :       TYPE(particle_type), DIMENSION(:), INTENT(IN), &
     169             :          POINTER                                         :: particle_set
     170             :       TYPE(cell_type), INTENT(IN), POINTER               :: cell
     171             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     172             :       REAL(KIND=dp), INTENT(OUT)                         :: energy
     173             :       TYPE(qs_force_type), DIMENSION(:), INTENT(INOUT), &
     174             :          OPTIONAL, POINTER                               :: force
     175             :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: atom_of_kind
     176             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
     177             :          OPTIONAL                                        :: virial
     178             : 
     179             :       MARK_USED(dispersion_env)
     180             :       MARK_USED(particle_set)
     181             :       MARK_USED(cell)
     182             :       MARK_USED(para_env)
     183             :       MARK_USED(energy)
     184             :       MARK_USED(force)
     185             :       MARK_USED(atom_of_kind)
     186             :       MARK_USED(virial)
     187             : 
     188             :       CPABORT("Build without DFTD4")
     189             : 
     190             :    END SUBROUTINE calculate_dispersion_d4_pairpot
     191             : 
     192             : #endif
     193             : 
     194             : END MODULE qs_dispersion_d4