|
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 MODULE hartree_local_methods 00007 00008 USE atomic_kind_types, ONLY: atomic_kind_type,& 00009 get_atomic_kind,& 00010 get_atomic_kind_set 00011 USE atprop_types, ONLY: atprop_array_init 00012 USE basis_set_types, ONLY: get_gto_basis_set,& 00013 gto_basis_set_type 00014 USE cell_types, ONLY: cell_type 00015 USE cp_control_types, ONLY: dft_control_type 00016 USE cp_para_types, ONLY: cp_para_env_type 00017 USE hartree_local_types, ONLY: allocate_ecoul_1center,& 00018 ecoul_1center_type,& 00019 hartree_local_type,& 00020 set_ecoul_1c 00021 USE input_constants, ONLY: tddfpt_singlet,& 00022 use_periodic 00023 USE kinds, ONLY: dp 00024 USE mathconstants, ONLY: fourpi,& 00025 pi 00026 USE message_passing, ONLY: mp_sum 00027 USE orbital_pointers, ONLY: indso,& 00028 nsoset 00029 USE pw_env_types, ONLY: pw_env_get,& 00030 pw_env_type 00031 USE pw_poisson_types, ONLY: pw_poisson_type 00032 USE qs_charges_types, ONLY: qs_charges_type 00033 USE qs_environment_types, ONLY: get_qs_env,& 00034 qs_environment_type 00035 USE qs_grid_atom, ONLY: grid_atom_type 00036 USE qs_harmonics_atom, ONLY: get_none0_cg_list,& 00037 harmonics_atom_type 00038 USE qs_local_rho_types, ONLY: rhoz_type 00039 USE qs_p_env_types, ONLY: qs_p_env_type 00040 USE qs_rho0_types, ONLY: get_rho0_mpole,& 00041 rho0_atom_type,& 00042 rho0_mpole_type 00043 USE qs_rho_atom_types, ONLY: get_rho_atom,& 00044 rho_atom_coeff,& 00045 rho_atom_type 00046 USE termination, ONLY: stop_program 00047 USE timings, ONLY: timeset,& 00048 timestop 00049 USE util, ONLY: get_limit 00050 #include "cp_common_uses.h" 00051 00052 IMPLICIT NONE 00053 00054 PRIVATE 00055 00056 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hartree_local_methods' 00057 00058 ! Public Subroutine 00059 00060 PUBLIC :: init_coulomb_local, calculate_Vh_1center, Vh_1c_gg_integrals 00061 00062 CONTAINS 00063 00064 ! ***************************************************************************** 00065 SUBROUTINE init_coulomb_local(hartree_local,natom,error) 00066 00067 TYPE(hartree_local_type), POINTER :: hartree_local 00068 INTEGER, INTENT(IN) :: natom 00069 TYPE(cp_error_type), INTENT(inout) :: error 00070 00071 CHARACTER(len=*), PARAMETER :: routineN = 'init_coulomb_local', 00072 routineP = moduleN//':'//routineN 00073 00074 INTEGER :: handle 00075 TYPE(ecoul_1center_type), DIMENSION(:), 00076 POINTER :: ecoul_1c 00077 00078 CALL timeset(routineN,handle) 00079 00080 NULLIFY(ecoul_1c) 00081 ! Allocate and Initialize 1-center Potentials and Integrals 00082 CALL allocate_ecoul_1center(ecoul_1c,natom, error) 00083 hartree_local%ecoul_1c => ecoul_1c 00084 00085 CALL timestop(handle) 00086 00087 END SUBROUTINE init_coulomb_local 00088 00089 ! ***************************************************************************** 00096 SUBROUTINE calculate_Vh_1center(vrad_h,vrad_s,rrad_h,rrad_s,rrad_0,rrad_z,grid_atom) 00097 00098 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: vrad_h, vrad_s 00099 TYPE(rho_atom_coeff), DIMENSION(:), 00100 INTENT(IN) :: rrad_h, rrad_s 00101 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rrad_0 00102 REAL(dp), DIMENSION(:), INTENT(IN) :: rrad_z 00103 TYPE(grid_atom_type), POINTER :: grid_atom 00104 00105 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_Vh_1center', 00106 routineP = moduleN//':'//routineN 00107 00108 INTEGER :: handle, ir, iso, ispin, 00109 l_ang, max_s_harm, nchannels, 00110 nr, nspins 00111 REAL(dp) :: I1_down, I1_up, I2_down, 00112 I2_up, prefactor 00113 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rho_1, rho_2 00114 REAL(dp), DIMENSION(:), POINTER :: wr 00115 REAL(dp), DIMENSION(:, :), POINTER :: oor2l, r2l 00116 00117 CALL timeset(routineN,handle) 00118 00119 nr = grid_atom%nr 00120 max_s_harm = SIZE(vrad_h,2) 00121 nspins = SIZE(rrad_h,1) 00122 nchannels = SIZE(rrad_0,2) 00123 00124 r2l => grid_atom%rad2l 00125 oor2l => grid_atom%oorad2l 00126 wr => grid_atom%wr 00127 00128 ALLOCATE(rho_1(nr,max_s_harm),rho_2(nr,max_s_harm)) 00129 rho_1 = 0.0_dp 00130 rho_2 = 0.0_dp 00131 00132 ! Case lm = 0 00133 rho_1(:,1) = rrad_z(:) 00134 rho_2(:,1) = rrad_0(:,1) 00135 00136 DO iso = 2,nchannels 00137 rho_2(:,iso) = rrad_0(:,iso) 00138 END DO 00139 00140 DO iso = 1,max_s_harm 00141 DO ispin = 1,nspins 00142 rho_1(:,iso) = rho_1(:,iso) + rrad_h(ispin)%r_coef(:,iso) 00143 rho_2(:,iso) = rho_2(:,iso) + rrad_s(ispin)%r_coef(:,iso) 00144 END DO 00145 00146 l_ang = indso(1,iso) 00147 prefactor = fourpi/(2._dp*l_ang+1._dp) 00148 00149 rho_1(:,iso) = rho_1(:,iso)*wr(:) 00150 rho_2(:,iso) = rho_2(:,iso)*wr(:) 00151 00152 I1_up = 0.0_dp 00153 I1_down = 0.0_dp 00154 I2_up = 0.0_dp 00155 I2_down = 0.0_dp 00156 00157 I1_up = r2l(nr,l_ang)*rho_1(nr,iso) 00158 I2_up = r2l(nr,l_ang)*rho_2(nr,iso) 00159 00160 DO ir = nr-1,1,-1 00161 I1_down = I1_down + oor2l(ir,l_ang+1)*rho_1(ir,iso) 00162 I2_down = I2_down + oor2l(ir,l_ang+1)*rho_2(ir,iso) 00163 END DO 00164 00165 vrad_h(nr,iso) = vrad_h(nr,iso) + prefactor*& 00166 (oor2l(nr,l_ang+1)*I1_up + r2l(nr,l_ang)*I1_down) 00167 vrad_s(nr,iso) = vrad_s(nr,iso) + prefactor*& 00168 (oor2l(nr,l_ang+1)*I2_up + r2l(nr,l_ang)*I2_down) 00169 00170 DO ir = nr-1,1,-1 00171 I1_up = I1_up + r2l(ir,l_ang)*rho_1(ir,iso) 00172 I1_down = I1_down -oor2l(ir,l_ang+1)*rho_1(ir,iso) 00173 I2_up = I2_up + r2l(ir,l_ang)*rho_2(ir,iso) 00174 I2_down = I2_down -oor2l(ir,l_ang+1)*rho_2(ir,iso) 00175 00176 vrad_h(ir,iso) = vrad_h(ir,iso) + prefactor*& 00177 (oor2l(ir,l_ang+1)*I1_up + r2l(ir,l_ang)*I1_down) 00178 vrad_s(ir,iso) = vrad_s(ir,iso) + prefactor*& 00179 (oor2l(ir,l_ang+1)*I2_up + r2l(ir,l_ang)*I2_down) 00180 00181 END DO 00182 00183 END DO 00184 00185 DEALLOCATE(rho_1,rho_2) 00186 00187 CALL timestop(handle) 00188 00189 END SUBROUTINE calculate_Vh_1center 00190 00191 ! ***************************************************************************** 00200 SUBROUTINE Vh_1c_gg_integrals(qs_env,energy_hartree_1c,tddft,do_triplet,p_env,error) 00201 00202 TYPE(qs_environment_type), POINTER :: qs_env 00203 REAL(kind=dp), INTENT(out) :: energy_hartree_1c 00204 LOGICAL, INTENT(IN), OPTIONAL :: tddft, do_triplet 00205 TYPE(qs_p_env_type), OPTIONAL, POINTER :: p_env 00206 TYPE(cp_error_type), INTENT(inout) :: error 00207 00208 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_gg_integrals', 00209 routineP = moduleN//':'//routineN 00210 00211 INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, 00212 istat, l_ang, llmax, lmax0, lmax_0, m1, max_iso, max_iso_not0, 00213 max_s_harm, maxl, maxso, mepos, n1, nat, natom, nchan_0, nkind, nr, 00214 nset, nsotot, nspins, num_pe 00215 INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list 00216 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list 00217 INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf 00218 LOGICAL :: atenergy, core_charge, 00219 failure, my_periodic, 00220 my_tddft, paw_atom 00221 REAL(dp) :: back_ch, factor 00222 REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr 00223 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_hh, aVh1b_ss, 00224 g0_h_w 00225 REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z 00226 REAL(dp), DIMENSION(:, :), POINTER :: g0_h, gsph, rrad_0, Vh1_h, 00227 Vh1_s, vrrad_0, zet 00228 REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg 00229 TYPE(atomic_kind_type), DIMENSION(:), 00230 POINTER :: atomic_kind_set 00231 TYPE(atomic_kind_type), POINTER :: atom_kind 00232 TYPE(cell_type), POINTER :: cell 00233 TYPE(cp_para_env_type), POINTER :: para_env 00234 TYPE(dft_control_type), POINTER :: dft_control 00235 TYPE(ecoul_1center_type), DIMENSION(:), 00236 POINTER :: ecoul_1c 00237 TYPE(grid_atom_type), POINTER :: grid_atom 00238 TYPE(gto_basis_set_type), POINTER :: orb_basis 00239 TYPE(harmonics_atom_type), POINTER :: harmonics 00240 TYPE(pw_env_type), POINTER :: pw_env 00241 TYPE(pw_poisson_type), POINTER :: poisson_env 00242 TYPE(qs_charges_type), POINTER :: qs_charges 00243 TYPE(rho0_atom_type), DIMENSION(:), 00244 POINTER :: rho0_atom_set 00245 TYPE(rho0_mpole_type), POINTER :: rho0_mpole 00246 TYPE(rho_atom_type), DIMENSION(:), 00247 POINTER :: rho_atom_set 00248 TYPE(rho_atom_type), POINTER :: rho_atom 00249 TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set 00250 00251 CALL timeset(routineN,handle) 00252 00253 failure = .FALSE. 00254 NULLIFY(cell,dft_control,para_env,poisson_env,pw_env,qs_charges) 00255 NULLIFY(atomic_kind_set, rho_atom_set, rho0_atom_set) 00256 NULLIFY(rho0_mpole,rhoz_set,ecoul_1c) 00257 NULLIFY(atom_kind,atom_list,grid_atom,harmonics) 00258 NULLIFY(orb_basis,lmin,lmax,npgf,zet) 00259 NULLIFY(gsph) 00260 00261 my_tddft = .FALSE. 00262 IF (PRESENT(tddft)) my_tddft = tddft 00263 core_charge = .NOT.my_tddft 00264 00265 CALL get_qs_env(qs_env=qs_env,& 00266 cell=cell,dft_control=dft_control,& 00267 para_env=para_env,& 00268 atomic_kind_set=atomic_kind_set,& 00269 pw_env=pw_env,qs_charges=qs_charges,error=error) 00270 00271 CALL pw_env_get(pw_env,poisson_env=poisson_env,error=error) 00272 my_periodic= (poisson_env%method==use_periodic) 00273 00274 back_ch = qs_charges%background*cell%deth 00275 00276 IF (my_tddft) THEN 00277 rho_atom_set => p_env%local_rho_set%rho_atom_set 00278 rho0_atom_set=> p_env%local_rho_set%rho0_atom_set 00279 rho0_mpole => p_env%local_rho_set%rho0_mpole 00280 ecoul_1c => p_env%hartree_local%ecoul_1c 00281 ELSE 00282 CALL get_qs_env(qs_env=qs_env, & 00283 rho_atom_set= rho_atom_set,& 00284 rho0_atom_set=rho0_atom_set, & 00285 rho0_mpole=rho0_mpole,& 00286 rhoz_set=rhoz_set,& 00287 ecoul_1c=ecoul_1c,error=error) 00288 END IF 00289 00290 nkind = SIZE(atomic_kind_set,1) 00291 nspins = dft_control%nspins 00292 00293 IF (my_tddft) THEN 00294 IF(PRESENT(do_triplet)) THEN 00295 IF (nspins==1.AND.do_triplet) RETURN 00296 ELSE 00297 IF (nspins==1.AND.dft_control%tddfpt_control%res_etype/=tddfpt_singlet) RETURN 00298 ENDIF 00299 END IF 00300 00301 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,maxg_iso_not0=max_iso) 00302 CALL get_rho0_mpole(rho0_mpole=rho0_mpole,lmax_0=lmax_0) 00303 00304 atenergy=.FALSE. 00305 IF (ASSOCIATED(qs_env%atprop)) THEN 00306 atenergy = qs_env%atprop%energy 00307 IF (atenergy) THEN 00308 CALL get_qs_env(qs_env=qs_env,natom=natom,error=error) 00309 CALL atprop_array_init(qs_env%atprop%ate1c,natom,error) 00310 END IF 00311 END IF 00312 00313 ! Put to 0 the local hartree energy contribution from 1 center integrals 00314 energy_hartree_1c = 0.0_dp 00315 00316 ! Here starts the loop over all the atoms 00317 DO ikind = 1,nkind 00318 00319 atom_kind => atomic_kind_set(ikind) 00320 CALL get_atomic_kind(atomic_kind=atom_kind, atom_list=atom_list,& 00321 orb_basis_set=orb_basis,& 00322 natom=nat,grid_atom=grid_atom,& 00323 harmonics=harmonics,ngrid_rad=nr,& 00324 max_iso_not0=max_iso_not0,paw_atom=paw_atom) 00325 IF(paw_atom) THEN 00326 !=========== PAW =============== 00327 CALL get_gto_basis_set(gto_basis_set=orb_basis,lmax=lmax,lmin=lmin,& 00328 maxso=maxso,npgf=npgf,maxl=maxl,& 00329 nset=nset,zet=zet) 00330 00331 max_s_harm = harmonics%max_s_harm 00332 llmax = harmonics%llmax 00333 00334 nsotot = maxso*nset 00335 ALLOCATE(gsph(nr,nsotot),STAT=istat) 00336 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00337 ALLOCATE(gexp(nr),STAT=istat) 00338 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00339 ALLOCATE(sqrtwr(nr),g0_h_w(nr,0:lmax_0),STAT=istat) 00340 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00341 00342 NULLIFY(Vh1_h,Vh1_s) 00343 ALLOCATE(Vh1_h(nr,max_iso_not0),STAT=istat) 00344 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00345 ALLOCATE(Vh1_s(nr,max_iso_not0),STAT=istat) 00346 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00347 00348 ALLOCATE(aVh1b_hh(nsotot,nsotot),STAT=istat) 00349 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00350 ALLOCATE(aVh1b_ss(nsotot,nsotot),STAT=istat) 00351 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00352 ALLOCATE(aVh1b_00(nsotot,nsotot),STAT=istat) 00353 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00354 ALLOCATE(cg_list(2,nsoset(maxl)**2,max_s_harm),cg_n_list(max_s_harm),STAT=istat) 00355 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00356 00357 NULLIFY(Qlm_gg,g0_h) 00358 CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, & 00359 l0_ikind=lmax0, & 00360 Qlm_gg=Qlm_gg, g0_h=g0_h) 00361 00362 nchan_0 = nsoset(lmax0) 00363 00364 IF(nchan_0 > max_iso_not0) CALL stop_program(routineN,moduleN,__LINE__,& 00365 "channels for rho0 > # max of spherical harmonics",para_env) 00366 00367 NULLIFY(rrad_z,my_CG) 00368 my_CG => harmonics%my_CG 00369 00370 ! set to zero temporary arrays 00371 sqrtwr=0.0_dp 00372 g0_h_w=0.0_dp 00373 gexp=0.0_dp 00374 gsph=0.0_dp 00375 00376 sqrtwr(1:nr) = SQRT(grid_atom%wr(1:nr)) 00377 DO l_ang = 0,lmax0 00378 g0_h_w(1:nr,l_ang) = g0_h(1:nr,l_ang)*grid_atom%wr(1:nr) 00379 END DO 00380 00381 m1 = 0 00382 DO iset1 = 1,nset 00383 n1 = nsoset(lmax(iset1)) 00384 DO ipgf1 = 1,npgf(iset1) 00385 gexp(1:nr) = EXP(-zet(ipgf1,iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr) 00386 DO is1 = nsoset(lmin(iset1)-1)+1,nsoset(lmax(iset1)) 00387 iso = is1 + (ipgf1-1)*n1 + m1 00388 l_ang = indso(1,is1) 00389 gsph(1:nr,iso) = grid_atom%rad2l(1:nr,l_ang)*gexp(1:nr) 00390 END DO ! is1 00391 END DO ! ipgf1 00392 m1 = m1 + maxso 00393 END DO ! iset1 00394 00395 ! Distribute the atoms of this kind 00396 num_pe = para_env%num_pe 00397 mepos = para_env%mepos 00398 bo = get_limit( nat, num_pe, mepos ) 00399 00400 DO iat = bo(1), bo(2) !1,nat 00401 iatom = atom_list(iat) 00402 rho_atom => rho_atom_set(iatom) 00403 00404 NULLIFY(rrad_z,vrrad_z,rrad_0,vrrad_0) 00405 IF(core_charge) THEN 00406 rrad_z => rhoz_set(ikind)%r_coef 00407 vrrad_z => rhoz_set(ikind)%vr_coef 00408 END IF 00409 rrad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef 00410 vrrad_0 => rho0_atom_set(iatom)%vrho0_rad_h%r_coef 00411 IF(my_periodic .AND. back_ch .GT. 1.E-3_dp) THEN 00412 factor = -2.0_dp*pi/3.0_dp*SQRT(fourpi)*qs_charges%background 00413 ELSE 00414 factor = 0._dp 00415 END IF 00416 00417 CALL Vh_1c_atom_potential(rho_atom,vrrad_0,& 00418 grid_atom,iatom,core_charge,vrrad_z,Vh1_h,Vh1_s,& 00419 nchan_0,nspins,max_iso_not0,factor,error) 00420 00421 CALL Vh_1c_atom_energy(energy_hartree_1c,ecoul_1c,rho_atom,rrad_0,& 00422 grid_atom,iatom,core_charge,rrad_z,Vh1_h,Vh1_s,& 00423 nchan_0,nspins,max_iso_not0,atenergy,qs_env%atprop%ate1c,error) 00424 00425 CALL Vh_1c_atom_integrals(rho_atom,& 00426 aVh1b_hh,aVh1b_ss,aVh1b_00,Vh1_h,Vh1_s,max_iso_not0,& 00427 max_s_harm,llmax,cg_list,cg_n_list,& 00428 nset,npgf,lmin,lmax,nsotot,maxso,nspins,nchan_0,gsph,g0_h_w,my_CG,Qlm_gg,error) 00429 00430 END DO ! iat 00431 00432 DEALLOCATE(aVh1b_hh,STAT=istat) 00433 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00434 DEALLOCATE(aVh1b_ss,STAT=istat) 00435 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00436 DEALLOCATE(aVh1b_00,STAT=istat) 00437 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00438 DEALLOCATE(Vh1_h,Vh1_s, STAT=istat) 00439 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00440 DEALLOCATE(cg_list,cg_n_list,STAT=istat) 00441 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00442 DEALLOCATE(gsph,STAT=istat) 00443 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00444 DEALLOCATE(gexp,STAT=istat) 00445 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00446 DEALLOCATE(sqrtwr,g0_h_w,STAT=istat) 00447 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00448 00449 ELSE 00450 !=========== NO PAW =============== 00451 ! This term is taken care of using the core density as in GPW 00452 CYCLE 00453 END IF ! paw 00454 END DO ! ikind 00455 00456 CALL mp_sum(energy_hartree_1c,para_env%group) 00457 00458 CALL timestop(handle) 00459 00460 END SUBROUTINE Vh_1c_gg_integrals 00461 00462 ! ***************************************************************************** 00463 00464 SUBROUTINE Vh_1c_atom_potential(rho_atom,vrrad_0,& 00465 grid_atom,iatom,core_charge,vrrad_z,Vh1_h,Vh1_s,& 00466 nchan_0,nspins,max_iso_not0,bfactor,error) 00467 00468 TYPE(rho_atom_type), POINTER :: rho_atom 00469 REAL(dp), DIMENSION(:, :), POINTER :: vrrad_0 00470 TYPE(grid_atom_type), POINTER :: grid_atom 00471 INTEGER, INTENT(IN) :: iatom 00472 LOGICAL, INTENT(IN) :: core_charge 00473 REAL(dp), DIMENSION(:), POINTER :: vrrad_z 00474 REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s 00475 INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0 00476 REAL(dp), INTENT(IN) :: bfactor 00477 TYPE(cp_error_type), INTENT(inout) :: error 00478 00479 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_atom_potential', 00480 routineP = moduleN//':'//routineN 00481 00482 INTEGER :: ir, iso, ispin, nr 00483 TYPE(rho_atom_coeff), DIMENSION(:), 00484 POINTER :: vr_h, vr_s 00485 00486 nr = grid_atom%nr 00487 00488 NULLIFY(vr_h,vr_s) 00489 CALL get_rho_atom(rho_atom=rho_atom,vrho_rad_h=vr_h,vrho_rad_s=vr_s) 00490 00491 Vh1_h = 0.0_dp 00492 Vh1_s = 0.0_dp 00493 00494 IF (core_charge) Vh1_h(:,1) = vrrad_z(:) 00495 00496 DO iso = 1,nchan_0 00497 Vh1_s(:,iso) = vrrad_0(:,iso) 00498 END DO 00499 00500 DO ispin = 1,nspins 00501 DO iso =1,max_iso_not0 00502 Vh1_h(:,iso) = Vh1_h(:,iso) + vr_h(ispin)%r_coef(:,iso) 00503 Vh1_s(:,iso) = Vh1_s(:,iso) + vr_s(ispin)%r_coef(:,iso) 00504 END DO 00505 END DO 00506 00507 IF(bfactor /= 0._dp) THEN 00508 DO ir = 1,nr 00509 Vh1_h(ir,1) = Vh1_h(ir,1) + bfactor * grid_atom%rad2(ir)*grid_atom%wr(ir) 00510 Vh1_s(ir,1) = Vh1_s(ir,1) + bfactor * grid_atom%rad2(ir)*grid_atom%wr(ir) 00511 END DO 00512 END IF 00513 00514 END SUBROUTINE Vh_1c_atom_potential 00515 00516 ! ***************************************************************************** 00517 00518 SUBROUTINE Vh_1c_atom_energy(energy_hartree_1c,ecoul_1c,rho_atom,rrad_0,& 00519 grid_atom,iatom,core_charge,rrad_z,Vh1_h,Vh1_s,& 00520 nchan_0,nspins,max_iso_not0,atenergy,ate1c,error) 00521 00522 REAL(dp), INTENT(INOUT) :: energy_hartree_1c 00523 TYPE(ecoul_1center_type), DIMENSION(:), 00524 POINTER :: ecoul_1c 00525 TYPE(rho_atom_type), POINTER :: rho_atom 00526 REAL(dp), DIMENSION(:, :), POINTER :: rrad_0 00527 TYPE(grid_atom_type), POINTER :: grid_atom 00528 INTEGER, INTENT(IN) :: iatom 00529 LOGICAL, INTENT(IN) :: core_charge 00530 REAL(dp), DIMENSION(:), POINTER :: rrad_z 00531 REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s 00532 INTEGER, INTENT(IN) :: nchan_0, nspins, max_iso_not0 00533 LOGICAL, INTENT(IN) :: atenergy 00534 REAL(dp), DIMENSION(:), POINTER :: ate1c 00535 TYPE(cp_error_type), INTENT(inout) :: error 00536 00537 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_atom_energy', 00538 routineP = moduleN//':'//routineN 00539 00540 INTEGER :: iso, ispin, nr 00541 REAL(dp) :: ecoul_1_0, ecoul_1_h, 00542 ecoul_1_s, ecoul_1_z 00543 TYPE(rho_atom_coeff), DIMENSION(:), 00544 POINTER :: r_h, r_s 00545 00546 nr = grid_atom%nr 00547 00548 NULLIFY(r_h,r_s) 00549 CALL get_rho_atom(rho_atom=rho_atom,rho_rad_h=r_h,rho_rad_s=r_s) 00550 00551 ! Calculate the contributions to Ecoul coming from Vh1_h*rhoz 00552 ecoul_1_z = 0.0_dp 00553 IF(core_charge) THEN 00554 ecoul_1_z = 0.5_dp*SUM(Vh1_h(:,1)*rrad_z(:)*grid_atom%wr(:)) 00555 END IF 00556 00557 ! Calculate the contributions to Ecoul coming from Vh1_s*rho0 00558 ecoul_1_0 = 0.0_dp 00559 DO iso = 1,nchan_0 00560 ecoul_1_0 = ecoul_1_0 + 0.5_dp*SUM(Vh1_s(:,iso)*rrad_0(:,iso)*grid_atom%wr(:)) 00561 END DO 00562 00563 ! Calculate the contributions to Ecoul coming from Vh1_h*rho1_h and Vh1_s*rho1_s 00564 ecoul_1_s = 0.0_dp 00565 ecoul_1_h = 0.0_dp 00566 DO ispin = 1,nspins 00567 DO iso = 1,max_iso_not0 00568 ecoul_1_s = ecoul_1_s + 0.5_dp*SUM(Vh1_s(:,iso)*r_s(ispin)%r_coef(:,iso)*grid_atom%wr(:)) 00569 ecoul_1_h = ecoul_1_h + 0.5_dp*SUM(Vh1_h(:,iso)*r_h(ispin)%r_coef(:,iso)*grid_atom%wr(:)) 00570 END DO 00571 END DO 00572 00573 CALL set_ecoul_1c(ecoul_1c,iatom,ecoul_1_z=ecoul_1_z,ecoul_1_0=ecoul_1_0) 00574 CALL set_ecoul_1c(ecoul_1c=ecoul_1c,iatom=iatom,ecoul_1_h=ecoul_1_h,ecoul_1_s=ecoul_1_s) 00575 00576 energy_hartree_1c = energy_hartree_1c + ecoul_1_z - ecoul_1_0 00577 energy_hartree_1c = energy_hartree_1c + ecoul_1_h - ecoul_1_s 00578 00579 ! atomic energy contribution 00580 IF (atenergy) THEN 00581 ate1c(iatom) = ate1c(iatom) + ecoul_1_z - ecoul_1_0 00582 END IF 00583 00584 END SUBROUTINE Vh_1c_atom_energy 00585 00586 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00587 00588 SUBROUTINE Vh_1c_atom_integrals(rho_atom,& 00589 aVh1b_hh,aVh1b_ss,aVh1b_00,Vh1_h,Vh1_s,max_iso_not0,& 00590 max_s_harm,llmax,cg_list,cg_n_list,& 00591 nset,npgf,lmin,lmax,nsotot,maxso,nspins,nchan_0,gsph,g0_h_w,my_CG,Qlm_gg,error) 00592 00593 TYPE(rho_atom_type), POINTER :: rho_atom 00594 REAL(dp), DIMENSION(:, :) :: aVh1b_hh, aVh1b_ss, aVh1b_00 00595 REAL(dp), DIMENSION(:, :), POINTER :: Vh1_h, Vh1_s 00596 INTEGER, INTENT(IN) :: max_iso_not0, max_s_harm, 00597 llmax 00598 INTEGER, DIMENSION(:, :, :) :: cg_list 00599 INTEGER, DIMENSION(:) :: cg_n_list 00600 INTEGER, INTENT(IN) :: nset 00601 INTEGER, DIMENSION(:), POINTER :: npgf, lmin, lmax 00602 INTEGER, INTENT(IN) :: nsotot, maxso, nspins, nchan_0 00603 REAL(dp), DIMENSION(:, :), POINTER :: gsph 00604 REAL(dp), DIMENSION(:, 0:) :: g0_h_w 00605 REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg 00606 TYPE(cp_error_type), INTENT(inout) :: error 00607 00608 CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_atom_integrals', 00609 routineP = moduleN//':'//routineN 00610 00611 INTEGER :: icg, ipgf1, ipgf2, ir, is1, is2, iset1, iset2, iso, iso1, 00612 iso2, ispin, l_ang, m1, m2, max_iso_not0_local, n1, n2, nr 00613 REAL(dp) :: gVg_0, gVg_h, gVg_s 00614 TYPE(rho_atom_coeff), DIMENSION(:), 00615 POINTER :: int_local_h, int_local_s 00616 00617 NULLIFY(int_local_h,int_local_s) 00618 CALL get_rho_atom(rho_atom=rho_atom,& 00619 ga_Vlocal_gb_h=int_local_h,& 00620 ga_Vlocal_gb_s=int_local_s) 00621 00622 ! Calculate the integrals of the potential with 2 primitives 00623 aVh1b_hh = 0.0_dp 00624 aVh1b_ss = 0.0_dp 00625 aVh1b_00 = 0.0_dp 00626 00627 nr = SIZE(gsph,1) 00628 00629 m1 = 0 00630 DO iset1 = 1,nset 00631 m2 = 0 00632 DO iset2 = 1,nset 00633 CALL get_none0_cg_list(my_CG,lmin(iset1),lmax(iset1),lmin(iset2),lmax(iset2),& 00634 max_s_harm,llmax,cg_list,cg_n_list,max_iso_not0_local,error) 00635 00636 n1 = nsoset(lmax(iset1)) 00637 DO ipgf1 = 1,npgf(iset1) 00638 n2 = nsoset(lmax(iset2)) 00639 DO ipgf2 = 1,npgf(iset2) 00640 ! with contributions to V1_s*rho0 00641 DO iso = 1,nchan_0 00642 l_ang = indso(1,iso) 00643 gVg_0 = SUM(Vh1_s(:,iso)*g0_h_w(:,l_ang)) 00644 DO icg = 1,cg_n_list(iso) 00645 is1 = cg_list(1,icg,iso) 00646 is2 = cg_list(2,icg,iso) 00647 00648 iso1 = is1 + n1*(ipgf1-1) + m1 00649 iso2 = is2 + n2*(ipgf2-1) + m2 00650 gVg_h = 0.0_dp 00651 gVg_s = 0.0_dp 00652 00653 DO ir = 1,nr 00654 gVg_h = gVg_h + gsph(ir,iso1)*gsph(ir,iso2)* Vh1_h(ir,iso) 00655 gVg_s = gVg_s + gsph(ir,iso1)*gsph(ir,iso2)* Vh1_s(ir,iso) 00656 END DO ! ir 00657 00658 aVh1b_hh(iso1,iso2) = aVh1b_hh(iso1,iso2) + gVg_h*my_CG(is1,is2,iso) 00659 aVh1b_ss(iso1,iso2) = aVh1b_ss(iso1,iso2) + gVg_s*my_CG(is1,is2,iso) 00660 aVh1b_00(iso1,iso2) = aVh1b_00(iso1,iso2) + gVg_0*Qlm_gg(iso1,iso2,iso) 00661 00662 END DO !icg 00663 END DO ! iso 00664 ! without contributions to V1_s*rho0 00665 DO iso = nchan_0+1,max_iso_not0 00666 DO icg = 1,cg_n_list(iso) 00667 is1 = cg_list(1,icg,iso) 00668 is2 = cg_list(2,icg,iso) 00669 00670 iso1 = is1 + n1*(ipgf1-1) + m1 00671 iso2 = is2 + n2*(ipgf2-1) + m2 00672 gVg_h = 0.0_dp 00673 gVg_s = 0.0_dp 00674 00675 DO ir = 1,nr 00676 gVg_h = gVg_h + gsph(ir,iso1)*gsph(ir,iso2)* Vh1_h(ir,iso) 00677 gVg_s = gVg_s + gsph(ir,iso1)*gsph(ir,iso2)* Vh1_s(ir,iso) 00678 END DO ! ir 00679 00680 aVh1b_hh(iso1,iso2) = aVh1b_hh(iso1,iso2) + gVg_h*my_CG(is1,is2,iso) 00681 aVh1b_ss(iso1,iso2) = aVh1b_ss(iso1,iso2) + gVg_s*my_CG(is1,is2,iso) 00682 00683 END DO !icg 00684 END DO ! iso 00685 END DO ! ipgf2 00686 END DO ! ipgf1 00687 m2 = m2 + maxso 00688 END DO ! iset2 00689 m1 = m1 + maxso 00690 END DO !iset1 00691 00692 DO ispin = 1,nspins 00693 CALL daxpy(nsotot*nsotot,1.0_dp,aVh1b_hh,1,int_local_h(ispin)%r_coef,1) 00694 CALL daxpy(nsotot*nsotot,1.0_dp,aVh1b_ss,1,int_local_s(ispin)%r_coef,1) 00695 CALL daxpy(nsotot*nsotot,-1.0_dp,aVh1b_00,1,int_local_h(ispin)%r_coef,1) 00696 CALL daxpy(nsotot*nsotot,-1.0_dp,aVh1b_00,1,int_local_s(ispin)%r_coef,1) 00697 END DO ! ispin 00698 00699 END SUBROUTINE Vh_1c_atom_integrals 00700 00701 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00702 00703 END MODULE hartree_local_methods 00704
1.7.3