LCOV - code coverage report
Current view: top level - src - qs_dispersion_pairpot.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 84.5 % 252 213
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 using pair potentials
      10              : !> \author JGH
      11              : ! **************************************************************************************************
      12              : MODULE qs_dispersion_pairpot
      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 bibliography,                    ONLY: Caldeweyher2017,&
      20              :                                               Caldeweyher2019,&
      21              :                                               Caldeweyher2020,&
      22              :                                               Goerigk2017,&
      23              :                                               cite_reference,&
      24              :                                               grimme2006,&
      25              :                                               grimme2010,&
      26              :                                               grimme2011
      27              :    USE cell_types,                      ONLY: cell_type
      28              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29              :                                               cp_logger_get_default_io_unit,&
      30              :                                               cp_logger_type
      31              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      32              :                                               cp_print_key_unit_nr
      33              :    USE cp_parser_methods,               ONLY: parser_get_next_line,&
      34              :                                               parser_get_object
      35              :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      36              :                                               parser_create,&
      37              :                                               parser_release
      38              :    USE eeq_input,                       ONLY: read_eeq_param
      39              :    USE input_constants,                 ONLY: vdw_pairpot_dftd2,&
      40              :                                               vdw_pairpot_dftd3,&
      41              :                                               vdw_pairpot_dftd3bj,&
      42              :                                               vdw_pairpot_dftd4,&
      43              :                                               xc_vdw_fun_pairpot
      44              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      45              :                                               section_vals_type,&
      46              :                                               section_vals_val_get
      47              :    USE kinds,                           ONLY: default_string_length,&
      48              :                                               dp
      49              :    USE message_passing,                 ONLY: mp_para_env_type
      50              :    USE physcon,                         ONLY: bohr,&
      51              :                                               kcalmol,&
      52              :                                               kjmol
      53              :    USE qs_dispersion_cnum,              ONLY: get_cn_radius,&
      54              :                                               setcn,&
      55              :                                               seten,&
      56              :                                               setr0ab,&
      57              :                                               setrcov
      58              :    USE qs_dispersion_d2,                ONLY: calculate_dispersion_d2_pairpot,&
      59              :                                               dftd2_param
      60              :    USE qs_dispersion_d3,                ONLY: calculate_dispersion_d3_pairpot,&
      61              :                                               dftd3_c6_param
      62              :    USE qs_dispersion_d4,                ONLY: calculate_dispersion_d4_pairpot
      63              :    USE qs_dispersion_types,             ONLY: dftd2_pp,&
      64              :                                               dftd3_pp,&
      65              :                                               dftd4_pp,&
      66              :                                               qs_atom_dispersion_type,&
      67              :                                               qs_dispersion_type
      68              :    USE qs_environment_types,            ONLY: get_qs_env,&
      69              :                                               qs_environment_type
      70              :    USE qs_force_types,                  ONLY: qs_force_type
      71              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      72              :                                               qs_kind_type,&
      73              :                                               set_qs_kind
      74              :    USE virial_types,                    ONLY: virial_type
      75              : #include "./base/base_uses.f90"
      76              : 
      77              :    IMPLICIT NONE
      78              : 
      79              :    PRIVATE
      80              : 
      81              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_pairpot'
      82              : 
      83              :    PUBLIC :: qs_dispersion_pairpot_init, calculate_dispersion_pairpot
      84              : 
      85              : ! **************************************************************************************************
      86              : 
      87              : CONTAINS
      88              : 
      89              : ! **************************************************************************************************
      90              : !> \brief ...
      91              : !> \param atomic_kind_set ...
      92              : !> \param qs_kind_set ...
      93              : !> \param dispersion_env ...
      94              : !> \param pp_section ...
      95              : !> \param para_env ...
      96              : ! **************************************************************************************************
      97         1080 :    SUBROUTINE qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
      98              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      99              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     100              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     101              :       TYPE(section_vals_type), OPTIONAL, POINTER         :: pp_section
     102              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     103              : 
     104              :       CHARACTER(len=*), PARAMETER :: routineN = 'qs_dispersion_pairpot_init'
     105              : 
     106              :       CHARACTER(LEN=2)                                   :: symbol
     107              :       CHARACTER(LEN=default_string_length)               :: aname, filename
     108              :       CHARACTER(LEN=default_string_length), &
     109         1080 :          DIMENSION(:), POINTER                           :: tmpstringlist
     110              :       INTEGER                                            :: elem, handle, i, ikind, j, max_elem, &
     111              :                                                             maxc, n_rep, nkind, nl, vdw_pp_type, &
     112              :                                                             vdw_type
     113         1080 :       INTEGER, DIMENSION(:), POINTER                     :: exlist
     114              :       LOGICAL                                            :: at_end, explicit, found, is_available
     115              :       REAL(KIND=dp)                                      :: dum
     116              :       TYPE(qs_atom_dispersion_type), POINTER             :: disp
     117              :       TYPE(section_vals_type), POINTER                   :: eeq_section
     118              : 
     119         1080 :       CALL timeset(routineN, handle)
     120              : 
     121         1080 :       nkind = SIZE(atomic_kind_set)
     122              : 
     123         1080 :       vdw_type = dispersion_env%type
     124         1080 :       SELECT CASE (vdw_type)
     125              :       CASE DEFAULT
     126              :          ! do nothing
     127              :       CASE (xc_vdw_fun_pairpot)
     128              :          ! setup information on pair potentials
     129         1080 :          vdw_pp_type = dispersion_env%type
     130         1080 :          SELECT CASE (dispersion_env%pp_type)
     131              :          CASE DEFAULT
     132              :             ! do nothing
     133              :          CASE (vdw_pairpot_dftd2)
     134           22 :             CALL cite_reference(Grimme2006)
     135           60 :             DO ikind = 1, nkind
     136           38 :                CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
     137           38 :                ALLOCATE (disp)
     138           38 :                disp%type = dftd2_pp
     139              :                ! get filename of parameter file
     140           38 :                filename = dispersion_env%parameter_file_name
     141              :                ! check for local parameters
     142           38 :                found = .FALSE.
     143           38 :                IF (PRESENT(pp_section)) THEN
     144           32 :                   CALL section_vals_val_get(pp_section, "ATOMPARM", n_rep_val=n_rep)
     145           32 :                   DO i = 1, n_rep
     146              :                      CALL section_vals_val_get(pp_section, "ATOMPARM", i_rep_val=i, &
     147            0 :                                                c_vals=tmpstringlist)
     148           32 :                      IF (TRIM(tmpstringlist(1)) == TRIM(symbol)) THEN
     149              :                         ! we assume the parameters are in atomic units!
     150            0 :                         READ (tmpstringlist(2), *) disp%c6
     151            0 :                         READ (tmpstringlist(3), *) disp%vdw_radii
     152            0 :                         found = .TRUE.
     153            0 :                         EXIT
     154              :                      END IF
     155              :                   END DO
     156              :                END IF
     157           38 :                IF (.NOT. found) THEN
     158              :                   ! check for internal parameters
     159           38 :                   CALL dftd2_param(elem, disp%c6, disp%vdw_radii, found)
     160              :                END IF
     161           38 :                IF (.NOT. found) THEN
     162              :                   ! check on file
     163            0 :                   INQUIRE (FILE=filename, EXIST=is_available)
     164            0 :                   IF (is_available) THEN
     165            0 :                      BLOCK
     166              :                         TYPE(cp_parser_type) :: parser
     167            0 :                         CALL parser_create(parser, filename, para_env=para_env)
     168              :                         DO
     169              :                            at_end = .FALSE.
     170            0 :                            CALL parser_get_next_line(parser, 1, at_end)
     171            0 :                            IF (at_end) EXIT
     172            0 :                            CALL parser_get_object(parser, aname)
     173            0 :                            IF (TRIM(aname) == TRIM(symbol)) THEN
     174            0 :                               CALL parser_get_object(parser, disp%c6)
     175              :                               ! we have to change the units J*nm^6*mol^-1 -> Hartree*Bohr^6
     176            0 :                               disp%c6 = disp%c6*1000._dp*bohr**6/kjmol
     177            0 :                               CALL parser_get_object(parser, disp%vdw_radii)
     178            0 :                               disp%vdw_radii = disp%vdw_radii*bohr
     179            0 :                               found = .TRUE.
     180            0 :                               EXIT
     181              :                            END IF
     182              :                         END DO
     183            0 :                         CALL parser_release(parser)
     184              :                      END BLOCK
     185              :                   END IF
     186              :                END IF
     187           38 :                IF (found) THEN
     188           38 :                   disp%defined = .TRUE.
     189              :                ELSE
     190            0 :                   disp%defined = .FALSE.
     191              :                END IF
     192              :                ! Check if the parameter is defined
     193           38 :                IF (.NOT. disp%defined) &
     194              :                   CALL cp_abort(__LOCATION__, &
     195              :                                 "Dispersion parameters for element ("//TRIM(symbol)//") are not defined! "// &
     196              :                                 "Please provide a valid set of parameters through the input section or "// &
     197            0 :                                 "through an external file! ")
     198           98 :                CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
     199              :             END DO
     200              :          CASE (vdw_pairpot_dftd3, vdw_pairpot_dftd3bj)
     201              :             !DFT-D3 Method initial setup
     202          392 :             CALL cite_reference(Grimme2010)
     203          392 :             CALL cite_reference(Grimme2011)
     204          392 :             CALL cite_reference(Goerigk2017)
     205          392 :             max_elem = 94
     206          392 :             maxc = 5
     207          392 :             dispersion_env%max_elem = max_elem
     208          392 :             dispersion_env%maxc = maxc
     209          392 :             ALLOCATE (dispersion_env%maxci(max_elem))
     210          392 :             ALLOCATE (dispersion_env%c6ab(max_elem, max_elem, maxc, maxc, 3))
     211          392 :             ALLOCATE (dispersion_env%r0ab(max_elem, max_elem))
     212          392 :             ALLOCATE (dispersion_env%rcov(max_elem))
     213          392 :             ALLOCATE (dispersion_env%eneg(max_elem))
     214          392 :             ALLOCATE (dispersion_env%r2r4(max_elem))
     215          392 :             ALLOCATE (dispersion_env%cn(max_elem))
     216              : 
     217              :             ! get filename of parameter file
     218          392 :             filename = dispersion_env%parameter_file_name
     219          392 :             CALL dftd3_c6_param(dispersion_env%c6ab, dispersion_env%maxci, filename, para_env)
     220          392 :             CALL setr0ab(dispersion_env%r0ab, dispersion_env%rcov, dispersion_env%r2r4)
     221              :             ! Electronegativity
     222          392 :             CALL seten(dispersion_env%eneg)
     223              :             ! the default coordination numbers
     224          392 :             CALL setcn(dispersion_env%cn)
     225              :             ! scale r4/r2 values of the atoms by sqrt(Z)
     226              :             ! sqrt is also globally close to optimum
     227              :             ! together with the factor 1/2 this yield reasonable
     228              :             ! c8 for he, ne and ar. for larger Z, C8 becomes too large
     229              :             ! which effectively mimics higher R^n terms neglected due
     230              :             ! to stability reasons
     231        37240 :             DO i = 1, max_elem
     232        36848 :                dum = 0.5_dp*dispersion_env%r2r4(i)*REAL(i, dp)**0.5_dp
     233              :                ! store it as sqrt because the geom. av. is taken
     234        37240 :                dispersion_env%r2r4(i) = SQRT(dum)
     235              :             END DO
     236              :             ! parameters
     237          392 :             dispersion_env%k1 = 16.0_dp
     238          392 :             dispersion_env%k2 = 4._dp/3._dp
     239              :             ! reasonable choices are between 3 and 5
     240              :             ! this gives smoth curves with maxima around the integer values
     241              :             ! k3=3 give for CN=0 a slightly smaller value than computed
     242              :             ! for the free atom. This also yields to larger CN for atoms
     243              :             ! in larger molecules but with the same chem. environment
     244              :             ! which is physically not right
     245              :             ! values >5 might lead to bumps in the potential
     246          392 :             dispersion_env%k3 = -4._dp
     247        37240 :             dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*bohr
     248              :             ! alpha default parameter
     249          392 :             dispersion_env%alp = 14._dp
     250              :             !
     251         1292 :             DO ikind = 1, nkind
     252          900 :                CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
     253          900 :                ALLOCATE (disp)
     254          900 :                disp%type = dftd3_pp
     255          900 :                IF (elem <= 94) THEN
     256          900 :                   disp%defined = .TRUE.
     257              :                ELSE
     258              :                   disp%defined = .FALSE.
     259              :                END IF
     260          900 :                IF (.NOT. disp%defined) &
     261              :                   CALL cp_abort(__LOCATION__, &
     262              :                                 "Dispersion parameters for element ("//TRIM(symbol)//") are not defined! "// &
     263              :                                 "Please provide a valid set of parameters through the input section or "// &
     264            0 :                                 "through an external file! ")
     265         2192 :                CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
     266              :             END DO
     267              : 
     268          392 :             IF (PRESENT(pp_section)) THEN
     269              :                ! Check for coordination numbers
     270           66 :                CALL section_vals_val_get(pp_section, "KIND_COORDINATION_NUMBERS", n_rep_val=n_rep)
     271           66 :                IF (n_rep > 0) THEN
     272           24 :                   ALLOCATE (dispersion_env%cnkind(n_rep))
     273           12 :                   DO i = 1, n_rep
     274              :                      CALL section_vals_val_get(pp_section, "KIND_COORDINATION_NUMBERS", i_rep_val=i, &
     275            6 :                                                c_vals=tmpstringlist)
     276            6 :                      READ (tmpstringlist(1), *) dispersion_env%cnkind(i)%cnum
     277           12 :                      READ (tmpstringlist(2), *) dispersion_env%cnkind(i)%kind
     278              :                   END DO
     279              :                END IF
     280           66 :                CALL section_vals_val_get(pp_section, "ATOM_COORDINATION_NUMBERS", n_rep_val=n_rep)
     281           66 :                IF (n_rep > 0) THEN
     282           10 :                   ALLOCATE (dispersion_env%cnlist(n_rep))
     283            6 :                   DO i = 1, n_rep
     284              :                      CALL section_vals_val_get(pp_section, "ATOM_COORDINATION_NUMBERS", i_rep_val=i, &
     285            4 :                                                c_vals=tmpstringlist)
     286            4 :                      nl = SIZE(tmpstringlist)
     287           12 :                      ALLOCATE (dispersion_env%cnlist(i)%atom(nl - 1))
     288            4 :                      dispersion_env%cnlist(i)%natom = nl - 1
     289            4 :                      READ (tmpstringlist(1), *) dispersion_env%cnlist(i)%cnum
     290           14 :                      DO j = 1, nl - 1
     291           12 :                         READ (tmpstringlist(j + 1), *) dispersion_env%cnlist(i)%atom(j)
     292              :                      END DO
     293              :                   END DO
     294              :                END IF
     295              :                ! Check for exclusion lists
     296           66 :                CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND", explicit=explicit)
     297           66 :                IF (explicit) THEN
     298            2 :                   CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND", i_vals=exlist)
     299            4 :                   DO j = 1, SIZE(exlist)
     300            2 :                      ikind = exlist(j)
     301            2 :                      CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp)
     302            4 :                      disp%defined = .FALSE.
     303              :                   END DO
     304              :                END IF
     305           66 :                CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", n_rep_val=n_rep)
     306           66 :                dispersion_env%nd3_exclude_pair = n_rep
     307           66 :                IF (n_rep > 0) THEN
     308            6 :                   ALLOCATE (dispersion_env%d3_exclude_pair(n_rep, 2))
     309            6 :                   DO i = 1, n_rep
     310              :                      CALL section_vals_val_get(pp_section, "D3_EXCLUDE_KIND_PAIR", i_rep_val=i, &
     311            4 :                                                i_vals=exlist)
     312           22 :                      dispersion_env%d3_exclude_pair(i, :) = exlist
     313              :                   END DO
     314              :                END IF
     315              :             END IF
     316              :          CASE (vdw_pairpot_dftd4)
     317              :             !most checks are done by the library
     318          666 :             CALL cite_reference(Caldeweyher2017)
     319          666 :             CALL cite_reference(Caldeweyher2019)
     320          666 :             CALL cite_reference(Caldeweyher2020)
     321         2074 :             DO ikind = 1, nkind
     322         1408 :                CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol, z=elem)
     323         1408 :                ALLOCATE (disp)
     324         1408 :                disp%type = dftd4_pp
     325         1408 :                disp%defined = .TRUE.
     326         3482 :                CALL set_qs_kind(qs_kind_set(ikind), dispersion=disp)
     327              :             END DO
     328              :             ! maybe needed in cnumber calculations
     329          666 :             max_elem = 94
     330          666 :             maxc = 5
     331          666 :             dispersion_env%max_elem = max_elem
     332          666 :             dispersion_env%maxc = maxc
     333          666 :             ALLOCATE (dispersion_env%maxci(max_elem))
     334          666 :             ALLOCATE (dispersion_env%rcov(max_elem))
     335          666 :             ALLOCATE (dispersion_env%eneg(max_elem))
     336          666 :             ALLOCATE (dispersion_env%cn(max_elem))
     337              :             ! the default covalent radii
     338          666 :             CALL setrcov(dispersion_env%rcov)
     339              :             ! the default coordination numbers
     340          666 :             CALL setcn(dispersion_env%cn)
     341              :             ! Electronegativity
     342          666 :             CALL seten(dispersion_env%eneg)
     343              :             ! parameters
     344          666 :             dispersion_env%k1 = 16.0_dp
     345          666 :             dispersion_env%k2 = 4._dp/3._dp
     346          666 :             dispersion_env%k3 = -4._dp
     347        63270 :             dispersion_env%rcov = dispersion_env%k2*dispersion_env%rcov*bohr
     348          666 :             dispersion_env%alp = 14._dp
     349              :             !
     350          666 :             dispersion_env%cnfun = 3
     351          666 :             dispersion_env%rc_cn = get_cn_radius(dispersion_env)
     352         1746 :             IF (PRESENT(pp_section)) THEN
     353           26 :                eeq_section => section_vals_get_subs_vals(pp_section, "EEQ")
     354           26 :                CALL read_eeq_param(eeq_section, dispersion_env%eeq_sparam)
     355              :             END IF
     356              :          END SELECT
     357              :       END SELECT
     358              : 
     359         1080 :       CALL timestop(handle)
     360              : 
     361         1080 :    END SUBROUTINE qs_dispersion_pairpot_init
     362              : 
     363              : ! **************************************************************************************************
     364              : !> \brief ...
     365              : !> \param qs_env ...
     366              : !> \param dispersion_env ...
     367              : !> \param energy ...
     368              : !> \param calculate_forces ...
     369              : !> \param atevdw ...
     370              : ! **************************************************************************************************
     371        21602 :    SUBROUTINE calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
     372              : 
     373              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     374              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     375              :       REAL(KIND=dp), INTENT(INOUT)                       :: energy
     376              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     377              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atevdw
     378              : 
     379              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_dispersion_pairpot'
     380              : 
     381              :       INTEGER                                            :: atom_a, handle, iatom, ikind, iw, natom, &
     382              :                                                             nkind, unit_nr
     383        21602 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     384              :       LOGICAL                                            :: atenergy, atex, debugall, use_virial
     385              :       REAL(KIND=dp)                                      :: evdw, gnorm
     386        21602 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: atomic_energy
     387              :       REAL(KIND=dp), DIMENSION(3)                        :: fdij
     388              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: dvirial, pv_loc, pv_virial_thread
     389        21602 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     390              :       TYPE(atprop_type), POINTER                         :: atprop
     391              :       TYPE(cell_type), POINTER                           :: cell
     392              :       TYPE(cp_logger_type), POINTER                      :: logger
     393              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     394        21602 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     395              :       TYPE(virial_type), POINTER                         :: virial
     396              : 
     397        21602 :       energy = 0._dp
     398              :       ! make valgrind happy
     399        21602 :       use_virial = .FALSE.
     400              : 
     401        21602 :       IF (dispersion_env%type .NE. xc_vdw_fun_pairpot) THEN
     402              :          RETURN
     403              :       END IF
     404              : 
     405         5676 :       CALL timeset(routineN, handle)
     406              : 
     407         5676 :       NULLIFY (atomic_kind_set)
     408              : 
     409              :       CALL get_qs_env(qs_env=qs_env, nkind=nkind, natom=natom, atomic_kind_set=atomic_kind_set, &
     410         5676 :                       cell=cell, virial=virial, para_env=para_env, atprop=atprop)
     411              : 
     412         5676 :       debugall = dispersion_env%verbose
     413              : 
     414         5676 :       NULLIFY (logger)
     415         5676 :       logger => cp_get_default_logger()
     416         5676 :       IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
     417              :          unit_nr = cp_print_key_unit_nr(logger, dispersion_env%dftd_section, "PRINT_DFTD", &
     418          262 :                                         extension=".dftd")
     419              :       ELSE
     420         5414 :          unit_nr = -1
     421              :       END IF
     422              : 
     423              :       ! atomic energy and stress arrays
     424         5676 :       atenergy = atprop%energy
     425              :       ! external atomic energy
     426         5676 :       atex = .FALSE.
     427         5676 :       IF (PRESENT(atevdw)) THEN
     428            2 :          atex = .TRUE.
     429              :       END IF
     430              : 
     431         5676 :       IF (unit_nr > 0) THEN
     432           16 :          WRITE (unit_nr, *)
     433           16 :          WRITE (unit_nr, *) " Pair potential vdW calculation"
     434           16 :          IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
     435            0 :             WRITE (unit_nr, *) " Dispersion potential type: DFT-D2"
     436            0 :             WRITE (unit_nr, *) " Scaling parameter (s6) ", dispersion_env%scaling
     437            0 :             WRITE (unit_nr, *) " Exponential prefactor  ", dispersion_env%exp_pre
     438           16 :          ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN
     439            4 :             WRITE (unit_nr, *) " Dispersion potential type: DFT-D3"
     440           12 :          ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
     441            4 :             WRITE (unit_nr, *) " Dispersion potential type: DFT-D3(BJ)"
     442            8 :          ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
     443            8 :             WRITE (unit_nr, *) " Dispersion potential type: DFT-D4"
     444              :          END IF
     445              :       END IF
     446              : 
     447         5676 :       CALL get_qs_env(qs_env=qs_env, force=force)
     448         5676 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     449         5676 :       IF (use_virial .AND. debugall) THEN
     450          104 :          dvirial = virial%pv_virial
     451              :       END IF
     452         5676 :       IF (use_virial) THEN
     453         8268 :          pv_loc = virial%pv_virial
     454              :       END IF
     455              : 
     456         5676 :       evdw = 0._dp
     457              :       pv_virial_thread(:, :) = 0._dp
     458              : 
     459         5676 :       CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
     460              : 
     461         5676 :       IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN
     462          116 :          CALL calculate_dispersion_d2_pairpot(qs_env, dispersion_env, evdw, calculate_forces, atevdw)
     463         5618 :       ELSEIF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
     464              :               dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
     465              :          CALL calculate_dispersion_d3_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
     466         9446 :                                               unit_nr, atevdw)
     467          894 :       ELSEIF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
     468          894 :          IF (dispersion_env%lrc) THEN
     469            0 :             CPABORT("Long range correction with DFTD4 not implemented")
     470              :          END IF
     471          894 :          IF (dispersion_env%srb) THEN
     472            0 :             CPABORT("Short range bond correction with DFTD4 not implemented")
     473              :          END IF
     474          894 :          IF (dispersion_env%domol) THEN
     475            0 :             CPABORT("Molecular approximation with DFTD4 not implemented")
     476              :          END IF
     477              :          !
     478          894 :          iw = -1
     479          894 :          IF (dispersion_env%verbose) iw = cp_logger_get_default_io_unit(logger)
     480              :          !
     481          894 :          IF (atenergy .OR. atex) THEN
     482            0 :             ALLOCATE (atomic_energy(natom))
     483              :             CALL calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, &
     484            0 :                                                  iw, atomic_energy=atomic_energy)
     485              :          ELSE
     486          894 :             CALL calculate_dispersion_d4_pairpot(qs_env, dispersion_env, evdw, calculate_forces, iw)
     487              :          END IF
     488              :          !
     489          894 :          IF (atex) THEN
     490            0 :             atevdw(1:natom) = atomic_energy(1:natom)
     491              :          END IF
     492          894 :          IF (atenergy) THEN
     493            0 :             CALL atprop_array_init(atprop%atevdw, natom)
     494            0 :             atprop%atevdw(1:natom) = atomic_energy(1:natom)
     495              :          END IF
     496          894 :          IF (atenergy .OR. atex) THEN
     497            0 :             DEALLOCATE (atomic_energy)
     498              :          END IF
     499              :       END IF
     500              : 
     501              :       ! set dispersion energy
     502         5676 :       CALL para_env%sum(evdw)
     503         5676 :       energy = evdw
     504         5676 :       IF (unit_nr > 0) THEN
     505           16 :          WRITE (unit_nr, *) " Total vdW energy [au]     :", evdw
     506           16 :          WRITE (unit_nr, *) " Total vdW energy [kcal]   :", evdw*kcalmol
     507           16 :          WRITE (unit_nr, *)
     508              :       END IF
     509         5676 :       IF (calculate_forces .AND. debugall) THEN
     510           30 :          IF (unit_nr > 0) THEN
     511            1 :             WRITE (unit_nr, *) " Dispersion Forces         "
     512            1 :             WRITE (unit_nr, *) " Atom   Kind                            Forces    "
     513              :          END IF
     514           30 :          gnorm = 0._dp
     515          166 :          DO iatom = 1, natom
     516          136 :             ikind = kind_of(iatom)
     517          136 :             atom_a = atom_of_kind(iatom)
     518          544 :             fdij(1:3) = force(ikind)%dispersion(:, atom_a)
     519          136 :             CALL para_env%sum(fdij)
     520          544 :             gnorm = gnorm + SUM(ABS(fdij))
     521          166 :             IF (unit_nr > 0) WRITE (unit_nr, "(i5,i7,3F20.14)") iatom, ikind, fdij
     522              :          END DO
     523           30 :          IF (unit_nr > 0) THEN
     524            1 :             WRITE (unit_nr, *)
     525            1 :             WRITE (unit_nr, *) "|G| = ", gnorm
     526            1 :             WRITE (unit_nr, *)
     527              :          END IF
     528           30 :          IF (use_virial) THEN
     529           78 :             dvirial = virial%pv_virial - dvirial
     530            6 :             CALL para_env%sum(dvirial)
     531            6 :             IF (unit_nr > 0) THEN
     532            0 :                WRITE (unit_nr, *) "Stress Tensor (dispersion)"
     533            0 :                WRITE (unit_nr, "(3G20.12)") dvirial
     534            0 :                WRITE (unit_nr, *) "  Tr(P)/3 :  ", (dvirial(1, 1) + dvirial(2, 2) + dvirial(3, 3))/3._dp
     535            0 :                WRITE (unit_nr, *)
     536              :             END IF
     537              :          END IF
     538              :       END IF
     539              : 
     540         5676 :       IF (calculate_forces .AND. use_virial) THEN
     541         3250 :          virial%pv_vdw = virial%pv_vdw + (virial%pv_virial - pv_loc)
     542              :       END IF
     543              : 
     544         5676 :       IF (ASSOCIATED(dispersion_env%dftd_section)) THEN
     545          262 :          CALL cp_print_key_finished_output(unit_nr, logger, dispersion_env%dftd_section, "PRINT_DFTD")
     546              :       END IF
     547              : 
     548         5676 :       CALL timestop(handle)
     549              : 
     550        27278 :    END SUBROUTINE calculate_dispersion_pairpot
     551              : 
     552              : END MODULE qs_dispersion_pairpot
        

Generated by: LCOV version 2.0-1