| 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 | !-------- |
|---|
| 171 | subroutine pseudo_write_xml(fname,p) |
|---|
| 172 | use flib_wxml |
|---|
| 173 | |
|---|
| 174 | character(len=*), intent(in) :: fname |
|---|
| 175 | type(pseudopotential_t) :: p |
|---|
| 176 | |
|---|
| 177 | integer :: i |
|---|
| 178 | type(xmlf_t) :: xf |
|---|
| 179 | |
|---|
| 180 | call xml_OpenFile(fname,xf,indent=.true.) |
|---|
| 181 | call xml_AddXMLDeclaration(xf) |
|---|
| 182 | call xml_NewElement(xf,"pseudo") |
|---|
| 183 | call xml_AddAttribute(xf,"version","0.5") |
|---|
| 184 | call xml_NewElement(xf,"header") |
|---|
| 185 | call xml_AddAttribute(xf,"symbol",trim(p%name)) |
|---|
| 186 | call xml_AddAttribute(xf,"zval",trim(str(p%zval))) |
|---|
| 187 | call xml_AddAttribute(xf,"creator",trim(p%creator)) |
|---|
| 188 | call xml_AddAttribute(xf,"date",trim(p%date)) |
|---|
| 189 | call xml_AddAttribute(xf,"flavor",trim(p%flavor)) |
|---|
| 190 | call 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)) |
|---|
| 204 | call xml_EndElement(xf,"header") |
|---|
| 205 | |
|---|
| 206 | call 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))) |
|---|
| 212 | call xml_EndElement(xf,"grid") |
|---|
| 213 | |
|---|
| 214 | call 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 | |
|---|
| 281 | end subroutine pseudo_write_xml |
|---|
| 282 | |
|---|
| 283 | ! |
|---|
| 284 | ! Experimental routine to extract information from "text" field. |
|---|
| 285 | ! and to set up more rational fields. |
|---|
| 286 | ! |
|---|
| 287 | subroutine pseudo_complete(p) |
|---|
| 288 | type(pseudopotential_t), intent(inout) :: p |
|---|
| 289 | |
|---|
| 290 | integer :: i, lmax, l, itext, n, status |
|---|
| 291 | real(dp) :: zup, zdown, ztot, rc_read |
|---|
| 292 | |
|---|
| 293 | p%creator = p%method(1) |
|---|
| 294 | p%date = p%method(2) |
|---|
| 295 | p%flavor = p%method(3) // p%method(4) // p%method(5) // p%method(6) |
|---|
| 296 | |
|---|
| 297 | lmax = 0 |
|---|
| 298 | do i = 1, p%npotd |
|---|
| 299 | lmax = max(lmax,p%ldown(i)) |
|---|
| 300 | enddo |
|---|
| 301 | p%lmax = lmax |
|---|
| 302 | allocate(p%principal_n(0:lmax)) |
|---|
| 303 | allocate(p%occupation(0:lmax)) |
|---|
| 304 | allocate(p%cutoff(0:lmax)) |
|---|
| 305 | ! |
|---|
| 306 | ! Decode text into useful information. Assumes l's are increasing from 0 |
|---|
| 307 | ! |
|---|
| 308 | if (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 |
|---|
| 321 | else |
|---|
| 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 | |
|---|
| 332 | endif |
|---|
| 333 | |
|---|
| 334 | end subroutine pseudo_complete |
|---|
| 335 | |
|---|
| 336 | |
|---|
| 337 | ! ---------------------------------------------------------------------- |
|---|
| 338 | subroutine get_unit(lun,iostat) |
|---|
| 339 | |
|---|
| 340 | ! Get an available Fortran unit number |
|---|
| 341 | |
|---|
| 342 | integer, intent(out) :: lun |
|---|
| 343 | integer, intent(out) :: iostat |
|---|
| 344 | |
|---|
| 345 | integer :: i |
|---|
| 346 | logical :: unit_used |
|---|
| 347 | |
|---|
| 348 | do 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 |
|---|
| 355 | enddo |
|---|
| 356 | iostat = -1 |
|---|
| 357 | lun = -1 |
|---|
| 358 | end subroutine get_unit |
|---|
| 359 | |
|---|
| 360 | end module m_pseudo_utils |
|---|
| 361 | |
|---|
| 362 | |
|---|
| 363 | |
|---|
