1
2! Copyright (C) 2014 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
3! This file is distributed under the terms of the GNU General Public License.
4! See the file COPYING for license details.
5
6subroutine genfspecies(zn,symb)
7use modmain
8use modmpi
9implicit none
10! arguments
11real(8), intent(in) :: zn
12character(*), intent(in) :: symb
13! local variables
14integer, parameter :: nit=4
15integer nst,ist,jst
16integer nmax,in,il,ik
17integer nrm,nr,ir,it
18integer n(maxstsp),l(maxstsp),k(maxstsp)
19integer idx(maxstsp),iv(maxstsp)
20real(8) rm,rmin,rmax
21real(8) mass,t1,t2,t3
22real(8) occ(maxstsp),eval(maxstsp),rv(maxstsp)
23character(256) name
24! allocatable arrays
25real(8), allocatable :: r(:),rho(:),vr(:),rwf(:,:,:)
26! external functions
27real(8), external :: massnucl
28name='Fractional species'
29! set up the initial occupancies
30occ(:)=0.d0
31t1=abs(zn)
32nmax=1
33ist=0
34do in=1,maxstsp
35  do il=0,in-1
36    do ik=max(il,1),il+1
37      t2=dble(2*ik)
38      t2=min(t2,t1)
39      ist=ist+1
40      n(ist)=in
41      l(ist)=il
42      k(ist)=ik
43      occ(ist)=t2
44      if (t2.gt.epsocc) nmax=in
45      t1=t1-t2
46      if (ist.eq.maxstsp) then
47        if (t1.gt.epsocc) then
48          write(*,*)
49          write(*,'("Error(genfspecies): too many states for fractional &
50           &species ",A)') trim(symb)
51          write(*,*)
52          stop
53        else
54          goto 10
55        end if
56      end if
57    end do
58  end do
59end do
6010 continue
61! minimum radius
62rmin=2.d-6/sqrt(abs(zn))
63! initial maximum radius
64rmax=100.d0
65! initial muffin-tin radius
66rm=2.d0
67! number of points to muffin-tin radius
68nrm=100*(nmax+1)
69! iterate the solution but not to self-consistency
70do it=1,nit
71! number of points to effective infinity
72  t1=log(rm/rmin)
73  t2=log(rmax/rmin)
74  t3=dble(nrm)*t2/t1
75  nr=int(t3)
76  allocate(r(nr),rho(nr),vr(nr),rwf(nr,2,maxstsp))
77! generate logarithmic radial mesh
78  t2=t1/dble(nrm-1)
79  do ir=1,nr
80    r(ir)=rmin*exp(dble(ir-1)*t2)
81  end do
82! solve the Kohn-Sham-Dirac equation for the atom
83  call atom(sol,.true.,zn,maxstsp,n,l,k,occ,3,0,nr,r,eval,rho,vr,rwf)
84! check for spurious eigenvalues
85  do ist=2,maxstsp
86    if (eval(ist).lt.eval(1)) eval(ist)=1.d6
87  end do
88! recompute the effective infinity
89  do ir=nr,1,-1
90    if (rho(ir).gt.1.d-20) then
91      rmax=1.75d0*r(ir)
92      exit
93    end if
94  end do
95! estimate the muffin-tin radius
96  do ir=nr,1,-1
97    if (rho(ir).gt.2.d-2) then
98      rm=r(ir)
99      exit
100    end if
101  end do
102  if (rm.lt.1.d0) rm=1.d0
103  if (rm.gt.3.2d0) rm=3.2d0
104! sort the eigenvalues
105  call sortidx(maxstsp,eval,idx)
106! recompute the occupancies
107  occ(:)=0.d0
108  t1=abs(zn)
109  do ist=1,maxstsp
110    jst=idx(ist)
111    ik=k(jst)
112    t2=dble(2*ik)
113    t2=min(t2,t1)
114    occ(jst)=t2
115    t1=t1-t2
116  end do
117  deallocate(r,rho,vr,rwf)
118end do
119! rearrange the arrays
120iv(:)=n(:)
121n(:)=iv(idx(:))
122iv(:)=l(:)
123l(:)=iv(idx(:))
124iv(:)=k(:)
125k(:)=iv(idx(:))
126rv(:)=occ(:)
127occ(:)=rv(idx(:))
128rv(:)=eval(:)
129eval(:)=rv(idx(:))
130! find the number of occupied states
131nst=0
132do ist=1,maxstsp
133  if (occ(ist).lt.epsocc) then
134    nst=ist
135    exit
136  end if
137end do
138! estimate the nuclear mass
139mass=massnucl(zn)
140! convert from 'atomic mass units' to atomic units
141mass=mass*amu
142! write the species file
143call writespecies(symb,name,zn,mass,rmin,rm,rmax,nrm,nst,n,l,k,occ,eval)
144if (mp_mpi) then
145  write(*,*)
146  write(*,'("Info(genfspecies): wrote fractional species file ",A,".in")') &
147   trim(symb)
148end if
149end subroutine
150
151