|
CP2K 2.5 (Revision 12981)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00007 MODULE xc_atom 00008 00009 USE cp_linked_list_xc_deriv, ONLY: cp_sll_xc_deriv_next,& 00010 cp_sll_xc_deriv_type 00011 USE f77_blas 00012 USE input_section_types, ONLY: section_vals_get_subs_vals,& 00013 section_vals_type 00014 USE kinds, ONLY: dp 00015 USE pw_pool_types, ONLY: pw_pool_type 00016 USE pw_types, ONLY: pw_p_type 00017 USE timings, ONLY: timeset,& 00018 timestop 00019 USE xc, ONLY: divide_by_norm_drho,& 00020 xc_calc_2nd_deriv 00021 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,& 00022 xc_dset_get_derivative 00023 USE xc_derivative_types, ONLY: xc_derivative_get,& 00024 xc_derivative_type 00025 USE xc_derivatives, ONLY: xc_functionals_eval 00026 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type 00027 USE xc_rho_set_types, ONLY: xc_rho_set_get,& 00028 xc_rho_set_type 00029 #include "cp_common_uses.h" 00030 00031 IMPLICIT NONE 00032 00033 PRIVATE 00034 00035 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_atom' 00036 00037 PUBLIC :: vxc_of_r_new, xc_rho_set_atom_update, xc_2nd_deriv_of_r, fill_rho_set 00038 00039 CONTAINS 00040 00041 ! ***************************************************************************** 00042 SUBROUTINE vxc_of_r_new(xc_fun_section,rho_set,deriv_set,deriv_order,needs,w,& 00043 lsd,na,nr,exc,vxc,vxg,vtau,energy_only,epr_xc,adiabatic_rescale_factor,error) 00044 00045 ! This routine updates rho_set by giving to it the rho and drho that are needed. 00046 ! Since for the local densities rho1_h and rho1_s local grids are used it is not possible 00047 ! to call xc_rho_set_update. 00048 ! As input of this routine one gets rho and drho on a one dimensional grid. 00049 ! The grid is the angular grid corresponding to a given point ir_pnt on the radial grid. 00050 ! The derivatives are calculated on this one dimensional grid, the results are stored in 00051 ! exc, vxc(1:na,ir_pnt,ispin), vxg(1:na,ir_pnt,ispin), vxg_cross(1:na,ir_pnt,ispin) 00052 ! Afterwords the arrays containing the derivatives are put to zero so that the routine 00053 ! can safely be called for the next radial point ir_pnt 00054 00055 TYPE(section_vals_type), POINTER :: xc_fun_section 00056 TYPE(xc_rho_set_type), POINTER :: rho_set 00057 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00058 INTEGER, INTENT(in) :: deriv_order 00059 TYPE(xc_rho_cflags_type), INTENT(in) :: needs 00060 REAL(dp), DIMENSION(:, :), POINTER :: w 00061 LOGICAL, INTENT(IN) :: lsd 00062 INTEGER, INTENT(in) :: na, nr 00063 REAL(dp) :: exc 00064 REAL(dp), DIMENSION(:, :, :), POINTER :: vxc 00065 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg 00066 REAL(dp), DIMENSION(:, :, :), POINTER :: vtau 00067 LOGICAL, INTENT(IN), OPTIONAL :: energy_only, epr_xc 00068 REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor 00069 TYPE(cp_error_type), INTENT(inout) :: error 00070 00071 CHARACTER(LEN=*), PARAMETER :: routineN = 'vxc_of_r_new', 00072 routineP = moduleN//':'//routineN 00073 00074 INTEGER :: handle, ia, idir, ir, 00075 my_deriv_order 00076 LOGICAL :: failure, gradient_f, 00077 my_epr_xc, my_only_energy 00078 REAL(dp) :: my_adiabatic_rescale_factor 00079 REAL(dp), DIMENSION(:, :, :), POINTER :: deriv_data 00080 REAL(KIND=dp) :: drho_cutoff 00081 TYPE(xc_derivative_type), POINTER :: deriv_att 00082 00083 CALL timeset(routineN,handle) 00084 failure = .FALSE. 00085 my_only_energy = .FALSE. 00086 IF(PRESENT(energy_only)) my_only_energy=energy_only 00087 00088 IF(PRESENT(adiabatic_rescale_factor)) THEN 00089 my_adiabatic_rescale_factor = adiabatic_rescale_factor 00090 ELSE 00091 my_adiabatic_rescale_factor = 1.0_dp 00092 END IF 00093 00094 ! needed for the epr routines 00095 my_epr_xc = .FALSE. 00096 IF(PRESENT(epr_xc)) my_epr_xc = epr_xc 00097 my_deriv_order = deriv_order 00098 IF(my_epr_xc) my_deriv_order = 2 00099 00100 gradient_f=(needs%drho_spin.OR.needs%norm_drho_spin.OR.& 00101 needs%drhoa_drhob.OR.needs%drho.OR.needs%norm_drho) 00102 00103 ! Calculate the derivatives 00104 CALL xc_functionals_eval(xc_fun_section, & 00105 lsd=lsd,& 00106 rho_set=rho_set, & 00107 deriv_set=deriv_set,& 00108 deriv_order=my_deriv_order, & 00109 error=error) 00110 00111 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, error=error) 00112 00113 NULLIFY (deriv_data) 00114 00115 IF (my_epr_xc) THEN 00116 ! nabla v_xc (using the vxg arrays) 00117 ! there's no point doing this when lsd = false 00118 IF(lsd) THEN 00119 deriv_att => xc_dset_get_derivative(deriv_set, "(rhoa)(rhoa)",error=error) 00120 IF (ASSOCIATED(deriv_att)) THEN 00121 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00122 DO ir = 1,nr 00123 DO ia = 1,na 00124 DO idir = 1,3 00125 vxg(idir,ia,ir,1) = rho_set%drhoa(idir)%array(ia,ir,1) * & 00126 deriv_data(ia,ir,1) 00127 END DO !idir 00128 END DO !ia 00129 END DO !ir 00130 NULLIFY(deriv_data) 00131 END IF 00132 deriv_att => xc_dset_get_derivative(deriv_set, "(rhob)(rhob)",error=error) 00133 IF (ASSOCIATED(deriv_att)) THEN 00134 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00135 DO ir = 1,nr 00136 DO ia = 1,na 00137 DO idir = 1,3 00138 vxg(idir,ia,ir,2) = rho_set%drhob(idir)%array(ia,ir,1) * & 00139 deriv_data(ia,ir,1) 00140 END DO !idir 00141 END DO !ia 00142 END DO !ir 00143 NULLIFY(deriv_data) 00144 END IF 00145 END IF 00146 ! EXC energy ! is that needed for epr? 00147 deriv_att => xc_dset_get_derivative(deriv_set,"", error=error) 00148 IF (ASSOCIATED(deriv_att)) THEN 00149 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00150 exc = 0.0_dp 00151 DO ir = 1, nr 00152 DO ia = 1, na 00153 exc = exc + deriv_data(ia,ir,1)*w(ia,ir) 00154 END DO 00155 END DO 00156 NULLIFY (deriv_data) 00157 END IF 00158 ELSE 00159 ! EXC energy 00160 deriv_att => xc_dset_get_derivative(deriv_set,"", error=error) 00161 IF (ASSOCIATED(deriv_att)) THEN 00162 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00163 exc = 0.0_dp 00164 DO ir = 1, nr 00165 DO ia = 1, na 00166 exc = exc + deriv_data(ia,ir,1)*w(ia,ir) 00167 END DO 00168 END DO 00169 NULLIFY (deriv_data) 00170 END IF 00171 ! Calculate the potential only if needed 00172 IF(.NOT.my_only_energy) THEN 00173 ! Derivative with respect to the density 00174 IF(lsd) THEN 00175 deriv_att => xc_dset_get_derivative(deriv_set, "(rhoa)",error=error) 00176 IF (ASSOCIATED(deriv_att)) THEN 00177 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00178 vxc(:,:,1) = deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00179 NULLIFY (deriv_data) 00180 END IF 00181 deriv_att => xc_dset_get_derivative(deriv_set, "(rhob)",error=error) 00182 IF (ASSOCIATED(deriv_att)) THEN 00183 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00184 vxc(:,:,2) = deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00185 NULLIFY(deriv_data) 00186 END IF 00187 deriv_att => xc_dset_get_derivative(deriv_set, "(rho)",error=error) 00188 IF(ASSOCIATED(deriv_att)) THEN 00189 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00190 vxc(:,:,1) = vxc(:,:,1) + deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00191 vxc(:,:,2) = vxc(:,:,2) + deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00192 NULLIFY (deriv_data) 00193 END IF 00194 ELSE 00195 deriv_att => xc_dset_get_derivative(deriv_set, "(rho)",error=error) 00196 IF (ASSOCIATED(deriv_att)) THEN 00197 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data, error=error) 00198 vxc(:,:,1) = deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00199 NULLIFY (deriv_data) 00200 END IF 00201 END IF ! lsd 00202 00203 ! Derivatives with respect to the gradient 00204 IF (lsd) THEN 00205 deriv_att => xc_dset_get_derivative(deriv_set, "(norm_drhoa)",error=error) 00206 IF (ASSOCIATED(deriv_att)) THEN 00207 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00208 DO ir = 1, nr 00209 DO ia = 1, na 00210 DO idir=1, 3 00211 IF (rho_set%norm_drhoa(ia,ir,1) > drho_cutoff) THEN 00212 vxg(idir,ia,ir,1) = rho_set%drhoa(idir)%array(ia,ir,1) * & 00213 deriv_data(ia,ir,1)*w(ia,ir) / & 00214 rho_set%norm_drhoa(ia,ir,1) * my_adiabatic_rescale_factor 00215 ELSE 00216 vxg(idir,ia,ir,1) = 0.0_dp 00217 END IF 00218 END DO 00219 END DO 00220 END DO 00221 NULLIFY(deriv_data) 00222 END IF 00223 deriv_att => xc_dset_get_derivative(deriv_set, "(norm_drhob)",error=error) 00224 IF (ASSOCIATED(deriv_att)) THEN 00225 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00226 DO ir = 1, nr 00227 DO ia = 1, na 00228 DO idir=1, 3 00229 IF (rho_set%norm_drhob(ia,ir,1) > drho_cutoff) THEN 00230 vxg(idir,ia,ir,2) = rho_set%drhob(idir)%array(ia,ir,1) * & 00231 deriv_data(ia,ir,1)*w(ia,ir) / & 00232 rho_set%norm_drhob(ia,ir,1) * my_adiabatic_rescale_factor 00233 ELSE 00234 vxg(idir,ia,ir,2) = 0.0_dp 00235 END IF 00236 END DO 00237 END DO 00238 END DO 00239 NULLIFY(deriv_data) 00240 END IF 00241 ! Cross Terms 00242 deriv_att => xc_dset_get_derivative(deriv_set, "(norm_drho)",error=error) 00243 IF(ASSOCIATED(deriv_att)) THEN 00244 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00245 DO ir = 1, nr 00246 DO ia = 1, na 00247 DO idir = 1, 3 00248 IF (rho_set%norm_drho(ia,ir,1) > drho_cutoff) THEN 00249 vxg(idir,ia,ir,1:2) = vxg(idir,ia,ir,1:2) + & 00250 (rho_set%drhoa(idir)%array(ia,ir,1) + rho_set%drhob(idir)%array(ia,ir,1)) * & 00251 deriv_data(ia,ir,1)*w(ia,ir)/rho_set%norm_drho(ia,ir,1) * my_adiabatic_rescale_factor 00252 END IF 00253 END DO 00254 END DO 00255 END DO 00256 NULLIFY (deriv_data) 00257 END IF 00258 ELSE 00259 deriv_att => xc_dset_get_derivative(deriv_set,"(norm_drho)",error=error) 00260 IF (ASSOCIATED(deriv_att)) THEN 00261 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00262 DO ir = 1, nr 00263 DO ia = 1, na 00264 IF (rho_set%norm_drho(ia,ir,1) > drho_cutoff) THEN 00265 DO idir=1, 3 00266 vxg(idir,ia,ir,1) = rho_set%drho(idir)%array(ia,ir,1) * & 00267 deriv_data(ia,ir,1)*w(ia,ir) / & 00268 rho_set%norm_drho(ia,ir,1) * my_adiabatic_rescale_factor 00269 END DO 00270 ELSE 00271 vxg(1:3,ia,ir,1) = 0.0_dp 00272 END IF 00273 END DO 00274 END DO 00275 NULLIFY (deriv_data) 00276 END IF 00277 END IF ! lsd 00278 ! Derivative with respect to tau 00279 IF (lsd) THEN 00280 deriv_att => xc_dset_get_derivative(deriv_set,"(taua)",error=error) 00281 IF (ASSOCIATED(deriv_att)) THEN 00282 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00283 vtau(:,:,1) = deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00284 NULLIFY (deriv_data) 00285 END IF 00286 deriv_att => xc_dset_get_derivative(deriv_set,"(taub)",error=error) 00287 IF (ASSOCIATED(deriv_att)) THEN 00288 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00289 vtau(:,:,2) = deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00290 NULLIFY (deriv_data) 00291 END IF 00292 ELSE 00293 deriv_att => xc_dset_get_derivative(deriv_set,"(tau)",error=error) 00294 IF (ASSOCIATED(deriv_att)) THEN 00295 CALL xc_derivative_get(deriv_att,deriv_data=deriv_data,error=error) 00296 vtau(:,:,1) = deriv_data(:,:,1)*w(:,:) * my_adiabatic_rescale_factor 00297 NULLIFY (deriv_data) 00298 END IF 00299 END IF ! lsd 00300 END IF ! only_energy 00301 END IF ! epr_xc 00302 00303 CALL timestop(handle) 00304 00305 END SUBROUTINE vxc_of_r_new 00306 00307 ! ***************************************************************************** 00308 SUBROUTINE xc_2nd_deriv_of_r(rho_set, rho1_set,xc_section,& 00309 deriv_set, needs, w, exc, vxc, vxg, error) 00310 00311 ! As input of this routine one gets rho and drho on a one dimensional grid. 00312 ! The grid is the angular grid corresponding to a given point ir on the radial grid. 00313 ! The derivatives are calculated on this one dimensional grid, the results are stored in 00314 ! exc, vxc(1:na,ir,ispin), vxg(1:na,ir,ispin), vxg_cross(1:na,ir,ispin) 00315 ! Afterwords the arrays containing the derivatives are put to zero so that the routine 00316 ! can safely be called for the next radial point ir 00317 00318 TYPE(xc_rho_set_type), POINTER :: rho_set, rho1_set 00319 TYPE(section_vals_type), POINTER :: xc_section 00320 TYPE(xc_derivative_set_type), POINTER :: deriv_set 00321 TYPE(xc_rho_cflags_type), INTENT(in) :: needs 00322 REAL(dp), DIMENSION(:, :), POINTER :: w 00323 REAL(dp) :: exc 00324 REAL(dp), DIMENSION(:, :, :), POINTER :: vxc 00325 REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg 00326 TYPE(cp_error_type), INTENT(inout) :: error 00327 00328 CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_2nd_deriv_of_r', 00329 routineP = moduleN//':'//routineN 00330 00331 INTEGER :: handle, ispin, nspins, stat 00332 LOGICAL :: failure, lsd 00333 REAL(dp) :: drho_cutoff 00334 TYPE(cp_sll_xc_deriv_type), POINTER :: pos 00335 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_pw 00336 TYPE(pw_pool_type), POINTER :: pw_pool 00337 TYPE(section_vals_type), POINTER :: xc_fun_section 00338 TYPE(xc_derivative_type), POINTER :: deriv_att 00339 00340 CALL timeset(routineN,handle) 00341 failure = .FALSE. 00342 00343 nspins = SIZE(vxc,3) 00344 lsd = (nspins==2) 00345 IF (ASSOCIATED(rho_set%rhoa)) THEN 00346 lsd = .TRUE. 00347 END IF 00348 CALL xc_rho_set_get(rho_set, drho_cutoff=drho_cutoff, error=error) 00349 xc_fun_section => section_vals_get_subs_vals(xc_section,& 00350 "XC_FUNCTIONAL",error=error) 00351 00352 ! Calculate the derivatives 00353 CALL xc_functionals_eval(xc_fun_section, & 00354 lsd=lsd,& 00355 rho_set=rho_set, & 00356 deriv_set=deriv_set,& 00357 deriv_order=2, & 00358 error=error) 00359 00360 CALL divide_by_norm_drho(deriv_set, rho_set, lsd, error) 00361 00362 ! multiply by -w 00363 pos => deriv_set%derivs 00364 DO WHILE (cp_sll_xc_deriv_next(pos,el_att=deriv_att,error=error)) 00365 !deriv_att%deriv_data(:,:,1) = -w(:,:)*deriv_att%deriv_data(:,:,1) 00366 deriv_att%deriv_data(:,:,1) = w(:,:)*deriv_att%deriv_data(:,:,1) 00367 END DO 00368 00369 NULLIFY(pw_pool) 00370 ALLOCATE(vxc_pw(nspins), stat=stat) 00371 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00372 DO ispin=1, nspins 00373 ALLOCATE(vxc_pw(ispin)%pw, stat=stat) 00374 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00375 vxc_pw(ispin)%pw%cr3d => vxc(:,:,ispin:ispin) 00376 END DO 00377 00378 CALL xc_calc_2nd_deriv(vxc_pw, deriv_set, rho_set, rho1_set, pw_pool, & 00379 xc_section,gapw=.TRUE., vxg=vxg, error=error) 00380 00381 DO ispin=1, nspins 00382 DEALLOCATE(vxc_pw(ispin)%pw, stat=stat) 00383 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00384 END DO 00385 DEALLOCATE(vxc_pw, stat=stat) 00386 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00387 00388 ! zero the derivative data for the next call 00389 pos => deriv_set%derivs 00390 DO WHILE (cp_sll_xc_deriv_next(pos,el_att=deriv_att,error=error)) 00391 deriv_att%deriv_data = 0.0_dp 00392 END DO 00393 00394 CALL timestop(handle) 00395 00396 END SUBROUTINE xc_2nd_deriv_of_r 00397 00398 ! ***************************************************************************** 00399 SUBROUTINE xc_rho_set_atom_update(rho_set,needs,nspins,bo) 00400 00401 ! This routine allocates the storage arrays for rho and drho 00402 ! In calculate_vxc_atom this is called once for each atomic_kind, 00403 ! After the loop over all the atoms of the kind and over all the points 00404 ! of the radial grid for each atom, rho_set is deallocated. 00405 ! Within the same kind, at each new point on the radial grid, the rho_set 00406 ! arrays rho and drho are overwritten. 00407 00408 TYPE(xc_rho_set_type), POINTER :: rho_set 00409 TYPE(xc_rho_cflags_type) :: needs 00410 INTEGER, INTENT(IN) :: nspins 00411 INTEGER, DIMENSION(2, 3), INTENT(IN) :: bo 00412 00413 CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_rho_set_atom_update', 00414 routineP = moduleN//':'//routineN 00415 00416 INTEGER :: idir 00417 00418 SELECT CASE(nspins) 00419 CASE(1) 00420 ! What is this for? 00421 IF (needs%rho_1_3) THEN 00422 NULLIFY(rho_set%rho_1_3) 00423 ALLOCATE(rho_set%rho_1_3(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00424 rho_set%owns%rho_1_3=.TRUE. 00425 rho_set%has%rho_1_3= .FALSE. 00426 END IF 00427 ! Allocate the storage space for the density 00428 IF (needs%rho) THEN 00429 NULLIFY(rho_set%rho) 00430 ALLOCATE(rho_set%rho(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00431 rho_set%owns%rho=.TRUE. 00432 rho_set%has%rho=.FALSE. 00433 END IF 00434 ! Allocate the storage space for the norm of the gradient of the density 00435 IF (needs%norm_drho) THEN 00436 NULLIFY(rho_set%norm_drho) 00437 ALLOCATE(rho_set%norm_drho(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00438 rho_set%owns%norm_drho=.TRUE. 00439 rho_set%has%norm_drho=.FALSE. 00440 END IF 00441 ! Allocate the storage space for the three components of the gradient of the density 00442 IF (needs%drho) THEN 00443 DO idir = 1,3 00444 NULLIFY (rho_set%drho(idir)%array) 00445 ALLOCATE(rho_set%drho(idir)%array(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00446 END DO 00447 rho_set%owns%drho=.TRUE. 00448 rho_set%has%drho=.FALSE. 00449 END IF 00450 CASE(2) 00451 ! Allocate the storage space for the total density 00452 IF (needs%rho) THEN 00453 ! this should never be the case unless you use LDA functionals with LSD 00454 NULLIFY(rho_set%rho) 00455 ALLOCATE(rho_set%rho(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00456 rho_set%owns%rho=.TRUE. 00457 rho_set%has%rho=.FALSE. 00458 END IF 00459 ! What is this for? 00460 IF (needs%rho_1_3) THEN 00461 NULLIFY(rho_set%rho_1_3) 00462 ALLOCATE(rho_set%rho_1_3(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00463 rho_set%owns%rho_1_3=.TRUE. 00464 rho_set%has%rho_1_3=.FALSE. 00465 END IF 00466 ! What is this for? 00467 IF (needs%rho_spin_1_3) THEN 00468 NULLIFY(rho_set%rhoa_1_3,rho_set%rhob_1_3) 00469 ALLOCATE(rho_set%rhoa_1_3(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00470 ALLOCATE(rho_set%rhob_1_3(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00471 rho_set%owns%rho_spin_1_3=.TRUE. 00472 rho_set%has%rho_spin_1_3=.FALSE. 00473 END IF 00474 ! Allocate the storage space for the spin densities rhoa and rhob 00475 IF (needs%rho_spin) THEN 00476 NULLIFY(rho_set%rhoa,rho_set%rhob) 00477 ALLOCATE(rho_set%rhoa(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00478 ALLOCATE(rho_set%rhob(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00479 rho_set%owns%rho_spin=.TRUE. 00480 rho_set%has%rho_spin=.FALSE. 00481 END IF 00482 ! Allocate the storage space for the norm of the gradient of the total density 00483 IF (needs%norm_drho) THEN 00484 NULLIFY(rho_set%norm_drho) 00485 ALLOCATE(rho_set%norm_drho(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00486 rho_set%owns%norm_drho=.TRUE. 00487 rho_set%has%norm_drho=.FALSE. 00488 END IF 00489 ! Allocate the storage space for the norm of the gradient of rhoa and of rhob separatedly 00490 IF (needs%norm_drho_spin) THEN 00491 NULLIFY(rho_set%norm_drhoa,rho_set%norm_drhob) 00492 ALLOCATE(rho_set%norm_drhoa(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00493 ALLOCATE(rho_set%norm_drhob(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00494 rho_set%owns%norm_drho_spin=.TRUE. 00495 rho_set%has%norm_drho_spin=.FALSE. 00496 END IF 00497 ! 00498 IF (needs%drhoa_drhob) THEN 00499 NULLIFY(rho_set%drhoa_drhob) 00500 ALLOCATE(rho_set%drhoa_drhob(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00501 rho_set%owns%drhoa_drhob=.TRUE. 00502 rho_set%has%drhoa_drhob=.FALSE. 00503 END IF 00504 ! Allocate the storage space for the components of the gradient for the total rho 00505 IF (needs%drho) THEN 00506 DO idir = 1,3 00507 NULLIFY(rho_set%drho(idir)%array) 00508 ALLOCATE(rho_set%drho(idir)%array(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00509 END DO 00510 rho_set%owns%drho=.TRUE. 00511 rho_set%has%drho=.FALSE. 00512 END IF 00513 ! Allocate the storage space for the components of the gradient for rhoa and rhob 00514 IF (needs%drho_spin) THEN 00515 DO idir = 1,3 00516 NULLIFY(rho_set%drhoa(idir)%array,rho_set%drhob(idir)%array) 00517 ALLOCATE(rho_set%drhoa(idir)%array(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00518 ALLOCATE(rho_set%drhob(idir)%array(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00519 END DO 00520 rho_set%owns%drho_spin=.TRUE. 00521 rho_set%has%drho_spin=.FALSE. 00522 END IF 00523 ! 00524 END SELECT 00525 00526 ! tau part 00527 IF (needs%tau) THEN 00528 NULLIFY(rho_set%tau) 00529 ALLOCATE(rho_set%tau(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00530 rho_set%owns%tau=.TRUE. 00531 END IF 00532 IF (needs%tau_spin) THEN 00533 NULLIFY(rho_set%tau_a,rho_set%tau_b) 00534 ALLOCATE(rho_set%tau_a(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00535 ALLOCATE(rho_set%tau_b(bo(1,1):bo(2,1),bo(1,2):bo(2,2),bo(1,3):bo(2,3))) 00536 rho_set%owns%tau_spin=.TRUE. 00537 rho_set%has%tau_spin=.FALSE. 00538 END IF 00539 00540 END SUBROUTINE xc_rho_set_atom_update 00541 00542 ! ***************************************************************************** 00543 SUBROUTINE fill_rho_set(rho_set,lsd,nspins,needs,rho,drho,tau,na,ir,error) 00544 00545 TYPE(xc_rho_set_type), POINTER :: rho_set 00546 LOGICAL, INTENT(IN) :: lsd 00547 INTEGER, INTENT(IN) :: nspins 00548 TYPE(xc_rho_cflags_type), INTENT(in) :: needs 00549 REAL(dp), DIMENSION(:, :), POINTER :: rho 00550 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho 00551 REAL(dp), DIMENSION(:, :), POINTER :: tau 00552 INTEGER, INTENT(IN) :: na, ir 00553 TYPE(cp_error_type), INTENT(inout) :: error 00554 00555 CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_rho_set', 00556 routineP = moduleN//':'//routineN 00557 REAL(KIND=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp) 00558 00559 INTEGER :: ia, idir, my_nspins 00560 LOGICAL :: failure, gradient_f, 00561 tddft_split 00562 00563 failure = .FALSE. 00564 my_nspins = nspins 00565 tddft_split = .FALSE. 00566 IF (lsd .AND. nspins==1) THEN 00567 my_nspins = 2 00568 tddft_split = .TRUE. 00569 END IF 00570 00571 ! some checks 00572 IF (lsd) THEN 00573 ELSE 00574 CPPrecondition(SIZE(rho,2)==1,cp_failure_level,routineP,error,failure) 00575 END IF 00576 SELECT CASE(my_nspins) 00577 CASE(1) 00578 CPPrecondition(.NOT.needs%rho_spin,cp_failure_level,routineP,error,failure) 00579 CPPrecondition(.NOT.needs%drho_spin,cp_failure_level,routineP,error,failure) 00580 CPPrecondition(.NOT.needs%norm_drho_spin,cp_failure_level,routineP,error,failure) 00581 CPPrecondition(.NOT.needs%drhoa_drhob,cp_failure_level,routineP,error,failure) 00582 CPPrecondition(.NOT.needs%rho_spin_1_3,cp_failure_level,routineP,error,failure) 00583 CASE(2) 00584 CASE default 00585 CPAssert(.FALSE.,cp_failure_level,routineP,error,failure) 00586 END SELECT 00587 00588 gradient_f=(needs%drho_spin.OR.needs%norm_drho_spin.OR.& 00589 needs%drhoa_drhob.OR.needs%drho.OR.needs%norm_drho) 00590 00591 SELECT CASE(my_nspins) 00592 CASE(1) 00593 ! Give rho to 1/3 00594 IF (needs%rho_1_3) THEN 00595 DO ia = 1,na 00596 rho_set%rho_1_3(ia,ir,1) = MAX(rho(ia,1),0.0_dp)**f13 00597 END DO 00598 rho_set%owns%rho_1_3=.TRUE. 00599 rho_set%has%rho_1_3=.TRUE. 00600 END IF 00601 ! Give the density 00602 IF (needs%rho) THEN 00603 DO ia = 1,na 00604 rho_set%rho(ia,ir,1) = rho(ia,1) 00605 END DO 00606 rho_set%owns%rho=.TRUE. 00607 rho_set%has%rho=.TRUE. 00608 END IF 00609 ! Give the norm of the gradient of the density 00610 IF (needs%norm_drho) THEN 00611 DO ia = 1,na 00612 rho_set%norm_drho(ia,ir,1) = drho(4,ia,ir,1) 00613 END DO 00614 rho_set%owns%norm_drho=.TRUE. 00615 rho_set%has%norm_drho=.TRUE. 00616 END IF 00617 ! Give the three components of the gradient of the density 00618 IF (needs%drho) THEN 00619 DO idir = 1,3 00620 DO ia = 1,na 00621 rho_set%drho(idir)%array(ia,ir,1) = drho(idir,ia,ir,1) 00622 END DO 00623 END DO 00624 rho_set%owns%drho=.TRUE. 00625 rho_set%has%drho=.TRUE. 00626 END IF 00627 CASE(2) 00628 ! Give the total density 00629 IF (needs%rho) THEN 00630 ! this should never be the case unless you use LDA functionals with LSD 00631 IF (.NOT.tddft_split) THEN 00632 DO ia = 1,na 00633 rho_set%rho(ia,ir,1) = rho(ia,1)+rho(ia,2) 00634 END DO 00635 ELSE 00636 DO ia = 1,na 00637 rho_set%rho(ia,ir,1) = rho(ia,1) 00638 END DO 00639 END IF 00640 rho_set%owns%rho=.TRUE. 00641 rho_set%has%rho=.TRUE. 00642 END IF 00643 ! Give the total density to 1/3 00644 IF (needs%rho_1_3) THEN 00645 IF (.NOT.tddft_split) THEN 00646 DO ia = 1,na 00647 rho_set%rho_1_3(ia,ir,1) = MAX(rho(ia,1)+rho(ia,2),0.0_dp)**f13 00648 END DO 00649 ELSE 00650 DO ia = 1,na 00651 rho_set%rho_1_3(ia,ir,1) = MAX(rho(ia,1),0.0_dp)**f13 00652 END DO 00653 END IF 00654 rho_set%owns%rho_1_3=.TRUE. 00655 rho_set%has%rho_1_3=.TRUE. 00656 END IF 00657 ! Give the spin densities to 1/3 00658 IF (needs%rho_spin_1_3) THEN 00659 IF (.NOT.tddft_split) THEN 00660 DO ia = 1,na 00661 rho_set%rhoa_1_3(ia,ir,1) = MAX(rho(ia,1),0.0_dp)**f13 00662 rho_set%rhob_1_3(ia,ir,1) = MAX(rho(ia,2),0.0_dp)**f13 00663 END DO 00664 ELSE 00665 DO ia = 1,na 00666 rho_set%rhoa_1_3(ia,ir,1) = MAX(0.5_dp*rho(ia,1),0.0_dp)**f13 00667 rho_set%rhob_1_3(ia,ir,1) = rho_set%rhoa_1_3(ia,ir,1) 00668 END DO 00669 END IF 00670 rho_set%owns%rho_spin_1_3=.TRUE. 00671 rho_set%has%rho_spin_1_3=.TRUE. 00672 END IF 00673 ! Give the spin densities rhoa and rhob 00674 IF (needs%rho_spin) THEN 00675 IF (.NOT.tddft_split) THEN 00676 DO ia = 1,na 00677 rho_set%rhoa(ia,ir,1) = rho(ia,1) 00678 rho_set%rhob(ia,ir,1) = rho(ia,2) 00679 END DO 00680 ELSE 00681 DO ia = 1,na 00682 rho_set%rhoa(ia,ir,1) = 0.5_dp*rho(ia,1) 00683 rho_set%rhob(ia,ir,1) = rho_set%rhoa(ia,ir,1) 00684 END DO 00685 END IF 00686 rho_set%owns%rho_spin=.TRUE. 00687 rho_set%has%rho_spin=.TRUE. 00688 END IF 00689 ! Give the norm of the gradient of the total density 00690 IF (needs%norm_drho) THEN 00691 IF (.NOT.tddft_split) THEN 00692 DO ia = 1,na 00693 rho_set%norm_drho(ia,ir,1) = SQRT(& 00694 (drho(1,ia,ir,1)+drho(1,ia,ir,2))**2+& 00695 (drho(2,ia,ir,1)+drho(2,ia,ir,2))**2+& 00696 (drho(3,ia,ir,1)+drho(3,ia,ir,2))**2) 00697 END DO 00698 ELSE 00699 DO ia = 1,na 00700 rho_set%norm_drho(ia,ir,1) = drho(4,ia,ir,1) 00701 !! rho_set%norm_drho(ia,ir,1) = SQRT(& 00702 !! (drho(1,ia,ir,1))**2+& 00703 !! (drho(2,ia,ir,1))**2+& 00704 !! (drho(3,ia,ir,1))**2) 00705 END DO 00706 END IF 00707 rho_set%owns%norm_drho=.TRUE. 00708 rho_set%has%norm_drho=.TRUE. 00709 END IF 00710 ! Give the norm of the gradient of rhoa and of rhob separatedly 00711 IF (needs%norm_drho_spin) THEN 00712 IF (.NOT.tddft_split) THEN 00713 DO ia = 1,na 00714 rho_set%norm_drhoa(ia,ir,1) = drho(4,ia,ir,1) 00715 rho_set%norm_drhob(ia,ir,1) = drho(4,ia,ir,2) 00716 END DO 00717 ELSE 00718 DO ia = 1,na 00719 rho_set%norm_drhoa(ia,ir,1) = 0.5_dp*drho(4,ia,ir,1) 00720 rho_set%norm_drhob(ia,ir,1) = rho_set%norm_drhoa(ia,ir,1) 00721 END DO 00722 END IF 00723 rho_set%owns%norm_drho_spin=.TRUE. 00724 rho_set%has%norm_drho_spin=.TRUE. 00725 END IF 00726 ! 00727 IF (needs%drhoa_drhob) THEN 00728 IF (.NOT.tddft_split) THEN 00729 DO ia =1,na 00730 rho_set%drhoa_drhob(ia,ir,1) =& 00731 (drho(1,ia,ir,1)*drho(1,ia,ir,2))+& 00732 (drho(2,ia,ir,1)*drho(2,ia,ir,2))+& 00733 (drho(3,ia,ir,1)*drho(3,ia,ir,2)) 00734 END DO 00735 ELSE 00736 DO ia =1,na 00737 rho_set%drhoa_drhob(ia,ir,1) = 0.25_dp * drho(4,ia,ir,1)**2 00738 !! rho_set%drhoa_drhob(ia,ir,1) =& 00739 !! (0.25_dp*drho(1,ia,ir,1)**2)+& 00740 !! (0.25_dp*drho(2,ia,ir,1)**2)+& 00741 !! (0.25_dp*drho(3,ia,ir,1)**2) 00742 END DO 00743 END IF 00744 rho_set%owns%drhoa_drhob=.TRUE. 00745 rho_set%has%drhoa_drhob=.TRUE. 00746 END IF 00747 ! Give the components of the gradient for the total rho 00748 IF (needs%drho) THEN 00749 IF (.NOT.tddft_split) THEN 00750 DO idir = 1,3 00751 DO ia = 1,na 00752 rho_set%drho(idir)%array(ia,ir,1) = drho(idir,ia,ir,1)+drho(idir,ia,ir,2) 00753 END DO 00754 END DO 00755 ELSE 00756 DO idir = 1,3 00757 DO ia = 1,na 00758 rho_set%drho(idir)%array(ia,ir,1) = drho(idir,ia,ir,1) 00759 END DO 00760 END DO 00761 END IF 00762 rho_set%owns%drho=.TRUE. 00763 rho_set%has%drho=.TRUE. 00764 END IF 00765 ! Give the components of the gradient for rhoa and rhob 00766 IF (needs%drho_spin) THEN 00767 IF (.NOT.tddft_split) THEN 00768 DO idir = 1,3 00769 DO ia = 1,na 00770 rho_set%drhoa(idir)%array(ia,ir,1) = drho(idir,ia,ir,1) 00771 rho_set%drhob(idir)%array(ia,ir,1) = drho(idir,ia,ir,2) 00772 END DO 00773 END DO 00774 ELSE 00775 DO idir = 1,3 00776 DO ia = 1,na 00777 rho_set%drhoa(idir)%array(ia,ir,1) = 0.5_dp*drho(idir,ia,ir,1) 00778 rho_set%drhob(idir)%array(ia,ir,1) = rho_set%drhoa(idir)%array(ia,ir,1) 00779 END DO 00780 END DO 00781 END IF 00782 rho_set%owns%drho_spin=.TRUE. 00783 rho_set%has%drho_spin=.TRUE. 00784 END IF 00785 ! 00786 END SELECT 00787 00788 ! tau part 00789 IF (needs%tau.OR.needs%tau_spin) THEN 00790 CPPrecondition(ASSOCIATED(tau),cp_failure_level,routineP,error,failure) 00791 CPPrecondition(SIZE(tau,2)==my_nspins,cp_failure_level,routineP,error,failure) 00792 END IF 00793 IF (needs%tau) THEN 00794 IF (my_nspins==2) THEN 00795 DO ia = 1,na 00796 rho_set%tau(ia,ir,1) = tau(ia,1)+tau(ia,2) 00797 END DO 00798 rho_set%owns%tau=.TRUE. 00799 rho_set%has%tau=.TRUE. 00800 ELSE 00801 DO ia = 1,na 00802 rho_set%tau(ia,ir,1) = tau(ia,1) 00803 END DO 00804 rho_set%owns%tau=.TRUE. 00805 rho_set%has%tau=.TRUE. 00806 END IF 00807 END IF 00808 IF (needs%tau_spin) THEN 00809 DO ia = 1,na 00810 rho_set%tau_a(ia,ir,1) = tau(ia,1) 00811 rho_set%tau_b(ia,ir,1) = tau(ia,2) 00812 END DO 00813 rho_set%owns%tau_spin=.TRUE. 00814 rho_set%has%tau_spin=.TRUE. 00815 END IF 00816 00817 END SUBROUTINE fill_rho_set 00818 00819 END MODULE xc_atom
1.7.3