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.2 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 real(8) t1,rn,vn,rho0,b
  26 ! allocatable arrays
  27 real(8), allocatable :: fr(:)
  28 real(8), allocatable :: gr(:)
  29 real(8), allocatable :: cf(:,:)
  30 ! initialise universal variables
  31 call init0
  32 ! read density and potentials from file
  33 call readstate
  34 ! allocate local arrays
  35 allocate(fr(nrmtmax))
  36 allocate(gr(nrmtmax))
  37 allocate(cf(3,nrmtmax))
  38 open(50,file='MOSSBAUER.OUT',action='WRITE',form='FORMATTED')
  39 do is=1,nspecies
  40 !--------------------------------!
  41 !     contact charge density     !
  42 !--------------------------------!
  43 ! determine approximate nuclear radius : 1.2*A^(1/3) fm
  44   t1=1.2d-15/0.5291772108d-10
  45   rn=t1*abs(spzn(is))**(1.d0/3.d0)
  46   do ir=1,nrmt(is)
  47     if (spr(ir,is).gt.rn) goto 10
  48   end do
  49   write(*,*)
  50   write(*,'("Error(mossbauer): nuclear radius too large : ",G18.10)') rn
  51   write(*,'(" for species ",I4)') is
  52   write(*,*)
  53   stop
  54 10 continue
  55   nr=ir
  56   rn=spr(nr,is)
  57 ! nuclear volume
  58   vn=(4.d0/3.d0)*pi*rn**3
  59   do ia=1,natoms(is)
  60     ias=idxas(ia,is)
  61 !--------------------------------!
  62 !     contact charge density     !
  63 !--------------------------------!
  64     fr(1:nr)=rhomt(1,1:nr,ias)*y00
  65     do ir=1,nr
  66       fr(ir)=(fourpi*spr(ir,is)**2)*fr(ir)
  67     end do
  68     call fderiv(-1,nr,spr(1,is),fr,gr,cf)
  69     rho0=gr(nr)/vn
  70     write(50,*)
  71     write(50,'("Species : ",I4," (",A,"), atom : ",I4)') is,trim(spsymb(is)),ia
  72     write(50,'(" approximate nuclear radius : ",G18.10)') rn
  73     write(50,'(" number of mesh points to nuclear radius : ",I6)') nr
  74     write(50,'(" contact charge density : ",G18.10)') rho0
  75 !------------------------------------------!
  76 !     contact magnetic hyperfine field     !
  77 !------------------------------------------!
  78     if (spinpol) then
  79       do ir=1,nr
  80         if (ndmag.eq.3) then
  81           t1=sqrt(magmt(1,ir,ias,1)**2+magmt(1,ir,ias,2)**2 &
  82            +magmt(1,ir,ias,3)**2)
  83         else
  84           t1=magmt(1,ir,ias,1)
  85         end if
  86         fr(ir)=t1*y00*fourpi*spr(ir,is)**2
  87       end do
  88       call fderiv(-1,nr,spr(1,is),fr,gr,cf)
  89       b=gr(nr)/vn
  90       write(50,'(" contact magnetic hyperfine field (mu_B) : ",G18.10)') b
  91     end if
  92   end do
  93 end do
  94 close(50)
  95 write(*,*)
  96 write(*,'("Info(mossbauer):")')
  97 write(*,'(" Mossbauer parameters written to MOSSBAUER.OUT")')
  98 write(*,*)
  99 deallocate(fr,gr,cf)
 100 return
 101 end subroutine
 102 !EOC