3
! Reads a psf or vps file and attempts to produce an ATOM inp file from it
5
use pseudopotential, only: pseudopotential_t
6
use pseudopotential, only: pseudo_read, get_ps_conf
10
integer, parameter :: dp = selected_real_kind(16,100)
12
character(len=20) :: label
13
type(pseudopotential_t) :: p
15
integer :: n_val_first, ncore, lmax, inp_lun, l, n
17
character(len=1) :: code_val_first
18
character(len=1) :: ispp
19
character(len=2) :: job
20
character(len=3) :: flavor = "tm2" ! for now
22
character(len=*), parameter :: &
23
ruler = "#23456789012345678901234567890123456789012345678901234567890 Ruler"
25
character(len=2), allocatable :: orb_arr(:)
26
real(dp), allocatable :: zdown_arr(:)
27
real(dp), allocatable :: zup_arr(:)
28
real(dp), allocatable :: rc_arr(:)
30
real(dp) :: aa, bb, rmax
33
call pseudo_read(label,p)
36
allocate(orb_arr(0:lmax))
37
allocate(zdown_arr(0:lmax))
38
allocate(zup_arr(0:lmax))
39
allocate(rc_arr(0:lmax))
41
call get_ps_conf(p%irel,lmax,p%text,p%gen_zval, &
42
orb_arr,zdown_arr,zup_arr,rc_arr)
44
read(orb_arr(0),"(i1)") n_val_first
45
read(orb_arr(0),"(1x,a1)") code_val_first
46
print *, "First valence orbital: ", orb_arr(0)
47
print *, "Pieces: ", n_val_first, code_val_first
48
if (code_val_first /= "s") call die("First valence pseudo is not s!")
50
call core_below_orb(lmax,orb_arr,ncore)
51
print *, "There are ", ncore, " core orbitals"
54
! call io_assign(inp_lun)
56
open(inp_lun, file=trim(label)//".inp", form="formatted", &
57
status="replace",action="write",position="rewind")
59
if (p%nicore(1:2) == "pc" ) then
64
write(inp_lun,"(3x,a2,a50)") job , " -- file generated from " // trim(label) // " ps file"
65
write(inp_lun,"(8x,a3)") flavor
74
write(inp_lun,"(3x,a2,3x,a2,a1)") p%name, p%icorr, ispp
78
write(inp_lun,"(6f10.3)") 0.0, 0.0, 0.0, aa, bb, rmax
79
write(inp_lun,"(2i5)") ncore, lmax+1
82
read(orb_arr(l),"(i1)") n
83
write(inp_lun,"(2i5,2f10.3,2x,a5)") n, l, zdown_arr(l), zup_arr(l), " #" // orb_arr(l)
86
write(inp_lun,"(f10.5)",advance="no") rc_arr(l)
89
write(inp_lun,"(f10.5)",advance="no") 0.0
93
write(inp_lun,"(2f10.5,a)") 0.0, -1.0, " FIX CORE RADIUS"
95
write(inp_lun,"(2f10.5,a)") 0.0, 0.0
99
! Final blank line and ruler
100
write(inp_lun,"(/,a)") ruler
102
! call io_close(inp_lun)
107
subroutine core_below_orb(lmax,orb_arr,ncore)
108
integer, intent(in) :: lmax
109
character(len=2), intent(in) :: orb_arr(0:lmax)
110
integer, intent(out) :: ncore
113
character(len=1) :: symbol
119
! 1s,2s,2p,3s,3p,3d,4s,4p,4d,5s,5p,4f,5d,6s,6p,
122
integer, parameter :: nshells = 19
123
integer, dimension(nshells) :: nc = (/ 1, 2, 2, 3, 3, 3, 4, 4, 4, &
124
5, 5, 4, 5, 6, 6, 7, 5, 6, 7 /)
125
integer, dimension(nshells) :: lc = (/ 0, 0, 1, 0, 1, 2, 0, 1, 2, &
126
0, 1, 3, 2, 0, 1, 0, 3, 2, 1 /)
127
character(len=1), dimension(4) :: sym = (/ "s", "p", "d", "f" /)
133
! Search for orbital among the valence ones
136
read(orb_arr(l),"(i1)") n
137
read(orb_arr(l),"(1x,a1)") symbol
138
if ((nc(i) == n .AND. sym(lc(i)+1) == symbol)) then
147
! print *, "-- ", nc(i), lc(i)
153
if (ncore == nshells) call die("Major snafu in ncore_below_orb")
154
!print *, "Number of core orbs: ", ncore
156
end subroutine core_below_orb
158
end program make_ps_input