src/species/species.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 program species
   7 implicit none
   8 ! local variables
   9 ! order of predictor-corrector polynomial
  10 integer, parameter :: np=4
  11 ! maximum number of states
  12 integer, parameter :: maxspst=40
  13 ! maximum angular momentum allowed
  14 integer, parameter :: lmax=50
  15 ! exchange-correlation type
  16 integer, parameter :: xctype=3
  17 integer, parameter :: xcgrad=0
  18 integer nz,spnst,spnr,nrmt
  19 integer nlx,nlorb,i,l,maxl
  20 integer ist,jst,ir,iostat
  21 real(8), parameter :: pi=3.1415926535897932385d0
  22 ! core-valence cut-off energy
  23 real(8), parameter :: ecvcut=-3.5d0
  24 ! semi-core cut-off energy
  25 real(8), parameter :: esccut=-0.35d0
  26 ! band offset energy
  27 real(8), parameter :: boe=0.15d0
  28 real(8) spmass,rmt,spzn,sprmin,sprmax,t1
  29 character(256) spsymb,spname
  30 ! automatic arrays
  31 logical spcore(maxspst)
  32 integer spn(maxspst),spl(maxspst),spk(maxspst)
  33 real(8) spocc(maxspst),eval(maxspst)
  34 ! allocatable arrays
  35 real(8), allocatable :: r(:),rho(:),vr(:),rwf(:,:,:)
  36 real(8), allocatable :: fr(:),gr(:),cf(:,:)
  37 open(40,file='species.dat',action='READ',status='OLD',form='FORMATTED')
  38 10 continue
  39 read(40,*,iostat=iostat) nz
  40 if (iostat.ne.0) stop
  41 read(40,*) spsymb
  42 read(40,*) spname
  43 read(40,*) spmass
  44 read(40,*) rmt
  45 read(40,*) spnst
  46 if (spnst.gt.maxspst) then
  47   write(*,*)
  48   write(*,'("Error(species): too many states for species ",A)') trim(spname)
  49   write(*,*)
  50   stop
  51 end if
  52 do ist=1,spnst
  53   read(40,*) spn(ist),spl(ist),spk(ist),i
  54   if (ist.ge.2) then
  55     if (spn(ist).lt.spn(ist-1)) then
  56       write(*,*)
  57       write(*,'("Error(species): states improperly ordered")')
  58       write(*,'(" for species ",A)') trim(spname)
  59       write(*,*)
  60       stop
  61     end if
  62   end if
  63   spocc(ist)=dble(i)
  64 end do
  65 read(40,*)
  66 write(*,'("Info(species): running Z = ",I4,", (",A,")")') nz,trim(spname)
  67 ! nuclear charge in units of e
  68 spzn=-dble(nz)
  69 ! minimum radial mesh point proportional to 1/sqrt(Z)
  70 sprmin=2.d-6/sqrt(abs(dble(spzn)))
  71 ! set the number of radial mesh points proportional to number of nodes
  72 nrmt=100*(spn(spnst)+1)
  73 ! find the optimal effective infinity
  74 sprmax=80.d0
  75 do i=1,2
  76   t1=log(sprmax/sprmin)/log(rmt/sprmin)
  77   spnr=int(t1*dble(nrmt))
  78   if (allocated(r)) deallocate(r)
  79   if (allocated(rho)) deallocate(rho)
  80   if (allocated(vr)) deallocate(vr)
  81   if (allocated(rwf)) deallocate(rwf)
  82   if (allocated(fr)) deallocate(fr)
  83   if (allocated(gr)) deallocate(gr)
  84   if (allocated(cf)) deallocate(cf)
  85   allocate(r(spnr))
  86   allocate(rho(spnr))
  87   allocate(vr(spnr))
  88   allocate(rwf(spnr,2,spnst))
  89   allocate(fr(spnr))
  90   allocate(gr(spnr))
  91   allocate(cf(3,spnr))
  92 ! generate the radial mesh
  93   call radmesh(spnr,nrmt,rmt,sprmin,r)
  94 ! solve the Kohn-Sham-Dirac equations for the atom
  95   call atom(.true.,spzn,spnst,spn,spl,spk,spocc,xctype,xcgrad,np,spnr,r,eval, &
  96    rho,vr,rwf)
  97   do ir=spnr,1,-1
  98     if (rho(ir).gt.1.d-20) then
  99       sprmax=1.5d0*r(ir)
 100       goto 20
 101     end if
 102   end do
 103 20 continue
 104 end do
 105 ! check total charge is correct
 106 do ir=1,spnr
 107   fr(ir)=4.d0*pi*rho(ir)*r(ir)**2
 108 end do
 109 call fderiv(-1,spnr,r,fr,gr,cf)
 110 if (abs(gr(spnr)+spzn).gt.1.d-5) then
 111   write(*,*)
 112   write(*,'("Error(species): charge mismatch")')
 113   write(*,*)
 114   stop
 115 end if
 116 ! find which states belong to core
 117 do ist=1,spnst
 118   if (eval(ist).lt.ecvcut) then
 119     spcore(ist)=.true.
 120   else
 121     spcore(ist)=.false.
 122   end if
 123 end do
 124 ! check that the state for same n and l but different k is also core
 125 do ist=1,spnst
 126   if (spcore(ist)) then
 127     do jst=1,spnst
 128       if ((spn(ist).eq.spn(jst)).and.(spl(ist).eq.spl(jst))) spcore(jst)=.true.
 129     end do
 130   end if
 131 end do
 132 ! find the total number of local orbitals
 133 nlorb=0
 134 maxl=0
 135 do ist=1,spnst
 136   if (.not.spcore(ist)) then
 137     if ((spl(ist).eq.0).or.(spl(ist).eq.spk(ist))) then
 138       if (eval(ist).lt.esccut) nlorb=nlorb+1
 139     end if
 140     if (spl(ist).gt.maxl) maxl=spl(ist)
 141   end if
 142 end do
 143 maxl=maxl+1
 144 if (maxl.gt.3) maxl=3
 145 nlorb=nlorb+maxl+1
 146 nlx=0
 147 ! open the atomic data file
 148 open(50,file=trim(spsymb)//'.in',action='WRITE',form='FORMATTED')
 149 write(50,'(" ''",A,"''",T45,": spsymb")') trim(spsymb)
 150 write(50,'(" ''",A,"''",T45,": spname")') trim(spname)
 151 write(50,'(G14.6,T45,": spzn")') spzn
 152 write(50,'(G18.10,T45,": spmass")') spmass
 153 write(50,'(G14.6,2F10.4,I6,T45,": sprmin, rmt, sprmax, nrmt")') sprmin,rmt, &
 154  sprmax,nrmt
 155 write(50,'(I4,T45,": spnst")') spnst
 156 write(50,'(3I4,G14.6,L1,T45,": spn, spl, spk, spocc, spcore")') spn(1),spl(1), &
 157  spk(1),spocc(1),spcore(1)
 158 do ist=2,spnst
 159   write(50,'(3I4,G14.6,L1)') spn(ist),spl(ist),spk(ist),spocc(ist),spcore(ist)
 160 end do
 161 write(50,'(I4,T45,": apword")') 1
 162 write(50,'(F8.4,I4,"  ",L1,T45,": apwe0, apwdm, apwve")') boe,0,.false.
 163 write(50,'(I4,T45,": nlx")') nlx
 164 write(50,'(I4,T45,": nlorb")') nlorb
 165 ! write the local-orbitals
 166 do l=0,maxl
 167   write(50,'(2I4,T45,": lorbl, lorbord")') l,2
 168   write(50,'(F8.4,I4,"  ",L1,T45,": lorbe0, lorbdm, lorbve")') boe,0,.false.
 169   write(50,'(F8.4,I4,"  ",L1)') boe,1,.false.
 170 end do
 171 do ist=1,spnst
 172   if (.not.spcore(ist)) then
 173     if ((spl(ist).eq.0).or.(spl(ist).eq.spk(ist))) then
 174       if (eval(ist).lt.esccut) then
 175         write(50,'(2I4,T45,": lorbl, lorbord")') spl(ist),3
 176         write(50,'(F8.4,I4,"  ",L1,T45,": lorbe0, lorbdm, lorbve")') boe,0, &
 177          .false.
 178         write(50,'(F8.4,I4,"  ",L1)') boe,1,.false.
 179         write(50,'(F8.4,I4,"  ",L1)') eval(ist)+0.5d0*boe,0,.true.
 180       end if
 181     end if
 182   end if
 183 end do
 184 close(50)
 185 ! read another element from file
 186 goto 10
 187 end program