|
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_force 00013 00014 USE atomic_kind_types, ONLY: atomic_kind_type,& 00015 get_atomic_kind_set 00016 USE cell_types, ONLY: cell_type 00017 USE cp_control_types, ONLY: dft_control_type 00018 USE cp_dbcsr_types, ONLY: cp_dbcsr_p_type 00019 USE cp_para_types, ONLY: cp_para_env_type 00020 USE f77_blas 00021 USE global_types, ONLY: global_environment_type 00022 USE harris_env_types, ONLY: harris_env_get,& 00023 harris_env_type 00024 USE harris_force_types, ONLY: harris_force_type 00025 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00026 section_vals_type 00027 USE kinds, ONLY: dp 00028 USE message_passing, ONLY: mp_sum 00029 USE particle_types, ONLY: particle_type 00030 USE pw_env_types, ONLY: pw_env_get,& 00031 pw_env_type 00032 USE pw_methods, ONLY: pw_axpy,& 00033 pw_copy,& 00034 pw_scale,& 00035 pw_transfer,& 00036 pw_zero 00037 USE pw_poisson_methods, ONLY: pw_poisson_solve 00038 USE pw_poisson_types, ONLY: pw_poisson_type 00039 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00040 pw_pool_give_back_pw,& 00041 pw_pool_type 00042 USE pw_types, ONLY: COMPLEXDATA1D,& 00043 REALDATA3D,& 00044 REALSPACE,& 00045 RECIPROCALSPACE,& 00046 pw_p_type,& 00047 pw_release 00048 USE qs_core_energies, ONLY: calculate_ecore_overlap,& 00049 calculate_ecore_self 00050 USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix 00051 USE qs_environment_types, ONLY: get_qs_env,& 00052 qs_environment_type,& 00053 set_qs_env 00054 USE qs_force_types, ONLY: deallocate_qs_force,& 00055 duplicate_qs_force,& 00056 qs_force_type,& 00057 zero_qs_force 00058 USE qs_integrate_potential, ONLY: integrate_v_core_rspace,& 00059 integrate_v_rspace 00060 USE qs_ks_types, ONLY: qs_ks_env_type 00061 USE qs_rho_methods, ONLY: diff_rho_type 00062 USE qs_rho_types, ONLY: qs_rho_get,& 00063 qs_rho_type 00064 USE qs_vxc, ONLY: qs_vxc_create 00065 USE timings, ONLY: timeset,& 00066 timestop 00067 USE xc, ONLY: xc_calc_2nd_deriv,& 00068 xc_prep_2nd_deriv 00069 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00070 xc_dset_release 00071 USE xc_rho_set_types, ONLY: xc_rho_set_release,& 00072 xc_rho_set_type 00073 #include "cp_common_uses.h" 00074 00075 IMPLICIT NONE 00076 PRIVATE 00077 00078 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'harris_force' 00079 00080 ! *** Public subroutines *** 00081 PUBLIC :: harris_force_correction, & 00082 harris_force_EVal, & 00083 harris_calc_nsc_force 00084 00085 !*** 00086 00087 CONTAINS 00088 00089 ! ***************************************************************************** 00103 SUBROUTINE harris_force_correction(qs_env, harris_env, globenv, error) 00104 00105 TYPE(qs_environment_type), POINTER :: qs_env 00106 TYPE(harris_env_type), POINTER :: harris_env 00107 TYPE(global_environment_type), POINTER :: globenv 00108 TYPE(cp_error_type), INTENT(INOUT) :: error 00109 00110 CHARACTER(len=*), PARAMETER :: routineN = 'harris_force_correction', 00111 routineP = moduleN//':'//routineN 00112 00113 INTEGER :: handle, i, iatom, ikind, 00114 ispin, natom, nkind, nspins, 00115 stat, unit_nr 00116 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 00117 natom_of_kind 00118 LOGICAL :: failure 00119 REAL(KIND=dp) :: Ehartree, Exc 00120 TYPE(atomic_kind_type), DIMENSION(:), 00121 POINTER :: atomic_kind_set 00122 TYPE(cell_type), POINTER :: cell 00123 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00124 POINTER :: matrix_ks, rho_ao 00125 TYPE(cp_logger_type), POINTER :: logger 00126 TYPE(cp_para_env_type), POINTER :: para_env 00127 TYPE(dft_control_type), POINTER :: dft_control 00128 TYPE(harris_force_type), POINTER :: harris_force 00129 TYPE(particle_type), DIMENSION(:), 00130 POINTER :: particle_set 00131 TYPE(pw_env_type), POINTER :: pw_env 00132 TYPE(pw_p_type) :: rho_tot_gspace, 00133 v_hartree_gspace, 00134 v_hartree_rspace 00135 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r, 00136 v_rspace_new, v_tau_rspace, 00137 v_xc 00138 TYPE(pw_p_type), POINTER :: rho_core 00139 TYPE(pw_poisson_type), POINTER :: poisson_env 00140 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00141 TYPE(qs_force_type), DIMENSION(:), 00142 POINTER :: force, force_tmp 00143 TYPE(qs_ks_env_type), POINTER :: ks_env 00144 TYPE(qs_rho_type), POINTER :: rho 00145 TYPE(section_vals_type), POINTER :: input, xc_section 00146 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00147 TYPE(xc_rho_set_type), POINTER :: rho_set 00148 00149 !E_overlap_core 00150 ! ------------------------------------------------------------------------ 00151 00152 CALL timeset(routineN,handle) 00153 00154 failure = .FALSE. 00155 NULLIFY(force, force_tmp, v_rspace_new, v_tau_rspace, rho, matrix_ks, & 00156 pw_env, auxbas_pw_pool, rho_core, cell, dft_control, xc_section, & 00157 atomic_kind_set, particle_set, deriv_set, rho_set, my_rho_r, & 00158 v_xc, harris_force, rho_r, ks_env, poisson_env, input, rho_ao) 00159 para_env=>qs_env%para_env 00160 logger => cp_error_get_logger(error) 00161 00162 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 00163 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 00164 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00165 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00166 CPPrecondition(ASSOCIATED(globenv), cp_failure_level, routineP, error, failure) 00167 CPPrecondition(globenv%ref_count>0, cp_failure_level, routineP, error, failure) 00168 00169 IF (.NOT. failure) THEN 00170 CALL harris_env_get(harris_env=harris_env, harris_force=harris_force, error=error) 00171 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force, & 00172 dft_control=dft_control, particle_set=particle_set, & 00173 cell=cell, input=input, error=error) 00174 00175 xc_section => section_vals_get_subs_vals(input, "DFT%XC", error=error) 00176 00177 natom = SIZE(particle_set) 00178 nkind = SIZE(force) 00179 nspins = dft_control%nspins 00180 00181 ALLOCATE (atom_of_kind(natom),STAT=stat) 00182 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00183 00184 ALLOCATE (kind_of(natom),STAT=stat) 00185 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00186 00187 ALLOCATE (natom_of_kind(nkind),STAT=stat) 00188 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00189 00190 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 00191 atom_of_kind=atom_of_kind, kind_of=kind_of, & 00192 natom_of_kind=natom_of_kind) 00193 00194 CALL duplicate_qs_force(qs_force_input=force, qs_force_output=force_tmp, & 00195 natom_of_kind=natom_of_kind) 00196 CALL zero_qs_force(force) 00197 00198 ! *** d/dR[Sum of eigenvalues] = d/dR[trace(H_in * rho_out)] *** ! 00199 CALL build_core_hamiltonian_matrix(qs_env=qs_env,calculate_forces=.TRUE.,error=error) 00200 CALL calculate_ecore_self(qs_env,error=error) 00201 CALL calculate_ecore_overlap(qs_env, para_env, & 00202 calculate_forces=.TRUE.,& 00203 error=error) 00204 00205 DO ikind = 1,SIZE(force) 00206 CALL mp_sum(force(ikind)%kinetic, para_env%group) 00207 CALL mp_sum(force(ikind)%overlap, para_env%group) 00208 CALL mp_sum(force(ikind)%gth_ppl, para_env%group) 00209 CALL mp_sum(force(ikind)%gth_nlcc, para_env%group) 00210 CALL mp_sum(force(ikind)%gth_ppnl, para_env%group) 00211 CALL mp_sum(force(ikind)%core_overlap, para_env%group) 00212 END DO 00213 00214 DO iatom = 1, natom 00215 ikind = kind_of(iatom) 00216 i = atom_of_kind(iatom) 00217 harris_force%f_kinetic(iatom, 1:3) = force(ikind)%kinetic(1:3,i) 00218 harris_force%f_overlap(iatom, 1:3) = force(ikind)%overlap(1:3,i) 00219 harris_force%f_gth_pp(iatom, 1:3) = force(ikind)%gth_ppl(1:3,i) & 00220 + force(ikind)%gth_ppnl(1:3,i) + force(ikind)%gth_nlcc(1:3,i) 00221 harris_force%f_ovrl(iatom, 1:3) = force(ikind)%core_overlap(1:3,i) 00222 harris_force%f_self(iatom, 1:3) = 0.0_dp 00223 END DO 00224 00225 CALL zero_qs_force(force) 00226 00227 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, rho_core=rho_core, & 00228 ks_env=ks_env, matrix_ks=matrix_ks, error=error) 00229 00230 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, & 00231 poisson_env=poisson_env, error=error) 00232 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 00233 use_data = COMPLEXDATA1D, & 00234 in_space = RECIPROCALSPACE, error=error) 00235 CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace%pw, & 00236 use_data=COMPLEXDATA1D, & 00237 in_space=RECIPROCALSPACE, error=error) 00238 00239 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace%pw, & 00240 use_data=REALDATA3D, in_space=REALSPACE,error=error) 00241 00242 ! *** Calculation of the forces due to the Hartree energy *** ! 00243 CALL pw_copy(rho_core%pw,rho_tot_gspace%pw,error=error) 00244 DO ispin = 1,nspins 00245 CALL pw_axpy(harris_env%rho%rho_g(ispin)%pw,rho_tot_gspace%pw,error=error) 00246 ! old code seems to be wrong for spin polarized case 00247 ! CALL pw_add(rho_core%pw, harris_env%rho%rho_g(ispin)%pw, rho_tot_gspace%pw) 00248 END DO 00249 00250 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree=Ehartree, & 00251 vhartree=v_hartree_gspace%pw,error=error) 00252 00253 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00254 00255 CALL qs_vxc_create(qs_env=qs_env, rho_struct=qs_env%rho, xc_section=xc_section, & 00256 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=Exc, & 00257 just_energy=.FALSE., error=error) 00258 00259 rho_ao => qs_env%rho%rho_ao 00260 00261 DO ispin = 1,nspins 00262 CALL pw_copy(v_rspace_new(ispin)%pw,v_rspace_new(ispin)%pw,error=error) 00263 CALL pw_axpy(v_hartree_rspace%pw,v_rspace_new(ispin)%pw,error=error) 00264 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 00265 00266 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), p=rho_ao(ispin), & 00267 h=matrix_ks(ispin), qs_env=qs_env, calculate_forces=.TRUE., & 00268 gapw=.FALSE., error=error) 00269 END DO 00270 00271 DO ikind = 1,SIZE(force) 00272 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 00273 END DO 00274 00275 DO iatom = 1, natom 00276 ikind = kind_of(iatom) 00277 i = atom_of_kind(iatom) 00278 harris_force%f_hartree(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 00279 END DO 00280 00281 CALL zero_qs_force(force) 00282 00283 ! *** Calculation of the forces of the xc-integral *** ! 00284 CALL qs_rho_get(rho_struct=rho, rho_r=rho_r, error=error) 00285 ALLOCATE(my_rho_r(SIZE(rho_r)), stat=stat) 00286 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00287 DO ispin=1, SIZE(rho_r) 00288 my_rho_r(ispin)%pw => harris_env%rho%rho_r(ispin)%pw 00289 END DO 00290 CALL xc_prep_2nd_deriv(deriv_set=deriv_set, rho_set=rho_set, & 00291 rho_r=my_rho_r, pw_pool=auxbas_pw_pool, & 00292 xc_section=xc_section, cell=cell, error=error) 00293 DEALLOCATE(my_rho_r, stat=stat) 00294 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00295 00296 ! ** v_xc ** ! 00297 ALLOCATE(v_xc(nspins), stat=stat) 00298 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00299 DO ispin=1, nspins 00300 NULLIFY(v_xc(ispin)%pw) 00301 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, use_data=REALDATA3D, & 00302 in_space=REALSPACE, error=error) 00303 CALL pw_zero(v_xc(ispin)%pw, error=error) 00304 END DO 00305 00306 CALL xc_calc_2nd_deriv(v_xc, deriv_set=deriv_set, rho_set=rho_set, & 00307 rho1_set=rho_set, pw_pool=auxbas_pw_pool, & 00308 xc_section=xc_section, gapw=.FALSE., error=error) 00309 00310 CALL xc_dset_release(deriv_set, error=error) 00311 00312 DO ispin = 1,nspins 00313 v_xc(ispin)%pw%cr3d = v_xc(ispin)%pw%cr3d * v_xc(ispin)%pw%pw_grid%dvol 00314 v_rspace_new(ispin)%pw%cr3d = (qs_env%rho%rho_r(ispin)%pw%cr3d & 00315 - harris_env%rho%rho_r(ispin)%pw%cr3d) & 00316 * v_xc(ispin)%pw%cr3d 00317 CALL pw_release(v_xc(ispin)%pw, error=error) 00318 END DO 00319 00320 DEALLOCATE(v_xc,stat=stat) 00321 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00322 00323 CALL xc_rho_set_release(rho_set,error=error) 00324 00325 ! + v_hartree[rho_out-rho_in] 00326 CALL pw_zero(rho_tot_gspace%pw, error=error) 00327 DO ispin = 1,nspins 00328 CALL pw_axpy(qs_env%rho%rho_g(ispin)%pw,rho_tot_gspace%pw,error=error) 00329 CALL pw_axpy(harris_env%rho%rho_g(ispin)%pw,rho_tot_gspace%pw,alpha=-1._dp,error=error) 00330 ! this original code seems to be wrong for spin polarized cases (overwrites rho_tot_gspace) 00331 ! CALL pw_subtract(qs_env%rho%rho_g(ispin)%pw, harris_env%rho%rho_g(ispin)%pw, rho_tot_gspace%pw) 00332 END DO 00333 00334 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree=Ehartree, & 00335 vhartree=v_hartree_gspace%pw,error=error) 00336 00337 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00338 00339 00340 rho_ao => harris_env%rho%rho_ao 00341 00342 DO ispin = 1,nspins 00343 CALL pw_copy(v_rspace_new(ispin)%pw,v_rspace_new(ispin)%pw,error=error) 00344 CALL pw_axpy(v_hartree_rspace%pw,v_rspace_new(ispin)%pw,error=error) 00345 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 00346 00347 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), p=rho_ao(ispin), & 00348 h=matrix_ks(ispin), qs_env=qs_env, calculate_forces=.TRUE., & 00349 gapw=.FALSE., error=error) 00350 END DO 00351 00352 IF (ASSOCIATED(v_rspace_new)) THEN 00353 DO ispin = 1,nspins 00354 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(ispin)%pw, error=error) 00355 END DO 00356 DEALLOCATE(v_rspace_new, stat=stat) 00357 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00358 END IF 00359 IF (ASSOCIATED(v_tau_rspace)) THEN 00360 DO ispin = 1,nspins 00361 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw, error=error) 00362 END DO 00363 DEALLOCATE(v_tau_rspace, stat=stat) 00364 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00365 END IF 00366 00367 DO ikind = 1,SIZE(force) 00368 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 00369 END DO 00370 00371 DO iatom = 1, natom 00372 ikind = kind_of(iatom) 00373 i = atom_of_kind(iatom) 00374 harris_force%f_integral_vxc(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 00375 END DO 00376 00377 DO iatom = 1, natom 00378 ikind = kind_of(iatom) 00379 i = atom_of_kind(iatom) 00380 harris_force%f_trace(iatom, 1:3) = harris_force%f_kinetic(iatom, 1:3) & 00381 + harris_force%f_gth_pp(iatom, 1:3) & 00382 + harris_force%f_overlap(iatom, 1:3) & 00383 + harris_force%f_hartree(iatom, 1:3) & 00384 + harris_force%f_integral_vxc(iatom, 1:3) 00385 END DO 00386 00387 DO iatom = 1,natom 00388 harris_force%f_hartree(iatom, 1:3) = 0.0_dp 00389 END DO 00390 00391 CALL zero_qs_force(force) 00392 00393 ! *** The forces due to n_core *** ! 00394 CALL pw_copy(rho_core%pw,rho_tot_gspace%pw,error=error) 00395 DO ispin = 1,nspins 00396 CALL pw_axpy(qs_env%rho%rho_g(ispin)%pw,rho_tot_gspace%pw,error=error) 00397 ! old code seems to be wrong for spin polarized case 00398 ! CALL pw_add(rho_core%pw, qs_env%rho%rho_g(ispin)%pw, rho_tot_gspace%pw) 00399 END DO 00400 00401 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree=Ehartree, & 00402 vhartree=v_hartree_gspace%pw,error=error) 00403 00404 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw, error=error) 00405 00406 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00407 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 00408 00409 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env,error=error) 00410 00411 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw, error=error) 00412 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw, error=error) 00413 00414 DO ikind = 1,SIZE(force) 00415 CALL mp_sum(force(ikind)%rho_core, para_env%group) 00416 END DO 00417 00418 DO iatom = 1, natom 00419 ikind = kind_of(iatom) 00420 i = atom_of_kind(iatom) 00421 harris_force%f_rho_core(iatom, 1:3) = force(ikind)%rho_core(1:3,i) 00422 END DO 00423 00424 CALL deallocate_qs_force(force) 00425 CALL set_qs_env(qs_env=qs_env, force=force_tmp, error=error) 00426 00427 CALL get_qs_env(qs_env=qs_env, force=force, error=error) 00428 harris_force%f_total(:) = 0.0_dp 00429 00430 DO iatom = 1, natom 00431 ikind = kind_of(iatom) 00432 i = atom_of_kind(iatom) 00433 00434 harris_force%f_EII(iatom, 1:3) = harris_force%f_ovrl(iatom, 1:3) & 00435 + harris_force%f_self(iatom, 1:3) & 00436 + harris_force%f_rho_core(iatom, 1:3) 00437 harris_force%f_harris(iatom, 1:3) = harris_force%f_trace(iatom, 1:3) & 00438 + harris_force%f_EII(iatom, 1:3) 00439 harris_force%f_total(1:3) = harris_force%f_total(1:3) & 00440 + harris_force%f_harris(iatom, 1:3) 00441 END DO 00442 00443 ! Output 00444 unit_nr=cp_logger_get_default_io_unit(logger) 00445 IF (unit_nr>0) THEN 00446 WRITE (unit_nr,*) ""; WRITE (unit_nr, *) "" 00447 WRITE (unit_nr,*) "The Harris functional force correction is performed!" 00448 WRITE (unit_nr,*) "" 00449 00450 WRITE (unit_nr, *) "F_Harris" 00451 DO iatom = 1, natom 00452 WRITE (unit_nr,*) harris_force%f_harris(iatom, 1:3) 00453 END DO 00454 WRITE (unit_nr, *) "F_Total" 00455 WRITE (unit_nr, *) harris_force%f_total(1:3) 00456 END IF 00457 00458 DEALLOCATE (atom_of_kind,STAT=stat) 00459 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00460 00461 DEALLOCATE (kind_of,STAT=stat) 00462 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00463 00464 DEALLOCATE (natom_of_kind,STAT=stat) 00465 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00466 00467 END IF 00468 00469 CALL timestop(handle) 00470 00471 END SUBROUTINE harris_force_correction 00472 00473 ! ***************************************************************************** 00487 SUBROUTINE harris_force_EVal(qs_env, harris_env, globenv, error) 00488 00489 TYPE(qs_environment_type), POINTER :: qs_env 00490 TYPE(harris_env_type), POINTER :: harris_env 00491 TYPE(global_environment_type), POINTER :: globenv 00492 TYPE(cp_error_type), INTENT(INOUT) :: error 00493 00494 CHARACTER(len=*), PARAMETER :: routineN = 'harris_force_EVal', 00495 routineP = moduleN//':'//routineN 00496 00497 INTEGER :: handle, i, iatom, ikind, 00498 ispin, natom, nkind, nspins, 00499 stat, unit_nr 00500 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 00501 natom_of_kind 00502 LOGICAL :: failure 00503 REAL(KIND=dp) :: E_overlap_core, Ehartree, Exc 00504 TYPE(atomic_kind_type), DIMENSION(:), 00505 POINTER :: atomic_kind_set 00506 TYPE(cell_type), POINTER :: cell 00507 TYPE(cp_dbcsr_p_type), DIMENSION(:), 00508 POINTER :: matrix_ks, rho_ao 00509 TYPE(cp_logger_type), POINTER :: logger 00510 TYPE(cp_para_env_type), POINTER :: para_env 00511 TYPE(dft_control_type), POINTER :: dft_control 00512 TYPE(harris_force_type), POINTER :: harris_force 00513 TYPE(particle_type), DIMENSION(:), 00514 POINTER :: particle_set 00515 TYPE(pw_env_type), POINTER :: pw_env 00516 TYPE(pw_p_type) :: rho_tot_gspace, 00517 v_hartree_gspace, 00518 v_hartree_rspace 00519 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r, 00520 v_rspace_new, v_tau_rspace, 00521 v_xc 00522 TYPE(pw_p_type), POINTER :: rho_core 00523 TYPE(pw_poisson_type), POINTER :: poisson_env 00524 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 00525 TYPE(qs_force_type), DIMENSION(:), 00526 POINTER :: force, force_tmp 00527 TYPE(qs_ks_env_type), POINTER :: ks_env 00528 TYPE(qs_rho_type), POINTER :: rho 00529 TYPE(section_vals_type), POINTER :: input, xc_section 00530 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00531 TYPE(xc_rho_set_type), POINTER :: rho_set 00532 00533 ! ------------------------------------------------------------------------ 00534 00535 CALL timeset(routineN,handle) 00536 para_env=>qs_env%para_env 00537 00538 failure = .FALSE. 00539 NULLIFY(force, force_tmp, v_rspace_new, v_tau_rspace, rho, matrix_ks, & 00540 pw_env, auxbas_pw_pool, rho_core, cell, dft_control, xc_section, & 00541 atomic_kind_set, particle_set, deriv_set, rho_set, my_rho_r, & 00542 v_xc, harris_force, rho_r, ks_env, poisson_env, input) 00543 logger => cp_error_get_logger(error) 00544 00545 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 00546 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 00547 CPPrecondition(ASSOCIATED(harris_env), cp_failure_level, routineP, error, failure) 00548 CPPrecondition(harris_env%ref_count>0, cp_failure_level, routineP, error, failure) 00549 CPPrecondition(ASSOCIATED(globenv), cp_failure_level, routineP, error, failure) 00550 CPPrecondition(globenv%ref_count>0, cp_failure_level, routineP, error, failure) 00551 00552 IF (.NOT. failure) THEN 00553 CALL harris_env_get(harris_env=harris_env, harris_force=harris_force, error=error) 00554 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force, & 00555 dft_control=dft_control, cell=cell, particle_set=particle_set, & 00556 input=input, error=error) 00557 00558 xc_section => section_vals_get_subs_vals(input, "DFT%XC", error=error) 00559 00560 natom = SIZE(particle_set) 00561 nkind = SIZE(force) 00562 nspins = dft_control%nspins 00563 00564 ALLOCATE (atom_of_kind(natom),STAT=stat) 00565 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00566 00567 ALLOCATE (kind_of(natom),STAT=stat) 00568 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00569 00570 ALLOCATE (natom_of_kind(nkind),STAT=stat) 00571 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00572 00573 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 00574 atom_of_kind=atom_of_kind, kind_of=kind_of, & 00575 natom_of_kind=natom_of_kind) 00576 00577 CALL duplicate_qs_force(qs_force_input=force, qs_force_output=force_tmp, & 00578 natom_of_kind=natom_of_kind) 00579 CALL zero_qs_force(force) 00580 00581 ! *** d/dR[Sum of eigenvalues] = d/dR[trace(H_in * rho_out)] *** ! 00582 CALL build_core_hamiltonian_matrix(qs_env=qs_env,calculate_forces=.TRUE.,error=error) 00583 00584 DO ikind = 1,SIZE(force) 00585 CALL mp_sum(force(ikind)%kinetic, para_env%group) 00586 CALL mp_sum(force(ikind)%overlap, para_env%group) 00587 CALL mp_sum(force(ikind)%gth_ppl, para_env%group) 00588 CALL mp_sum(force(ikind)%gth_nlcc, para_env%group) 00589 CALL mp_sum(force(ikind)%gth_ppnl, para_env%group) 00590 END DO 00591 00592 DO iatom = 1, natom 00593 ikind = kind_of(iatom) 00594 i = atom_of_kind(iatom) 00595 harris_force%f_kinetic(iatom, 1:3) = force(ikind)%kinetic(1:3,i) 00596 harris_force%f_overlap(iatom, 1:3) = force(ikind)%overlap(1:3,i) 00597 harris_force%f_gth_pp(iatom, 1:3) = force(ikind)%gth_ppl(1:3,i) & 00598 + force(ikind)%gth_ppnl(1:3,i) + force(ikind)%gth_nlcc(1:3,i) 00599 END DO 00600 00601 CALL zero_qs_force(force) 00602 00603 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, rho_core=rho_core, & 00604 ks_env=ks_env, matrix_ks=matrix_ks, & 00605 error=error) 00606 00607 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, & 00608 poisson_env=poisson_env, error=error) 00609 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 00610 use_data = COMPLEXDATA1D, & 00611 in_space = RECIPROCALSPACE, error=error) 00612 CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace%pw, & 00613 use_data=COMPLEXDATA1D, & 00614 in_space=RECIPROCALSPACE, error=error) 00615 00616 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace%pw, & 00617 use_data=REALDATA3D, in_space=REALSPACE,error=error) 00618 00619 ! *** Calculation of the forces due to the Hartree energy *** ! 00620 CALL pw_copy(rho_core%pw,rho_tot_gspace%pw,error=error) 00621 DO ispin = 1,nspins 00622 CALL pw_axpy(harris_env%rho%rho_g(ispin)%pw,rho_tot_gspace%pw,error=error) 00623 ! old code seems to be wrong for spin polarized case 00624 ! CALL pw_add(rho_core%pw, harris_env%rho%rho_g(ispin)%pw, rho_tot_gspace%pw) 00625 END DO 00626 00627 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree=Ehartree, & 00628 vhartree=v_hartree_gspace%pw,error=error) 00629 00630 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw, error=error) 00631 00632 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00633 00634 CALL qs_vxc_create(qs_env=qs_env, rho_struct=qs_env%rho, xc_section=xc_section, & 00635 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=Exc, & 00636 just_energy=.FALSE., error=error) 00637 00638 rho_ao => qs_env%rho%rho_ao 00639 00640 DO ispin = 1,nspins 00641 CALL pw_copy(v_rspace_new(ispin)%pw,v_rspace_new(ispin)%pw,error=error) 00642 CALL pw_axpy(v_hartree_rspace%pw,v_rspace_new(ispin)%pw,error=error) 00643 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 00644 00645 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), p=rho_ao(ispin), & 00646 h=matrix_ks(ispin), qs_env=qs_env, calculate_forces=.TRUE., & 00647 gapw=.FALSE., error=error) 00648 END DO 00649 00650 DO ikind = 1,SIZE(force) 00651 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 00652 END DO 00653 00654 DO iatom = 1, natom 00655 ikind = kind_of(iatom) 00656 i = atom_of_kind(iatom) 00657 harris_force%f_hartree(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 00658 END DO 00659 00660 CALL zero_qs_force(force) 00661 00662 ! ** Force computation of f_V ** ! 00663 DO ispin = 1,nspins 00664 CALL pw_poisson_solve(poisson_env, qs_env%rho%rho_g(ispin)%pw, ehartree=Ehartree, & 00665 vhartree=v_hartree_gspace%pw,error=error) 00666 END DO 00667 00668 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00669 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 00670 00671 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env,error=error) 00672 00673 DO ikind = 1,SIZE(force) 00674 CALL mp_sum(force(ikind)%rho_core, para_env%group) 00675 END DO 00676 00677 DO iatom = 1, natom 00678 ikind = kind_of(iatom) 00679 i = atom_of_kind(iatom) 00680 harris_force%f_V(iatom, 1:3) = force(ikind)%rho_core(1:3,i) 00681 END DO 00682 00683 CALL zero_qs_force(force) 00684 00685 ! *** Calculation of the forces of the xc-integral *** ! 00686 CALL qs_rho_get(rho_struct=rho, rho_r=rho_r, error=error) 00687 ALLOCATE(my_rho_r(SIZE(rho_r)), stat=stat) 00688 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00689 DO ispin=1, SIZE(rho_r) 00690 my_rho_r(ispin)%pw => harris_env%rho%rho_r(ispin)%pw 00691 END DO 00692 CALL xc_prep_2nd_deriv(deriv_set=deriv_set, rho_set=rho_set, & 00693 rho_r=my_rho_r, pw_pool=auxbas_pw_pool, & 00694 xc_section=xc_section, cell=cell, error=error) 00695 DEALLOCATE(my_rho_r, stat=stat) 00696 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00697 00698 ! ** v_xc ** ! 00699 ALLOCATE(v_xc(nspins), stat=stat) 00700 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00701 DO ispin=1, nspins 00702 NULLIFY(v_xc(ispin)%pw) 00703 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, use_data=REALDATA3D, & 00704 in_space=REALSPACE, error=error) 00705 CALL pw_zero(v_xc(ispin)%pw, error=error) 00706 END DO 00707 00708 CALL xc_calc_2nd_deriv(v_xc, deriv_set=deriv_set, rho_set=rho_set, & 00709 rho1_set=rho_set, pw_pool=auxbas_pw_pool, & 00710 xc_section=xc_section, gapw=.FALSE., error=error) 00711 00712 CALL xc_dset_release(deriv_set, error=error) 00713 00714 DO ispin = 1,nspins 00715 v_xc(ispin)%pw%cr3d = v_xc(ispin)%pw%cr3d * v_xc(ispin)%pw%pw_grid%dvol 00716 v_rspace_new(ispin)%pw%cr3d = rho_r(ispin)%pw%cr3d * v_xc(ispin)%pw%cr3d 00717 00718 CALL pw_release(v_xc(ispin)%pw, error=error) 00719 END DO 00720 00721 DEALLOCATE(v_xc,stat=stat) 00722 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00723 00724 CALL xc_rho_set_release(rho_set,error=error) 00725 00726 ! + v_hartree[rho_out] 00727 DO ispin = 1,nspins 00728 CALL pw_poisson_solve(poisson_env, qs_env%rho%rho_g(ispin)%pw, ehartree=Ehartree, & 00729 vhartree=v_hartree_gspace%pw,error=error) 00730 END DO 00731 00732 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00733 00734 00735 rho_ao => harris_env%rho%rho_ao 00736 00737 DO ispin = 1,nspins 00738 CALL pw_copy(v_rspace_new(ispin)%pw,v_rspace_new(ispin)%pw,error=error) 00739 CALL pw_axpy(v_hartree_rspace%pw,v_rspace_new(ispin)%pw,error=error) 00740 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 00741 00742 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), p=rho_ao(ispin), & 00743 h=matrix_ks(ispin), qs_env=qs_env, calculate_forces=.TRUE., & 00744 gapw=.FALSE., error=error) 00745 END DO 00746 00747 IF (ASSOCIATED(v_rspace_new)) THEN 00748 DO ispin = 1,nspins 00749 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_rspace_new(ispin)%pw,error=error) 00750 END DO 00751 DEALLOCATE(v_rspace_new,stat=stat) 00752 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00753 END IF 00754 IF (ASSOCIATED(v_tau_rspace)) THEN 00755 DO ispin = 1,nspins 00756 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_tau_rspace(ispin)%pw,error=error) 00757 END DO 00758 DEALLOCATE(v_tau_rspace,stat=stat) 00759 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00760 END IF 00761 00762 DO ikind = 1,SIZE(force) 00763 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 00764 END DO 00765 00766 DO iatom = 1, natom 00767 ikind = kind_of(iatom) 00768 i = atom_of_kind(iatom) 00769 harris_force%f_integral_vxc(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 00770 END DO 00771 00772 DO iatom = 1, natom 00773 ikind = kind_of(iatom) 00774 i = atom_of_kind(iatom) 00775 harris_force%f_trace(iatom, 1:3) = harris_force%f_kinetic(iatom, 1:3) & 00776 + harris_force%f_gth_pp(iatom, 1:3) & 00777 + harris_force%f_overlap(iatom, 1:3) & 00778 + harris_force%f_V(iatom, 1:3) & 00779 + harris_force%f_hartree(iatom, 1:3) & 00780 + harris_force%f_integral_vxc(iatom, 1:3) 00781 END DO 00782 00783 DO iatom = 1,natom 00784 harris_force%f_integral_vxc(iatom, 1:3) = 0.0_dp 00785 END DO 00786 00787 CALL zero_qs_force(force) 00788 ! *** End of force calculation due to the sum_of_eigenvalues *** ! 00789 00790 CALL harris_env_get(harris_env=harris_env, rho=rho, error=error) 00791 00792 ! *** The forces due to the core-core repulsion *** ! 00793 CALL pw_poisson_solve(poisson_env, rho_core%pw, ehartree=Ehartree, & 00794 vhartree=v_hartree_gspace%pw,error=error) 00795 00796 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00797 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 00798 00799 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env,error=error) 00800 00801 DO ikind = 1,SIZE(force) 00802 CALL mp_sum(force(ikind)%rho_core, para_env%group) 00803 END DO 00804 00805 DO iatom = 1, natom 00806 ikind = kind_of(iatom) 00807 i = atom_of_kind(iatom) 00808 harris_force%f_rho_core(iatom, 1:3) = force(ikind)%rho_core(1:3,i) 00809 END DO 00810 00811 CALL zero_qs_force(force) 00812 00813 ! *** The forces due to exchange and correlation *** ! 00814 CALL qs_vxc_create(qs_env=qs_env, rho_struct=qs_env%rho, xc_section=xc_section, & 00815 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=Exc, & 00816 just_energy=.FALSE., error=error) 00817 00818 ! *** Calculation of the forces of the xc-integral *** ! 00819 CALL qs_rho_get(rho, rho_r=rho_r, error=error) 00820 ALLOCATE(my_rho_r(SIZE(rho_r)), stat=stat) 00821 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00822 DO ispin=1, SIZE(rho_r) 00823 my_rho_r(ispin)%pw => harris_env%rho%rho_r(ispin)%pw 00824 END DO 00825 CALL xc_prep_2nd_deriv(deriv_set=deriv_set, rho_set=rho_set, & 00826 rho_r=my_rho_r, pw_pool=auxbas_pw_pool, & 00827 xc_section=xc_section, cell=cell, error=error) 00828 DEALLOCATE(my_rho_r, stat=stat) 00829 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00830 00831 ! ** v_xc ** ! 00832 nspins = dft_control%nspins 00833 ALLOCATE(v_xc(nspins), stat=stat) 00834 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 00835 DO ispin=1, nspins 00836 NULLIFY(v_xc(ispin)%pw) 00837 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, use_data=REALDATA3D, & 00838 in_space=REALSPACE, error=error) 00839 CALL pw_zero(v_xc(ispin)%pw, error=error) 00840 END DO 00841 00842 CALL xc_calc_2nd_deriv(v_xc, deriv_set=deriv_set, rho_set=rho_set, & 00843 rho1_set=rho_set, pw_pool=auxbas_pw_pool, & 00844 xc_section=xc_section, gapw=.FALSE., error=error) 00845 00846 CALL xc_dset_release(deriv_set, error=error) 00847 00848 DO ispin = 1,nspins 00849 v_xc(ispin)%pw%cr3d = v_xc(ispin)%pw%cr3d * v_xc(ispin)%pw%pw_grid%dvol 00850 v_rspace_new(ispin)%pw%cr3d = rho_r(ispin)%pw%cr3d * v_xc(ispin)%pw%cr3d 00851 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 00852 00853 CALL pw_release(v_xc(ispin)%pw, error=error) 00854 END DO 00855 00856 DEALLOCATE(v_xc,stat=stat) 00857 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00858 00859 CALL xc_rho_set_release(rho_set,error=error) 00860 00861 DO ispin = 1,nspins 00862 CALL pw_poisson_solve(poisson_env, rho%rho_g(ispin)%pw, ehartree=Ehartree, & 00863 vhartree=v_hartree_gspace%pw,error=error) 00864 END DO 00865 00866 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 00867 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 00868 00869 rho_ao => harris_env%rho%rho_ao 00870 DO ispin = 1,nspins 00871 CALL pw_copy(v_rspace_new(ispin)%pw,v_rspace_new(ispin)%pw,error=error) 00872 CALL pw_axpy(v_hartree_rspace%pw,v_rspace_new(ispin)%pw,error=error) 00873 00874 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), p=rho_ao(ispin), & 00875 h=matrix_ks(ispin), qs_env=qs_env, calculate_forces=.TRUE., & 00876 gapw=.FALSE., error=error) 00877 END DO 00878 00879 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw, error=error) 00880 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw, error=error) 00881 00882 IF (ASSOCIATED(v_rspace_new)) THEN 00883 DO ispin = 1,nspins 00884 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_rspace_new(ispin)%pw,error=error) 00885 END DO 00886 DEALLOCATE(v_rspace_new,stat=stat) 00887 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00888 END IF 00889 IF (ASSOCIATED(v_tau_rspace)) THEN 00890 DO ispin = 1,nspins 00891 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_tau_rspace(ispin)%pw,error=error) 00892 END DO 00893 DEALLOCATE(v_tau_rspace,stat=stat) 00894 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00895 END IF 00896 00897 DO ikind = 1,SIZE(force) 00898 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 00899 END DO 00900 00901 DO iatom = 1, natom 00902 ikind = kind_of(iatom) 00903 i = atom_of_kind(iatom) 00904 harris_force%f_integral_vxc(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 00905 END DO 00906 00907 CALL zero_qs_force(force) 00908 00909 CALL calculate_ecore_overlap(qs_env=qs_env, para_env=para_env, & 00910 calculate_forces=.TRUE., E_overlap_core=E_overlap_core,& 00911 error=error) 00912 00913 DO iatom = 1, natom 00914 ikind = kind_of(iatom) 00915 i = atom_of_kind(iatom) 00916 harris_force%f_ovrl(iatom, 1:3) = force(ikind)%core_overlap(1:3,i) 00917 END DO 00918 00919 CALL deallocate_qs_force(force) 00920 CALL set_qs_env(qs_env=qs_env, force=force_tmp, error=error) 00921 00922 CALL get_qs_env(qs_env=qs_env, force=force, error=error) 00923 harris_force%f_total(:) = 0.0_dp 00924 00925 DO ikind = 1,SIZE(force) 00926 CALL mp_sum(force(ikind)%core_overlap, para_env%group) 00927 END DO 00928 00929 DO iatom = 1, natom 00930 ikind = kind_of(iatom) 00931 i = atom_of_kind(iatom) 00932 00933 !harris_force%f_ovrl(iatom, 1:3) = force(ikind)%core_overlap(1:3,i) 00934 harris_force%f_self(iatom, 1:3) = 0.0_dp 00935 harris_force%f_EII(iatom, 1:3) = harris_force%f_ovrl(iatom, 1:3) & 00936 + harris_force%f_self(iatom, 1:3) & 00937 + harris_force%f_rho_core(iatom, 1:3) 00938 harris_force%f_harris(iatom, 1:3) = harris_force%f_trace(iatom, 1:3) & 00939 - harris_force%f_integral_vxc(iatom, 1:3) & 00940 + harris_force%f_EII(iatom, 1:3) 00941 harris_force%f_total(1:3) = harris_force%f_total(1:3) & 00942 + harris_force%f_harris(iatom, 1:3) 00943 END DO 00944 00945 ! Output 00946 unit_nr=cp_logger_get_default_io_unit(logger) 00947 IF (unit_nr>0) THEN 00948 WRITE (unit_nr,*) ""; WRITE (unit_nr, *) "" 00949 WRITE (unit_nr,*) "The Harris functional force correction is performed!" 00950 WRITE (unit_nr,*) "" 00951 00952 WRITE (unit_nr, *) "F_Harris" 00953 DO iatom = 1, natom 00954 WRITE (unit_nr,*) harris_force%f_harris(iatom, 1:3) 00955 END DO 00956 WRITE (unit_nr, *) "F_Total" 00957 WRITE (unit_nr, *) harris_force%f_total(1:3) 00958 END IF 00959 00960 DEALLOCATE (atom_of_kind,STAT=stat) 00961 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00962 00963 DEALLOCATE (kind_of,STAT=stat) 00964 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00965 00966 DEALLOCATE (natom_of_kind,STAT=stat) 00967 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 00968 00969 END IF 00970 00971 CALL timestop(handle) 00972 00973 END SUBROUTINE harris_force_EVal 00974 00975 ! ***************************************************************************** 00984 SUBROUTINE harris_calc_nsc_force(qs_env, error) 00985 TYPE(qs_environment_type), POINTER :: qs_env 00986 TYPE(cp_error_type), INTENT(INOUT) :: error 00987 00988 CHARACTER(len=*), PARAMETER :: routineN = 'harris_calc_nsc_force', 00989 routineP = moduleN//':'//routineN 00990 00991 INTEGER :: handle, i, iatom, ikind, 00992 ispin, natom, nkind, nspins, 00993 stat 00994 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 00995 natom_of_kind 00996 LOGICAL :: failure 00997 REAL(KIND=dp) :: Ehartree, Exc 00998 TYPE(atomic_kind_type), DIMENSION(:), 00999 POINTER :: atomic_kind_set 01000 TYPE(cell_type), POINTER :: cell 01001 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01002 POINTER :: matrix_ks, rho_ao 01003 TYPE(cp_para_env_type), POINTER :: para_env 01004 TYPE(dft_control_type), POINTER :: dft_control 01005 TYPE(harris_env_type), POINTER :: harris_env 01006 TYPE(harris_force_type), POINTER :: harris_force 01007 TYPE(particle_type), DIMENSION(:), 01008 POINTER :: particle_set 01009 TYPE(pw_env_type), POINTER :: pw_env 01010 TYPE(pw_p_type) :: v_hartree_gspace, 01011 v_hartree_rspace 01012 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_rho_r, rho_r, 01013 v_rspace_new, v_tau_rspace, 01014 v_xc 01015 TYPE(pw_poisson_type), POINTER :: poisson_env 01016 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01017 TYPE(qs_force_type), DIMENSION(:), 01018 POINTER :: force, force_tmp 01019 TYPE(section_vals_type), POINTER :: input, xc_section 01020 TYPE(xc_derivative_set_type), POINTER :: deriv_set 01021 TYPE(xc_rho_set_type), POINTER :: rho_set 01022 01023 ! ------------------------------------------------------------------------ 01024 01025 CALL timeset(routineN,handle) 01026 01027 failure = .FALSE. 01028 para_env=>qs_env%para_env 01029 NULLIFY(harris_env, harris_force, force, force_tmp, particle_set, & 01030 v_rspace_new, v_tau_rspace, rho_r, my_rho_r, v_xc, dft_control, & 01031 cell, deriv_set, rho_set, auxbas_pw_pool, matrix_ks, pw_env, & 01032 poisson_env, xc_section, input) 01033 01034 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 01035 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 01036 01037 IF (.NOT. failure) THEN 01038 CALL get_qs_env(qs_env=qs_env, harris_env=harris_env, dft_control=dft_control, & 01039 atomic_kind_set=atomic_kind_set, force=force, pw_env=pw_env, & 01040 particle_set=particle_set, matrix_ks=matrix_ks, cell=cell, & 01041 error=error) 01042 01043 xc_section => section_vals_get_subs_vals(input, "DFT%XC", error=error) 01044 01045 CALL harris_env_get(harris_env=harris_env, harris_force=harris_force, error=error) 01046 01047 natom = SIZE(particle_set) 01048 nkind = SIZE(force) 01049 nspins = dft_control%nspins 01050 01051 ALLOCATE (atom_of_kind(natom),STAT=stat) 01052 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01053 01054 ALLOCATE (kind_of(natom),STAT=stat) 01055 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01056 01057 ALLOCATE (natom_of_kind(nkind),STAT=stat) 01058 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01059 01060 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 01061 atom_of_kind=atom_of_kind, kind_of=kind_of, & 01062 natom_of_kind=natom_of_kind) 01063 01064 CALL diff_rho_type(rho_input1=qs_env%rho, rho_input2=harris_env%rho, & 01065 rho_output=harris_env%rho_diff, qs_env=qs_env, error=error) 01066 01067 CALL duplicate_qs_force(qs_force_input=force, qs_force_output=force_tmp, & 01068 natom_of_kind=natom_of_kind) 01069 CALL zero_qs_force(force) 01070 01071 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, & 01072 poisson_env=poisson_env, error=error) 01073 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 01074 use_data = COMPLEXDATA1D, & 01075 in_space = RECIPROCALSPACE, error=error) 01076 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace%pw, & 01077 use_data=REALDATA3D, in_space=REALSPACE,error=error) 01078 01079 ! *** Calculation of the forces due to the Hartree energy *** ! 01080 DO ispin = 1,nspins 01081 CALL pw_poisson_solve(poisson_env, harris_env%rho_diff%rho_g(ispin)%pw, & 01082 ehartree=Ehartree, vhartree=v_hartree_gspace%pw,error=error) 01083 END DO 01084 01085 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 01086 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 01087 01088 rho_ao => harris_env%rho_diff%rho_ao 01089 DO ispin = 1,nspins 01090 CALL integrate_v_rspace(v_rspace=v_hartree_rspace, & 01091 p=rho_ao(ispin), & 01092 h=matrix_ks(ispin), qs_env=qs_env, & 01093 calculate_forces=.TRUE., gapw=.FALSE., & 01094 error=error) 01095 END DO 01096 01097 DO ikind = 1,SIZE(force) 01098 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 01099 END DO 01100 01101 DO iatom = 1, natom 01102 ikind = kind_of(iatom) 01103 i = atom_of_kind(iatom) 01104 harris_force%f_nsc(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 01105 END DO 01106 01107 CALL zero_qs_force(force) 01108 01109 DO ispin = 1,nspins 01110 CALL pw_poisson_solve(poisson_env, harris_env%rho_diff%rho_g(ispin)%pw, & 01111 ehartree=Ehartree, vhartree=v_hartree_gspace%pw,error=error) 01112 END DO 01113 01114 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 01115 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 01116 01117 01118 rho_ao => harris_env%rho_diff%rho_ao 01119 DO ispin = 1,nspins 01120 CALL integrate_v_rspace(v_rspace=v_hartree_rspace, & 01121 p=rho_ao(ispin), & 01122 h=matrix_ks(ispin), qs_env=qs_env, & 01123 calculate_forces=.TRUE., gapw=.FALSE., & 01124 error=error) 01125 END DO 01126 01127 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw, error=error) 01128 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw, error=error) 01129 01130 DO ikind = 1,SIZE(force) 01131 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 01132 END DO 01133 01134 DO iatom = 1, natom 01135 ikind = kind_of(iatom) 01136 i = atom_of_kind(iatom) 01137 harris_force%f_nsc(iatom, 1:3) = harris_force%f_nsc(iatom, 1:3) & 01138 + force(ikind)%rho_elec(1:3,i) 01139 END DO 01140 01141 CALL zero_qs_force(force) 01142 01143 ! *** Calculation of the forces of the xc-integral *** ! 01144 01145 CALL qs_vxc_create(qs_env=qs_env, rho_struct=harris_env%rho_diff, xc_section=xc_section, & 01146 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=Exc, & 01147 just_energy=.FALSE., error=error) 01148 01149 rho_ao => harris_env%rho_diff%rho_ao 01150 01151 DO ispin = 1,nspins 01152 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 01153 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 01154 p=rho_ao(ispin), & 01155 h=matrix_ks(ispin), qs_env=qs_env, & 01156 calculate_forces=.TRUE., gapw=.FALSE., & 01157 error=error) 01158 END DO 01159 01160 DO ikind = 1,SIZE(force) 01161 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 01162 END DO 01163 01164 DO iatom = 1, natom 01165 ikind = kind_of(iatom) 01166 i = atom_of_kind(iatom) 01167 !harris_force%f_nsc(iatom, 1:3) = force(ikind)%rho_elec(1:3,i) 01168 END DO 01169 01170 CALL zero_qs_force(force) 01171 01172 CALL qs_rho_get(rho_struct=harris_env%rho_diff, rho_r=rho_r, error=error) 01173 !CALL qs_rho_get(rho_struct=harris_env%rho_diff, rho_r=rho_r, error=error) 01174 ALLOCATE(my_rho_r(SIZE(rho_r)), stat=stat) 01175 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01176 DO ispin=1, SIZE(rho_r) 01177 my_rho_r(ispin)%pw => rho_r(ispin)%pw 01178 END DO 01179 CALL xc_prep_2nd_deriv(deriv_set=deriv_set, rho_set=rho_set, & 01180 rho_r=my_rho_r, pw_pool=auxbas_pw_pool, & 01181 xc_section=xc_section, cell=cell, error=error) 01182 DEALLOCATE(my_rho_r, stat=stat) 01183 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01184 01185 ! ** v_xc ** ! 01186 ALLOCATE(v_xc(nspins), stat=stat) 01187 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01188 DO ispin=1, nspins 01189 NULLIFY(v_xc(ispin)%pw) 01190 CALL pw_pool_create_pw(auxbas_pw_pool, v_xc(ispin)%pw, use_data=REALDATA3D, & 01191 in_space=REALSPACE, error=error) 01192 CALL pw_zero(v_xc(ispin)%pw, error=error) 01193 END DO 01194 01195 CALL xc_calc_2nd_deriv(v_xc, deriv_set=deriv_set, rho_set=rho_set, & 01196 rho1_set=rho_set, pw_pool=auxbas_pw_pool, & 01197 xc_section=xc_section, gapw=.FALSE., error=error) 01198 01199 CALL xc_dset_release(deriv_set, error=error) 01200 01201 DO ispin = 1,nspins 01202 v_xc(ispin)%pw%cr3d = v_xc(ispin)%pw%cr3d * v_xc(ispin)%pw%pw_grid%dvol 01203 v_rspace_new(ispin)%pw%cr3d = rho_r(ispin)%pw%cr3d * v_xc(ispin)%pw%cr3d 01204 CALL pw_scale(v_rspace_new(ispin)%pw, v_rspace_new(ispin)%pw%pw_grid%dvol, error=error) 01205 01206 CALL pw_release(v_xc(ispin)%pw, error=error) 01207 END DO 01208 01209 DEALLOCATE(v_xc,stat=stat) 01210 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01211 01212 CALL xc_rho_set_release(rho_set,error=error) 01213 01214 rho_ao => harris_env%rho_diff%rho_ao 01215 01216 DO ispin = 1,nspins 01217 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), & 01218 p=rho_ao(ispin), & 01219 h=matrix_ks(ispin), qs_env=qs_env, & 01220 calculate_forces=.TRUE., gapw=.FALSE., & 01221 error=error) 01222 END DO 01223 01224 IF (ASSOCIATED(v_rspace_new)) THEN 01225 DO ispin = 1,nspins 01226 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_rspace_new(ispin)%pw,error=error) 01227 END DO 01228 DEALLOCATE(v_rspace_new,stat=stat) 01229 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01230 END IF 01231 IF (ASSOCIATED(v_tau_rspace)) THEN 01232 DO ispin = 1,nspins 01233 CALL pw_pool_give_back_pw(auxbas_pw_pool,v_tau_rspace(ispin)%pw,error=error) 01234 END DO 01235 DEALLOCATE(v_tau_rspace,stat=stat) 01236 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01237 END IF 01238 01239 DO ikind = 1,SIZE(force) 01240 CALL mp_sum(force(ikind)%rho_elec, para_env%group) 01241 END DO 01242 01243 DO iatom = 1, natom 01244 ikind = kind_of(iatom) 01245 i = atom_of_kind(iatom) 01246 !harris_force%f_nsc(iatom, 1:3) = harris_force%f_nsc(iatom,1:3) + force(ikind)%rho_elec(1:3,i) 01247 END DO 01248 01249 END IF 01250 01251 CALL timestop(handle) 01252 01253 END SUBROUTINE harris_calc_nsc_force 01254 01255 ! ***************************************************************************** 01269 SUBROUTINE harris_force_test_rho_core(qs_env, harris_force, error) 01270 TYPE(qs_environment_type), POINTER :: qs_env 01271 TYPE(harris_force_type), POINTER :: harris_force 01272 TYPE(cp_error_type), INTENT(INOUT) :: error 01273 01274 CHARACTER(len=*), PARAMETER :: routineN = 'harris_force_test_rho_core', 01275 routineP = moduleN//':'//routineN 01276 01277 INTEGER :: i, iatom, ikind, natom, 01278 nkind, stat, unit_nr 01279 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 01280 natom_of_kind 01281 LOGICAL :: failure 01282 TYPE(atomic_kind_type), DIMENSION(:), 01283 POINTER :: atomic_kind_set 01284 TYPE(cp_logger_type), POINTER :: logger 01285 TYPE(particle_type), DIMENSION(:), 01286 POINTER :: particle_set 01287 TYPE(qs_force_type), DIMENSION(:), 01288 POINTER :: force 01289 01290 ! ------------------------------------------------------------------------ 01291 01292 failure = .FALSE. 01293 logger => cp_error_get_logger(error) 01294 01295 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 01296 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 01297 CPPrecondition(ASSOCIATED(harris_force), cp_failure_level, routineP, error, failure) 01298 CPPrecondition(harris_force%ref_count>0, cp_failure_level, routineP, error, failure) 01299 01300 unit_nr=cp_logger_get_default_io_unit(logger) 01301 IF (unit_nr>0) THEN 01302 WRITE (unit_nr,*) ""; WRITE (unit_nr, *) "" 01303 WRITE (unit_nr,*) "The Harris force correction rho_core test is performed!" 01304 WRITE (unit_nr,*) "" 01305 END IF 01306 01307 IF (.NOT. failure) THEN 01308 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force, & 01309 particle_set=particle_set, error=error) 01310 01311 natom = SIZE(particle_set) 01312 nkind = SIZE(force) 01313 01314 ALLOCATE (atom_of_kind(natom),STAT=stat) 01315 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01316 01317 ALLOCATE (kind_of(natom),STAT=stat) 01318 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01319 01320 ALLOCATE (natom_of_kind(nkind),STAT=stat) 01321 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01322 01323 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 01324 atom_of_kind=atom_of_kind, kind_of=kind_of, & 01325 natom_of_kind=natom_of_kind) 01326 01327 IF (unit_nr>0) THEN 01328 DO iatom = 1, natom 01329 ikind = kind_of(iatom) 01330 i = atom_of_kind(iatom) 01331 WRITE (unit_nr,*) force(ikind)%rho_core(1:3,i) 01332 WRITE (unit_nr,*) harris_force%f_rho_core(iatom,1:3) + & 01333 harris_force%f_cross_integrate_v_core(iatom,1:3) 01334 END DO 01335 END IF 01336 01337 DEALLOCATE (atom_of_kind,STAT=stat) 01338 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01339 01340 DEALLOCATE (kind_of,STAT=stat) 01341 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01342 01343 DEALLOCATE (natom_of_kind,STAT=stat) 01344 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01345 01346 END IF 01347 01348 END SUBROUTINE harris_force_test_rho_core 01349 01350 ! ***************************************************************************** 01364 SUBROUTINE harris_force_test_rho_elec(qs_env, harris_force, error) 01365 TYPE(qs_environment_type), POINTER :: qs_env 01366 TYPE(harris_force_type), POINTER :: harris_force 01367 TYPE(cp_error_type), INTENT(INOUT) :: error 01368 01369 CHARACTER(len=*), PARAMETER :: routineN = 'harris_force_test_rho_elec', 01370 routineP = moduleN//':'//routineN 01371 01372 INTEGER :: i, iatom, ikind, natom, 01373 nkind, stat, unit_nr 01374 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 01375 natom_of_kind 01376 LOGICAL :: failure 01377 TYPE(atomic_kind_type), DIMENSION(:), 01378 POINTER :: atomic_kind_set 01379 TYPE(cp_logger_type), POINTER :: logger 01380 TYPE(particle_type), DIMENSION(:), 01381 POINTER :: particle_set 01382 TYPE(qs_force_type), DIMENSION(:), 01383 POINTER :: force 01384 01385 ! ------------------------------------------------------------------------ 01386 01387 failure = .FALSE. 01388 logger => cp_error_get_logger(error) 01389 01390 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 01391 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 01392 CPPrecondition(ASSOCIATED(harris_force), cp_failure_level, routineP, error, failure) 01393 CPPrecondition(harris_force%ref_count>0, cp_failure_level, routineP, error, failure) 01394 01395 unit_nr=cp_logger_get_default_io_unit(logger) 01396 IF (unit_nr>0) THEN 01397 WRITE (unit_nr,*) ""; WRITE (unit_nr, *) "" 01398 WRITE (unit_nr,*) "The Harris force correction rho_elec test is performed!" 01399 WRITE (unit_nr,*) "" 01400 END IF 01401 01402 IF (.NOT. failure) THEN 01403 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force, & 01404 particle_set=particle_set, error=error) 01405 01406 natom = SIZE(particle_set) 01407 nkind = SIZE(force) 01408 01409 ALLOCATE (atom_of_kind(natom),STAT=stat) 01410 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01411 01412 ALLOCATE (kind_of(natom),STAT=stat) 01413 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01414 01415 ALLOCATE (natom_of_kind(nkind),STAT=stat) 01416 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01417 01418 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 01419 atom_of_kind=atom_of_kind, kind_of=kind_of, & 01420 natom_of_kind=natom_of_kind) 01421 01422 IF (unit_nr>0) THEN 01423 DO iatom = 1, natom 01424 ikind = kind_of(iatom) 01425 i = atom_of_kind(iatom) 01426 WRITE (unit_nr,*) force(ikind)%rho_elec(1:3,i) 01427 WRITE (unit_nr,*) harris_force%f_hartree(iatom,1:3) + & 01428 harris_force%f_xc(iatom,1:3) + & 01429 harris_force%f_cross_integrate_v(iatom,1:3) 01430 END DO 01431 END IF 01432 01433 DEALLOCATE (atom_of_kind,STAT=stat) 01434 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01435 01436 DEALLOCATE (kind_of,STAT=stat) 01437 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01438 01439 DEALLOCATE (natom_of_kind,STAT=stat) 01440 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01441 01442 END IF 01443 01444 END SUBROUTINE harris_force_test_rho_elec 01445 01446 ! ***************************************************************************** 01458 SUBROUTINE harris_force_test_integral_vxc(qs_env, harris_force, error) 01459 TYPE(qs_environment_type), POINTER :: qs_env 01460 TYPE(harris_force_type), POINTER :: harris_force 01461 TYPE(cp_error_type), INTENT(INOUT) :: error 01462 01463 CHARACTER(len=*), PARAMETER :: 01464 routineN = 'harris_force_test_integral_vxc', 01465 routineP = moduleN//':'//routineN 01466 01467 INTEGER :: i, iatom, ikind, ispin, 01468 natom, nkind, nspins, stat, 01469 unit_nr 01470 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, 01471 natom_of_kind 01472 LOGICAL :: failure 01473 REAL(KIND=dp) :: Ehartree, Exc 01474 TYPE(atomic_kind_type), DIMENSION(:), 01475 POINTER :: atomic_kind_set 01476 TYPE(cp_dbcsr_p_type), DIMENSION(:), 01477 POINTER :: matrix_ks, rho_ao 01478 TYPE(cp_logger_type), POINTER :: logger 01479 TYPE(dft_control_type), POINTER :: dft_control 01480 TYPE(particle_type), DIMENSION(:), 01481 POINTER :: particle_set 01482 TYPE(pw_env_type), POINTER :: pw_env 01483 TYPE(pw_p_type) :: rho_tot_gspace, 01484 v_hartree_gspace, 01485 v_hartree_rspace 01486 TYPE(pw_p_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace 01487 TYPE(pw_p_type), POINTER :: rho_core 01488 TYPE(pw_poisson_type), POINTER :: poisson_env 01489 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 01490 TYPE(qs_force_type), DIMENSION(:), 01491 POINTER :: force, force_tmp 01492 TYPE(qs_rho_type), POINTER :: rho 01493 TYPE(section_vals_type), POINTER :: input, xc_section 01494 01495 ! ------------------------------------------------------------------------ 01496 01497 failure = .FALSE. 01498 NULLIFY(rho, v_rspace_new, matrix_ks, auxbas_pw_pool, pw_env, & 01499 rho_core, force_tmp, v_tau_rspace, poisson_env, dft_control, input,& 01500 xc_section) 01501 logger => cp_error_get_logger(error) 01502 01503 CPPrecondition(ASSOCIATED(qs_env), cp_failure_level, routineP, error, failure) 01504 CPPrecondition(qs_env%ref_count>0, cp_failure_level, routineP, error, failure) 01505 CPPrecondition(ASSOCIATED(harris_force), cp_failure_level, routineP, error, failure) 01506 CPPrecondition(harris_force%ref_count>0, cp_failure_level, routineP, error, failure) 01507 01508 unit_nr=cp_logger_get_default_io_unit(logger) 01509 IF (unit_nr>0) THEN 01510 WRITE (unit_nr,*) ""; WRITE (unit_nr, *) "" 01511 WRITE (unit_nr,*) "The Harris force correction integral_vxc test is performed!" 01512 WRITE (unit_nr,*) "" 01513 END IF 01514 01515 IF (.NOT. failure) THEN 01516 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, force=force, & 01517 dft_control=dft_control, input=input, particle_set=particle_set, error=error) 01518 01519 xc_section => section_vals_get_subs_vals(input, "DFT%XC", error=error) 01520 01521 natom = SIZE(particle_set) 01522 nkind = SIZE(force) 01523 nspins = dft_control%nspins 01524 01525 ALLOCATE (atom_of_kind(natom),STAT=stat) 01526 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01527 01528 ALLOCATE (kind_of(natom),STAT=stat) 01529 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01530 01531 ALLOCATE (natom_of_kind(nkind),STAT=stat) 01532 CPPostcondition(stat==0, cp_failure_level, routineP, error, failure) 01533 01534 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 01535 atom_of_kind=atom_of_kind, kind_of=kind_of, & 01536 natom_of_kind=natom_of_kind) 01537 01538 ! *** Debug Code: Calculate the forces of the sum of eigenvalues *** ! 01539 CALL duplicate_qs_force(qs_force_input=force, qs_force_output=force_tmp, & 01540 natom_of_kind=natom_of_kind) 01541 CALL zero_qs_force(force) 01542 01543 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, rho_core=rho_core, & 01544 matrix_ks=matrix_ks, error=error) 01545 01546 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool, & 01547 poisson_env=poisson_env, error=error) 01548 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 01549 use_data = COMPLEXDATA1D,& 01550 in_space = RECIPROCALSPACE, error=error) 01551 CALL pw_pool_create_pw(auxbas_pw_pool, rho_tot_gspace%pw, & 01552 use_data=COMPLEXDATA1D, & 01553 in_space=RECIPROCALSPACE, error=error) 01554 CALL pw_copy(rho_core%pw,rho_tot_gspace%pw, error=error) 01555 DO ispin = 1,nspins 01556 CALL pw_axpy(rho%rho_g(ispin)%pw,rho_tot_gspace%pw,error=error) 01557 ! old code seems to be wrong for spin polarized case 01558 ! CALL pw_add(rho_core%pw, rho%rho_g(ispin)%pw, rho_tot_gspace%pw) 01559 END DO 01560 01561 CALL pw_poisson_solve(poisson_env, rho_tot_gspace%pw, ehartree=Ehartree, & 01562 vhartree=v_hartree_gspace%pw,error=error) 01563 01564 CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw, & 01565 error=error) 01566 01567 CALL pw_pool_create_pw(auxbas_pw_pool, v_hartree_rspace%pw, & 01568 use_data=REALDATA3D, in_space=REALSPACE,error=error) 01569 CALL pw_transfer(v_hartree_gspace%pw, v_hartree_rspace%pw, error=error) 01570 01571 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_gspace%pw, & 01572 error=error) 01573 01574 NULLIFY(v_rspace_new) 01575 CALL qs_vxc_create(qs_env=qs_env, rho_struct=qs_env%rho, xc_section=xc_section, & 01576 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=Exc, & 01577 just_energy = .FALSE., error=error) 01578 01579 ! sum-up: v_hartree_rspace = v_hartree_rspace + v_rspace_new(ispin) 01580 DO ispin = 1,nspins 01581 CALL pw_axpy(v_rspace_new(ispin)%pw, v_hartree_rspace%pw,error=error) 01582 CALL pw_scale(v_hartree_rspace%pw, v_hartree_rspace%pw%pw_grid%dvol, error=error) 01583 END DO 01584 01585 IF (ASSOCIATED(v_rspace_new)) THEN 01586 DO ispin = 1,nspins 01587 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_rspace_new(ispin)%pw, error=error) 01588 END DO 01589 DEALLOCATE(v_rspace_new,stat=stat) 01590 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01591 END IF 01592 IF (ASSOCIATED(v_tau_rspace)) THEN 01593 DO ispin = 1,nspins 01594 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_tau_rspace(ispin)%pw, error=error) 01595 END DO 01596 DEALLOCATE(v_tau_rspace,stat=stat) 01597 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01598 END IF 01599 01600 rho_ao => rho%rho_ao 01601 DO ispin = 1,nspins 01602 CALL integrate_v_rspace(v_rspace=v_hartree_rspace, p=rho_ao(ispin), & 01603 h=matrix_ks(ispin), qs_env=qs_env, & 01604 calculate_forces=.TRUE., gapw=.FALSE.,error=error) 01605 END DO 01606 01607 CALL pw_pool_give_back_pw(auxbas_pw_pool, v_hartree_rspace%pw, error=error) 01608 01609 DO iatom = 1, natom 01610 ikind = kind_of(iatom) 01611 i = atom_of_kind(iatom) 01612 WRITE (unit_nr,*) force(ikind)%rho_elec(1:3,i) 01613 END DO 01614 01615 CALL deallocate_qs_force(force) 01616 CALL set_qs_env(qs_env=qs_env, force=force_tmp,error=error) 01617 CALL get_qs_env(qs_env=qs_env, force=force, error=error) 01618 01619 IF (unit_nr>0) THEN 01620 WRITE (unit_nr, *) "" 01621 DO iatom = 1, natom 01622 ikind = kind_of(iatom) 01623 i = atom_of_kind(iatom) 01624 WRITE (unit_nr, *) force(ikind)%rho_elec(1:3,i) 01625 WRITE (unit_nr, *) harris_force%f_cross_integrate_v(iatom,1:3) & 01626 + harris_force%f_hartree(iatom,1:3) & 01627 + harris_force%f_xc(iatom,1:3) 01628 END DO 01629 WRITE (unit_nr, *) "" 01630 DO iatom = 1, natom 01631 ikind = kind_of(iatom) 01632 i = atom_of_kind(iatom) 01633 WRITE (unit_nr, *) harris_force%f_cross_integrate_v(iatom,1:3) 01634 WRITE (unit_nr, *) harris_force%f_cross_integrate_v_core(iatom,1:3) 01635 WRITE (unit_nr, *) harris_force%f_V(iatom,1:3) 01636 END DO 01637 WRITE (unit_nr, *) "" 01638 DO iatom = 1, natom 01639 ikind = kind_of(iatom) 01640 i = atom_of_kind(iatom) 01641 WRITE (unit_nr,*) harris_force%f_trace(iatom,1:3) 01642 WRITE (unit_nr,*) force(ikind)%kinetic(1:3,i) & 01643 + force(ikind)%rho_elec(1:3,i) & 01644 + harris_force%f_delta_integral_vxc(iatom,1:3) & 01645 + harris_force%f_hartree(iatom,1:3) & 01646 + harris_force%f_cross_integrate_v_core(iatom,1:3) & 01647 + force(ikind)%overlap(1:3,i) & 01648 + force(ikind)%gth_ppl(1:3,i) & 01649 + force(ikind)%gth_nlcc(1:3,i) & 01650 + force(ikind)%gth_ppnl(1:3,i) 01651 END DO 01652 01653 END IF 01654 01655 DEALLOCATE (atom_of_kind,STAT=stat) 01656 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01657 01658 DEALLOCATE (kind_of,STAT=stat) 01659 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01660 01661 DEALLOCATE (natom_of_kind,STAT=stat) 01662 CPPostconditionNoFail(stat==0, cp_warning_level, routineP, error) 01663 01664 END IF 01665 01666 END SUBROUTINE harris_force_test_integral_vxc 01667 01668 END MODULE harris_force
1.7.3