|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE harris_functional 00013 USE cp_control_types, ONLY: dft_control_type 00014 USE cp_dbcsr_interface, ONLY: cp_dbcsr_trace 00015 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00016 USE cp_para_types, ONLY: cp_para_env_type 00017 USE f77_blas 00018 USE harris_energy_types, ONLY: harris_energy_type 00019 USE harris_env_types, ONLY: harris_env_get,& 00020 harris_env_type 00021 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00022 section_vals_type 00023 USE kinds, ONLY: dp 00024 USE pw_env_types, ONLY: pw_env_get,& 00025 pw_env_type 00026 USE pw_methods, ONLY: pw_integral_ab 00027 USE pw_poisson_methods, ONLY: pw_poisson_solve 00028 USE pw_poisson_types, ONLY: pw_poisson_type 00029 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00030 pw_pool_give_back_pw,& 00031 pw_pool_type 00032 USE pw_types, ONLY: COMPLEXDATA1D,& 00033 RECIPROCALSPACE,& 00034 pw_p_type 00035 USE qs_core_energies, ONLY: calculate_ecore_overlap,& 00036 calculate_ecore_self 00037 USE qs_environment_types, ONLY: get_qs_env,& 00038 qs_environment_type 00039 USE qs_mo_types, ONLY: get_mo_set,& 00040 mo_set_p_type 00041 USE qs_rho_methods, ONLY: duplicate_rho_type 00042 USE qs_rho_types, ONLY: qs_rho_get,& 00043 qs_rho_type 00044 USE qs_scf_types, ONLY: qs_scf_env_type 00045 USE qs_vxc, ONLY: qs_vxc_create 00046 USE timings, ONLY: timeset,& 00047 timestop 00048 #include "cp_common_uses.h" 00049 00050 IMPLICIT NONE 00051 PRIVATE 00052 00053 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'harris_functional' 00054 00055 ! *** Public subroutines *** 00056 PUBLIC :: harris_energy_correction, & 00057 harris_eigenvalue_summation, & 00058 harris_eigenvalue_calculation, & 00059 harris_eigenvalue_trace_KS_Pmix, & 00060 harris_postprocessing 00061 00062 !*** 00063 00064 CONTAINS 00065 00066 ! ***************************************************************************** 00083 SUBROUTINE harris_energy_correction(qs_env, harris_env, para_env, EII_necessary, fast, error) 00084 00085 TYPE(qs_environment_type), POINTER :: qs_env 00086 TYPE(harris_env_type), POINTER :: harris_env 00087 TYPE(cp_para_env_type), POINTER :: para_env 00088 LOGICAL, INTENT(IN), OPTIONAL :: EII_necessary, fast 00089 TYPE(cp_error_type), INTENT(INOUT) :: error 00090 00091 CHARACTER(len=*), PARAMETER :: routineN = 'harris_energy_correction', 00092 routineP = moduleN//':'//routineN 00093 00094 INTEGER :: handle, ispin, nspins, stat 00095 LOGICAL :: failure, fast_flag, 00096 my_EII_necessary 00097 TYPE(dft_control_type), POINTER :: dft_control 00098 TYPE(harris_energy_type), POINTER :: harris_energy 00099 TYPE(pw_env_type), POINTER :: pw_env 00100 TYPE(pw_p_type) :: v_hartree_gspace 00101 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace 00102 TYPE(pw_p_type), POINTER :: rho_core 00103 TYPE(pw_poisson_type), POINTER :: poisson_env 00104 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00105 TYPE(qs_rho_type), POINTER :: rho 00106 TYPE(section_vals_type), POINTER :: input, xc_section 00107 00108 CALL timeset(routineN,handle) 00109 00110 failure = .FALSE. 00111 NULLIFY(rho, pw_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, & 00112 rho_core, rho, harris_energy, poisson_env, dft_control, input, xc_section) 00113 00114 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 00115 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 00116 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00117 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00118 00119 my_EII_necessary = .TRUE. 00120 IF (PRESENT(EII_necessary)) my_EII_necessary=EII_necessary 00121 00122 IF (PRESENT(fast)) THEN 00123 fast_flag = fast 00124 ELSE 00125 fast_flag = .FALSE. 00126 END IF 00127 00128 IF (.NOT. failure) THEN 00129 CALL harris_env_get(harris_env=harris_env, harris_energy=harris_energy, error=error) 00130 CALL get_qs_env(qs_env=qs_env, input=input, pw_env=pw_env, rho=rho, rho_core=rho_core, & 00131 dft_control=dft_control, error=error) 00132 00133 xc_section => section_vals_get_subs_vals(input, "DFT%XC", error=error) 00134 00135 nspins = dft_control%nspins 00136 00137 CALL duplicate_rho_type(rho_input=rho, rho_output=harris_env%rho, & 00138 qs_env=qs_env, error=error) 00139 00140 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env, & 00141 error=error) 00142 00143 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 00144 use_data=COMPLEXDATA1D, & 00145 in_space=RECIPROCALSPACE, error=error) 00146 00147 DO ispin = 1,nspins 00148 CALL pw_poisson_solve(poisson_env, rho%rho_g(ispin)%pw, & 00149 ehartree=harris_energy%Ehartree_elec, & 00150 vhartree=v_hartree_gspace%pw,error=error) 00151 END DO 00152 00153 IF (my_EII_necessary) THEN 00154 CALL pw_poisson_solve(poisson_env, rho_core%pw, ehartree=harris_energy%Ehartree_core, & 00155 vhartree=v_hartree_gspace%pw,error=error) 00156 END IF 00157 00158 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 00159 error=error) 00160 00161 CALL qs_vxc_create(qs_env=qs_env, rho_struct=harris_env%rho, xc_section=xc_section, & 00162 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=harris_energy%Exc, & 00163 just_energy=.FALSE., error=error) 00164 00165 DO ispin = 1,nspins 00166 harris_energy%integral_vxc = pw_integral_ab(v_rspace_new(ispin)%pw, & 00167 rho%rho_r(ispin)%pw,error=error) 00168 END DO 00169 00170 IF (ASSOCIATED(v_rspace_new)) THEN 00171 DO ispin = 1,nspins 00172 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(ispin)%pw, error=error) 00173 END DO 00174 DEALLOCATE(v_rspace_new,stat=stat) 00175 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00176 END IF 00177 IF (ASSOCIATED(v_tau_rspace)) THEN 00178 DO ispin = 1,nspins 00179 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw, error=error) 00180 END DO 00181 DEALLOCATE(v_tau_rspace,stat=stat) 00182 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00183 END IF 00184 00185 IF (my_EII_necessary) THEN 00186 IF (.NOT. fast_flag) THEN 00187 qs_env%energy%core_overlap = 0.0_dp 00188 qs_env%energy%core_self = 0.0_dp 00189 00190 CALL calculate_ecore_overlap(qs_env=qs_env, para_env=para_env, & 00191 calculate_forces=.FALSE., & 00192 E_overlap_core=harris_energy%Ecore_overlap,& 00193 error=error) 00194 CALL calculate_ecore_self(qs_env, E_self_core=harris_energy%Ecore_self,error=error) 00195 ELSE 00196 harris_energy%Ecore_overlap = qs_env%energy%core_overlap 00197 harris_energy%Ecore_self = qs_env%energy%core_self 00198 END IF 00199 00200 harris_energy%EII = harris_energy%Ecore_overlap + harris_energy%Ecore_self + & 00201 harris_energy%Ehartree_core 00202 !IF (PRESENT(EII_necessary)) EII_necessary=.FALSE. 00203 END IF 00204 00205 harris_energy%Eharris_correction = - harris_energy%Ehartree_elec & 00206 - harris_energy%integral_vxc & 00207 + harris_energy%Exc + harris_energy%EII 00208 END IF 00209 00210 CALL timestop(handle) 00211 END SUBROUTINE harris_energy_correction 00212 00213 ! ***************************************************************************** 00225 SUBROUTINE harris_eigenvalue_summation(qs_env, harris_env, error) 00226 00227 TYPE(qs_environment_type), POINTER :: qs_env 00228 TYPE(harris_env_type), POINTER :: harris_env 00229 TYPE(cp_error_type), INTENT(INOUT) :: error 00230 00231 CHARACTER(len=*), PARAMETER :: routineN = 'harris_eigenvalue_summation', 00232 routineP = moduleN//':'//routineN 00233 00234 INTEGER :: handle, homo, ispin, 00235 iterator, nspins 00236 LOGICAL :: failure 00237 REAL(KIND=dp) :: sum_of_eigenvalues 00238 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, 00239 occupation_numbers 00240 TYPE(dft_control_type), POINTER :: dft_control 00241 TYPE(harris_energy_type), POINTER :: harris_energy 00242 TYPE(mo_set_p_type), DIMENSION(:), 00243 POINTER :: mo_array 00244 00245 ! ------------------------------------------------------------------------ 00246 00247 CALL timeset(routineN,handle) 00248 00249 failure = .FALSE. 00250 NULLIFY(eigenvalues, occupation_numbers, mo_array, harris_energy, dft_control) 00251 00252 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 00253 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 00254 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00255 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00256 00257 sum_of_eigenvalues = 0.0_dp 00258 00259 IF (.NOT. failure) THEN 00260 CALL harris_env_get(harris_env=harris_env, harris_energy=harris_energy, error=error) 00261 00262 CALL get_qs_env(qs_env, dft_control=dft_control, mos=mo_array, error=error) 00263 00264 nspins = dft_control%nspins 00265 DO ispin = 1,nspins 00266 CALL get_mo_set(mo_array(ispin)%mo_set, homo=homo, eigenvalues=eigenvalues, & 00267 occupation_numbers = occupation_numbers) 00268 END DO 00269 00270 DO iterator=1,homo 00271 sum_of_eigenvalues = sum_of_eigenvalues + occupation_numbers(iterator) & 00272 * eigenvalues(iterator) 00273 END DO 00274 00275 ! Write the sum of eigenvalues back to the harris_energy_type 00276 harris_energy%sum_of_eigenvalues = sum_of_eigenvalues 00277 00278 harris_energy%Eharris = harris_energy%sum_of_eigenvalues + & 00279 harris_energy%Eharris_correction 00280 00281 END IF 00282 00283 CALL timestop(handle) 00284 00285 END SUBROUTINE harris_eigenvalue_summation 00286 00287 ! ***************************************************************************** 00300 SUBROUTINE harris_eigenvalue_calculation(qs_env, harris_env, error) 00301 00302 TYPE(qs_environment_type), POINTER :: qs_env 00303 TYPE(harris_env_type), POINTER :: harris_env 00304 TYPE(cp_error_type), INTENT(INOUT) :: error 00305 00306 CHARACTER(len=*), PARAMETER :: 00307 routineN = 'harris_eigenvalue_calculation', 00308 routineP = moduleN//':'//routineN 00309 00310 INTEGER :: handle, ispin, nspins 00311 LOGICAL :: failure 00312 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00313 POINTER :: matrix_ks, rho_ao 00314 TYPE(cp_para_env_type), POINTER :: para_env 00315 TYPE(dft_control_type), POINTER :: dft_control 00316 TYPE(harris_energy_type), POINTER :: harris_energy 00317 TYPE(qs_rho_type), POINTER :: rho 00318 00319 ! ------------------------------------------------------------------------ 00320 00321 CALL timeset(routineN,handle) 00322 00323 failure = .FALSE. 00324 NULLIFY(rho, rho_ao, matrix_ks, harris_energy, dft_control) 00325 00326 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 00327 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 00328 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00329 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00330 00331 IF (.NOT. failure) THEN 00332 CALL harris_env_get(harris_env=harris_env, harris_energy=harris_energy, error=error) 00333 00334 CALL get_qs_env(qs_env, rho=rho, dft_control=dft_control, matrix_ks=matrix_ks, & 00335 para_env=para_env, error=error) 00336 CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, error=error) 00337 00338 nspins = dft_control%nspins 00339 DO ispin = 1,nspins 00340 CALL cp_dbcsr_trace(matrix_ks(ispin)%matrix, rho_ao(ispin)%matrix, & 00341 trace=harris_energy%sum_of_eigenvalues, error=error) 00342 END DO 00343 00344 harris_energy%Eharris = harris_energy%sum_of_eigenvalues & 00345 + harris_energy%Eharris_correction 00346 END IF 00347 00348 CALL timestop(handle) 00349 00350 END SUBROUTINE harris_eigenvalue_calculation 00351 00352 ! ***************************************************************************** 00366 SUBROUTINE harris_eigenvalue_trace_KS_Pmix(scf_env, qs_env, harris_env, error) 00367 00368 TYPE(qs_scf_env_type), POINTER :: scf_env 00369 TYPE(qs_environment_type), POINTER :: qs_env 00370 TYPE(harris_env_type), POINTER :: harris_env 00371 TYPE(cp_error_type), INTENT(INOUT) :: error 00372 00373 CHARACTER(len=*), PARAMETER :: 00374 routineN = 'harris_eigenvalue_trace_KS_Pmix', 00375 routineP = moduleN//':'//routineN 00376 00377 INTEGER :: handle, ispin, nspins 00378 LOGICAL :: failure 00379 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00380 POINTER :: matrix_ks 00381 TYPE(cp_para_env_type), POINTER :: para_env 00382 TYPE(dft_control_type), POINTER :: dft_control 00383 TYPE(harris_energy_type), POINTER :: harris_energy 00384 00385 ! ------------------------------------------------------------------------ 00386 00387 CALL timeset(routineN,handle) 00388 00389 failure = .FALSE. 00390 NULLIFY(matrix_ks, harris_energy, dft_control) 00391 00392 CPPrecondition(ASSOCIATED(scf_env),cp_failure_level, routineP, error, failure) 00393 CPPrecondition(scf_env%ref_count>0, CP_failure_level, routineP, error, failure) 00394 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 00395 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 00396 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00397 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00398 00399 IF (.NOT. failure) THEN 00400 CALL harris_env_get(harris_env=harris_env, harris_energy=harris_energy, error=error) 00401 00402 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks, & 00403 para_env=para_env, error=error) 00404 00405 nspins = dft_control%nspins 00406 DO ispin = 1,nspins 00407 CALL cp_dbcsr_trace(matrix_ks(ispin)%matrix, scf_env%p_mix_new(ispin)%matrix, & 00408 trace=harris_energy%sum_of_eigenvalues, & 00409 error=error) 00410 END DO 00411 00412 harris_energy%Eharris = harris_energy%sum_of_eigenvalues + & 00413 harris_energy%Eharris_correction 00414 END IF 00415 00416 CALL timestop(handle) 00417 00418 END SUBROUTINE harris_eigenvalue_trace_KS_Pmix 00419 00420 ! ***************************************************************************** 00429 SUBROUTINE harris_postprocessing(harris_env, error) 00430 00431 TYPE(harris_env_type), POINTER :: harris_env 00432 TYPE(cp_error_type), INTENT(INOUT) :: error 00433 00434 CHARACTER(len=*), PARAMETER :: routineN = 'harris_postprocessing', 00435 routineP = moduleN//':'//routineN 00436 00437 INTEGER :: unit_nr 00438 LOGICAL :: failure 00439 TYPE(cp_logger_type), POINTER :: logger 00440 TYPE(harris_energy_type), POINTER :: harris_energy 00441 00442 ! ------------------------------------------------------------------------ 00443 00444 failure = .FALSE. 00445 NULLIFY(harris_energy) 00446 00447 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00448 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00449 logger => cp_error_get_logger(error) 00450 00451 IF (.NOT. failure) THEN 00452 CALL harris_env_get(harris_env=harris_env, harris_energy=harris_energy, error=error) 00453 00454 ! Output 00455 unit_nr=cp_logger_get_default_io_unit(logger) 00456 IF (unit_nr>0) THEN 00457 WRITE (unit_nr,*) "" 00458 WRITE (unit_nr,*) "The Harris functional energy correction is performed!" 00459 WRITE (unit_nr,*) "" 00460 WRITE (unit_nr,*) "Ehartree n_elec =", harris_energy%Ehartree_elec 00461 WRITE (unit_nr,*) "Ehartree n_core =", harris_energy%Ehartree_core 00462 00463 WRITE (unit_nr,*) "Exc =", harris_energy%Exc 00464 00465 WRITE (unit_nr,*) "The XC potential integral =", harris_energy%integral_vxc 00466 00467 WRITE (unit_nr,*) "Ecore_overlap =", harris_energy%Ecore_overlap 00468 WRITE (unit_nr,*) "Ecore_self =", harris_energy%Ecore_self 00469 00470 WRITE (unit_nr,*) "EII =", harris_energy%EII 00471 00472 WRITE (unit_nr,*) "Eharris correction energy =", harris_energy%Eharris_correction 00473 00474 WRITE (unit_nr,*) "The sum of the eigenvalues =", harris_energy%sum_of_eigenvalues 00475 WRITE (unit_nr,*) "Eharris =", harris_energy%Eharris 00476 WRITE (unit_nr,*) "" 00477 END IF 00478 END IF 00479 00480 END SUBROUTINE harris_postprocessing 00481 00482 END MODULE harris_functional
1.7.3