src/mossbauer.f90

   1
   2 ! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: mossbauer
   8 ! !INTERFACE:
   9 subroutine mossbauer
  10 ! !USES:
  11 use modmain
  12 ! !DESCRIPTION:
  13 !   Computes the contact charge density and contact magnetic hyperfine field for
  14 !   each atom and outputs the data to the file {\tt MOSSBAUER.OUT}. The nuclear
  15 !   radius used for the contact quantities is approximated by the empirical
  16 !   formula $R_{\rm N}=1.25 Z^{1/3}$ fm, where $Z$ is the atomic number.
  17 !
  18 ! !REVISION HISTORY:
  19 !   Created May 2004 (JKD)
  20 !EOP
  21 !BOC
  22 implicit none
  23 ! local variables
  24 integer is,ia,ias,ir,nr
  25 ! nuclear radius constant in Bohr
  26 real(8), parameter :: r0=1.25d-15/0.52917720859d-10
  27 real(8) rn,vn,rho0,b,t1
  28 ! allocatable arrays
  29 real(8), allocatable :: fr(:)
  30 real(8), allocatable :: gr(:)
  31 real(8), allocatable :: cf(:,:)
  32 ! initialise universal variables
  33 call init0
  34 ! read density and potentials from file
  35 call readstate
  36 ! allocate local arrays
  37 allocate(fr(nrmtmax))
  38 allocate(gr(nrmtmax))
  39 allocate(cf(3,nrmtmax))
  40 open(50,file='MOSSBAUER.OUT',action='WRITE',form='FORMATTED')
  41 do is=1,nspecies
  42 !--------------------------------!
  43 !     contact charge density     !
  44 !--------------------------------!
  45 ! approximate nuclear radius : r0*A^(1/3)
  46   rn=r0*abs(spzn(is))**(1.d0/3.d0)
  47   do ir=1,nrmt(is)
  48     if (spr(ir,is).gt.rn) goto 10
  49   end do
  50   write(*,*)
  51   write(*,'("Error(mossbauer): nuclear radius too large : ",G18.10)') rn
  52   write(*,'(" for species ",I4)') is
  53   write(*,*)
  54   stop
  55 10 continue
  56   nr=ir
  57   rn=spr(nr,is)
  58 ! nuclear volume
  59   vn=(4.d0/3.d0)*pi*rn**3
  60   do ia=1,natoms(is)
  61     ias=idxas(ia,is)
  62 !--------------------------------!
  63 !     contact charge density     !
  64 !--------------------------------!
  65     fr(1:nr)=rhomt(1,1:nr,ias)*y00
  66     do ir=1,nr
  67       fr(ir)=(fourpi*spr(ir,is)**2)*fr(ir)
  68     end do
  69     call fderiv(-1,nr,spr(:,is),fr,gr,cf)
  70     rho0=gr(nr)/vn
  71     write(50,*)
  72     write(50,'("Species : ",I4," (",A,"), atom : ",I4)') is,trim(spsymb(is)),ia
  73     write(50,'(" approximate nuclear radius : ",G18.10)') rn
  74     write(50,'(" number of mesh points to nuclear radius : ",I6)') nr
  75     write(50,'(" contact charge density : ",G18.10)') rho0
  76 !------------------------------------------!
  77 !     contact magnetic hyperfine field     !
  78 !------------------------------------------!
  79     if (spinpol) then
  80       do ir=1,nr
  81         if (ncmag) then
  82 ! non-collinear
  83           t1=sqrt(magmt(1,ir,ias,1)**2+magmt(1,ir,ias,2)**2 &
  84            +magmt(1,ir,ias,3)**2)
  85         else
  86 ! collinear
  87           t1=magmt(1,ir,ias,1)
  88         end if
  89         fr(ir)=t1*y00*fourpi*spr(ir,is)**2
  90       end do
  91       call fderiv(-1,nr,spr(:,is),fr,gr,cf)
  92       b=gr(nr)/vn
  93       write(50,'(" contact magnetic hyperfine field (mu_B) : ",G18.10)') b
  94     end if
  95   end do
  96 end do
  97 close(50)
  98 write(*,*)
  99 write(*,'("Info(mossbauer):")')
 100 write(*,'(" Mossbauer parameters written to MOSSBAUER.OUT")')
 101 write(*,*)
 102 deallocate(fr,gr,cf)
 103 return
 104 end subroutine
 105 !EOC