1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Set disperson types for DFT calculations 8!> \author JGH (04.2014) 9! ************************************************************************************************** 10MODULE qs_dispersion_utils 11 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind 14 USE cp_log_handling, ONLY: cp_get_default_logger,& 15 cp_logger_type 16 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 17 cp_print_key_unit_nr 18 USE input_constants, ONLY: vdw_nl_DRSLL,& 19 vdw_nl_LMKLL,& 20 vdw_nl_RVV10,& 21 vdw_pairpot_dftd2,& 22 vdw_pairpot_dftd3,& 23 vdw_pairpot_dftd3bj,& 24 xc_vdw_fun_nonloc,& 25 xc_vdw_fun_pairpot 26 USE input_section_types, ONLY: section_vals_get_subs_vals,& 27 section_vals_type,& 28 section_vals_val_get 29 USE kinds, ONLY: dp 30 USE physcon, ONLY: bohr,& 31 kjmol 32 USE qs_dispersion_pairpot, ONLY: qs_scaling_dftd3,& 33 qs_scaling_dftd3bj,& 34 qs_scaling_init 35 USE qs_dispersion_types, ONLY: qs_atom_dispersion_type,& 36 qs_dispersion_type 37 USE qs_environment_types, ONLY: get_qs_env,& 38 qs_environment_type 39 USE qs_kind_types, ONLY: get_qs_kind,& 40 qs_kind_type 41#include "./base/base_uses.f90" 42 43 IMPLICIT NONE 44 45 PRIVATE 46 47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dispersion_utils' 48 49 PUBLIC :: qs_dispersion_env_set, qs_write_dispersion 50 51! ************************************************************************************************** 52CONTAINS 53! ************************************************************************************************** 54!> \brief ... 55!> \param dispersion_env ... 56!> \param xc_section ... 57! ************************************************************************************************** 58 SUBROUTINE qs_dispersion_env_set(dispersion_env, xc_section) 59 TYPE(qs_dispersion_type), POINTER :: dispersion_env 60 TYPE(section_vals_type), POINTER :: xc_section 61 62 CHARACTER(len=*), PARAMETER :: routineN = 'qs_dispersion_env_set', & 63 routineP = moduleN//':'//routineN 64 65 LOGICAL :: exfun, explicit 66 REAL(dp), POINTER :: params(:), scal(:) 67 TYPE(section_vals_type), POINTER :: nl_section, pp_section, vdw_section, & 68 xc_fun_section 69 70 CPASSERT(ASSOCIATED(dispersion_env)) 71 72 ! set general defaults 73 dispersion_env%doabc = .FALSE. 74 dispersion_env%c9cnst = .FALSE. 75 dispersion_env%lrc = .FALSE. 76 dispersion_env%srb = .FALSE. 77 dispersion_env%verbose = .FALSE. 78 NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, & 79 dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist) 80 NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, & 81 dispersion_env%d2y_dx2) 82 NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn) 83 NULLIFY (dispersion_env%dftd_section) 84 NULLIFY (vdw_section, xc_fun_section) 85 vdw_section => section_vals_get_subs_vals(xc_section, "vdw_potential") 86 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL") 87 CALL section_vals_val_get(vdw_section, "POTENTIAL_TYPE", i_val=dispersion_env%type) 88 IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN 89 NULLIFY (pp_section) 90 pp_section => section_vals_get_subs_vals(vdw_section, "PAIR_POTENTIAL") 91 CALL section_vals_val_get(pp_section, "VERBOSE_OUTPUT", l_val=dispersion_env%verbose) 92 CALL section_vals_val_get(pp_section, "TYPE", i_val=dispersion_env%pp_type) 93 IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN 94 ! functional parameters for Grimme D2 type 95 CALL section_vals_val_get(pp_section, "EXP_PRE", r_val=dispersion_env%exp_pre) 96 CALL section_vals_val_get(pp_section, "SCALING", explicit=explicit) 97 IF (.NOT. explicit) THEN 98 CALL section_vals_val_get(pp_section, "REFERENCE_FUNCTIONAL", explicit=exfun) 99 CPASSERT(exfun) 100 CALL qs_scaling_init(dispersion_env%scaling, vdw_section) 101 ELSE 102 CALL section_vals_val_get(pp_section, "SCALING", r_val=dispersion_env%scaling) 103 END IF 104 ELSE 105 dispersion_env%exp_pre = 0._dp 106 dispersion_env%scaling = 0._dp 107 END IF 108 IF (dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. & 109 dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN 110 ! functional parameters for Grimme DFT-D3 type 111 CALL section_vals_val_get(pp_section, "EPS_CN", r_val=dispersion_env%eps_cn) 112 CALL section_vals_val_get(pp_section, "CALCULATE_C9_TERM", l_val=dispersion_env%doabc) 113 CALL section_vals_val_get(pp_section, "REFERENCE_C9_TERM", l_val=dispersion_env%c9cnst) 114 CALL section_vals_val_get(pp_section, "LONG_RANGE_CORRECTION", l_val=dispersion_env%lrc) 115 CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION", l_val=dispersion_env%srb) 116 CALL section_vals_val_get(pp_section, "SHORT_RANGE_CORRECTION_PARAMETERS", r_vals=params) 117 dispersion_env%srb_params(1:4) = params(1:4) 118 ! KG corrections 119 CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION", l_val=dispersion_env%domol) 120 CALL section_vals_val_get(pp_section, "MOLECULE_CORRECTION_C8", r_val=dispersion_env%kgc8) 121 IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN 122 CALL section_vals_val_get(pp_section, "D3_SCALING", explicit=explicit) 123 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN 124 CALL section_vals_val_get(pp_section, "D3BJ_SCALING", explicit=explicit) 125 END IF 126 IF (.NOT. explicit) THEN 127 CALL section_vals_val_get(pp_section, "REFERENCE_FUNCTIONAL", explicit=exfun) 128 CPASSERT(exfun) 129 IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN 130 CALL qs_scaling_dftd3(dispersion_env%s6, dispersion_env%sr6, dispersion_env%s8, vdw_section) 131 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN 132 CALL qs_scaling_dftd3bj(dispersion_env%s6, dispersion_env%a1, dispersion_env%s8, & 133 dispersion_env%a2, vdw_section) 134 END IF 135 ELSE 136 IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN 137 ! zero damping 138 CALL section_vals_val_get(pp_section, "D3_SCALING", r_vals=scal) 139 dispersion_env%s6 = scal(1) 140 dispersion_env%sr6 = scal(2) 141 dispersion_env%s8 = scal(3) 142 dispersion_env%a1 = 0.0_dp 143 dispersion_env%a2 = 0.0_dp 144 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN 145 ! BJ damping 146 CALL section_vals_val_get(pp_section, "D3BJ_SCALING", r_vals=scal) 147 dispersion_env%s6 = scal(1) 148 dispersion_env%a1 = scal(2) 149 dispersion_env%s8 = scal(3) 150 dispersion_env%a2 = scal(4) 151 dispersion_env%sr6 = 0.0_dp 152 END IF 153 END IF 154 ELSE 155 dispersion_env%s6 = 0._dp 156 dispersion_env%sr6 = 0._dp 157 dispersion_env%s8 = 0._dp 158 dispersion_env%a1 = 0._dp 159 dispersion_env%a2 = 0._dp 160 dispersion_env%eps_cn = 0._dp 161 END IF 162 CALL section_vals_val_get(pp_section, "R_CUTOFF", r_val=dispersion_env%rc_disp) 163 CALL section_vals_val_get(pp_section, "PARAMETER_FILE_NAME", & 164 c_val=dispersion_env%parameter_file_name) 165 ! set DFTD section for output handling 166 dispersion_env%dftd_section => pp_section 167 ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN 168 NULLIFY (nl_section) 169 nl_section => section_vals_get_subs_vals(vdw_section, "NON_LOCAL") 170 CALL section_vals_val_get(nl_section, "VERBOSE_OUTPUT", l_val=dispersion_env%verbose) 171 CALL section_vals_val_get(nl_section, "KERNEL_FILE_NAME", & 172 c_val=dispersion_env%kernel_file_name) 173 CALL section_vals_val_get(nl_section, "TYPE", i_val=dispersion_env%nl_type) 174 CALL section_vals_val_get(nl_section, "CUTOFF", r_val=dispersion_env%pw_cutoff) 175 CALL section_vals_val_get(nl_section, "PARAMETERS", r_vals=params) 176 dispersion_env%b_value = params(1) 177 dispersion_env%c_value = params(2) 178 END IF 179 END SUBROUTINE qs_dispersion_env_set 180 181! ************************************************************************************************** 182 183! ************************************************************************************************** 184!> \brief ... 185!> \param qs_env ... 186!> \param dispersion_env ... 187!> \param ounit ... 188! ************************************************************************************************** 189 SUBROUTINE qs_write_dispersion(qs_env, dispersion_env, ounit) 190 TYPE(qs_environment_type), POINTER :: qs_env 191 TYPE(qs_dispersion_type), POINTER :: dispersion_env 192 INTEGER, INTENT(in), OPTIONAL :: ounit 193 194 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_write_dispersion', & 195 routineP = moduleN//':'//routineN 196 197 CHARACTER(LEN=2) :: symbol 198 INTEGER :: i, ikind, nkind, output_unit 199 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 200 TYPE(cp_logger_type), POINTER :: logger 201 TYPE(qs_atom_dispersion_type), POINTER :: disp 202 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 203 TYPE(section_vals_type), POINTER :: dft_section 204 205 IF (PRESENT(ounit)) THEN 206 output_unit = ounit 207 ELSE 208 NULLIFY (logger) 209 logger => cp_get_default_logger() 210 211 dft_section => section_vals_get_subs_vals(qs_env%input, "DFT") 212 output_unit = cp_print_key_unit_nr(logger, dft_section, & 213 "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log") 214 END IF 215 216 IF (output_unit > 0) THEN 217 ! vdW type specific output 218 IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN 219 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T67,'Pair Potential')") 220 ! Pair potentials 221 IF (dispersion_env%pp_type == vdw_pairpot_dftd2) THEN 222 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Potential Form: S. Grimme, JCC 27: 1787 (2006)')") 223 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Cutoff Radius [Bohr]:',T73,F8.2)") dispersion_env%rc_disp 224 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Scaling Factor:',T73,F8.4)") dispersion_env%scaling 225 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T35,'Exp Prefactor for Damping:',T73,F8.1)") dispersion_env%exp_pre 226 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set) 227 nkind = SIZE(atomic_kind_set) 228 DO ikind = 1, nkind 229 CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=symbol) 230 CALL get_qs_kind(qs_kind_set(ikind), dispersion=disp) 231 IF (disp%defined) THEN 232 WRITE (output_unit, fmt="(' vdW PARAMETER| ',T18,'Atom=',A2,"// & 233 "T28,'C6[J*nm^6*mol^-1]=',F8.4,T63,'r(vdW)[A]=',F8.4)") & 234 symbol, disp%c6/(1000._dp*bohr**6/kjmol), disp%vdw_radii/bohr 235 ELSE 236 WRITE (output_unit, fmt="(' vdW PARAMETER| ',T20,'Atom=',A2,T70,'not defined')") 237 END IF 238 END DO 239 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3) THEN 240 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'DFT-D3 (Version 3.1)')") 241 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Potential Form: S. Grimme et al, JCP 132: 154104 (2010)')") 242 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Zero Damping')") 243 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff Radius [Bohr]:',T73,F8.2)") dispersion_env%rc_disp 244 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s6 Scaling Factor:',T73,F8.4)") dispersion_env%s6 245 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'sr6 Scaling Factor:',T73,F8.4)") dispersion_env%sr6 246 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s8 Scaling Factor:',T73,F8.4)") dispersion_env%s8 247 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff for CN calculation:',T69,E12.4)") dispersion_env%eps_cn 248 ELSE IF (dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN 249 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'DFT-D3 (Version 3.1)')") 250 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Potential Form: S. Grimme et al, JCP 132: 154104 (2010)')") 251 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'BJ Damping: S. Grimme et al, JCC 32: 1456 (2011)')") 252 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff Radius [Bohr]:',T73,F8.2)") dispersion_env%rc_disp 253 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s6 Scaling Factor:',T73,F8.4)") dispersion_env%s6 254 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'a1 Damping Factor:',T73,F8.4)") dispersion_env%a1 255 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'s8 Scaling Factor:',T73,F8.4)") dispersion_env%s8 256 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'a2 Damping Factor:',T73,F8.4)") dispersion_env%a2 257 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T26,'Cutoff for CN calculation:',T69,E12.4)") dispersion_env%eps_cn 258 END IF 259 ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN 260 WRITE (output_unit, fmt="(' vdW POTENTIAL| ',T61,'Non-local Functional')") 261 WRITE (output_unit, & 262 fmt="(' vdW POTENTIAL| ','Implementation: G. Roman-Perez, J. Soler, PRL 103: 096102 (2009)')") 263 WRITE (output_unit, & 264 fmt="(' vdW POTENTIAL| ',T38,' T. Thonhauser et al, PRB 76: 125112 (2007)')") 265 WRITE (output_unit, & 266 fmt="(' vdW POTENTIAL| ',T22,' R. Sabatini et al, J.Phys:Condens Matter 24: 424209 (2012)')") 267 WRITE (output_unit, & 268 fmt="(' vdW POTENTIAL| ',T16,' Based on QE implementation by Brian Kolb, Timo Thonhauser (2009)')") 269 SELECT CASE (dispersion_env%nl_type) 270 CASE DEFAULT 271 ! unknown functional 272 CPABORT("") 273 CASE (vdw_nl_DRSLL) 274 WRITE (output_unit, & 275 fmt="(' vdW POTENTIAL| ','DRSLL Functional: M. Dion et al, PRL 92: 246401 (2004)')") 276 CASE (vdw_nl_LMKLL) 277 WRITE (output_unit, & 278 fmt="(' vdW POTENTIAL| ','LMKLL Functional: K. Lee et al, PRB 82: 081101 (2010)')") 279 CASE (vdw_nl_RVV10) 280 WRITE (output_unit, & 281 fmt="(' vdW POTENTIAL| ','RVV10 Functional: R. Sabatini et al, PRB 87: 041108(R) (2013)')") 282 END SELECT 283 IF (dispersion_env%verbose) THEN 284 WRITE (output_unit, & 285 fmt="(' vdW POTENTIAL| ',' Carrying out vdW-DF run using the following parameters:')") 286 WRITE (output_unit, fmt="(' vdW POTENTIAL| ','Nqs =',I8,' Nr_points =',I8,' r_max =',F10.3)") & 287 dispersion_env%nqs, dispersion_env%nr_points, dispersion_env%r_max 288 WRITE (output_unit, fmt="(' vdW POTENTIAL| ','q_mesh =')") 289 WRITE (output_unit, fmt="(8X,4F18.8)") (dispersion_env%q_mesh(i), i=1, dispersion_env%nqs) 290 WRITE (output_unit, & 291 fmt="(' vdW POTENTIAL| ','Density cutoff for convolution [a.u.]:',T71,F10.1)") & 292 dispersion_env%pw_cutoff 293 END IF 294 END IF 295 END IF 296 IF (.NOT. PRESENT(ounit)) THEN 297 CALL cp_print_key_finished_output(output_unit, logger, dft_section, & 298 "PRINT%DFT_CONTROL_PARAMETERS") 299 END IF 300 301 END SUBROUTINE qs_write_dispersion 302 303! ************************************************************************************************** 304 305END MODULE qs_dispersion_utils 306 307