~jose-soler/siesta/unfolding

« back to all changes in this revision

Viewing changes to Pseudo/atom/Util/Ps2Inp/make_ps_input.f90

  • Committer: Alberto Garcia
  • Date: 2016-01-25 16:00:16 UTC
  • mto: (483.3.1 rel-4.0)
  • mto: This revision was merged to the branch mainline in revision 485.
  • Revision ID: albertog@icmab.es-20160125160016-c1qivg1zw06kgyfd
Prepare GPL release

* Include proper headers

* Add Docs/Contributors.txt and NOTICE.txt files.

* Update READMEs and LICENSE files in several directories.

* Remove Pseudo/atom, Util/test-xml

* Remove DOM files from Src/xmlparser

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
program make_ps_input
2
 
!
3
 
! Reads a psf or vps file and attempts to produce an ATOM inp file from it
4
 
!
5
 
  use pseudopotential, only: pseudopotential_t
6
 
  use pseudopotential, only: pseudo_read, get_ps_conf
7
 
 
8
 
  implicit none
9
 
 
10
 
      integer, parameter :: dp = selected_real_kind(16,100)
11
 
 
12
 
  character(len=20) :: label
13
 
  type(pseudopotential_t) :: p
14
 
 
15
 
  integer   :: n_val_first, ncore, lmax, inp_lun, l, n
16
 
 
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
21
 
 
22
 
  character(len=*), parameter ::   &
23
 
       ruler = "#23456789012345678901234567890123456789012345678901234567890      Ruler"
24
 
 
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(:)
29
 
 
30
 
  real(dp) :: aa, bb, rmax
31
 
 
32
 
  read(5,*) label
33
 
  call pseudo_read(label,p)
34
 
 
35
 
  lmax = p%npotd-1
36
 
  allocate(orb_arr(0:lmax))
37
 
  allocate(zdown_arr(0:lmax))
38
 
  allocate(zup_arr(0:lmax))
39
 
  allocate(rc_arr(0:lmax))
40
 
 
41
 
  call get_ps_conf(p%irel,lmax,p%text,p%gen_zval, &
42
 
                   orb_arr,zdown_arr,zup_arr,rc_arr)
43
 
 
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!")
49
 
 
50
 
  call core_below_orb(lmax,orb_arr,ncore)
51
 
  print *, "There are ", ncore, " core orbitals"
52
 
 
53
 
!
54
 
!  call io_assign(inp_lun)
55
 
  inp_lun = 2
56
 
  open(inp_lun, file=trim(label)//".inp", form="formatted",  &
57
 
      status="replace",action="write",position="rewind")
58
 
 
59
 
  if (p%nicore(1:2) == "pc" ) then
60
 
     job = "pe"
61
 
  else
62
 
     job = "pg"
63
 
  endif
64
 
  write(inp_lun,"(3x,a2,a50)") job , " -- file generated from " // trim(label) // " ps file"
65
 
  write(inp_lun,"(8x,a3)") flavor
66
 
  select case (p%irel)
67
 
     case ("isp") 
68
 
        ispp = "s"
69
 
     case ("rel")
70
 
        ispp = "r"
71
 
     case default
72
 
        ispp = " "
73
 
  end select
74
 
  write(inp_lun,"(3x,a2,3x,a2,a1)") p%name, p%icorr, ispp
75
 
  aa = 0.0
76
 
  bb = 0.0
77
 
  rmax = 0.0
78
 
  write(inp_lun,"(6f10.3)") 0.0, 0.0, 0.0, aa, bb, rmax
79
 
  write(inp_lun,"(2i5)") ncore, lmax+1
80
 
 
81
 
  do l = 0, lmax
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)
84
 
  end do
85
 
  do l = 0, lmax
86
 
     write(inp_lun,"(f10.5)",advance="no") rc_arr(l)
87
 
  enddo
88
 
  do l = lmax + 1, 3
89
 
     write(inp_lun,"(f10.5)",advance="no") 0.0
90
 
  end do
91
 
 
92
 
  if (job == "pe") then
93
 
     write(inp_lun,"(2f10.5,a)") 0.0, -1.0, " FIX CORE RADIUS"
94
 
  else
95
 
     write(inp_lun,"(2f10.5,a)") 0.0, 0.0
96
 
  endif
97
 
 
98
 
!
99
 
! Final blank line and ruler
100
 
  write(inp_lun,"(/,a)") ruler
101
 
 
102
 
!  call io_close(inp_lun)
103
 
  close(inp_lun)
104
 
 
105
 
CONTAINS
106
 
 
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
111
 
 
112
 
integer :: l, i
113
 
character(len=1) :: symbol
114
 
logical  :: found
115
 
 
116
 
 
117
 
!  data for orbitals:                                                                                
118
 
!                                                                                                    
119
 
!              1s,2s,2p,3s,3p,3d,4s,4p,4d,5s,5p,4f,5d,6s,6p,                                         
120
 
!                     7s, 5f, 6d, 7p                                                                 
121
 
!                                                                                                    
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" /)
128
 
 
129
 
 
130
 
ncore = 0
131
 
do i = 1, nshells
132
 
 
133
 
! Search for orbital among the valence ones
134
 
  found = .false.
135
 
  do l = 0, lmax
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
139
 
        found = .true.
140
 
        exit
141
 
     endif
142
 
  enddo
143
 
 
144
 
  if (found) then
145
 
     exit
146
 
  else
147
 
!     print *, "-- ", nc(i), lc(i)
148
 
     ncore = ncore + 1
149
 
  endif
150
 
 
151
 
enddo
152
 
 
153
 
if (ncore == nshells) call die("Major snafu in ncore_below_orb")
154
 
!print *, "Number of core orbs: ", ncore
155
 
 
156
 
end subroutine core_below_orb
157
 
  
158
 
end program make_ps_input