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