1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief 8!> 9!> 10!> \par History 11!> refactoring 03-2011 [MI] 12!> \author MI 13! ************************************************************************************************** 14MODULE xc_adiabatic_utils 15 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_para_types, ONLY: cp_para_env_type 18 USE dbcsr_api, ONLY: dbcsr_p_type 19 USE hfx_communication, ONLY: scale_and_add_fock_to_ks_matrix 20 USE hfx_derivatives, ONLY: derivatives_four_center 21 USE input_constants, ONLY: do_adiabatic_hybrid_mcy3,& 22 do_adiabatic_model_pade 23 USE input_section_types, ONLY: section_vals_get,& 24 section_vals_get_subs_vals,& 25 section_vals_type,& 26 section_vals_val_get 27 USE kinds, ONLY: dp 28 USE pw_types, ONLY: pw_p_type 29 USE qs_energy_types, ONLY: qs_energy_type 30 USE qs_environment_types, ONLY: get_qs_env,& 31 qs_environment_type 32 USE qs_ks_types, ONLY: qs_ks_env_type 33 USE qs_rho_types, ONLY: qs_rho_get,& 34 qs_rho_type 35 USE qs_vxc, ONLY: qs_vxc_create 36 USE qs_vxc_atom, ONLY: calculate_vxc_atom 37 USE xc_adiabatic_methods, ONLY: rescale_MCY3_pade 38#include "./base/base_uses.f90" 39 40 IMPLICIT NONE 41 42 PRIVATE 43 44 ! *** Public subroutines *** 45 PUBLIC :: rescale_xc_potential 46 47 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_adiabatic_utils' 48 49CONTAINS 50 51! ************************************************************************************************** 52!> \brief 53!> 54!> \param qs_env ... 55!> \param ks_matrix ... 56!> \param rho ... 57!> \param energy ... 58!> \param v_rspace_new ... 59!> \param v_tau_rspace ... 60!> \param hf_energy ... 61!> \param just_energy ... 62!> \param calculate_forces ... 63!> \param use_virial ... 64! ************************************************************************************************** 65 SUBROUTINE rescale_xc_potential(qs_env, ks_matrix, rho, energy, v_rspace_new, v_tau_rspace, & 66 hf_energy, just_energy, calculate_forces, use_virial) 67 68 TYPE(qs_environment_type), POINTER :: qs_env 69 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix 70 TYPE(qs_rho_type), POINTER :: rho 71 TYPE(qs_energy_type), POINTER :: energy 72 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace 73 REAL(dp), DIMENSION(:) :: hf_energy 74 LOGICAL, INTENT(in) :: just_energy, calculate_forces, use_virial 75 76 CHARACTER(LEN=*), PARAMETER :: routineN = 'rescale_xc_potential', & 77 routineP = moduleN//':'//routineN 78 79 INTEGER :: adiabatic_functional, adiabatic_model, & 80 handle, n_rep_hf, nimages 81 LOGICAL :: do_adiabatic_rescaling, do_hfx, gapw, & 82 gapw_xc 83 REAL(dp) :: adiabatic_lambda, adiabatic_omega, & 84 scale_dDFA, scale_ddW0, scale_dEx1, & 85 scale_dEx2, total_energy_xc 86 TYPE(cp_para_env_type), POINTER :: para_env 87 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao 88 TYPE(dft_control_type), POINTER :: dft_control 89 TYPE(qs_ks_env_type), POINTER :: ks_env 90 TYPE(qs_rho_type), POINTER :: rho_xc 91 TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, & 92 hfx_sections, input, xc_section 93 94 CALL timeset(routineN, handle) 95 NULLIFY (para_env, dft_control, adiabatic_rescaling_section, hfx_sections, & 96 input, xc_section, rho_xc, ks_env, rho_ao) 97 98 CALL get_qs_env(qs_env, & 99 dft_control=dft_control, & 100 para_env=para_env, & 101 input=input, & 102 rho_xc=rho_xc, & 103 ks_env=ks_env) 104 105 nimages = dft_control%nimages 106 CPASSERT(nimages == 1) 107 108 CALL qs_rho_get(rho, rho_ao_kp=rho_ao) 109 110 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING") 111 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling) 112 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF") 113 CALL section_vals_get(hfx_sections, explicit=do_hfx) 114 CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf) 115 116 gapw = dft_control%qs_control%gapw 117 gapw_xc = dft_control%qs_control%gapw_xc 118 119 CALL section_vals_val_get(adiabatic_rescaling_section, "FUNCTIONAL_TYPE", & 120 i_val=adiabatic_functional) 121 CALL section_vals_val_get(adiabatic_rescaling_section, "FUNCTIONAL_MODEL", & 122 i_val=adiabatic_model) 123 CALL section_vals_val_get(adiabatic_rescaling_section, "LAMBDA", & 124 r_val=adiabatic_lambda) 125 CALL section_vals_val_get(adiabatic_rescaling_section, "OMEGA", & 126 r_val=adiabatic_omega) 127 SELECT CASE (adiabatic_functional) 128 CASE (do_adiabatic_hybrid_mcy3) 129 SELECT CASE (adiabatic_model) 130 CASE (do_adiabatic_model_pade) 131 IF (n_rep_hf /= 2) & 132 CALL cp_abort(__LOCATION__, & 133 " For this kind of adiababatic hybrid functional 2 HF sections have to be provided. "// & 134 " Please check your input file.") 135 CALL rescale_MCY3_pade(qs_env, hf_energy, energy, adiabatic_lambda, & 136 adiabatic_omega, scale_dEx1, scale_ddW0, scale_dDFA, & 137 scale_dEx2, total_energy_xc) 138 139 !! Scale and add Fock matrix to KS matrix 140 IF (do_hfx) THEN 141 CALL scale_and_add_fock_to_ks_matrix(para_env, qs_env, ks_matrix, 1, & 142 scale_dEx1) 143 CALL scale_and_add_fock_to_ks_matrix(para_env, qs_env, ks_matrix, 2, & 144 scale_dEx2) 145 END IF 146 IF (calculate_forces) THEN 147 CPASSERT(.NOT. use_virial) 148 !! we also have to scale the forces!!!! 149 CALL derivatives_four_center(qs_env, rho_ao, hfx_sections, para_env, 1, use_virial, & 150 adiabatic_rescale_factor=scale_dEx1) 151 CALL derivatives_four_center(qs_env, rho_ao, hfx_sections, para_env, 2, use_virial, & 152 adiabatic_rescale_factor=scale_dEx2) 153 END IF 154 155 ! Calculate vxc and rescale it 156 xc_section => section_vals_get_subs_vals(input, "DFT%XC") 157 IF (gapw_xc) THEN 158 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_xc, xc_section=xc_section, & 159 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, & 160 just_energy=just_energy, adiabatic_rescale_factor=scale_dDFA) 161 ELSE 162 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, & 163 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, & 164 just_energy=just_energy, adiabatic_rescale_factor=scale_dDFA) 165 END IF 166 167 ! Calculate vxc and rescale it 168 IF (gapw .OR. gapw_xc) THEN 169 CALL calculate_vxc_atom(qs_env, just_energy, adiabatic_rescale_factor=scale_dDFA) 170 END IF 171 !! Hack for the total energy expression 172 energy%ex = 0.0_dp 173 energy%exc1 = 0.0_dp 174 energy%exc = total_energy_xc 175 176 END SELECT 177 END SELECT 178 CALL timestop(handle) 179 180 END SUBROUTINE rescale_xc_potential 181 182END MODULE xc_adiabatic_utils 183 184