1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for calculating local energy and stress tensor 8!> \author JGH 9!> \par History 10!> - 07.2019 created 11! ************************************************************************************************** 12MODULE qs_local_properties 13 USE bibliography, ONLY: Cohen2000,& 14 Filippetti2000,& 15 Rogers2002,& 16 cite_reference 17 USE cp_control_types, ONLY: dft_control_type 18 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set 19 USE cp_log_handling, ONLY: cp_get_default_logger,& 20 cp_logger_get_default_io_unit,& 21 cp_logger_type 22 USE dbcsr_api, ONLY: dbcsr_copy,& 23 dbcsr_p_type,& 24 dbcsr_set,& 25 dbcsr_type 26 USE input_section_types, ONLY: section_vals_get_subs_vals,& 27 section_vals_type 28 USE kinds, ONLY: dp 29 USE mathlib, ONLY: det_3x3 30 USE pw_env_types, ONLY: pw_env_get,& 31 pw_env_type 32 USE pw_methods, ONLY: pw_axpy,& 33 pw_copy,& 34 pw_derive,& 35 pw_integrate_function,& 36 pw_multiply,& 37 pw_transfer,& 38 pw_zero 39 USE pw_pool_types, ONLY: pw_pool_create_pw,& 40 pw_pool_give_back_pw,& 41 pw_pool_type 42 USE pw_types, ONLY: COMPLEXDATA1D,& 43 REALDATA3D,& 44 REALSPACE,& 45 RECIPROCALSPACE,& 46 pw_p_type 47 USE qs_collocate_density, ONLY: calculate_rho_elec 48 USE qs_core_energies, ONLY: calculate_ptrace 49 USE qs_energy_matrix_w, ONLY: qs_energies_compute_w 50 USE qs_energy_types, ONLY: qs_energy_type 51 USE qs_environment_types, ONLY: get_qs_env,& 52 qs_environment_type 53 USE qs_ks_methods, ONLY: calc_rho_tot_gspace 54 USE qs_ks_types, ONLY: qs_ks_env_type,& 55 set_ks_env 56 USE qs_rho_types, ONLY: qs_rho_get,& 57 qs_rho_type 58 USE qs_vxc, ONLY: qs_xc_density 59 USE virial_types, ONLY: virial_type 60#include "./base/base_uses.f90" 61 62 IMPLICIT NONE 63 64 PRIVATE 65 66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_local_properties' 67 68 PUBLIC :: qs_local_energy, qs_local_stress 69 70! ************************************************************************************************** 71 72CONTAINS 73 74! ************************************************************************************************** 75!> \brief Routine to calculate the local energy 76!> \param qs_env the qs_env to update 77!> \param energy_density ... 78!> \par History 79!> 07.2019 created 80!> \author JGH 81! ************************************************************************************************** 82 SUBROUTINE qs_local_energy(qs_env, energy_density) 83 TYPE(qs_environment_type), POINTER :: qs_env 84 TYPE(pw_p_type), INTENT(INOUT) :: energy_density 85 86 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_local_energy', & 87 routineP = moduleN//':'//routineN 88 89 INTEGER :: handle, img, iounit, ispin, nimages, & 90 nkind, nspins 91 LOGICAL :: gapw, gapw_xc 92 REAL(KIND=dp) :: eban, eband, eh, exc, ovol 93 REAL(KIND=dp), DIMENSION(2) :: band_energy 94 TYPE(cp_logger_type), POINTER :: logger 95 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao 96 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s, matrix_w, rho_ao_kp 97 TYPE(dbcsr_type), POINTER :: matrix 98 TYPE(dft_control_type), POINTER :: dft_control 99 TYPE(pw_env_type), POINTER :: pw_env 100 TYPE(pw_p_type) :: band_density, edens_g, edens_r, & 101 hartree_density, rho_tot_gspace, & 102 rho_tot_rspace, v_hartree_rspace, & 103 xc_density 104 TYPE(pw_p_type), DIMENSION(:), POINTER :: rho_r 105 TYPE(pw_p_type), POINTER :: rho_core 106 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 107 TYPE(qs_energy_type), POINTER :: energy 108 TYPE(qs_ks_env_type), POINTER :: ks_env 109 TYPE(qs_rho_type), POINTER :: rho, rho_struct 110 TYPE(section_vals_type), POINTER :: input, xc_section 111 112 CALL timeset(routineN, handle) 113 114 CALL cite_reference(Cohen2000) 115 116 CPASSERT(ASSOCIATED(qs_env)) 117 logger => cp_get_default_logger() 118 iounit = cp_logger_get_default_io_unit() 119 120 ! Check for GAPW method : additional terms for local densities 121 CALL get_qs_env(qs_env, nkind=nkind, dft_control=dft_control) 122 gapw = dft_control%qs_control%gapw 123 gapw_xc = dft_control%qs_control%gapw_xc 124 125 nimages = dft_control%nimages 126 nspins = dft_control%nspins 127 128 ! get working arrays 129 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 130 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 131 CALL pw_pool_create_pw(auxbas_pw_pool, band_density%pw, use_data=REALDATA3D, in_space=REALSPACE) 132 CALL pw_pool_create_pw(auxbas_pw_pool, hartree_density%pw, use_data=REALDATA3D, in_space=REALSPACE) 133 CALL pw_pool_create_pw(auxbas_pw_pool, xc_density%pw, use_data=REALDATA3D, in_space=REALSPACE) 134 135 ! w matrix 136 CALL get_qs_env(qs_env, matrix_w_kp=matrix_w) 137 IF (.NOT. ASSOCIATED(matrix_w)) THEN 138 CALL get_qs_env(qs_env, & 139 ks_env=ks_env, & 140 matrix_s_kp=matrix_s) 141 matrix => matrix_s(1, 1)%matrix 142 CALL dbcsr_allocate_matrix_set(matrix_w, nspins, nimages) 143 DO ispin = 1, nspins 144 DO img = 1, nimages 145 ALLOCATE (matrix_w(ispin, img)%matrix) 146 CALL dbcsr_copy(matrix_w(ispin, img)%matrix, matrix, name="W MATRIX") 147 CALL dbcsr_set(matrix_w(ispin, img)%matrix, 0.0_dp) 148 END DO 149 END DO 150 CALL set_ks_env(ks_env, matrix_w_kp=matrix_w) 151 END IF 152 ! band structure energy density 153 CALL qs_energies_compute_w(qs_env, .TRUE.) 154 band_energy = 0.0_dp 155 CALL get_qs_env(qs_env, ks_env=ks_env, matrix_w_kp=matrix_w) 156 CALL pw_pool_create_pw(auxbas_pw_pool, & 157 edens_r%pw, & 158 use_data=REALDATA3D, & 159 in_space=REALSPACE) 160 CALL pw_pool_create_pw(auxbas_pw_pool, & 161 edens_g%pw, & 162 use_data=COMPLEXDATA1D, & 163 in_space=RECIPROCALSPACE) 164 CALL pw_zero(band_density%pw) 165 DO ispin = 1, nspins 166 rho_ao => matrix_w(ispin, :) 167 CALL calculate_rho_elec(matrix_p_kp=rho_ao, & 168 rho=edens_r, & 169 rho_gspace=edens_g, & 170 total_rho=band_energy(ispin), & 171 ks_env=ks_env, soft_valid=(gapw .OR. gapw_xc)) 172 CALL pw_axpy(edens_r%pw, band_density%pw) 173 END DO 174 CALL pw_pool_give_back_pw(auxbas_pw_pool, edens_r%pw) 175 CALL pw_pool_give_back_pw(auxbas_pw_pool, edens_g%pw) 176 177 ! Hartree energy density correction = -0.5 * V_H(r) * [rho(r) - rho_core(r)] 178 CALL pw_pool_create_pw(auxbas_pw_pool, & 179 rho_tot_gspace%pw, & 180 use_data=COMPLEXDATA1D, & 181 in_space=RECIPROCALSPACE) 182 CALL pw_pool_create_pw(auxbas_pw_pool, & 183 rho_tot_rspace%pw, & 184 use_data=REALDATA3D, & 185 in_space=REALSPACE) 186 CALL get_qs_env(qs_env, & 187 v_hartree_rspace=v_hartree_rspace%pw, & 188 rho_core=rho_core, rho=rho) 189 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho) 190 CALL qs_rho_get(rho, rho_r=rho_r) 191 CALL pw_transfer(rho_core%pw, rho_tot_rspace%pw) 192 DO ispin = 1, nspins 193 CALL pw_axpy(rho_r(ispin)%pw, rho_tot_rspace%pw, alpha=-1.0_dp) 194 END DO 195 CALL pw_zero(hartree_density%pw) 196 ovol = 0.5_dp/hartree_density%pw%pw_grid%dvol 197 CALL pw_multiply(hartree_density%pw, v_hartree_rspace%pw, rho_tot_rspace%pw, alpha=ovol) 198 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw) 199 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_rspace%pw) 200 201 IF (dft_control%do_admm) THEN 202 CALL cp_warn(__LOCATION__, "ADMM not supported for local energy calculation") 203 END IF 204 IF (gapw_xc .OR. gapw) THEN 205 CALL cp_warn(__LOCATION__, "GAPW/GAPW_XC not supported for local energy calculation") 206 END IF 207 ! XC energy density correction = E_xc(r) - V_xc(r)*rho(r) 208 CALL get_qs_env(qs_env, input=input) 209 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 210 CALL get_qs_env(qs_env=qs_env, rho=rho_struct) 211 ! 212 CALL qs_xc_density(ks_env, rho_struct, xc_section, xc_density) 213 ! 214 ! energies 215 CALL get_qs_env(qs_env=qs_env, energy=energy) 216 eban = pw_integrate_function(band_density%pw) 217 eh = pw_integrate_function(hartree_density%pw) 218 exc = pw_integrate_function(xc_density%pw) 219 220 ! band energy 221 CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks) 222 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp) 223 CALL calculate_ptrace(matrix_ks, rho_ao_kp, eband, nspins) 224 225 ! get full density 226 CALL pw_copy(band_density%pw, energy_density%pw) 227 CALL pw_axpy(hartree_density%pw, energy_density%pw) 228 CALL pw_axpy(xc_density%pw, energy_density%pw) 229 230 IF (iounit > 0) THEN 231 WRITE (UNIT=iounit, FMT="(/,T3,A)") REPEAT("=", 78) 232 WRITE (UNIT=iounit, FMT="(T4,A,T52,A,T75,A)") "Local Energy Calculation", "GPW/GAPW", "Local" 233 WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8,T65,F15.8)") "Band Energy", eband, eban 234 WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "Hartree Energy Correction", eh 235 WRITE (UNIT=iounit, FMT="(T4,A,T65,F15.8)") "XC Energy Correction", exc 236 WRITE (UNIT=iounit, FMT="(T4,A,T45,F15.8,T65,F15.8)") "Total Energy", & 237 energy%total, eban + eh + exc + energy%core_overlap + energy%core_self + energy%dispersion 238 WRITE (UNIT=iounit, FMT="(T3,A)") REPEAT("=", 78) 239 END IF 240 241 ! return temp arrays 242 CALL pw_pool_give_back_pw(auxbas_pw_pool, band_density%pw) 243 CALL pw_pool_give_back_pw(auxbas_pw_pool, hartree_density%pw) 244 CALL pw_pool_give_back_pw(auxbas_pw_pool, xc_density%pw) 245 246 CALL timestop(handle) 247 248 END SUBROUTINE qs_local_energy 249 250! ************************************************************************************************** 251!> \brief Routine to calculate the local stress 252!> \param qs_env the qs_env to update 253!> \param stress_tensor ... 254!> \param pressure ... 255!> \param beta ... 256!> \par History 257!> 07.2019 created 258!> \author JGH 259! ************************************************************************************************** 260 SUBROUTINE qs_local_stress(qs_env, stress_tensor, pressure, beta) 261 TYPE(qs_environment_type), POINTER :: qs_env 262 TYPE(pw_p_type), DIMENSION(:, :), INTENT(INOUT), & 263 OPTIONAL :: stress_tensor 264 TYPE(pw_p_type), INTENT(INOUT), OPTIONAL :: pressure 265 REAL(KIND=dp), INTENT(IN), OPTIONAL :: beta 266 267 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_local_stress', & 268 routineP = moduleN//':'//routineN 269 270 INTEGER :: handle, i, iounit, j, nimages, nkind, & 271 nspins 272 LOGICAL :: do_pressure, do_stress, gapw, gapw_xc, & 273 use_virial 274 REAL(KIND=dp) :: my_beta 275 REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc 276 TYPE(cp_logger_type), POINTER :: logger 277 TYPE(dft_control_type), POINTER :: dft_control 278 TYPE(pw_env_type), POINTER :: pw_env 279 TYPE(pw_p_type) :: v_hartree_gspace, v_hartree_rspace, & 280 xc_density 281 TYPE(pw_p_type), DIMENSION(3) :: efield 282 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 283 TYPE(qs_ks_env_type), POINTER :: ks_env 284 TYPE(qs_rho_type), POINTER :: rho_struct 285 TYPE(section_vals_type), POINTER :: input, xc_section 286 TYPE(virial_type), POINTER :: virial 287 288 CALL cp_warn(__LOCATION__, "Local Stress Tensor code is not working, skipping") 289 RETURN 290 291 CALL timeset(routineN, handle) 292 293 CALL cite_reference(Filippetti2000) 294 CALL cite_reference(Rogers2002) 295 296 CPASSERT(ASSOCIATED(qs_env)) 297 298 IF (PRESENT(stress_tensor)) THEN 299 do_stress = .TRUE. 300 ELSE 301 do_stress = .FALSE. 302 END IF 303 IF (PRESENT(pressure)) THEN 304 do_pressure = .TRUE. 305 ELSE 306 do_pressure = .FALSE. 307 END IF 308 IF (PRESENT(beta)) THEN 309 my_beta = beta 310 ELSE 311 my_beta = 0.0_dp 312 END IF 313 314 logger => cp_get_default_logger() 315 iounit = cp_logger_get_default_io_unit() 316 317 !!!!!! 318 CALL cp_warn(__LOCATION__, "Local Stress Tensor code is not tested") 319 !!!!!! 320 321 ! Check for GAPW method : additional terms for local densities 322 CALL get_qs_env(qs_env, nkind=nkind, dft_control=dft_control) 323 gapw = dft_control%qs_control%gapw 324 gapw_xc = dft_control%qs_control%gapw_xc 325 326 nimages = dft_control%nimages 327 nspins = dft_control%nspins 328 329 ! get working arrays 330 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env) 331 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 332 CALL pw_pool_create_pw(auxbas_pw_pool, xc_density%pw, use_data=REALDATA3D, in_space=REALSPACE) 333 334 ! init local stress tensor 335 IF (do_stress) THEN 336 DO i = 1, 3 337 DO j = 1, 3 338 CALL pw_zero(stress_tensor(i, j)%pw) 339 END DO 340 END DO 341 END IF 342 343 IF (dft_control%do_admm) THEN 344 CALL cp_warn(__LOCATION__, "ADMM not supported for local energy calculation") 345 END IF 346 IF (gapw_xc .OR. gapw) THEN 347 CALL cp_warn(__LOCATION__, "GAPW/GAPW_XC not supported for local energy calculation") 348 END IF 349 ! XC energy density correction = E_xc(r) - V_xc(r)*rho(r) 350 CALL get_qs_env(qs_env, ks_env=ks_env, input=input, rho=rho_struct) 351 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 352 ! 353 CALL qs_xc_density(ks_env, rho_struct, xc_section, xc_density) 354 355 ! Electrical field terms 356 CALL get_qs_env(qs_env, v_hartree_rspace=v_hartree_rspace%pw) 357 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 358 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 359 CALL pw_transfer(v_hartree_rspace%pw, v_hartree_gspace%pw) 360 DO i = 1, 3 361 NULLIFY (efield(i)%pw) 362 CALL pw_pool_create_pw(auxbas_pw_pool, efield(i)%pw, & 363 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 364 CALL pw_copy(v_hartree_gspace%pw, efield(i)%pw) 365 END DO 366 CALL pw_derive(efield(1)%pw, (/1, 0, 0/)) 367 CALL pw_derive(efield(2)%pw, (/0, 1, 0/)) 368 CALL pw_derive(efield(3)%pw, (/0, 0, 1/)) 369 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw) 370 DO i = 1, 3 371 CALL pw_pool_give_back_pw(auxbas_pw_pool, efield(i)%pw) 372 END DO 373 374 pv_loc = 0.0_dp 375 376 CALL get_qs_env(qs_env=qs_env, virial=virial) 377 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 378 IF (.NOT. use_virial) THEN 379 CALL cp_warn(__LOCATION__, "Local stress should be used with standard stress calculation.") 380 END IF 381 IF (iounit > 0 .AND. use_virial) THEN 382 WRITE (UNIT=iounit, FMT="(/,T3,A)") REPEAT("=", 78) 383 WRITE (UNIT=iounit, FMT="(T4,A)") "Local Stress Calculation" 384 WRITE (UNIT=iounit, FMT="(T42,A,T64,A)") " 1/3 Trace", " Determinant" 385 WRITE (UNIT=iounit, FMT="(T4,A,T42,F16.8,T64,F16.8)") "Total Stress", & 386 (pv_loc(1, 1) + pv_loc(2, 2) + pv_loc(3, 3))/3.0_dp, det_3x3(pv_loc) 387 WRITE (UNIT=iounit, FMT="(T3,A)") REPEAT("=", 78) 388 END IF 389 390 CALL timestop(handle) 391 392 END SUBROUTINE qs_local_stress 393 394! ************************************************************************************************** 395 396END MODULE qs_local_properties 397