TABLE OF CONTENTS
ABINIT/random_stopping_power [ Functions ]
NAME
random_stopping_power
FUNCTION
Calculate the electronic random stopping power
COPYRIGHT
Copyright (C) 2001-2018 ABINIT group (AS,FB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
screening
CHILDREN
get_bz_item,spline,splint,wrtout
SOURCE
33 #if defined HAVE_CONFIG_H 34 #include "config.h" 35 #endif 36 37 #include "abi_common.h" 38 39 subroutine random_stopping_power(iqibz,npvel,pvelmax,Ep,Gsph_epsG0,Qmesh,Vcp,Cryst,Dtfil,epsm1,rspower) 40 41 use defs_basis 42 use defs_datatypes 43 use defs_abitypes 44 use m_errors 45 use m_splines 46 47 use m_io_tools, only : open_file 48 use m_bz_mesh, only : get_BZ_item, kmesh_t 49 use m_gsphere, only : gsphere_t 50 use m_crystal, only : crystal_t 51 use m_gwdefs, only : GW_TOLQ0, em1params_t 52 use m_vcoul, only : vcoul_t 53 54 !This section has been created automatically by the script Abilint (TD). 55 !Do not modify the following lines by hand. 56 #undef ABI_FUNC 57 #define ABI_FUNC 'random_stopping_power' 58 use interfaces_14_hidewrite 59 !End of the abilint section 60 61 implicit none 62 63 !Arguments ------------------------------------ 64 !scalars 65 integer,intent(in) :: iqibz,npvel 66 real(dp),intent(in) :: pvelmax(3) 67 68 type(em1params_t),intent(in) :: Ep 69 type(gsphere_t),intent(in) :: Gsph_epsG0 70 type(kmesh_t),intent(in) :: Qmesh 71 type(vcoul_t),intent(in) :: Vcp 72 type(crystal_t),intent(in) :: Cryst 73 type(Datafiles_type),intent(in) :: Dtfil 74 75 complex(gwpc),intent(in) :: epsm1(Ep%npwe,Ep%npwe,Ep%nomega) 76 77 real(dp),intent(inout) :: rspower(npvel) 78 79 !Local variables ------------------------------ 80 integer :: ipvel,ig 81 integer :: iq_bz,iq_ibz,isym_q,itim_q 82 integer :: iomega,iomegap,nomega_re 83 integer :: unt_rsp 84 integer,allocatable :: iomega_re(:) 85 86 real(dp),parameter :: zp=1.0_dp ! Hard-coded charge of the impinging particle 87 real(dp) :: omega_p 88 real(dp) :: im_epsm1_int(1) 89 real(dp) :: qbz(3),qpgcart(3),qpgred(3) 90 real(dp) :: pvel(3,npvel),pvel_norm(npvel) 91 real(dp) :: ypp_i(Ep%nomega) 92 real(dp) :: vcoul(Ep%npwe) 93 real(dp),allocatable :: im_epsm1_diag_qbz(:,:),tmp_data(:) 94 real(dp),allocatable :: omega_re(:) 95 96 character(len=500) :: msg 97 character(len=fnlen+4) :: fname 98 99 !************************************************************************ 100 101 ! 102 ! First set up the velocities array from the input variables npvel and pvelmax(3) 103 ! Remember pvelmax is in Cartesian coordinates and so is pvel 104 do ipvel=1,npvel 105 pvel(:,ipvel) = REAL(ipvel,dp) / REAL(npvel,dp) * pvelmax(:) 106 pvel_norm(ipvel) = SQRT( SUM( pvel(:,ipvel)**2 ) ) 107 enddo 108 ! 109 ! Select the purely real frequency in Ep%omega 110 nomega_re=0 111 do iomega=1,Ep%nomega 112 if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then 113 nomega_re=nomega_re+1 114 endif 115 enddo 116 ABI_ALLOCATE(omega_re,(nomega_re)) 117 ABI_ALLOCATE(iomega_re,(nomega_re)) 118 ABI_ALLOCATE(im_epsm1_diag_qbz,(Ep%npwe,Ep%nomega)) 119 ABI_ALLOCATE(tmp_data,(Ep%nomega)) 120 121 iomegap=0 122 do iomega=1,Ep%nomega 123 if( AIMAG(Ep%omega(iomega)) < 1.0e-4_dp ) then 124 iomegap=iomegap+1 125 iomega_re(iomegap)=iomega 126 omega_re(iomegap)=REAL(Ep%omega(iomega),dp) 127 endif 128 enddo 129 130 ! 131 ! Loop over all the q-points in the full Brillouin zone and select only the 132 ! ones that corresponds to the correct q-point in the irreducible wedge we are 133 ! currently treating (index iqibz) 134 do iq_bz=1,Qmesh%nbz 135 136 ! Perform the check and obtain the symmetry information 137 call get_BZ_item(Qmesh,iq_bz,qbz,iq_ibz,isym_q,itim_q) 138 if( iqibz /= iq_ibz ) cycle 139 140 ! Apply the symmetry operation to the diagonal of epsm1 141 do iomega=1,nomega_re 142 do ig=1,Ep%npwe 143 im_epsm1_diag_qbz(Gsph_epsG0%rottb(ig,itim_q,isym_q),iomega)= AIMAG( epsm1(ig,ig,iomega_re(iomega)) ) 144 enddo 145 enddo 146 ! Apply the symmetry operation to the Coulomb interaction 147 do ig=1,Ep%npwe 148 vcoul(Gsph_epsG0%rottb(ig,itim_q,isym_q))=Vcp%vc_sqrt(ig,iqibz)**2 149 enddo 150 151 ! 152 ! Sum over G vectors 153 do ig=1,Ep%npwe 154 ! 155 ! Loop over velocities 156 do ipvel=1,npvel 157 158 qpgred(:) = qbz(:) + Gsph_epsG0%gvec(:,ig) 159 ! Transform q + G from reduced to cartesian with the symmetry operation 160 qpgcart(:) = two_pi * Cryst%gprimd(:,1) * qpgred(1) & 161 & + two_pi * Cryst%gprimd(:,2) * qpgred(2) & 162 & + two_pi * Cryst%gprimd(:,3) * qpgred(3) 163 164 ! omega_p = ( q + G ) . v 165 omega_p = DOT_PRODUCT( qpgcart(:) , pvel(:,ipvel) ) 166 167 ! Check that the calculated frequency omega_p is within the omega 168 ! range of epsm1 and thus that the interpolation will go fine 169 if ( ABS(omega_p) > omega_re(nomega_re) ) then 170 write(msg,'(a,e16.4,2a,e16.4)') ' freqremax is currently ',omega_re(nomega_re),ch10,& 171 & ' increase it to at least ',omega_p 172 MSG_WARNING(msg) 173 endif 174 175 ! Perform the spline interpolation to obtain epsm1 at the desired omega = omega_p 176 tmp_data = im_epsm1_diag_qbz(ig,:) 177 178 call spline( omega_re, tmp_data, nomega_re, 1.0e+32_dp, 1.0e+32_dp, ypp_i) 179 call splint( nomega_re, omega_re, tmp_data, ypp_i, 1, (/ ABS(omega_p) /), im_epsm1_int ) 180 181 ! 182 ! Apply the odd parity of Im epsm1 in omega to recover the causal response function 183 if (omega_p<zero) im_epsm1_int(1)=-im_epsm1_int(1) 184 185 ! 186 ! Calculate 4 * pi / |q+G|**2 * omega_p * Im{ epsm1_GG(q,omega_p) } 187 ! 188 im_epsm1_int(1) = omega_p * vcoul(ig) * im_epsm1_int(1) 189 190 ! Accumulate the final result without the prefactor 191 ! (It will be included at the very end) 192 rspower(ipvel) = rspower(ipvel) + im_epsm1_int(1) 193 194 end do ! end of velocity loop 195 end do ! end G-loop 196 197 enddo ! end of q loop in the full BZ 198 199 ! If it is the last q, write down the result in the main output file and in a 200 ! separate file _RSP (for Random Stopping Power) 201 if (iqibz == Qmesh%nibz ) then 202 203 ! Multiply by the prefactors 204 ! Note that this expression differs from Eq. (3.11) in Campillo PRB 58, 10309 (1998). 205 ! A factor one half is missing in the paper. 206 rspower(:) = - zp**2 / ( Cryst%ucvol * Qmesh%nbz * pvel_norm(:) ) * rspower(:) 207 208 write(msg,'(2a)') ch10,' ==== Random stopping power along Cartesian direction === ' 209 call wrtout(std_out,msg,'COLL') 210 call wrtout(ab_out,msg,'COLL') 211 write(msg,'(a,3(f12.4,2x),a)') ' ==== ',pvelmax(:),'====' 212 call wrtout(std_out,msg,'COLL') 213 call wrtout(ab_out,msg,'COLL') 214 write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) ' 215 call wrtout(std_out,msg,'COLL') 216 call wrtout(ab_out,msg,'COLL') 217 do ipvel=1,npvel 218 write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel) 219 call wrtout(std_out,msg,'COLL') 220 call wrtout(ab_out,msg,'COLL') 221 enddo 222 write(msg,'(2a)') ' ========================================================= ',ch10 223 call wrtout(std_out,msg,'COLL') 224 call wrtout(ab_out,msg,'COLL') 225 226 fname=TRIM(Dtfil%filnam_ds(4))//'_RSP' 227 if (open_file(fname,msg,newunit=unt_rsp,status='unknown',form='formatted') /= 0) then 228 MSG_ERROR(msg) 229 end if 230 231 write(msg,'(a)') '# ==== Random stopping power along Cartesian direction === ' 232 call wrtout(unt_rsp,msg,'COLL') 233 write(msg,'(a,3(f12.4,2x))') '# ==== ',pvelmax(:) 234 call wrtout(unt_rsp,msg,'COLL') 235 write(msg,'(a)') '# |v| (a.u.) , RSP (a.u.) ' 236 call wrtout(unt_rsp,msg,'COLL') 237 do ipvel=1,npvel 238 write(msg,'(f16.8,4x,f16.8)') pvel_norm(ipvel),rspower(ipvel) 239 call wrtout(unt_rsp,msg,'COLL') 240 enddo 241 close(unt_rsp) 242 end if 243 244 ABI_DEALLOCATE(omega_re) 245 ABI_DEALLOCATE(iomega_re) 246 ABI_DEALLOCATE(im_epsm1_diag_qbz) 247 ABI_DEALLOCATE(tmp_data) 248 249 end subroutine random_stopping_power