TABLE OF CONTENTS


ABINIT/vso_realspace_local [ Functions ]

[ Top ] [ 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