|
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 ! ***************************************************************************** 00010 MODULE atom_xc 00011 USE atom_types, ONLY: GTH_PSEUDO,& 00012 atom_type,& 00013 opmat_type 00014 USE atom_utils, ONLY: atom_core_density,& 00015 atom_density,& 00016 integrate_grid,& 00017 numpot_matrix 00018 USE input_constants, ONLY: xc_none 00019 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00020 section_vals_type,& 00021 section_vals_val_get 00022 USE kinds, ONLY: dp 00023 USE mathconstants, ONLY: fourpi 00024 USE timings, ONLY: timeset,& 00025 timestop 00026 USE xc_atom, ONLY: xc_rho_set_atom_update 00027 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00028 xc_dset_create,& 00029 xc_dset_get_derivative,& 00030 xc_dset_release,& 00031 xc_dset_zero_all 00032 USE xc_derivative_types, ONLY: xc_derivative_get,& 00033 xc_derivative_type 00034 USE xc_derivatives, ONLY: xc_functionals_eval,& 00035 xc_functionals_get_needs 00036 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00037 USE xc_rho_set_types, ONLY: xc_rho_set_create,& 00038 xc_rho_set_release,& 00039 xc_rho_set_type 00040 #include "cp_common_uses.h" 00041 00042 IMPLICIT NONE 00043 00044 PRIVATE 00045 00046 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_xc' 00047 00048 PUBLIC :: calculate_atom_vxc_lda, calculate_atom_vxc_lsd 00049 00050 ! ***************************************************************************** 00051 00052 CONTAINS 00053 00054 ! ***************************************************************************** 00055 SUBROUTINE calculate_atom_vxc_lda(xcmat,atom,xc_section,error) 00056 TYPE(opmat_type), POINTER :: xcmat 00057 TYPE(atom_type), INTENT(INOUT) :: atom 00058 TYPE(section_vals_type), POINTER :: xc_section 00059 TYPE(cp_error_type), INTENT(inout) :: error 00060 00061 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atom_vxc_lda', 00062 routineP = moduleN//':'//routineN 00063 00064 INTEGER :: deriv_order, handle, i, ierr, 00065 l, myfun, n1, n2, n3, nr, 00066 nspins 00067 INTEGER, DIMENSION(2, 3) :: bounds 00068 LOGICAL :: failure = .FALSE., lsd, nlcc 00069 REAL(KIND=dp) :: density_cut, gradient_cut, 00070 tau_cut 00071 REAL(KIND=dp), DIMENSION(:), POINTER :: exc, vxc 00072 REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, rho, tau 00073 REAL(KIND=dp), DIMENSION(:, :, :), 00074 POINTER :: taumat, xcpot 00075 TYPE(section_vals_type), POINTER :: xc_fun_section 00076 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00077 TYPE(xc_derivative_type), POINTER :: deriv 00078 TYPE(xc_rho_cflags_type) :: needs 00079 TYPE(xc_rho_set_type), POINTER :: rho_set 00080 00081 ! ------------------------------------------------------------------------- 00082 00083 CALL timeset(routineN,handle) 00084 00085 nlcc=.FALSE. 00086 IF ( atom%potential%ppot_type==GTH_PSEUDO ) THEN 00087 nlcc = atom%potential%gth_pot%nlcc 00088 END IF 00089 00090 IF ( ASSOCIATED(xc_section) ) THEN 00091 NULLIFY(rho_set) 00092 00093 xc_fun_section => section_vals_get_subs_vals(xc_section,"XC_FUNCTIONAL",error=error) 00094 CALL section_vals_val_get(xc_fun_section,"_SECTION_PARAMETERS_",i_val=myfun,error=error) 00095 00096 IF(myfun == xc_none) THEN 00097 atom%energy%exc = 0._dp 00098 ELSE 00099 CALL section_vals_val_get(xc_section,"DENSITY_CUTOFF",r_val=density_cut,error=error) 00100 CALL section_vals_val_get(xc_section,"GRADIENT_CUTOFF",r_val=gradient_cut,error=error) 00101 CALL section_vals_val_get(xc_section,"TAU_CUTOFF",r_val=tau_cut,error=error) 00102 00103 lsd = .FALSE. 00104 nspins = 1 00105 needs = xc_functionals_get_needs(xc_fun_section,lsd=lsd,add_basic_components=.FALSE.,error=error) 00106 00107 ! Prepare the structures needed to calculate and store the xc derivatives 00108 00109 ! Array dimension: here anly one dimensional arrays are used, 00110 ! i.e. only the first column of deriv_data is read. 00111 ! The other to dimensions are set to size equal 1 00112 nr = atom%basis%grid%nr 00113 bounds(1:2,1:3) = 1 00114 bounds(2,1) = nr 00115 00116 ! create a place where to put the derivatives 00117 NULLIFY(deriv_set) 00118 CALL xc_dset_create(deriv_set, local_bounds=bounds, error=error) 00119 ! create the place where to store the argument for the functionals 00120 CALL xc_rho_set_create(rho_set,bounds,rho_cutoff=density_cut,& 00121 drho_cutoff=gradient_cut,tau_cutoff=tau_cut,error=error) 00122 ! allocate the required 3d arrays where to store rho and drho 00123 CALL xc_rho_set_atom_update(rho_set,needs,nspins,bounds) 00124 00125 NULLIFY(rho,drho,tau) 00126 IF ( needs%rho ) THEN 00127 ALLOCATE(rho(nr,1),STAT=ierr) 00128 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00129 CALL atom_density(rho(:,1),atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00130 IF ( nlcc ) THEN 00131 CALL atom_core_density(rho(:,1),atom%potential%gth_pot,typ="RHO",rr=atom%basis%grid%rad,error=error) 00132 END IF 00133 END IF 00134 IF ( needs%norm_drho ) THEN 00135 ALLOCATE(drho(nr,1),STAT=ierr) 00136 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00137 CALL atom_density(drho(:,1),atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,typ="DER",error=error) 00138 IF ( nlcc ) THEN 00139 CALL atom_core_density(drho(:,1),atom%potential%gth_pot,typ="DER",rr=atom%basis%grid%rad,error=error) 00140 END IF 00141 END IF 00142 IF ( needs%tau ) THEN 00143 ALLOCATE(tau(nr,1),STAT=ierr) 00144 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00145 CALL atom_density(tau(:,1),atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,& 00146 typ="KIN",rr=atom%basis%grid%rad2,error=error) 00147 END IF 00148 00149 CALL fill_rho_set(rho_set,nspins,needs,rho,drho,tau,nr,error=error) 00150 00151 CALL xc_dset_zero_all(deriv_set, error) 00152 00153 deriv_order = 1 00154 CALL xc_functionals_eval(xc_fun_section,lsd=lsd,rho_set=rho_set,deriv_set=deriv_set,& 00155 deriv_order=deriv_order,error=error) 00156 00157 ! Integration to get the matrix elements and energy 00158 deriv => xc_dset_get_derivative(deriv_set,"",allocate_deriv=.FALSE., error=error) 00159 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00160 atom%energy%exc = fourpi*integrate_grid(xcpot(:,1,1),atom%basis%grid) 00161 00162 ! dump grid density and xcpot (xc energy?) 00163 IF (.FALSE.) THEN 00164 OPEN(UNIT=17,FILE="atom.dat") 00165 DO i=1,SIZE(xcpot(:,1,1)) 00166 WRITE(17,*) atom%basis%grid%rad(i),rho(i,1),xcpot(i,1,1) 00167 ENDDO 00168 CLOSE(UNIT=17) 00169 ENDIF 00170 00171 IF ( needs%rho ) THEN 00172 deriv => xc_dset_get_derivative(deriv_set,"(rho)",allocate_deriv=.FALSE.,error=error) 00173 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00174 CALL numpot_matrix(xcmat%op,xcpot(:,1,1),atom%basis,0,error) 00175 DEALLOCATE(rho,STAT=ierr) 00176 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00177 END IF 00178 IF ( needs%norm_drho ) THEN 00179 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",allocate_deriv=.FALSE.,error=error) 00180 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00181 CALL numpot_matrix(xcmat%op,xcpot(:,1,1),atom%basis,1,error) 00182 DEALLOCATE(drho,STAT=ierr) 00183 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00184 END IF 00185 IF ( needs%tau ) THEN 00186 deriv => xc_dset_get_derivative(deriv_set,"(tau)",allocate_deriv=.FALSE.,error=error) 00187 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00188 n1 = SIZE(xcmat%op,1) 00189 n2 = SIZE(xcmat%op,2) 00190 n3 = SIZE(xcmat%op,3) 00191 ALLOCATE(taumat(n1,n2,0:n3-1),STAT=ierr) 00192 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00193 taumat = 0._dp 00194 00195 xcpot(:,1,1) = 0.5_dp * xcpot(:,1,1) 00196 CALL numpot_matrix(xcmat%op,xcpot(:,1,1),atom%basis,2,error) 00197 xcpot(:,1,1) = xcpot(:,1,1)/atom%basis%grid%rad2(:) 00198 CALL numpot_matrix(taumat,xcpot(:,1,1),atom%basis,0,error) 00199 DO l=0,3 00200 xcmat%op(:,:,l) = xcmat%op(:,:,l) + REAL(l*(l+1),dp)*taumat(:,:,l) 00201 END DO 00202 00203 DEALLOCATE(tau,STAT=ierr) 00204 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00205 DEALLOCATE(taumat,STAT=ierr) 00206 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00207 END IF 00208 00209 ! Release the xc structure used to store the xc derivatives 00210 CALL xc_dset_release(deriv_set, error=error) 00211 CALL xc_rho_set_release(rho_set,error=error) 00212 00213 END IF !xc_none 00214 00215 ELSE 00216 00217 ! we don't have an xc_section, use a default setup 00218 nr = atom%basis%grid%nr 00219 ALLOCATE(rho(nr,1),exc(nr),vxc(nr),STAT=ierr) 00220 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00221 00222 CALL atom_density(rho(:,1),atom%orbitals%pmat,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00223 IF ( nlcc ) THEN 00224 CALL atom_core_density(rho(:,1),atom%potential%gth_pot,typ="RHO",rr=atom%basis%grid%rad,error=error) 00225 END IF 00226 CALL lda_pade(rho(:,1),exc,vxc,error) 00227 00228 atom%energy%exc = fourpi*integrate_grid(exc,atom%basis%grid) 00229 CALL numpot_matrix(xcmat%op,vxc,atom%basis,0,error) 00230 00231 DEALLOCATE(rho,exc,vxc,STAT=ierr) 00232 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00233 00234 END IF 00235 00236 CALL timestop(handle) 00237 00238 END SUBROUTINE calculate_atom_vxc_lda 00239 00240 ! ***************************************************************************** 00241 SUBROUTINE calculate_atom_vxc_lsd(xcmata,xcmatb,atom,xc_section,error) 00242 TYPE(opmat_type), POINTER :: xcmata, xcmatb 00243 TYPE(atom_type), INTENT(INOUT) :: atom 00244 TYPE(section_vals_type), POINTER :: xc_section 00245 TYPE(cp_error_type), INTENT(inout) :: error 00246 00247 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atom_vxc_lsd', 00248 routineP = moduleN//':'//routineN 00249 00250 INTEGER :: deriv_order, handle, ierr, l, 00251 myfun, n1, n2, n3, nr, nspins 00252 INTEGER, DIMENSION(2, 3) :: bounds 00253 LOGICAL :: failure = .FALSE., lsd 00254 REAL(KIND=dp) :: density_cut, gradient_cut, 00255 tau_cut 00256 REAL(KIND=dp), DIMENSION(:), POINTER :: exc, vxca, vxcb 00257 REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, rho, tau 00258 REAL(KIND=dp), DIMENSION(:, :, :), 00259 POINTER :: taumat, xcpot 00260 TYPE(section_vals_type), POINTER :: xc_fun_section 00261 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00262 TYPE(xc_derivative_type), POINTER :: deriv 00263 TYPE(xc_rho_cflags_type) :: needs 00264 TYPE(xc_rho_set_type), POINTER :: rho_set 00265 00266 ! ------------------------------------------------------------------------- 00267 00268 CALL timeset(routineN,handle) 00269 00270 IF ( ASSOCIATED(xc_section) ) THEN 00271 NULLIFY(rho_set) 00272 00273 xc_fun_section => section_vals_get_subs_vals(xc_section,"XC_FUNCTIONAL",error=error) 00274 CALL section_vals_val_get(xc_fun_section,"_SECTION_PARAMETERS_",i_val=myfun,error=error) 00275 00276 IF(myfun == xc_none) THEN 00277 atom%energy%exc = 0._dp 00278 ELSE 00279 CALL section_vals_val_get(xc_section,"DENSITY_CUTOFF",r_val=density_cut,error=error) 00280 CALL section_vals_val_get(xc_section,"GRADIENT_CUTOFF",r_val=gradient_cut,error=error) 00281 CALL section_vals_val_get(xc_section,"TAU_CUTOFF",r_val=tau_cut,error=error) 00282 00283 lsd = .TRUE. 00284 nspins = 2 00285 needs = xc_functionals_get_needs(xc_fun_section,lsd=lsd,add_basic_components=.FALSE.,error=error) 00286 00287 ! Prepare the structures needed to calculate and store the xc derivatives 00288 00289 ! Array dimension: here anly one dimensional arrays are used, 00290 ! i.e. only the first column of deriv_data is read. 00291 ! The other to dimensions are set to size equal 1 00292 nr = atom%basis%grid%nr 00293 bounds(1:2,1:3) = 1 00294 bounds(2,1) = nr 00295 00296 ! create a place where to put the derivatives 00297 NULLIFY(deriv_set) 00298 CALL xc_dset_create(deriv_set, local_bounds=bounds, error=error) 00299 ! create the place where to store the argument for the functionals 00300 CALL xc_rho_set_create(rho_set,bounds,rho_cutoff=density_cut,& 00301 drho_cutoff=gradient_cut,tau_cutoff=tau_cut,error=error) 00302 ! allocate the required 3d arrays where to store rho and drho 00303 CALL xc_rho_set_atom_update(rho_set,needs,nspins,bounds) 00304 00305 NULLIFY(rho,drho,tau) 00306 IF ( needs%rho_spin ) THEN 00307 ALLOCATE(rho(nr,2),STAT=ierr) 00308 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00309 CALL atom_density(rho(:,1),atom%orbitals%pmata,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00310 CALL atom_density(rho(:,2),atom%orbitals%pmatb,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00311 END IF 00312 IF ( needs%norm_drho_spin ) THEN 00313 ALLOCATE(drho(nr,2),STAT=ierr) 00314 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00315 CALL atom_density(drho(:,1),atom%orbitals%pmata,atom%basis,atom%state%maxl_occ,typ="DER",error=error) 00316 CALL atom_density(drho(:,2),atom%orbitals%pmatb,atom%basis,atom%state%maxl_occ,typ="DER",error=error) 00317 END IF 00318 IF ( needs%tau_spin ) THEN 00319 ALLOCATE(tau(nr,2),STAT=ierr) 00320 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00321 CALL atom_density(tau(:,1),atom%orbitals%pmata,atom%basis,atom%state%maxl_occ,& 00322 typ="KIN",rr=atom%basis%grid%rad2,error=error) 00323 CALL atom_density(tau(:,2),atom%orbitals%pmatb,atom%basis,atom%state%maxl_occ,& 00324 typ="KIN",rr=atom%basis%grid%rad2,error=error) 00325 END IF 00326 00327 CALL fill_rho_set(rho_set,nspins,needs,rho,drho,tau,nr,error=error) 00328 00329 CALL xc_dset_zero_all(deriv_set, error) 00330 00331 deriv_order = 1 00332 CALL xc_functionals_eval(xc_fun_section,lsd=lsd,rho_set=rho_set,deriv_set=deriv_set,& 00333 deriv_order=deriv_order,error=error) 00334 00335 ! Integration to get the matrix elements and energy 00336 deriv => xc_dset_get_derivative(deriv_set,"",allocate_deriv=.FALSE., error=error) 00337 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00338 atom%energy%exc = fourpi*integrate_grid(xcpot(:,1,1),atom%basis%grid) 00339 00340 IF ( needs%rho_spin ) THEN 00341 deriv => xc_dset_get_derivative(deriv_set,"(rhoa)",allocate_deriv=.FALSE.,error=error) 00342 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00343 CALL numpot_matrix(xcmata%op,xcpot(:,1,1),atom%basis,0,error) 00344 deriv => xc_dset_get_derivative(deriv_set,"(rhob)",allocate_deriv=.FALSE.,error=error) 00345 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00346 CALL numpot_matrix(xcmatb%op,xcpot(:,1,1),atom%basis,0,error) 00347 DEALLOCATE(rho,STAT=ierr) 00348 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00349 END IF 00350 IF ( needs%norm_drho_spin ) THEN 00351 ! drhoa 00352 NULLIFY(deriv) 00353 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhoa)",allocate_deriv=.FALSE.,error=error) 00354 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00355 CALL numpot_matrix(xcmata%op,xcpot(:,1,1),atom%basis,1,error) 00356 ! drhob 00357 NULLIFY(deriv) 00358 deriv => xc_dset_get_derivative(deriv_set,"(norm_drhob)",allocate_deriv=.FALSE.,error=error) 00359 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00360 CALL numpot_matrix(xcmatb%op,xcpot(:,1,1),atom%basis,1,error) 00361 ! Cross Terms 00362 NULLIFY(deriv) 00363 deriv => xc_dset_get_derivative(deriv_set,"(norm_drho)",error=error) 00364 IF(ASSOCIATED(deriv)) THEN 00365 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00366 CALL numpot_matrix(xcmata%op,xcpot(:,1,1),atom%basis,1,error) 00367 CALL numpot_matrix(xcmatb%op,xcpot(:,1,1),atom%basis,1,error) 00368 END IF 00369 DEALLOCATE(drho,STAT=ierr) 00370 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00371 END IF 00372 IF ( needs%tau_spin ) THEN 00373 n1 = SIZE(xcmata%op,1) 00374 n2 = SIZE(xcmata%op,2) 00375 n3 = SIZE(xcmata%op,3) 00376 ALLOCATE(taumat(n1,n2,0:n3-1),STAT=ierr) 00377 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00378 00379 deriv => xc_dset_get_derivative(deriv_set,"(tau_a)",allocate_deriv=.FALSE.,error=error) 00380 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00381 taumat = 0._dp 00382 xcpot(:,1,1) = 0.5_dp * xcpot(:,1,1) 00383 CALL numpot_matrix(xcmata%op,xcpot(:,1,1),atom%basis,2,error) 00384 xcpot(:,1,1) = xcpot(:,1,1)/atom%basis%grid%rad2(:) 00385 CALL numpot_matrix(taumat,xcpot(:,1,1),atom%basis,0,error) 00386 DO l=0,3 00387 xcmata%op(:,:,l) = xcmata%op(:,:,l) + REAL(l*(l+1),dp)*taumat(:,:,l) 00388 END DO 00389 00390 deriv => xc_dset_get_derivative(deriv_set,"(tau_b)",allocate_deriv=.FALSE.,error=error) 00391 CALL xc_derivative_get(deriv,deriv_data=xcpot,error=error) 00392 taumat = 0._dp 00393 xcpot(:,1,1) = 0.5_dp * xcpot(:,1,1) 00394 CALL numpot_matrix(xcmatb%op,xcpot(:,1,1),atom%basis,2,error) 00395 xcpot(:,1,1) = xcpot(:,1,1)/atom%basis%grid%rad2(:) 00396 CALL numpot_matrix(taumat,xcpot(:,1,1),atom%basis,0,error) 00397 DO l=0,3 00398 xcmatb%op(:,:,l) = xcmatb%op(:,:,l) + REAL(l*(l+1),dp)*taumat(:,:,l) 00399 END DO 00400 00401 DEALLOCATE(tau,STAT=ierr) 00402 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00403 DEALLOCATE(taumat,STAT=ierr) 00404 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00405 END IF 00406 00407 ! Release the xc structure used to store the xc derivatives 00408 CALL xc_dset_release(deriv_set, error=error) 00409 CALL xc_rho_set_release(rho_set,error=error) 00410 00411 END IF !xc_none 00412 00413 ELSE 00414 00415 ! we don't have an xc_section, use a default setup 00416 nr = atom%basis%grid%nr 00417 ALLOCATE(rho(nr,2),exc(nr),vxca(nr),vxcb(nr),STAT=ierr) 00418 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00419 00420 CALL atom_density(rho(:,1),atom%orbitals%pmata,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00421 CALL atom_density(rho(:,2),atom%orbitals%pmatb,atom%basis,atom%state%maxl_occ,typ="RHO",error=error) 00422 CALL lsd_pade(rho(:,1),rho(:,2),exc,vxca,vxcb,error) 00423 00424 atom%energy%exc = fourpi*integrate_grid(exc,atom%basis%grid) 00425 CALL numpot_matrix(xcmata%op,vxca,atom%basis,0,error) 00426 CALL numpot_matrix(xcmatb%op,vxcb,atom%basis,0,error) 00427 00428 DEALLOCATE(rho,exc,vxca,vxcb,STAT=ierr) 00429 CPPostcondition(ierr==0, cp_failure_level, routineP, error, failure) 00430 00431 END IF 00432 00433 CALL timestop(handle) 00434 00435 END SUBROUTINE calculate_atom_vxc_lsd 00436 00437 ! ***************************************************************************** 00438 SUBROUTINE fill_rho_set(rho_set,nspins,needs,rho,drho,tau,na,error) 00439 00440 TYPE(xc_rho_set_type), POINTER :: rho_set 00441 INTEGER, INTENT(IN) :: nspins 00442 TYPE(xc_rho_cflags_type), INTENT(in) :: needs 00443 REAL(dp), DIMENSION(:, :), POINTER :: rho, drho, tau 00444 INTEGER, INTENT(IN) :: na 00445 TYPE(cp_error_type), INTENT(inout) :: error 00446 00447 CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_rho_set', 00448 routineP = moduleN//':'//routineN 00449 REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp) 00450 00451 INTEGER :: ia 00452 LOGICAL :: failure 00453 00454 failure = .FALSE. 00455 00456 SELECT CASE(nspins) 00457 CASE(1) 00458 CPPrecondition(.NOT.needs%rho_spin,cp_failure_level,routineP,error,failure) 00459 CPPrecondition(.NOT.needs%drho_spin,cp_failure_level,routineP,error,failure) 00460 CPPrecondition(.NOT.needs%norm_drho_spin,cp_failure_level,routineP,error,failure) 00461 CPPrecondition(.NOT.needs%drhoa_drhob,cp_failure_level,routineP,error,failure) 00462 CPPrecondition(.NOT.needs%rho_spin_1_3,cp_failure_level,routineP,error,failure) 00463 CPPrecondition(.NOT.needs%tau_spin,cp_failure_level,routineP,error,failure) 00464 CPPrecondition(.NOT.needs%drho,cp_failure_level,routineP,error,failure) 00465 ! Give rho to 1/3 00466 IF (needs%rho_1_3) THEN 00467 DO ia = 1,na 00468 rho_set%rho_1_3(ia,1,1) = MAX(rho(ia,1),0.0_dp)**f13 00469 END DO 00470 rho_set%owns%rho_1_3=.TRUE. 00471 rho_set%has%rho_1_3=.TRUE. 00472 END IF 00473 ! Give the density 00474 IF (needs%rho) THEN 00475 DO ia = 1,na 00476 rho_set%rho(ia,1,1) = rho(ia,1) 00477 END DO 00478 rho_set%owns%rho=.TRUE. 00479 rho_set%has%rho=.TRUE. 00480 END IF 00481 ! Give the norm of the gradient of the density 00482 IF (needs%norm_drho) THEN 00483 DO ia = 1,na 00484 rho_set%norm_drho(ia,1,1) = drho(ia,1) 00485 END DO 00486 rho_set%owns%norm_drho=.TRUE. 00487 rho_set%has%norm_drho=.TRUE. 00488 END IF 00489 CASE(2) 00490 CPPrecondition(.NOT.needs%drho,cp_failure_level,routineP,error,failure) 00491 CPPrecondition(.NOT.needs%drho_spin,cp_failure_level,routineP,error,failure) 00492 CPPrecondition(.NOT.needs%drhoa_drhob,cp_failure_level,routineP,error,failure) 00493 ! Give the total density 00494 IF (needs%rho) THEN 00495 DO ia = 1,na 00496 rho_set%rho(ia,1,1) = rho(ia,1) + rho(ia,2) 00497 END DO 00498 rho_set%owns%rho=.TRUE. 00499 rho_set%has%rho=.TRUE. 00500 END IF 00501 ! Give the norm of the total gradient of the density 00502 IF (needs%norm_drho) THEN 00503 DO ia = 1,na 00504 rho_set%norm_drho(ia,1,1) = drho(ia,1) + drho(ia,2) 00505 END DO 00506 rho_set%owns%norm_drho=.TRUE. 00507 rho_set%has%norm_drho=.TRUE. 00508 END IF 00509 ! Give rho_spin 00510 IF (needs%rho_spin) THEN 00511 DO ia = 1,na 00512 rho_set%rhoa(ia,1,1) = rho(ia,1) 00513 rho_set%rhob(ia,1,1) = rho(ia,2) 00514 END DO 00515 rho_set%owns%rho_spin=.TRUE. 00516 rho_set%has%rho_spin=.TRUE. 00517 END IF 00518 ! Give rho_spin to 1/3 00519 IF (needs%rho_spin_1_3) THEN 00520 DO ia = 1,na 00521 rho_set%rhoa_1_3(ia,1,1) = MAX(rho(ia,1),0.0_dp)**f13 00522 rho_set%rhob_1_3(ia,1,1) = MAX(rho(ia,2),0.0_dp)**f13 00523 END DO 00524 rho_set%owns%rho_1_3=.TRUE. 00525 rho_set%has%rho_1_3=.TRUE. 00526 END IF 00527 ! Give the norm of the gradient of rhoa and of rhob separatedly 00528 IF (needs%norm_drho_spin) THEN 00529 DO ia = 1,na 00530 rho_set%norm_drhoa(ia,1,1) = drho(ia,1) 00531 rho_set%norm_drhob(ia,1,1) = drho(ia,2) 00532 END DO 00533 rho_set%owns%norm_drho_spin=.TRUE. 00534 rho_set%has%norm_drho_spin=.TRUE. 00535 END IF 00536 ! 00537 END SELECT 00538 00539 ! tau part 00540 IF (needs%tau) THEN 00541 IF (nspins==2) THEN 00542 DO ia = 1,na 00543 rho_set%tau(ia,1,1) = tau(ia,1)+tau(ia,2) 00544 END DO 00545 rho_set%owns%tau=.TRUE. 00546 rho_set%has%tau=.TRUE. 00547 ELSE 00548 DO ia = 1,na 00549 rho_set%tau(ia,1,1) = tau(ia,1) 00550 END DO 00551 rho_set%owns%tau=.TRUE. 00552 rho_set%has%tau=.TRUE. 00553 END IF 00554 END IF 00555 IF (needs%tau_spin) THEN 00556 CPPrecondition(nspins==2,cp_failure_level,routineP,error,failure) 00557 DO ia = 1,na 00558 rho_set%tau_a(ia,1,1) = tau(ia,1) 00559 rho_set%tau_b(ia,1,1) = tau(ia,2) 00560 END DO 00561 rho_set%owns%tau_spin=.TRUE. 00562 rho_set%has%tau_spin=.TRUE. 00563 END IF 00564 00565 END SUBROUTINE fill_rho_set 00566 ! ***************************************************************************** 00567 00568 SUBROUTINE lda_pade(rho,exc,vxc,error) 00569 00570 REAL(dp), DIMENSION(:) :: rho, exc, vxc 00571 TYPE(cp_error_type), INTENT(inout) :: error 00572 00573 CHARACTER(LEN=*), PARAMETER :: routineN = 'lda_pade', 00574 routineP = moduleN//':'//routineN 00575 REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429E+0_dp, 00576 a1 = 0.2217058676663745E+1_dp, a2 = 0.7405551735357053E+0_dp, 00577 a3 = 0.1968227878617998E-1_dp, b1 = 1.0000000000000000E+0_dp, 00578 b2 = 0.4504130959426697E+1_dp, b3 = 0.1110667363742916E+1_dp, 00579 b4 = 0.2359291751427506E-1_dp, f13 = (1.0_dp/3.0_dp), 00580 rsfac = 0.6203504908994000166680065_dp 00581 00582 INTEGER :: i, n 00583 LOGICAL :: failure 00584 REAL(KIND=dp) :: depade, dpv, dq, epade, p, q, 00585 rs 00586 00587 failure = .FALSE. 00588 00589 n = SIZE(rho) 00590 exc(1:n) = 0._dp 00591 vxc(1:n) = 0._dp 00592 00593 DO i=1,n 00594 IF ( rho(i) > 1.e-20_dp ) THEN 00595 rs = rsfac * rho(i)**(-f13) 00596 p = a0 + (a1 + (a2 + a3*rs)*rs)*rs 00597 q = (b1 + (b2 + (b3 + b4*rs)*rs)*rs)*rs 00598 epade = -p/q 00599 00600 dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs)*rs 00601 dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs)*rs)*rs 00602 depade = f13 * rs * (dpv*q - p*dq) / (q*q) 00603 00604 exc(i) = epade * rho(i) 00605 vxc(i) = epade + depade 00606 END IF 00607 END DO 00608 00609 END SUBROUTINE lda_pade 00610 ! ***************************************************************************** 00611 00612 SUBROUTINE lsd_pade ( rhoa, rhob, exc, vxca, vxcb, error ) 00613 00614 REAL(dp), DIMENSION(:) :: rhoa, rhob, exc, vxca, vxcb 00615 TYPE(cp_error_type), INTENT(inout) :: error 00616 00617 CHARACTER(LEN=*), PARAMETER :: routineN = 'lsd_pade', 00618 routineP = moduleN//':'//routineN 00619 REAL(KIND=dp), PARAMETER :: a0 = 0.4581652932831429E+0_dp, 00620 a1 = 0.2217058676663745E+1_dp, a2 = 0.7405551735357053E+0_dp, 00621 a3 = 0.1968227878617998E-1_dp, b1 = 1.0000000000000000E+0_dp, 00622 b2 = 0.4504130959426697E+1_dp, b3 = 0.1110667363742916E+1_dp, 00623 b4 = 0.2359291751427506E-1_dp, da0 = 0.119086804055547E+0_dp, 00624 da1 = 0.6157402568883345E+0_dp, da2 = 0.1574201515892867E+0_dp, 00625 da3 = 0.3532336663397157E-2_dp, db1 = 0.0000000000000000E+0_dp, 00626 db2 = 0.2673612973836267E+0_dp, db3 = 0.2052004607777787E+0_dp, 00627 db4 = 0.4200005045691381E-2_dp, f13 = (1.0_dp/3.0_dp), 00628 f43 = (4.0_dp/3.0_dp), fxfac = 1.923661050931536319759455_dp , 00629 rsfac = 0.6203504908994000166680065_dp 00630 00631 INTEGER :: i, n 00632 LOGICAL :: failure 00633 REAL(KIND=dp) :: dc, dpv, dq, dr, dx, fa0, 00634 fa1, fa2, fa3, fb1, fb2, fb3, 00635 fb4, fx1, fx2, p, q, rhoab, 00636 rs, x, xp, xq 00637 00638 ! 1/(2^(4/3) - 2) 00639 00640 failure = .FALSE. 00641 00642 n = SIZE(rhoa) 00643 exc(1:n) = 0._dp 00644 vxca(1:n) = 0._dp 00645 vxcb(1:n) = 0._dp 00646 00647 DO i=1,n 00648 rhoab = rhoa(i) + rhob(i) 00649 IF ( rhoab > 1.e-20_dp ) THEN 00650 rs = rsfac * rhoab**(-f13) 00651 00652 x = (rhoa(i) - rhob(i)) / rhoab 00653 IF ( x < -1.0_dp ) THEN 00654 fx1 = 1.0_dp 00655 fx2 = -f43*fxfac*2.0_dp**f13 00656 ELSE IF ( x > 1.0_dp ) THEN 00657 fx1 = 1.0_dp 00658 fx2 = f43*fxfac*2.0_dp**f13 00659 ELSE 00660 fx1 = ( (1.0_dp+x)**f43 + (1.0_dp-x)**f43 - 2.0_dp ) * fxfac 00661 fx2 = ( (1.0_dp+x)**f13 - (1.0_dp-x)**f13 ) * fxfac * f43 00662 END IF 00663 00664 fa0 = a0 + fx1*da0 00665 fa1 = a1 + fx1*da1 00666 fa2 = a2 + fx1*da2 00667 fa3 = a3 + fx1*da3 00668 fb1 = b1 + fx1*db1 00669 fb2 = b2 + fx1*db2 00670 fb3 = b3 + fx1*db3 00671 fb4 = b4 + fx1*db4 00672 00673 p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs 00674 q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs 00675 dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs 00676 dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + & 00677 4.0_dp*fb4*rs)*rs)*rs 00678 xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs 00679 xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs 00680 00681 dr = (dpv*q - p*dq)/(q*q) 00682 dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx2/rhoab 00683 dc = f13*rs*dr - p/q 00684 00685 exc(i) = -p/q*rhoab 00686 vxca(i) = dc - dx*rhob(i) 00687 vxcb(i) = dc + dx*rhoa(i) 00688 END IF 00689 END DO 00690 00691 END SUBROUTINE lsd_pade 00692 00693 ! ***************************************************************************** 00694 00695 END MODULE atom_xc
1.7.3