1
2      module m_pseudo_utils
3!
4!     Main module for ps input and output, based on
5!     a derived type representation closely resembling
6!     the Froyen data structures.
7!
8!
9!     The radial coordinate is reparametrized to allow a more
10!     precise sampling of the area near the origin.
11
12!      r: 0->...
13!      x: 0->...
14!      i: 1->...
15
16!           r(x) = grid_scale *  [ exp( grid_step*x) - 1 ]
17!
18!     **** WARNING *****
19!     In SIESTA, grid_scale = b, grid_step = a
20!     In ATOM,   grid_scale = a, grid_step = b
21!     ******************
22!
23!     pseudo%nr and pseudo%nrval are identical
24!     (for ATOM and SIESTA use)
25!
26!     Working precision should be double precision
27!     for backwards binary compatibility
28!
29      private
30
31      public :: pseudo_read_formatted
32      public :: pseudo_header_print
33      public :: pseudo_write_xml
34      public :: pseudo_complete
35      private :: get_unit
36
37      integer, private, parameter :: dp = selected_real_kind(14)
38
39      type, public  ::  pseudopotential_t
40        character(len=2)        :: name
41        integer                 :: nr
42        integer                 :: nrval
43        real(dp)                :: zval
44        logical                 :: relativistic
45        character(len=2)        :: icorr
46        character(len=3)        :: irel
47        character(len=4)        :: nicore
48        real(dp)                :: grid_scale
49        real(dp)                :: grid_step
50        character(len=10)       :: method(6)
51        character(len=70)       :: text
52        integer                 :: npotu
53        integer                 :: npotd
54        real(dp), pointer       :: r(:)
55        real(dp), pointer       :: chcore(:)
56        real(dp), pointer       :: chval(:)
57        real(dp), pointer       :: vdown(:,:)
58        real(dp), pointer       :: vup(:,:)
59        integer, pointer        :: ldown(:)
60        integer, pointer        :: lup(:)
61!
62!       Extra fields for more functionality
63!
64        character(len=10)       :: creator
65        character(len=10)       :: date
66        character(len=40)       :: flavor
67        integer                 :: lmax
68        integer, pointer        :: principal_n(:)
69        real(dp), pointer       :: occupation(:)
70        real(dp), pointer       :: cutoff(:)
71      end type pseudopotential_t
72!
73!     These determine the format for ASCII files
74!
75      character(len=*), parameter, private ::  &
76               fmt_int = "(tr1,i2)" ,                  &
77               fmt_nam = "(tr1,a2,tr1,a2,tr1,a3,tr1,a4)", &
78               fmt_met = "(tr1,6a10,/,tr1,a70)" ,       &
79               fmt_pot= "(tr1,2i3,i5,3es20.12)" ,      &
80               fmt_rad = "(4(es20.12))"        ,      &
81               fmt_txt = "(tr1,a)"
82
83        CONTAINS
84
85!----
86        subroutine pseudo_read_formatted(fname,p)
87        character(len=*), intent(in) :: fname
88        type(pseudopotential_t)                    :: p
89
90        integer  :: io_ps, i, j, status
91        character(len=70) :: dummy
92        real(kind=dp)  :: r2
93
94        call get_unit(io_ps,status)
95        if (status /= 0) stop "cannot get unit number"
96        open(unit=io_ps,file=fname,form="formatted",status="old",&
97             action="read",position="rewind")
98        write(6,"(3a)") "Reading pseudopotential information ", &
99            "in formatted form from ", trim(fname)
100
101
102        read(io_ps,fmt=fmt_nam) p%name, p%icorr, p%irel, p%nicore
103        read(io_ps,fmt_met) (p%method(i),i=1,6), p%text
104        read(io_ps,fmt=fmt_pot) p%npotd, p%npotu, p%nr, &
105                                p%grid_scale, p%grid_step, p%zval
106
107        p%nrval = p%nr + 1
108        p%nr = p%nr + 1
109        allocate(p%r(1:p%nrval))
110        read(io_ps,fmt=fmt_txt) dummy
111        read(io_ps,fmt=fmt_rad) (p%r(j),j=2,p%nrval)
112        p%r(1) = 0.d0
113        r2=p%r(2)/(p%r(3)-p%r(2))
114
115        if (p%npotd.gt.0) then
116           allocate(p%vdown(1:p%npotd,1:p%nrval))
117           allocate(p%ldown(1:p%npotd))
118        endif
119        do i=1,p%npotd
120           read(io_ps,fmt=fmt_txt) dummy
121           read(io_ps,fmt=fmt_int) p%ldown(i)
122           read(io_ps,fmt=fmt_rad) (p%vdown(i,j), j=2,p%nrval)
123           p%vdown(i,1) = p%vdown(i,2) - r2*(p%vdown(i,3)-p%vdown(i,2))
124        enddo
125
126        if (p%npotu.gt.0) then
127           allocate(p%vup(1:p%npotu,1:p%nrval))
128           allocate(p%lup(1:p%npotu))
129        endif
130        do i=1,p%npotu
131           read(io_ps,fmt=fmt_txt) dummy
132           read(io_ps,fmt=fmt_int) p%lup(i)
133           read(io_ps,fmt=fmt_rad) (p%vup(i,j), j=2,p%nrval)
134           p%vup(i,1) = p%vup(i,2) - r2*(p%vup(i,3)-p%vup(i,2))
135        enddo
136
137        allocate(p%chcore(1:p%nrval))
138        allocate(p%chval(1:p%nrval))
139
140        read(io_ps,fmt=fmt_txt) dummy
141        read(io_ps,fmt=fmt_rad) (p%chcore(j),j=2,p%nrval)
142        p%chcore(1) = p%chcore(2) - r2*(p%chcore(3)-p%chcore(2))
143
144        read(io_ps,fmt=fmt_txt) dummy
145        read(io_ps,fmt=fmt_rad) (p%chval(j),j=2,p%nrval)
146        p%chval(1) = p%chval(2) - r2*(p%chval(3)-p%chval(2))
147
148        close(io_ps)
149        end subroutine pseudo_read_formatted
150!------
151
152        subroutine vps_init(p)
153        type(pseudopotential_t)  :: p
154        nullify(p%lup,p%ldown,p%r,p%chcore,p%chval,p%vdown,p%vup)
155        end subroutine vps_init
156
157!-------
158        subroutine pseudo_header_print(lun,p)
159        integer, intent(in) :: lun
160        type(pseudopotential_t)  :: p
161
162        integer :: i
163
164        write(lun,"(a)") "<pseudopotential_header>"
165        write(lun,fmt=fmt_nam) p%name, p%icorr, p%irel, p%nicore
166        write(lun,fmt_met) (p%method(i),i=1,6), p%text
167        write(lun,"(a)") "</pseudopotential_header>"
168
169        end subroutine pseudo_header_print
170!--------
171subroutine pseudo_write_xml(fname,p)
172use xmlf90_wxml
173
174character(len=*), intent(in) :: fname
175type(pseudopotential_t)      :: p
176
177integer  :: i
178type(xmlf_t) :: xf
179
180call xml_OpenFile(fname,xf,indent=.true.)
181call xml_AddXMLDeclaration(xf)
182call xml_NewElement(xf,"pseudo")
183call xml_AddAttribute(xf,"version","0.5")
184call xml_NewElement(xf,"header")
185call xml_AddAttribute(xf,"symbol",trim(p%name))
186!call xml_AddAttribute(xf,"zval",trim(str(p%zval)))
187call xml_AddAttribute(xf,"zval",p%zval) ! Overloaded
188call xml_AddAttribute(xf,"creator",trim(p%creator))
189call xml_AddAttribute(xf,"date",trim(p%date))
190call xml_AddAttribute(xf,"flavor",trim(p%flavor))
191call xml_AddAttribute(xf,"correlation",trim(p%icorr))
192
193      select case (trim(p%irel))
194         case ("isp")
195            call xml_AddAttribute(xf,"relativistic","no")
196            call xml_AddAttribute(xf,"polarized","yes")
197         case ("nrl")
198            call xml_AddAttribute(xf,"relativistic","no")
199            call xml_AddAttribute(xf,"polarized","no")
200         case ("rel")
201            call xml_AddAttribute(xf,"relativistic","yes")
202            call xml_AddAttribute(xf,"polarized","no")
203      end select
204      call xml_AddAttribute(xf,"core-corrections",trim(p%nicore))
205call xml_EndElement(xf,"header")
206
207call xml_NewElement(xf,"grid")
208      call xml_AddAttribute(xf,"type","log")
209      call xml_AddAttribute(xf,"units","bohr")
210      call xml_AddAttribute(xf,"scale",p%grid_scale)
211      call xml_AddAttribute(xf,"step",p%grid_step)
212      call xml_AddAttribute(xf,"npts",p%nr-1)
213      call xml_EndElement(xf,"grid")
214
215call xml_NewElement(xf,"semilocal")
216
217      call xml_AddAttribute(xf,"units","rydberg")
218      call xml_AddAttribute(xf,"format","r*V")
219      call xml_AddAttribute(xf,"npots-down",p%npotd)
220      call xml_AddAttribute(xf,"npots-up",p%npotu)
221
222      do i=1,p%npotd
223         call xml_NewElement(xf,"vps")
224         call xml_AddAttribute(xf,"principal-n", &
225                p%principal_n(p%ldown(i)))
226         call xml_AddAttribute(xf,"l",p%ldown(i))
227         call xml_AddAttribute(xf,"cutoff", &
228                p%cutoff(p%ldown(i)))
229         call xml_AddAttribute(xf,"occupation", &
230                p%occupation(p%ldown(i)))
231         call xml_AddAttribute(xf,"spin","-1")
232
233         call xml_NewElement(xf,"radfunc")
234         call xml_NewElement(xf,"data")
235         call xml_AddArray(xf,p%vdown(i,2:p%nr),fmt_rad)
236         call xml_EndElement(xf,"data")
237         call xml_EndElement(xf,"radfunc")
238         call xml_EndElement(xf,"vps")
239      enddo
240
241      do i=1,p%npotu
242         call xml_NewElement(xf,"vps")
243         call xml_AddAttribute(xf,"principal-n", &
244                trim(str(p%principal_n(p%lup(i)))))
245         call xml_AddAttribute(xf,"l",trim(str(p%lup(i))))
246         call xml_AddAttribute(xf,"cutoff", &
247                trim(str(p%cutoff(p%lup(i)))))
248         call xml_AddAttribute(xf,"occupation", &
249                trim(str(p%occupation(p%lup(i)))))
250         call xml_AddAttribute(xf,"spin","-1")
251
252         call xml_NewElement(xf,"radfunc")
253         call xml_NewElement(xf,"data")
254         call xml_AddArray(xf,p%vup(i,2:p%nr),fmt_rad)
255         call xml_EndElement(xf,"data")
256         call xml_EndElement(xf,"radfunc")
257         call xml_EndElement(xf,"vps")
258      enddo
259
260      call xml_EndElement(xf,"semilocal")
261
262      call xml_NewElement(xf,"valence-charge")
263         call xml_NewElement(xf,"radfunc")
264         call xml_NewElement(xf,"data")
265         call xml_AddArray(xf,p%chval(2:p%nr),fmt_rad)
266         call xml_EndElement(xf,"data")
267         call xml_EndElement(xf,"radfunc")
268         call xml_EndElement(xf,"valence-charge")
269
270      call xml_NewElement(xf,"pseudocore-charge")
271         call xml_NewElement(xf,"radfunc")
272         call xml_NewElement(xf,"data")
273         call xml_AddArray(xf,p%chcore(2:p%nr),fmt_rad)
274         call xml_EndElement(xf,"data")
275         call xml_EndElement(xf,"radfunc")
276         call xml_EndElement(xf,"pseudocore-charge")
277
278   call xml_EndElement(xf,"pseudo")
279
280   call xml_Close(xf)
281
282end subroutine pseudo_write_xml
283
284!
285!  Experimental routine to extract information from "text" field.
286!  and to set up more rational fields.
287!
288subroutine pseudo_complete(p)
289type(pseudopotential_t), intent(inout) :: p
290
291integer  :: i, lmax, l, itext, n, status
292real(dp) :: zup, zdown, ztot, rc_read
293
294p%creator = p%method(1)
295p%date = p%method(2)
296p%flavor = p%method(3) // p%method(4) // p%method(5) // p%method(6)
297
298lmax = 0
299do i = 1, p%npotd
300   lmax = max(lmax,p%ldown(i))
301enddo
302p%lmax = lmax
303allocate(p%principal_n(0:lmax))
304allocate(p%occupation(0:lmax))
305allocate(p%cutoff(0:lmax))
306!
307! Decode text into useful information. Assumes l's are increasing from 0
308!
309if (p%irel=="isp") then
310      print *, "Polarized........*************"
311      print *, "|", p%text, "|"
312      do l=0,min(lmax,3)
313         itext=l*17
314         read(p%text(itext+1:),iostat=status, &
315                     fmt="(i1,tr1,f4.2,tr1,f4.2,tr1,f4.2)") &
316                     n, zdown, zup, rc_read
317         if (status /=0) STOP "fallo text"
318         p%principal_n(l) = n
319         p%occupation(l) = zdown+zup
320         p%cutoff(l) = rc_read
321      enddo
322else
323      do l=0,min(lmax,3)
324         itext=l*17
325         read(p%text(itext+1:),iostat=status,fmt="(i1,tr1,f5.2,tr4,f5.2)") &
326                      n, ztot, rc_read
327         if (status /=0) STOP "fallo text"
328         p%principal_n(l) = n
329         p%occupation(l) = ztot
330         p%cutoff(l) = rc_read
331      enddo
332
333endif
334
335end subroutine pseudo_complete
336
337
338! ----------------------------------------------------------------------
339subroutine get_unit(lun,iostat)
340
341! Get an available Fortran unit number
342
343integer, intent(out)  :: lun
344integer, intent(out)  :: iostat
345
346integer :: i
347logical :: unit_used
348
349do i = 10, 99
350   lun = i
351   inquire(unit=lun,opened=unit_used)
352   if (.not. unit_used) then
353      iostat = 0
354      return
355   endif
356enddo
357iostat = -1
358lun = -1
359end subroutine get_unit
360
361end module m_pseudo_utils
362
363
364
365