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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculation of dispersion in DFTB
      10              : !> \author JGH
      11              : ! **************************************************************************************************
      12              : MODULE qs_dftb_dispersion
      13              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14              :                                               get_atomic_kind_set
      15              :    USE atprop_types,                    ONLY: atprop_array_init,&
      16              :                                               atprop_type
      17              :    USE cp_control_types,                ONLY: dft_control_type,&
      18              :                                               dftb_control_type
      19              :    USE input_constants,                 ONLY: dispersion_d2,&
      20              :                                               dispersion_d3,&
      21              :                                               dispersion_d3bj,&
      22              :                                               dispersion_uff
      23              :    USE kinds,                           ONLY: dp
      24              :    USE message_passing,                 ONLY: mp_para_env_type
      25              :    USE particle_types,                  ONLY: particle_type
      26              :    USE qs_dftb_types,                   ONLY: qs_dftb_atom_type,&
      27              :                                               qs_dftb_pairpot_type
      28              :    USE qs_dftb_utils,                   ONLY: get_dftb_atom_param
      29              :    USE qs_dispersion_pairpot,           ONLY: calculate_dispersion_pairpot
      30              :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      31              :    USE qs_energy_types,                 ONLY: qs_energy_type
      32              :    USE qs_environment_types,            ONLY: get_qs_env,&
      33              :                                               qs_environment_type
      34              :    USE qs_force_types,                  ONLY: qs_force_type
      35              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      36              :                                               qs_kind_type
      37              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      38              :                                               neighbor_list_iterate,&
      39              :                                               neighbor_list_iterator_create,&
      40              :                                               neighbor_list_iterator_p_type,&
      41              :                                               neighbor_list_iterator_release,&
      42              :                                               neighbor_list_set_p_type
      43              :    USE virial_methods,                  ONLY: virial_pair_force
      44              :    USE virial_types,                    ONLY: virial_type
      45              : #include "./base/base_uses.f90"
      46              : 
      47              :    IMPLICIT NONE
      48              : 
      49              :    PRIVATE
      50              : 
      51              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dftb_dispersion'
      52              : 
      53              :    PUBLIC :: calculate_dftb_dispersion
      54              : 
      55              : CONTAINS
      56              : 
      57              : ! **************************************************************************************************
      58              : !> \brief ...
      59              : !> \param qs_env ...
      60              : !> \param para_env ...
      61              : !> \param calculate_forces ...
      62              : ! **************************************************************************************************
      63         2822 :    SUBROUTINE calculate_dftb_dispersion(qs_env, para_env, calculate_forces)
      64              : 
      65              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      66              :       TYPE(mp_para_env_type), POINTER                    :: para_env
      67              :       LOGICAL, INTENT(IN)                                :: calculate_forces
      68              : 
      69              :       TYPE(dft_control_type), POINTER                    :: dft_control
      70              :       TYPE(dftb_control_type), POINTER                   :: dftb_control
      71              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
      72              :       TYPE(qs_energy_type), POINTER                      :: energy
      73              : 
      74              :       CALL get_qs_env(qs_env=qs_env, &
      75              :                       energy=energy, &
      76         2822 :                       dft_control=dft_control)
      77              : 
      78         2822 :       energy%dispersion = 0._dp
      79              : 
      80         2822 :       dftb_control => dft_control%qs_control%dftb_control
      81         2822 :       IF (dftb_control%dispersion) THEN
      82         2926 :          SELECT CASE (dftb_control%dispersion_type)
      83              :          CASE (dispersion_uff)
      84         1376 :             CALL calculate_dispersion_uff(qs_env, para_env, calculate_forces)
      85              :          CASE (dispersion_d3, dispersion_d3bj, dispersion_d2)
      86          174 :             CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
      87              :             CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
      88          174 :                                               energy%dispersion, calculate_forces)
      89              :          CASE DEFAULT
      90         1550 :             CPABORT("")
      91              :          END SELECT
      92              :       END IF
      93              : 
      94         2822 :    END SUBROUTINE calculate_dftb_dispersion
      95              : 
      96              : ! **************************************************************************************************
      97              : !> \brief ...
      98              : !> \param qs_env ...
      99              : !> \param para_env ...
     100              : !> \param calculate_forces ...
     101              : ! **************************************************************************************************
     102         1376 :    SUBROUTINE calculate_dispersion_uff(qs_env, para_env, calculate_forces)
     103              : 
     104              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     105              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     106              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     107              : 
     108              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_uff'
     109              : 
     110              :       INTEGER                                            :: atom_a, atom_b, handle, iatom, ikind, &
     111              :                                                             jatom, jkind, nkind
     112         1376 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind
     113              :       LOGICAL                                            :: use_virial
     114         1376 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: define_kind
     115         1376 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: rc_kind
     116              :       REAL(KIND=dp)                                      :: a, b, c, devdw, dij, dr, eij, evdw, fac, &
     117              :                                                             rc, x0ij, xij, xp
     118              :       REAL(KIND=dp), DIMENSION(3)                        :: fdij, rij
     119         1376 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     120              :       TYPE(atprop_type), POINTER                         :: atprop
     121              :       TYPE(dft_control_type), POINTER                    :: dft_control
     122              :       TYPE(dftb_control_type), POINTER                   :: dftb_control
     123              :       TYPE(neighbor_list_iterator_p_type), &
     124         1376 :          DIMENSION(:), POINTER                           :: nl_iterator
     125              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     126         1376 :          POINTER                                         :: sab_vdw
     127         1376 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     128              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_kind_a
     129              :       TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
     130         1376 :          POINTER                                         :: dftb_potential
     131              :       TYPE(qs_dftb_pairpot_type), POINTER                :: dftb_param_ij
     132              :       TYPE(qs_energy_type), POINTER                      :: energy
     133         1376 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     134         1376 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     135              :       TYPE(virial_type), POINTER                         :: virial
     136              : 
     137         1376 :       CALL timeset(routineN, handle)
     138              : 
     139         1376 :       NULLIFY (atomic_kind_set, sab_vdw, atprop)
     140              : 
     141              :       CALL get_qs_env(qs_env=qs_env, &
     142              :                       energy=energy, &
     143              :                       atomic_kind_set=atomic_kind_set, &
     144              :                       qs_kind_set=qs_kind_set, &
     145              :                       virial=virial, atprop=atprop, &
     146         1376 :                       dft_control=dft_control)
     147              : 
     148         1376 :       energy%dispersion = 0._dp
     149              : 
     150         1376 :       dftb_control => dft_control%qs_control%dftb_control
     151              : 
     152         1376 :       IF (dftb_control%dispersion) THEN
     153              : 
     154         1376 :          NULLIFY (dftb_potential)
     155         1376 :          CALL get_qs_env(qs_env=qs_env, dftb_potential=dftb_potential)
     156         1376 :          IF (calculate_forces) THEN
     157          450 :             NULLIFY (force, particle_set)
     158              :             CALL get_qs_env(qs_env=qs_env, &
     159              :                             particle_set=particle_set, &
     160          450 :                             force=force)
     161          450 :             CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     162          450 :             use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     163              :          ELSE
     164              :             use_virial = .FALSE.
     165              :          END IF
     166         1376 :          nkind = SIZE(atomic_kind_set)
     167         6880 :          ALLOCATE (define_kind(nkind), rc_kind(nkind))
     168         4156 :          DO ikind = 1, nkind
     169         2780 :             CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_kind_a)
     170         4156 :             CALL get_dftb_atom_param(dftb_kind_a, defined=define_kind(ikind), rcdisp=rc_kind(ikind))
     171              :          END DO
     172              : 
     173         1376 :          evdw = 0._dp
     174              : 
     175         1376 :          IF (atprop%energy) THEN
     176           58 :             CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
     177           58 :             CALL atprop_array_init(atprop%atevdw, natom=SIZE(particle_set))
     178              :          END IF
     179              : 
     180         1376 :          CALL get_qs_env(qs_env=qs_env, sab_vdw=sab_vdw)
     181         1376 :          CALL neighbor_list_iterator_create(nl_iterator, sab_vdw)
     182      8663828 :          DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     183      8662452 :             CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
     184      8662452 :             IF ((.NOT. define_kind(ikind)) .OR. (.NOT. define_kind(jkind))) CYCLE
     185      8662452 :             rc = rc_kind(ikind) + rc_kind(jkind)
     186              :             ! vdW potential
     187     34649808 :             dr = SQRT(SUM(rij(:)**2))
     188      8663828 :             IF (dr <= rc .AND. dr > 0.001_dp) THEN
     189      8655844 :                fac = 1._dp
     190      8655844 :                IF (iatom == jatom) fac = 0.5_dp
     191              :                ! retrieve information on potential
     192      8655844 :                dftb_param_ij => dftb_potential(ikind, jkind)
     193              :                ! vdW parameter
     194      8655844 :                xij = dftb_param_ij%xij
     195      8655844 :                dij = dftb_param_ij%dij
     196      8655844 :                x0ij = dftb_param_ij%x0ij
     197      8655844 :                a = dftb_param_ij%a
     198      8655844 :                b = dftb_param_ij%b
     199      8655844 :                c = dftb_param_ij%c
     200      8655844 :                IF (dr > x0ij) THEN
     201              :                   ! This is the standard London contribution.
     202              :                   ! UFF1 - Eq. 20 (long-range)
     203      8640421 :                   xp = xij/dr
     204      8640421 :                   eij = dij*(-2._dp*xp**6 + xp**12)*fac
     205      8640421 :                   evdw = evdw + eij
     206      8640421 :                   IF (calculate_forces .AND. (dr > 0.001_dp)) THEN
     207      2058704 :                      devdw = dij*12._dp*(xp**6 - xp**12)/dr*fac
     208      2058704 :                      atom_a = atom_of_kind(iatom)
     209      2058704 :                      atom_b = atom_of_kind(jatom)
     210      8234816 :                      fdij(:) = devdw*rij(:)/dr
     211              :                      force(ikind)%dispersion(:, atom_a) = &
     212      8234816 :                         force(ikind)%dispersion(:, atom_a) - fdij(:)
     213              :                      force(jkind)%dispersion(:, atom_b) = &
     214      8234816 :                         force(jkind)%dispersion(:, atom_b) + fdij(:)
     215              :                   END IF
     216              :                ELSE
     217              :                   ! Shorter distance.
     218              :                   ! London contribution should converge to a finite value.
     219              :                   ! Using a parabola of the form (y = A - Bx**5 -Cx**10).
     220              :                   ! Analytic parameters by forcing energy, first and second
     221              :                   ! derivatives to be continuous.
     222        15423 :                   eij = (A - B*dr**5 - C*dr**10)*fac
     223        15423 :                   evdw = evdw + eij
     224        15423 :                   IF (calculate_forces .AND. (dr > 0.001_dp)) THEN
     225         6139 :                      atom_a = atom_of_kind(iatom)
     226         6139 :                      atom_b = atom_of_kind(jatom)
     227         6139 :                      devdw = (-5*B*dr**4 - 10*C*dr**9)*fac
     228        24556 :                      fdij(:) = devdw*rij(:)/dr
     229              :                      force(ikind)%dispersion(:, atom_a) = &
     230        24556 :                         force(ikind)%dispersion(:, atom_a) - fdij(:)
     231              :                      force(jkind)%dispersion(:, atom_b) = &
     232        24556 :                         force(jkind)%dispersion(:, atom_b) + fdij(:)
     233              :                   END IF
     234              :                END IF
     235      8655844 :                IF (atprop%energy) THEN
     236      1368540 :                   atprop%atevdw(iatom) = atprop%atevdw(iatom) + 0.5_dp*eij
     237      1368540 :                   atprop%atevdw(jatom) = atprop%atevdw(jatom) + 0.5_dp*eij
     238              :                END IF
     239      8655844 :                IF (calculate_forces .AND. (dr > 0.001_dp) .AND. use_virial) THEN
     240      1094831 :                   CALL virial_pair_force(virial%pv_virial, -1._dp, fdij, rij)
     241              :                END IF
     242              :             END IF
     243              :          END DO
     244         1376 :          CALL neighbor_list_iterator_release(nl_iterator)
     245              : 
     246              :          ! set dispersion energy
     247         1376 :          CALL para_env%sum(evdw)
     248         1376 :          energy%dispersion = evdw
     249              : 
     250              :       END IF
     251              : 
     252         1376 :       CALL timestop(handle)
     253              : 
     254         2752 :    END SUBROUTINE calculate_dispersion_uff
     255              : 
     256              : END MODULE qs_dftb_dispersion
     257              : 
        

Generated by: LCOV version 2.0-1