source: XMLF90/doc/Examples/wxml/i.m_pseudo_utils.f90 @ 6

Last change on this file since 6 was 6, checked in by ymipsl, 13 years ago

Import des sources XMLF90

File size: 11.3 KB
Line 
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 flib_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))
186call xml_AddAttribute(xf,"zval",trim(str(p%zval)))
187call xml_AddAttribute(xf,"creator",trim(p%creator))
188call xml_AddAttribute(xf,"date",trim(p%date))
189call xml_AddAttribute(xf,"flavor",trim(p%flavor))
190call xml_AddAttribute(xf,"correlation",trim(p%icorr))
191
192      select case (trim(p%irel))
193         case ("isp") 
194            call xml_AddAttribute(xf,"relativistic","no")
195            call xml_AddAttribute(xf,"polarized","yes")
196         case ("nrl")
197            call xml_AddAttribute(xf,"relativistic","no")
198            call xml_AddAttribute(xf,"polarized","no")
199         case ("rel")
200            call xml_AddAttribute(xf,"relativistic","yes")
201            call xml_AddAttribute(xf,"polarized","no")
202      end select
203      call xml_AddAttribute(xf,"core-corrections",trim(p%nicore))
204call xml_EndElement(xf,"header")
205
206call xml_NewElement(xf,"grid")
207      call xml_AddAttribute(xf,"type","log")
208      call xml_AddAttribute(xf,"units","bohr")
209      call xml_AddAttribute(xf,"scale",trim(str(p%grid_scale)))
210      call xml_AddAttribute(xf,"step",trim(str(p%grid_step)))
211      call xml_AddAttribute(xf,"npts",trim(str(p%nr-1)))
212call xml_EndElement(xf,"grid")
213
214call xml_NewElement(xf,"semilocal")
215
216      call xml_AddAttribute(xf,"units","rydberg")
217      call xml_AddAttribute(xf,"format","r*V")
218      call xml_AddAttribute(xf,"npots-down",trim(str(p%npotd)))
219      call xml_AddAttribute(xf,"npots-up",trim(str(p%npotu)))
220
221      do i=1,p%npotd
222         call xml_NewElement(xf,"vps")
223         call xml_AddAttribute(xf,"principal-n", &
224                trim(str(p%principal_n(p%ldown(i)))))
225         call xml_AddAttribute(xf,"l",trim(str(p%ldown(i))))
226         call xml_AddAttribute(xf,"cutoff", &
227                trim(str(p%cutoff(p%ldown(i)))))
228         call xml_AddAttribute(xf,"occupation", &
229                trim(str(p%occupation(p%ldown(i)))))
230         call xml_AddAttribute(xf,"spin","-1")
231
232         call xml_NewElement(xf,"radfunc")
233         call xml_NewElement(xf,"data")
234         call xml_AddArray(xf,p%vdown(i,2:p%nr),fmt_rad)
235         call xml_EndElement(xf,"data")
236         call xml_EndElement(xf,"radfunc")
237         call xml_EndElement(xf,"vps")
238      enddo
239
240      do i=1,p%npotu
241         call xml_NewElement(xf,"vps")
242         call xml_AddAttribute(xf,"principal-n", &
243                trim(str(p%principal_n(p%lup(i)))))
244         call xml_AddAttribute(xf,"l",trim(str(p%lup(i))))
245         call xml_AddAttribute(xf,"cutoff", &
246                trim(str(p%cutoff(p%lup(i)))))
247         call xml_AddAttribute(xf,"occupation", &
248                trim(str(p%occupation(p%lup(i)))))
249         call xml_AddAttribute(xf,"spin","-1")
250
251         call xml_NewElement(xf,"radfunc")
252         call xml_NewElement(xf,"data")
253         call xml_AddArray(xf,p%vup(i,2:p%nr),fmt_rad)
254         call xml_EndElement(xf,"data")
255         call xml_EndElement(xf,"radfunc")
256         call xml_EndElement(xf,"vps")
257      enddo
258
259      call xml_EndElement(xf,"semilocal")
260     
261      call xml_NewElement(xf,"valence-charge")
262         call xml_NewElement(xf,"radfunc")
263         call xml_NewElement(xf,"data")
264         call xml_AddArray(xf,p%chval(2:p%nr),fmt_rad)
265         call xml_EndElement(xf,"data")
266         call xml_EndElement(xf,"radfunc")
267         call xml_EndElement(xf,"valence-charge")
268
269      call xml_NewElement(xf,"pseudocore-charge")
270         call xml_NewElement(xf,"radfunc")
271         call xml_NewElement(xf,"data")
272         call xml_AddArray(xf,p%chcore(2:p%nr),fmt_rad)
273         call xml_EndElement(xf,"data")
274         call xml_EndElement(xf,"radfunc")
275         call xml_EndElement(xf,"pseudocore-charge")
276
277   call xml_EndElement(xf,"pseudo")
278
279   call xml_Close(xf)
280
281end subroutine pseudo_write_xml
282
283!
284!  Experimental routine to extract information from "text" field.
285!  and to set up more rational fields.
286!
287subroutine pseudo_complete(p)
288type(pseudopotential_t), intent(inout) :: p
289
290integer  :: i, lmax, l, itext, n, status
291real(dp) :: zup, zdown, ztot, rc_read
292
293p%creator = p%method(1)
294p%date = p%method(2)
295p%flavor = p%method(3) // p%method(4) // p%method(5) // p%method(6)
296
297lmax = 0
298do i = 1, p%npotd
299   lmax = max(lmax,p%ldown(i))
300enddo
301p%lmax = lmax
302allocate(p%principal_n(0:lmax))
303allocate(p%occupation(0:lmax))
304allocate(p%cutoff(0:lmax))
305!
306! Decode text into useful information. Assumes l's are increasing from 0
307!
308if (p%irel=="isp") then
309      print *, "Polarized........*************"
310      print *, "|", p%text, "|"
311      do l=0,min(lmax,3)
312         itext=l*17
313         read(p%text(itext+1:),iostat=status, &
314                     fmt="(i1,tr1,f4.2,tr1,f4.2,tr1,f4.2)") &
315                     n, zdown, zup, rc_read
316         if (status /=0) STOP "fallo text"
317         p%principal_n(l) = n
318         p%occupation(l) = zdown+zup
319         p%cutoff(l) = rc_read
320      enddo
321else
322      do l=0,min(lmax,3)
323         itext=l*17
324         read(p%text(itext+1:),iostat=status,fmt="(i1,tr1,f5.2,tr4,f5.2)") &
325                      n, ztot, rc_read
326         if (status /=0) STOP "fallo text"
327         p%principal_n(l) = n
328         p%occupation(l) = ztot
329         p%cutoff(l) = rc_read
330      enddo
331
332endif
333
334end subroutine pseudo_complete
335
336
337! ----------------------------------------------------------------------
338subroutine get_unit(lun,iostat)
339
340! Get an available Fortran unit number
341
342integer, intent(out)  :: lun
343integer, intent(out)  :: iostat
344
345integer :: i
346logical :: unit_used
347
348do i = 10, 99
349   lun = i
350   inquire(unit=lun,opened=unit_used)
351   if (.not. unit_used) then
352      iostat = 0
353      return
354   endif
355enddo
356iostat = -1
357lun = -1
358end subroutine get_unit
359
360end module m_pseudo_utils
361
362
363
Note: See TracBrowser for help on using the repository browser.