|
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 ! ***************************************************************************** 00014 MODULE qs_vxc 00015 00016 USE cell_types, ONLY: cell_type 00017 USE cp_control_types, ONLY: dft_control_type 00018 USE input_constants, ONLY: sic_ad,& 00019 sic_eo,& 00020 sic_mauri_spz,& 00021 sic_mauri_us,& 00022 sic_none,& 00023 xc_none,& 00024 xc_vdw_fun_nonloc 00025 USE input_section_types, ONLY: section_vals_type,& 00026 section_vals_val_get 00027 USE kinds, ONLY: dp 00028 USE pw_env_types, ONLY: pw_env_get,& 00029 pw_env_type 00030 USE pw_grids, ONLY: pw_grid_compare 00031 USE pw_methods, ONLY: pw_axpy,& 00032 pw_copy,& 00033 pw_scale,& 00034 pw_transfer,& 00035 pw_zero 00036 USE pw_pool_types, ONLY: pw_pool_create_pw,& 00037 pw_pool_give_back_pw,& 00038 pw_pool_type 00039 USE pw_types, ONLY: COMPLEXDATA1D,& 00040 REALDATA3D,& 00041 REALSPACE,& 00042 RECIPROCALSPACE,& 00043 pw_p_type,& 00044 pw_release,& 00045 pw_type 00046 USE qs_dispersion_nonloc, ONLY: calculate_dispersion_nonloc 00047 USE qs_environment_types, ONLY: get_qs_env,& 00048 qs_environment_type 00049 USE qs_mo_types, ONLY: get_mo_set,& 00050 mo_set_p_type 00051 USE qs_rho_types, ONLY: qs_rho_type 00052 USE timings, ONLY: timeset,& 00053 timestop 00054 USE virial_types, ONLY: virial_type 00055 USE xc, ONLY: xc_exc_calc,& 00056 xc_vxc_pw_create1 00057 #include "cp_common_uses.h" 00058 00059 IMPLICIT NONE 00060 00061 PRIVATE 00062 00063 ! *** Public subroutines *** 00064 PUBLIC :: qs_vxc_create 00065 00066 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc' 00067 00068 CONTAINS 00069 00070 ! ***************************************************************************** 00095 SUBROUTINE qs_vxc_create( qs_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc,& 00096 just_energy, adiabatic_rescale_factor, error) 00097 00098 TYPE(qs_environment_type), POINTER :: qs_env 00099 TYPE(qs_rho_type), POINTER :: rho_struct 00100 TYPE(section_vals_type), POINTER :: xc_section 00101 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau 00102 REAL(KIND=dp) :: exc 00103 LOGICAL, INTENT(in), OPTIONAL :: just_energy 00104 REAL(KIND=dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor 00105 TYPE(cp_error_type), INTENT(inout) :: error 00106 00107 CHARACTER(len=*), PARAMETER :: routineN = 'qs_vxc_create', 00108 routineP = moduleN//':'//routineN 00109 00110 INTEGER :: handle, ispin, myfun, stat, 00111 vdw 00112 LOGICAL :: do_adiabatic_rescaling, 00113 failure, my_just_energy, 00114 sic_scaling_b_zero, uf_grid, 00115 vdW_nl 00116 REAL(KIND=dp) :: exc_m, factor, 00117 my_adiabatic_rescale_factor, 00118 my_scaling, nelec_s_inv, 00119 nelec_spin(2) 00120 TYPE(cell_type), POINTER :: cell 00121 TYPE(dft_control_type), POINTER :: dft_control 00122 TYPE(mo_set_p_type), DIMENSION(:), 00123 POINTER :: mo_array 00124 TYPE(pw_env_type), POINTER :: pw_env 00125 TYPE(pw_p_type), DIMENSION(:), POINTER :: my_vxc_rho, my_vxc_tau, 00126 rho_g, rho_m_gspace, 00127 rho_m_rspace, rho_r, tau 00128 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, vdw_pw_pool, 00129 xc_pw_pool 00130 TYPE(pw_type), POINTER :: tmp_g, tmp_g2, tmp_pw 00131 TYPE(virial_type), POINTER :: virial 00132 00133 CALL timeset( routineN ,handle) 00134 00135 failure=.FALSE. 00136 CPPrecondition(.NOT.ASSOCIATED(vxc_rho),cp_failure_level,routineP,error,failure) 00137 CPPrecondition(.NOT.ASSOCIATED(vxc_tau),cp_failure_level,routineP,error,failure) 00138 NULLIFY(dft_control, pw_env, auxbas_pw_pool, xc_pw_pool, vdw_pw_pool, cell, my_vxc_rho, & 00139 tmp_pw, tmp_g, tmp_g2, my_vxc_tau, rho_g, rho_r, tau, rho_m_rspace, & 00140 rho_m_gspace) 00141 00142 my_just_energy=.FALSE. 00143 IF (PRESENT(just_energy)) my_just_energy=just_energy 00144 my_adiabatic_rescale_factor = 1.0_dp 00145 do_adiabatic_rescaling = .FALSE. 00146 IF( PRESENT(adiabatic_rescale_factor)) THEN 00147 my_adiabatic_rescale_factor = adiabatic_rescale_factor 00148 do_adiabatic_rescaling = .TRUE. 00149 END IF 00150 00151 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 00152 pw_env=pw_env, cell=cell, virial=virial, error=error) 00153 00154 CALL section_vals_val_get(xc_section,"XC_FUNCTIONAL%_SECTION_PARAMETERS_",& 00155 i_val=myfun,error=error) 00156 CALL section_vals_val_get(xc_section,"VDW_POTENTIAL%POTENTIAL_TYPE",& 00157 i_val=vdw,error=error) 00158 vdW_nl = (vdw==xc_vdw_fun_nonloc) 00159 00160 IF (myfun/=xc_none .OR. vdW_nl) THEN 00161 00162 ! test if the real space density is available 00163 CPPrecondition(ASSOCIATED(rho_struct),cp_failure_level,routineP,error,failure) 00164 CPPrecondition(rho_struct%ref_count>0,cp_failure_level,routineP,error,failure) 00165 CPPrecondition(rho_struct%rho_r_valid,cp_failure_level,routineP,error,failure) 00166 CALL cp_assert( dft_control%nspins == 1 .OR. dft_control%nspins == 2,& 00167 cp_failure_level,cp_assertion_failed,routineP,& 00168 "nspins must be 1 or 2",error,failure) 00169 ! there are some options related to SIC here. 00170 ! Normal DFT computes E(rho_alpha,rho_beta) (or its variant E(2*rho_alpha) for non-LSD) 00171 ! SIC can E(rho_alpha,rho_beta)-b*(E(rho_alpha,rho_beta)-E(rho_beta,rho_beta)) 00172 ! or compute E(rho_alpha,rho_beta)-b*E(rho_alpha-rho_beta,0) 00173 00174 ! my_scaling is the scaling needed of the standard E(rho_alpha,rho_beta) term 00175 my_scaling=1.0_dp 00176 SELECT CASE (dft_control%sic_method_id) 00177 CASE ( sic_none ) 00178 ! all fine 00179 CASE ( sic_mauri_spz, sic_ad ) 00180 ! no idea yet what to do here in that case 00181 CPPrecondition(.NOT.rho_struct%tau_r_valid ,cp_failure_level,routineP,error,failure) 00182 CASE ( sic_mauri_us ) 00183 my_scaling=1.0_dp-dft_control%sic_scaling_b 00184 ! no idea yet what to do here in that case 00185 CPPrecondition(.NOT.rho_struct%tau_r_valid ,cp_failure_level,routineP,error,failure) 00186 CASE ( sic_eo ) 00187 ! NOTHING TO BE DONE 00188 CASE DEFAULT 00189 ! this case has not yet been treated here 00190 CALL cp_assert(.FALSE., cp_failure_level,cp_assertion_failed,routineP,"NYI",error,failure) 00191 END SELECT 00192 00193 IF (dft_control%sic_scaling_b .EQ. 0.0_dp) THEN 00194 sic_scaling_b_zero = .TRUE. 00195 ELSE 00196 sic_scaling_b_zero = .FALSE. 00197 ENDIF 00198 00199 IF ( .NOT. failure ) THEN 00200 CALL pw_env_get(pw_env,xc_pw_pool=xc_pw_pool,auxbas_pw_pool=auxbas_pw_pool,& 00201 error=error) 00202 uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid,xc_pw_pool%pw_grid) 00203 00204 ALLOCATE(rho_r(dft_control%nspins),stat=stat) 00205 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00206 IF (.not.uf_grid) THEN 00207 DO ispin=1,dft_control%nspins 00208 rho_r(ispin)%pw => rho_struct%rho_r(ispin)%pw 00209 END DO 00210 00211 IF (rho_struct%tau_r_valid) THEN 00212 ALLOCATE(tau(dft_control%nspins),stat=stat) 00213 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00214 DO ispin=1,dft_control%nspins 00215 tau(ispin)%pw => rho_struct%tau_r(ispin)%pw 00216 END DO 00217 END IF 00218 00219 ! for gradient corrected functional the density in g space might 00220 ! be useful so if we have it, we pass it in 00221 IF ( rho_struct%rho_g_valid ) THEN 00222 ALLOCATE(rho_g(dft_control%nspins),stat=stat) 00223 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00224 DO ispin=1,dft_control%nspins 00225 rho_g(ispin)%pw => rho_struct%rho_g(ispin)%pw 00226 END DO 00227 END IF 00228 ELSE 00229 CPPrecondition(rho_struct%rho_g_valid,cp_failure_level,routineP,error,failure) 00230 ALLOCATE(rho_g(dft_control%nspins),stat=stat) 00231 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00232 DO ispin=1,dft_control%nspins 00233 CALL pw_pool_create_pw(xc_pw_pool,rho_g(ispin)%pw,& 00234 in_space=RECIPROCALSPACE, use_data=COMPLEXDATA1D,error=error) 00235 CALL pw_transfer(rho_struct%rho_g(ispin)%pw,rho_g(ispin)%pw, error=error) 00236 END DO 00237 DO ispin=1,dft_control%nspins 00238 CALL pw_pool_create_pw(xc_pw_pool,rho_r(ispin)%pw,& 00239 in_space=REALSPACE, use_data=REALDATA3D,error=error) 00240 CALL pw_transfer(rho_g(ispin)%pw,rho_r(ispin)%pw, error=error) 00241 END DO 00242 IF (rho_struct%tau_r_valid) THEN 00243 ! tau with finer grids is not implemented (at least not correctly), which this asserts 00244 CALL cp_unimplemented_error(fromWhere=routineP, & 00245 message="tau with finer grids", & 00246 error=error, error_level=cp_failure_level) 00247 ! ALLOCATE(tau(dft_control%nspins),stat=stat) 00248 ! CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00249 ! DO ispin=1,dft_control%nspins 00250 ! CALL pw_pool_create_pw(xc_pw_pool,tau(ispin)%pw,& 00251 ! in_space=REALSPACE, use_data=REALDATA3D,error=error) 00252 ! 00253 ! CALL pw_pool_create_pw(xc_pw_pool,tmp_g,& 00254 ! in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D,error=error) 00255 ! CALL pw_pool_create_pw(auxbas_pw_pool,tmp_g2,& 00256 ! in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D,error=error) 00257 ! CALL pw_transfer(tau(ispin)%pw,tmp_g, error=error) 00258 ! CALL pw_transfer(tmp_g,tmp_g2, error=error) 00259 ! CALL pw_transfer(tmp_g2,tmp_pw, error=error) 00260 ! CALL pw_pool_give_back_pw(auxbas_pw_pool,tmp_g2,error=error) 00261 ! CALL pw_pool_give_back_pw(xc_pw_pool,tmp_g,error=error) 00262 ! END DO 00263 END IF 00264 END IF 00265 00266 ! add the nlcc densities 00267 IF (ASSOCIATED(qs_env%rho_nlcc)) THEN 00268 factor=1.0_dp 00269 DO ispin=1,dft_control%nspins 00270 CALL pw_axpy ( qs_env%rho_nlcc%pw, rho_r(ispin)%pw, factor, error) 00271 CALL pw_axpy ( qs_env%rho_nlcc_g%pw, rho_g(ispin)%pw, factor, error) 00272 ENDDO 00273 ENDIF 00274 00275 ! 00276 ! here the rho_r, rho_g, tau is what it should be 00277 ! we get back the right my_vxc_rho and my_vxc_tau as required 00278 ! 00279 IF (my_just_energy) THEN 00280 exc=xc_exc_calc(rho_r=rho_r,tau=tau,& 00281 rho_g=rho_g, xc_section=xc_section,& 00282 cell=cell, pw_pool=xc_pw_pool,& 00283 error=error) 00284 00285 ELSE 00286 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho,vxc_tau=my_vxc_tau, rho_r=rho_r,& 00287 rho_g=rho_g,tau=tau,exc=exc,& 00288 xc_section=xc_section,& 00289 cell=cell, pw_pool=xc_pw_pool,& 00290 virial=virial,& 00291 error=error) 00292 END IF 00293 00294 ! remove the nlcc densities (keep stuff in original state) 00295 IF (ASSOCIATED(qs_env%rho_nlcc)) THEN 00296 factor=-1.0_dp 00297 DO ispin=1,dft_control%nspins 00298 CALL pw_axpy ( qs_env%rho_nlcc%pw, rho_r(ispin)%pw, factor, error) 00299 CALL pw_axpy ( qs_env%rho_nlcc_g%pw, rho_g(ispin)%pw, factor, error) 00300 ENDDO 00301 ENDIF 00302 00303 ! calclulate non-local vdW functional 00304 ! only if this XC_SECTION has it 00305 ! if yes, we use the dispersion_env from qs_env 00306 ! this is dangerous, as it assumes a special connection xc_section -> qs_env 00307 IF (vdW_nl) THEN 00308 ! no spin polarization allowed 00309 CPPrecondition(dft_control%nspins==1,cp_failure_level,routineP,error,failure) 00310 CALL pw_env_get(pw_env,vdw_pw_pool=vdw_pw_pool,error=error) 00311 CALL calculate_dispersion_nonloc(qs_env,my_vxc_rho,rho_r,rho_g,my_just_energy,& 00312 cell,vdw_pw_pool,xc_pw_pool,virial=virial,error=error) 00313 END IF 00314 00315 !! Apply rescaling to the potential if requested 00316 IF(.NOT. my_just_energy) THEN 00317 IF(do_adiabatic_rescaling) THEN 00318 IF( ASSOCIATED(my_vxc_rho)) THEN 00319 DO ispin=1,SIZE(my_vxc_rho) 00320 my_vxc_rho(ispin)%pw%cr3d=my_vxc_rho(ispin)%pw%cr3d*my_adiabatic_rescale_factor 00321 END DO 00322 END IF 00323 END IF 00324 END IF 00325 00326 IF (my_scaling .NE. 1.0_dp) THEN 00327 exc=exc * my_scaling 00328 IF (ASSOCIATED(my_vxc_rho)) THEN 00329 DO ispin=1,SIZE(my_vxc_rho) 00330 my_vxc_rho(ispin)%pw%cr3d=my_vxc_rho(ispin)%pw%cr3d*my_scaling 00331 ENDDO 00332 ENDIF 00333 IF (ASSOCIATED(my_vxc_tau)) THEN 00334 DO ispin=1,SIZE(my_vxc_tau) 00335 my_vxc_tau(ispin)%pw%cr3d=my_vxc_tau(ispin)%pw%cr3d*my_scaling 00336 ENDDO 00337 ENDIF 00338 ENDIF 00339 00340 ! we have pw data for the xc, qs_ks requests coeff structure, here we transfer 00341 ! pw -> coeff 00342 IF (ASSOCIATED(my_vxc_rho)) THEN 00343 ALLOCATE(vxc_rho(dft_control%nspins),stat=stat) 00344 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00345 DO ispin=1,dft_control%nspins 00346 vxc_rho(ispin)%pw => my_vxc_rho(ispin)%pw 00347 END DO 00348 DEALLOCATE(my_vxc_rho,stat=stat) 00349 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00350 END IF 00351 IF (ASSOCIATED(my_vxc_tau)) THEN 00352 ALLOCATE(vxc_tau(dft_control%nspins),stat=stat) 00353 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00354 DO ispin=1,dft_control%nspins 00355 vxc_tau(ispin)%pw => my_vxc_tau(ispin)%pw 00356 END DO 00357 DEALLOCATE(my_vxc_tau,stat=stat) 00358 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00359 END IF 00360 00361 ! compute again the xc but now for Exc(m,o) and the opposite sign 00362 IF (dft_control%sic_method_id .EQ. sic_mauri_spz .AND. .NOT. sic_scaling_b_zero) THEN 00363 ALLOCATE(rho_m_rspace(2),rho_m_gspace(2)) 00364 CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(1)%pw,& 00365 use_data = COMPLEXDATA1D,& 00366 in_space = RECIPROCALSPACE, error=error) 00367 CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(1)%pw,& 00368 use_data = REALDATA3D,& 00369 in_space = REALSPACE, error=error) 00370 CALL pw_copy(rho_struct%rho_r(1)%pw,rho_m_rspace(1)%pw, error=error) 00371 CALL pw_axpy(rho_struct%rho_r(2)%pw,rho_m_rspace(1)%pw,alpha=-1._dp, error=error) 00372 CALL pw_copy(rho_struct%rho_g(1)%pw,rho_m_gspace(1)%pw, error=error) 00373 CALL pw_axpy(rho_struct%rho_g(2)%pw,rho_m_gspace(1)%pw,alpha=-1._dp, error=error) 00374 ! bit sad, these will be just zero... 00375 CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(2)%pw,& 00376 use_data = COMPLEXDATA1D,& 00377 in_space = RECIPROCALSPACE, error=error) 00378 CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(2)%pw,& 00379 use_data = REALDATA3D,& 00380 in_space = REALSPACE, error=error) 00381 CALL pw_zero(rho_m_rspace(2)%pw, error=error) 00382 CALL pw_zero(rho_m_gspace(2)%pw, error=error) 00383 00384 rho_g(1)%pw => rho_m_gspace(1)%pw 00385 rho_g(2)%pw => rho_m_gspace(2)%pw 00386 rho_r(1)%pw => rho_m_rspace(1)%pw 00387 rho_r(2)%pw => rho_m_rspace(2)%pw 00388 00389 IF (my_just_energy) THEN 00390 exc_m=xc_exc_calc(rho_r=rho_r,tau=tau,& 00391 rho_g=rho_g, xc_section=xc_section,& 00392 cell=cell, pw_pool=xc_pw_pool,& 00393 error=error) 00394 ELSE 00395 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho,vxc_tau=my_vxc_tau, rho_r=rho_r,& 00396 rho_g=rho_g,tau=tau,exc=exc_m,& 00397 xc_section=xc_section,& 00398 cell=cell, pw_pool=xc_pw_pool,& 00399 error=error) 00400 END IF 00401 00402 exc = exc - dft_control%sic_scaling_b * exc_m 00403 00404 ! and take care of the potential only vxc_rho is taken into account 00405 IF (.NOT. my_just_energy) THEN 00406 vxc_rho(1)%pw%cr3d=vxc_rho(1)%pw%cr3d-dft_control%sic_scaling_b *& 00407 my_vxc_rho(1)%pw%cr3d 00408 vxc_rho(2)%pw%cr3d=vxc_rho(2)%pw%cr3d+dft_control%sic_scaling_b *& 00409 my_vxc_rho(1)%pw%cr3d ! 1=m 00410 CALL pw_release(my_vxc_rho(1)%pw,error=error) 00411 CALL pw_release(my_vxc_rho(2)%pw,error=error) 00412 DEALLOCATE(my_vxc_rho,stat=stat) 00413 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00414 ENDIF 00415 00416 DO ispin=1,2 00417 CALL pw_pool_give_back_pw(xc_pw_pool,rho_m_rspace(ispin)%pw,& 00418 error=error) 00419 CALL pw_pool_give_back_pw(xc_pw_pool,rho_m_gspace(ispin)%pw,& 00420 error=error) 00421 ENDDO 00422 DEALLOCATE(rho_m_rspace) 00423 DEALLOCATE(rho_m_gspace) 00424 00425 ENDIF 00426 00427 ! now we have - sum_s N_s * Exc(rho_s/N_s,0) 00428 IF ( dft_control%sic_method_id .EQ. sic_ad .AND. .NOT. sic_scaling_b_zero ) THEN 00429 00430 ! find out how many elecs we have 00431 CALL get_qs_env(qs_env,mos=mo_array,error=error) 00432 CALL get_mo_set(mo_set=mo_array(1)%mo_set,n_el_f=nelec_spin(1)) 00433 CALL get_mo_set(mo_set=mo_array(2)%mo_set,n_el_f=nelec_spin(2)) 00434 00435 ALLOCATE(rho_m_rspace(2),rho_m_gspace(2)) 00436 DO ispin=1,2 00437 CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(ispin)%pw,& 00438 use_data = COMPLEXDATA1D,& 00439 in_space = RECIPROCALSPACE, error=error) 00440 CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(ispin)%pw,& 00441 use_data = REALDATA3D,& 00442 in_space = REALSPACE, error=error) 00443 ENDDO 00444 00445 rho_g(1)%pw => rho_m_gspace(1)%pw 00446 rho_g(2)%pw => rho_m_gspace(2)%pw 00447 rho_r(1)%pw => rho_m_rspace(1)%pw 00448 rho_r(2)%pw => rho_m_rspace(2)%pw 00449 00450 DO ispin=1,2 00451 IF (nelec_spin(ispin) .GT. 0.0_dp ) THEN 00452 nelec_s_inv=1.0_dp/nelec_spin(ispin) 00453 ELSE 00454 ! does it matter if there are no electrons with this spin (H) ? 00455 nelec_s_inv=0.0_dp 00456 ENDIF 00457 CALL pw_copy(rho_struct%rho_r(ispin)%pw,rho_m_rspace(1)%pw, error=error) 00458 CALL pw_copy(rho_struct%rho_g(ispin)%pw,rho_m_gspace(1)%pw, error=error) 00459 CALL pw_scale(rho_m_rspace(1)%pw,nelec_s_inv, error=error) 00460 CALL pw_scale(rho_m_gspace(1)%pw,nelec_s_inv, error=error) 00461 CALL pw_zero(rho_m_rspace(2)%pw, error=error) 00462 CALL pw_zero(rho_m_gspace(2)%pw, error=error) 00463 00464 IF (my_just_energy) THEN 00465 exc_m=xc_exc_calc(rho_r=rho_r,tau=tau,& 00466 rho_g=rho_g, xc_section=xc_section,& 00467 cell=cell, pw_pool=xc_pw_pool,& 00468 error=error) 00469 ELSE 00470 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho,vxc_tau=my_vxc_tau, rho_r=rho_r,& 00471 rho_g=rho_g,tau=tau,exc=exc_m,& 00472 xc_section=xc_section,& 00473 cell=cell, pw_pool=xc_pw_pool,& 00474 error=error) 00475 END IF 00476 00477 exc = exc - dft_control%sic_scaling_b * nelec_spin(ispin) * exc_m 00478 00479 ! and take care of the potential only vxc_rho is taken into account 00480 IF (.NOT. my_just_energy) THEN 00481 vxc_rho(ispin)%pw%cr3d=vxc_rho(ispin)%pw%cr3d-dft_control%sic_scaling_b *& 00482 my_vxc_rho(1)%pw%cr3d 00483 CALL pw_release(my_vxc_rho(1)%pw,error=error) 00484 CALL pw_release(my_vxc_rho(2)%pw,error=error) 00485 DEALLOCATE(my_vxc_rho,stat=stat) 00486 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00487 ENDIF 00488 ENDDO 00489 00490 DO ispin=1,2 00491 CALL pw_pool_give_back_pw(xc_pw_pool,rho_m_rspace(ispin)%pw,& 00492 error=error) 00493 CALL pw_pool_give_back_pw(xc_pw_pool,rho_m_gspace(ispin)%pw,& 00494 error=error) 00495 ENDDO 00496 DEALLOCATE(rho_m_rspace) 00497 DEALLOCATE(rho_m_gspace) 00498 00499 ENDIF 00500 00501 ! compute again the xc but now for Exc(n_down,n_down) 00502 IF (dft_control%sic_method_id .EQ. sic_mauri_us .AND. .NOT. sic_scaling_b_zero ) THEN 00503 rho_r(1)%pw => rho_struct%rho_r(2)%pw 00504 rho_r(2)%pw => rho_struct%rho_r(2)%pw 00505 IF ( rho_struct%rho_g_valid ) THEN 00506 rho_g(1)%pw => rho_struct%rho_g(2)%pw 00507 rho_g(2)%pw => rho_struct%rho_g(2)%pw 00508 ENDIF 00509 00510 IF (my_just_energy) THEN 00511 exc_m=xc_exc_calc(rho_r=rho_r,tau=tau,& 00512 rho_g=rho_g, xc_section=xc_section,& 00513 cell=cell, pw_pool=xc_pw_pool,& 00514 error=error) 00515 ELSE 00516 CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho,vxc_tau=my_vxc_tau, rho_r=rho_r,& 00517 rho_g=rho_g,tau=tau,exc=exc_m,& 00518 xc_section=xc_section,& 00519 cell=cell, pw_pool=xc_pw_pool,& 00520 error=error) 00521 END IF 00522 00523 exc = exc + dft_control%sic_scaling_b * exc_m 00524 00525 ! and take care of the potential 00526 IF (.NOT. my_just_energy) THEN 00527 ! both go to minority spin 00528 vxc_rho(2)%pw%cr3d = vxc_rho(2)%pw%cr3d + & 00529 2.0_dp * dft_control%sic_scaling_b * my_vxc_rho(1)%pw%cr3d 00530 CALL pw_release(my_vxc_rho(1)%pw,error=error) 00531 CALL pw_release(my_vxc_rho(2)%pw,error=error) 00532 DEALLOCATE(my_vxc_rho) 00533 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00534 ENDIF 00535 00536 ENDIF 00537 00538 ! 00539 ! cleanups 00540 ! 00541 IF (uf_grid) THEN 00542 DO ispin=1,SIZE(rho_r) 00543 CALL pw_pool_give_back_pw(xc_pw_pool,rho_r(ispin)%pw,error=error) 00544 END DO 00545 IF (ASSOCIATED(vxc_rho)) THEN 00546 DO ispin=1,SIZE(vxc_rho) 00547 CALL pw_pool_create_pw(auxbas_pw_pool,tmp_pw,& 00548 in_space=REALSPACE,use_data=REALDATA3D,error=error) 00549 00550 CALL pw_pool_create_pw(xc_pw_pool,tmp_g,& 00551 in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D,error=error) 00552 CALL pw_pool_create_pw(auxbas_pw_pool,tmp_g2,& 00553 in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D,error=error) 00554 CALL pw_transfer(vxc_rho(ispin)%pw,tmp_g, error=error) 00555 CALL pw_transfer(tmp_g,tmp_g2, error=error) 00556 CALL pw_transfer(tmp_g2,tmp_pw, error=error) 00557 CALL pw_pool_give_back_pw(auxbas_pw_pool,tmp_g2,error=error) 00558 CALL pw_pool_give_back_pw(xc_pw_pool,tmp_g,error=error) 00559 !FM CALL pw_zero(tmp_pw,error=error) 00560 !FM CALL pw_restrict_s3(vxc_rho(ispin)%pw,tmp_pw,& 00561 !FM auxbas_pw_pool,param_section=interp_section,error=error) 00562 CALL pw_pool_give_back_pw(xc_pw_pool,vxc_rho(ispin)%pw,error=error) 00563 vxc_rho(ispin)%pw => tmp_pw 00564 NULLIFY(tmp_pw) 00565 END DO 00566 END IF 00567 IF (ASSOCIATED(vxc_tau)) THEN 00568 DO ispin=1,SIZE(vxc_tau) 00569 CALL pw_pool_create_pw(auxbas_pw_pool,tmp_pw,& 00570 in_space=REALSPACE,use_data=REALDATA3D,error=error) 00571 00572 CALL pw_pool_create_pw(xc_pw_pool,tmp_g,& 00573 in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D,error=error) 00574 CALL pw_pool_create_pw(auxbas_pw_pool,tmp_g2,& 00575 in_space=RECIPROCALSPACE,use_data=COMPLEXDATA1D,error=error) 00576 CALL pw_transfer(vxc_tau(ispin)%pw,tmp_g, error=error) 00577 CALL pw_transfer(tmp_g,tmp_g2, error=error) 00578 CALL pw_transfer(tmp_g2,tmp_pw, error=error) 00579 CALL pw_pool_give_back_pw(auxbas_pw_pool,tmp_g2,error=error) 00580 CALL pw_pool_give_back_pw(xc_pw_pool,tmp_g,error=error) 00581 !FM CALL pw_zero(tmp_pw,error=error) 00582 !FM CALL pw_restrict_s3(vxc_rho(ispin)%pw,tmp_pw,& 00583 !FM auxbas_pw_pool,param_section=interp_section,error=error) 00584 CALL pw_pool_give_back_pw(xc_pw_pool,vxc_tau(ispin)%pw,error=error) 00585 vxc_tau(ispin)%pw => tmp_pw 00586 NULLIFY(tmp_pw) 00587 END DO 00588 END IF 00589 00590 END IF 00591 DEALLOCATE(rho_r,stat=stat) 00592 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00593 IF (ASSOCIATED(rho_g)) THEN 00594 IF (uf_grid) THEN 00595 DO ispin=1,SIZE(rho_g) 00596 CALL pw_pool_give_back_pw(xc_pw_pool,rho_g(ispin)%pw,error=error) 00597 END DO 00598 END IF 00599 DEALLOCATE(rho_g,stat=stat) 00600 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00601 END IF 00602 IF (ASSOCIATED(tau)) THEN 00603 IF (uf_grid) THEN 00604 DO ispin=1,SIZE(tau) 00605 CALL pw_pool_give_back_pw(xc_pw_pool,tau(ispin)%pw,error=error) 00606 END DO 00607 END IF 00608 DEALLOCATE(tau,stat=stat) 00609 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 00610 END IF 00611 00612 END IF 00613 END IF 00614 00615 CALL timestop(handle) 00616 00617 END SUBROUTINE qs_vxc_create 00618 00619 END MODULE qs_vxc
1.7.3