TABLE OF CONTENTS
ABINIT/vso_realspace_local [ Functions ]
NAME
vso_realspace_local
FUNCTION
Calculate real space (local - (r,r)) values of the SO part of the pseudopotential. Reconstructed explicitly in the HGH/GTH case.
COPYRIGHT
Copyright (C) 2005-2018 ABINIT group (Mver) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
spin_current
CHILDREN
gamma_function,spline,splint,xred2xcart
SOURCE
33 #if defined HAVE_CONFIG_H 34 #include "config.h" 35 #endif 36 37 #include "abi_common.h" 38 39 subroutine vso_realspace_local(dtset,hdr,position_op,psps,vso_realspace) 40 41 use m_profiling_abi 42 43 use defs_basis 44 use defs_datatypes 45 use defs_abitypes 46 use m_splines 47 48 use m_geometry, only : xred2xcart 49 50 !This section has been created automatically by the script Abilint (TD). 51 !Do not modify the following lines by hand. 52 #undef ABI_FUNC 53 #define ABI_FUNC 'vso_realspace_local' 54 use interfaces_28_numeric_noabirule 55 !End of the abilint section 56 57 implicit none 58 59 !Arguments ------------------------------- 60 61 type(hdr_type),intent(inout) :: hdr 62 type(dataset_type),intent(in) :: dtset 63 type(pseudopotential_type),intent(in) :: psps 64 65 real(dp),intent(in) :: position_op(3,dtset%ngfft(1),dtset%ngfft(2),dtset%ngfft(3)) 66 67 real(dp),intent(out) :: vso_realspace(2,dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3),& 68 & dtset%nspinor,dtset%nspinor,3) 69 70 !Local variables ------------------------- 71 72 integer :: i,j,l, lmax,ipsp,iatom, ir1,ir2,ir3 73 integer :: rcexponent,irealsp 74 integer :: nradgrid,iradgrid 75 76 real(dp) :: gammai, gammaj, relative_position(3), radial_cutoff, norm_rel_pos 77 real(dp) :: expfact,lfact, vso_interpol, x,y,z 78 real(dp) :: xcart(3,dtset%natom),splint_x(1),splint_y(1) 79 80 real(dp), allocatable :: radial_grid(:) 81 real(dp), allocatable :: prefact_ijl(:,:,:,:),tmpvso(:),tmpvso_pp(:) 82 real(dp), allocatable :: vso_radial(:,:),vso_radial_pp(:,:),tmp_spline(:) 83 real(dp), allocatable :: offdiag_l_fact(:,:,:),kpar_matrix(:,:) 84 85 ! ********************************************************************* 86 87 !recalculate xcart (option = 1) 88 call xred2xcart(dtset%natom,hdr%rprimd,xcart,hdr%xred) 89 90 91 lmax = psps%mpsang-1 92 93 !content of gth pseudo type: 94 !These are {rloc, C(1...4)} coefficients for psppar(0, :, :) indices, 95 !Followed by the h coefficients for psppar(1:2, 1:, :) indices. 96 !size (0:2, 0:4, npsp) 97 !potential radius r_l is in psppar(l+1,0,ipsp) 98 !real(dp), pointer :: psppar(:, :, :) 99 !The covalence radii for each pseudo (?) size (npsp) 100 !real(dp), pointer :: radii_cov(:) 101 !Cut-off radii for core part and long-range part. 102 !radii_cf(:, 1) is for the long-range cut-off and 103 !radii_cf(:, 2) is for the core cut-off. 104 !size (npsp, 2) 105 !real(dp), pointer :: radii_cf(:, :) 106 !Spin orbit coefficients in HGH/GTH formats: k11p 107 !etc... see psp3ini.F90 108 !dimension = num l channels, 3 coeffs, num psp = 109 !(1:lmax+1,1:3,npsp) 110 !real(dp), pointer :: psp_k_par(:, :, :) 111 112 !v_SO^l (r,r') = sum_i sum_j sum_m Y_{lm} (\hat{r}) p_i^l (r) k_{ij}^l p_j^l(r') Y^{*}_lm (\hat{r'}) 113 ! 114 !v_SO^l (r,r) = sum_ij p_i^l (r) k_{ij}^l p_j^l(r) sum_m Y_{lm} (\hat{r}) Y^{*}_lm (\hat{r}) 115 != (2l+1)/4\pi sum_ij p_i^l (r) k_{ij}^l p_j^l(r) (eq B.17 Patrick Rinke thesis) 116 !p are gaussian projectors (from HGH paper prb 58 3641) 117 !sum_l v_SO^l (r,r) is a purely radial quantity (function of |r|), so spline it 118 119 !maximum distance needed in unit cell 120 radial_cutoff = four * maxval(psps%gth_params%psppar(:, 0, :)) 121 122 !setup radial grid; Should we use a logarithmic grid? The spline functions can 123 !take it... 124 nradgrid = 201 ! this is heuristic 125 ABI_ALLOCATE(radial_grid,(nradgrid)) 126 do iradgrid=1,nradgrid 127 radial_grid(iradgrid) = (iradgrid-1)*radial_cutoff/(nradgrid-1) 128 end do 129 130 !calculate prefactors independent of r 131 ABI_ALLOCATE(prefact_ijl,(3,3,0:lmax,psps%npsp)) 132 ABI_ALLOCATE(offdiag_l_fact,(3,3,0:lmax)) 133 ABI_ALLOCATE(kpar_matrix,(3,3)) 134 135 !these factors complete the full 3x3 matrix of k (or h) parameters for the 136 !HGH pseudos 137 offdiag_l_fact = zero 138 !l=0 139 offdiag_l_fact(1,2,0) = -half*sqrt(three/five) 140 offdiag_l_fact(1,3,0) = half*sqrt(five/21._dp) 141 offdiag_l_fact(2,3,0) = -half*sqrt(100._dp/63._dp) 142 !l=1 143 offdiag_l_fact(1,2,1) = -half*sqrt(five/seven) 144 offdiag_l_fact(1,3,1) = sixth*sqrt(35._dp/11._dp) 145 offdiag_l_fact(2,3,1) = -sixth*14._dp/sqrt(11._dp) 146 !l=2 147 if (lmax >= 2) then 148 offdiag_l_fact(1,2,2) = -half*sqrt(seven/nine) 149 offdiag_l_fact(1,3,2) = half*sqrt(63._dp/143._dp) 150 offdiag_l_fact(2,3,2) = -half*18._dp /sqrt(143._dp) 151 end if 152 !l=3 153 if (lmax >= 3) then 154 offdiag_l_fact(1,2,3) = zero 155 offdiag_l_fact(1,3,3) = zero 156 offdiag_l_fact(2,3,3) = zero 157 end if 158 !get prefactors for evaluation of V_SO: terms that do not depend on r 159 prefact_ijl = zero 160 do l=0,lmax 161 ! first the diagonal i=j term 162 do i=1,3 163 call gamma_function(l+(4._dp*i-1._dp)*0.5_dp, gammai) 164 gammai = sqrt(gammai) 165 rcexponent = 2*l+2*i+2*i-1 166 do ipsp=1,psps%npsp 167 prefact_ijl(i,i,l,ipsp) = psps%gth_params%psp_k_par(l+1,i,ipsp) & 168 & / ( (psps%gth_params%psppar(l+1,0,ipsp))**(rcexponent) & 169 & * gammai * gammai) 170 end do 171 end do 172 ! now the off diagonal elements 173 do ipsp=1,psps%npsp 174 kpar_matrix(1,2) = offdiag_l_fact (1,2,l)* psps%gth_params%psp_k_par(l+1,2,ipsp) 175 kpar_matrix(2,1) = kpar_matrix(1,2) 176 kpar_matrix(1,3) = offdiag_l_fact (1,3,l)* psps%gth_params%psp_k_par(l+1,3,ipsp) 177 kpar_matrix(3,1) = kpar_matrix(1,3) 178 kpar_matrix(2,3) = offdiag_l_fact (2,3,l)* psps%gth_params%psp_k_par(l+1,3,ipsp) 179 kpar_matrix(3,2) = kpar_matrix(2,3) 180 end do 181 182 ! for the f case only the 1,1 matrix element is non 0 - it is done above and 183 ! all these terms are actually 0 184 if (l > 2) cycle 185 186 do i=1,3 187 call gamma_function(l+(4._dp*i-1._dp)*0.5_dp, gammai) 188 gammai = sqrt(gammai) 189 do j=1,3 190 if (j==i) cycle 191 rcexponent = 2*l+2*i+2*j-1 192 call gamma_function(l+(4._dp*j-1._dp)*0.5_dp,gammaj) 193 gammaj = sqrt(gammaj) 194 do ipsp=1,psps%npsp 195 prefact_ijl(i,j,l,ipsp) = kpar_matrix(i,j) & 196 & / ( (psps%gth_params%psppar(l+1,0,ipsp))**rcexponent & 197 & * gammai * gammaj ) 198 end do 199 end do 200 end do 201 end do 202 203 ABI_DEALLOCATE(kpar_matrix) 204 ABI_DEALLOCATE(offdiag_l_fact) 205 206 prefact_ijl = prefact_ijl * two 207 208 !calculate v_SO on radial grid 209 ! MGNAG Runtime Error: *** Arithmetic exception: Floating invalid operation - aborting 210 ABI_ALLOCATE(vso_radial,(nradgrid,psps%npsp)) 211 vso_radial = zero 212 do l=0,lmax 213 lfact=(2._dp*l+1._dp)/four/pi 214 do iradgrid=1,nradgrid 215 norm_rel_pos = radial_grid(iradgrid) 216 do ipsp=1,psps%npsp 217 expfact = exp(-norm_rel_pos**2 / & 218 & (psps%gth_params%psppar(l+1,0,ipsp))**2) 219 220 do i=1,3 221 do j=1,3 222 rcexponent = 2*l +2*i+2*j-4 223 if(prefact_ijl(i,j,l,ipsp)/=0) then !vz_d 0**0 224 vso_radial(iradgrid,ipsp) = vso_radial(iradgrid,ipsp) + & 225 & prefact_ijl(i,j,l,ipsp)*(norm_rel_pos**rcexponent) * expfact 226 end if !vz_d 227 end do ! j 228 end do ! i 229 end do ! ipsp 230 end do ! iradgrid 231 end do ! lmax 232 233 !spline v_SO(radial coord): get second derivative coefficients 234 ABI_ALLOCATE(vso_radial_pp,(nradgrid,psps%npsp)) 235 236 ABI_ALLOCATE(tmp_spline,(nradgrid)) 237 ABI_ALLOCATE(tmpvso,(nradgrid)) 238 ABI_ALLOCATE(tmpvso_pp,(nradgrid)) 239 do ipsp=1,psps%npsp 240 tmpvso = vso_radial(:,ipsp) 241 call spline( radial_grid, tmpvso, nradgrid, zero, radial_grid(nradgrid), tmpvso_pp ) 242 vso_radial_pp(:,ipsp) = tmpvso_pp 243 end do 244 ABI_DEALLOCATE(tmp_spline) 245 ABI_DEALLOCATE(tmpvso) 246 ABI_DEALLOCATE(tmpvso_pp) 247 248 !to optimize this I should precalculate the distances which are actually needed by 249 !symmetry, or only sum over irreducible points in space and use weights 250 251 !for each physical atom present in unit cell 252 vso_realspace = zero 253 do iatom=1,dtset%natom 254 ! atom type will be dtset%typat(iatom) 255 256 ! for each point on grid 257 do ir3=1,dtset%ngfft(3) 258 do ir2=1,dtset%ngfft(2) 259 do ir1=1,dtset%ngfft(1) 260 irealsp = ir1 + (ir2-1)*dtset%ngfft(1) + (ir3-1)*dtset%ngfft(2)*dtset%ngfft(1) 261 262 ! relative position from atom to point 263 relative_position = position_op(:,ir1,ir2,ir3) - xcart(:,iatom) 264 x=relative_position(1) 265 y=relative_position(2) 266 z=relative_position(3) 267 268 ! calculate norm^2 269 norm_rel_pos = relative_position(1)**2+relative_position(2)**2+relative_position(3)**2 270 271 ! if norm^2 is too large, skip this point 272 if (norm_rel_pos > radial_cutoff*radial_cutoff) cycle 273 274 ! calculate norm 275 splint_x(1) = sqrt(norm_rel_pos) 276 277 ! spline interpolate vso only depends on position (through pos - atomic position) 278 call splint (nradgrid,radial_grid,vso_radial(:,dtset%typat(iatom)),& 279 & vso_radial_pp(:,dtset%typat(iatom)),1,splint_x,splint_y) 280 vso_interpol=splint_y(1) 281 282 ! multiply by vectorial spin factor (S x r) 283 ! NOTE: this r is taken relative to atom center. It could be that the r operator should 284 ! applied in an absolute way wrt the origin... 285 ! 286 ! Is this correct: accumulated sum over atoms ? 287 vso_realspace(1,irealsp,:,:,1) = vso_realspace(1,irealsp,:,:,1) + & 288 & vso_interpol * reshape((/y, zero,zero,-y/),(/2,2/)) 289 vso_realspace(2,irealsp,:,:,1) = vso_realspace(2,irealsp,:,:,1) + & 290 & vso_interpol * reshape((/zero,z, -z, zero/),(/2,2/)) 291 292 vso_realspace(1,irealsp,:,:,2) = vso_realspace(1,irealsp,:,:,2) + & 293 & vso_interpol * reshape((/-x, z, z, x/),(/2,2/)) 294 vso_realspace(2,irealsp,:,:,2) = vso_realspace(2,irealsp,:,:,2) + & 295 & vso_interpol * reshape((/zero,zero,zero,zero/),(/2,2/)) 296 297 vso_realspace(1,irealsp,:,:,3) = vso_realspace(1,irealsp,:,:,3) + & 298 & vso_interpol * reshape((/zero,-y, -y, zero/),(/2,2/)) 299 vso_realspace(2,irealsp,:,:,3) = vso_realspace(2,irealsp,:,:,3) + & 300 & vso_interpol * reshape((/zero,-x, -x, zero/),(/2,2/)) 301 302 303 end do ! ir3 304 end do ! ir2 305 end do ! ir1 306 end do ! iatom 307 308 ABI_DEALLOCATE(prefact_ijl) 309 ABI_DEALLOCATE(vso_radial) 310 ABI_DEALLOCATE(vso_radial_pp) 311 ABI_DEALLOCATE(radial_grid) 312 313 end subroutine vso_realspace_local