1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8#:include 'common.fypp' 9 10!> Contains type for representing the data stored in the old SK-file format and subroutines to read 11!> that data from file. 12module dftbp_oldskdata 13 use dftbp_assert 14 use dftbp_accuracy 15 use dftbp_constants 16 use dftbp_repspline, only : TRepSplineIn 17 use dftbp_reppoly, only : TRepPolyIn 18 use dftbp_message 19 use dftbp_rangeseparated, only : TRangeSepSKTag 20 implicit none 21 private 22 23 public :: TOldSKData, readFromFile, readSplineRep 24 25 26 !> Represents the Slater-Koster data in an SK file. 27 type TOldSKData 28 29 !> Grid separation 30 real(dp) :: dist 31 32 !> Nr. of grid points 33 integer :: nGrid 34 35 !> Atomic eigenvalues. 36 real(dp) :: skSelf(4) 37 38 !> Hubbard Us 39 real(dp) :: skHubbU(4) 40 41 !> Occupations 42 real(dp) :: skOcc(4) 43 44 !> Mass of the atom 45 real(dp) :: mass 46 47 !> Table for H 48 real(dp), allocatable :: skHam(:,:) 49 50 !> Table for S 51 real(dp), allocatable :: skOver(:,:) 52 end type TOldSKData 53 54 55 !> Reads the data from an SK-file. 56 interface readFromFile 57 module procedure OldSKData_readFromFile 58 end interface readFromFile 59 60 61 !> Reads the repulsive from an open file 62 interface readSplineRep 63 module procedure OldSKData_readSplineRep 64 end interface readSplineRep 65 66 67 !> nr. of interactions in the SK-file 68 integer, parameter :: nSKInter = 20 69 70 !> nr. of ints. in old SK-file 71 integer, parameter :: nSKInterOld = 10 72 73 74 !> Mapping between old (spd) and new (spdf) interactions in the SK-table 75 integer, parameter :: iSKInterOld(nSKInterOld) & 76 & = (/8, 9, 10, 13, 14, 15, 16, 18, 19, 20/) 77 78contains 79 80 81 !> Reads the data from an SK-file. 82 subroutine OldSKData_readFromFile(skData, fileName, homo, iSp1, iSp2, repSplineIn, repPolyIn,& 83 & rangeSepSK) 84 85 !> Contains the content of the SK-file on exit 86 type(TOldSKData), intent(out) :: skData 87 88 !> Name of the file to read the data from 89 character(len=*), intent(in) :: fileName 90 91 !> Is it a homonuclear SK-file? 92 logical, intent(in) :: homo 93 94 !> Index of 1st interacting species (for error messages only) 95 integer, intent(in), optional :: iSp1 96 97 !> Index of 2nd interacting species (for error messages only) 98 integer, intent(in), optional :: iSp2 99 100 !> Repulsive spline part of the SK-file. 101 type(TRepSplineIn), intent(out), optional :: repSplineIn 102 103 !> Repulsive polynomial part of the SK-file. 104 type(TRepPolyIn), intent(out), optional :: repPolyIn 105 106 !> Reads rangeseparation parameter from SK file 107 type(TRangeSepSKTag), intent(inout), optional :: rangeSepSK 108 109 integer :: file 110 character(lc) :: chDummy 111 112 !> extended format for f orbitals 113 logical :: tExtended 114 115 integer :: nShell 116 integer :: ii, iGrid 117 real(dp) :: rDummy 118 real(dp) :: coeffs(2:9), polyCutoff 119 integer :: iostat 120 121 @:ASSERT(present(repSplineIn) .eqv. present(iSp1)) 122 @:ASSERT(present(iSp1) .eqv. present(iSp2)) 123 124 open(newunit=file, file=fileName, status="old", action="read", iostat=iostat) 125 call checkIoError(iostat, fileName, "Unable to open file") 126 rewind(file) 127 128 read (file, '(A1)', iostat=iostat) chDummy 129 call checkIoError(iostat, fileName, "Unable to read 1st line") 130 if (chDummy == '@') then 131 tExtended = .true. 132 nShell = 4 133 else 134 tExtended = .false. 135 nShell = 3 136 rewind(file) 137 end if 138 139 read (file,*, iostat=iostat) skData%dist, skData%nGrid 140 call checkIoError(iostat, fileName, "Unable to read 1st data line") 141 skData%nGrid = skData%nGrid - 1 142 if (homo) then 143 skData%skSelf(nShell+1:) = 0.0_dp; 144 skData%skHubbU(nShell+1:) = 0.0_dp; 145 skData%skOcc(nShell+1:) = 0.0_dp; 146 read (file,*, iostat=iostat) (skData%skSelf(ii), ii = nShell, 1, -1), & 147 &rDummy, & 148 &(skData%skHubbU(ii), ii = nShell, 1, -1),& 149 &(skData%skOcc(ii), ii = nShell, 1, -1) 150 call checkIoError(iostat, fileName, "Unable to read 2nd data line") 151 read (file,*, iostat=iostat) skData%mass, (coeffs(ii), ii = 2, 9), & 152 &polyCutoff, rDummy, (rDummy, ii = 12, 20) 153 call checkIoError(iostat, fileName, "Unable to read 3rd data line") 154 skData%mass = skData%mass * amu__au !convert to atomic units 155 else 156 read (file,*, iostat=iostat) rDummy, (coeffs(ii), ii = 2, 9),& 157 & polyCutoff, (rDummy, ii = 11, 20) 158 call checkIoError(iostat, fileName, "Unable to read 1st data line") 159 end if 160 161 if (present(repPolyIn)) then 162 repPolyIn%polyCoeffs(:) = coeffs(:) 163 repPolyIn%cutoff = polyCutoff 164 end if 165 166 allocate(skData%skHam(skData%nGrid, nSKInter)) 167 allocate(skData%skOver(skData%nGrid, nSKInter)) 168 skData%skHam(:,:) = 0.0_dp 169 skData%skOver(:,:) = 0.0_dp 170 do iGrid = 1, skData%nGrid 171 if (tExtended) then 172 read (file, *, iostat=iostat) & 173 &(skData%skHam(iGrid, ii), ii = 1, nSKInter), & 174 &(skData%skOver(iGrid, ii), ii = 1, nSKInter) 175 call checkIoError(iostat, fileName, "Reading error for integrals") 176 else 177 read(file,*, iostat=iostat) & 178 &(skData%skHam(iGrid, iSKInterOld(ii)), ii = 1, nSKInterOld), & 179 &(skData%skOver(iGrid, iSKInterOld(ii)), ii = 1, nSKInterOld) 180 call checkIoError(iostat, fileName, "Reading error for integrals") 181 end if 182 end do 183 184 if (.not. present(repSplineIn)) then 185 close(file) 186 return 187 end if 188 189 call readSplineRep(file, fileName, repSplineIn, iSp1, iSp2) 190 191 !> Read range separation parameter 192 if (present(rangeSepSK)) then 193 call readRangeSep(file, fileName, rangeSepSK) 194 end if 195 196 197 close(file) 198 199 end subroutine OldSKData_readFromFile 200 201 202 !> Reads the repulsive from an open file. 203 subroutine OldSKData_readsplinerep(fp, fname, repsplinein, isp1, isp2) 204 !! File identifier. 205 integer, intent(in) :: fp 206 !! Name of the file (for printing help messages). 207 character(*), intent(in) :: fname 208 !! Input structure for the spline repulsives 209 type(trepsplinein), intent(inout) :: repsplinein 210 !! Index of the first species in the repulsive (for messages) 211 integer, intent(in), optional :: isp1 212 !! Index of the second species in the repulsive (for messsages) 213 integer, intent(in), optional :: isp2 214 215 integer :: iostat 216 integer :: nint, ii, jj 217 character(lc) :: chdummy 218 logical :: hasspline 219 real(dp), allocatable :: xend(:) 220 221 ! Look for spline 222 do 223 read(fp, '(A)', iostat=iostat) chdummy 224 if (iostat /= 0) then 225 hasspline = .false. 226 exit 227 elseif (chdummy == "Spline") then 228 hasspline = .true. 229 exit 230 end if 231 end do 232 233 if (.not. hasspline) then 234 write(chdummy, "(A,A,A)") "No spline repulsive found in file '",& 235 & trim(fname), "'" 236 call error(chdummy) 237 end if 238 239 read(fp, *, iostat=iostat) nint, repsplinein%cutoff 240 call checkioerror(iostat, fname, "Error in reading nint and cutoff") 241 read(fp, *, iostat=iostat) (repsplinein%expcoeffs(ii), ii = 1, 3) 242 call checkioerror(iostat, fname, "Error in reading exponential coeffs") 243 allocate(repsplinein%xstart(nint)) 244 allocate(repsplinein%spcoeffs(4, nint - 1)) 245 allocate(xend(nint)) 246 247 do jj = 1, nint - 1 248 read(fp, *, iostat=iostat) repsplinein%xstart(jj), xend(jj),& 249 & (repsplinein%spcoeffs(ii,jj), ii = 1, 4) 250 call checkioerror(iostat, fname, "Error in reading spline coeffs") 251 end do 252 read(fp, *, iostat=iostat) repsplinein%xstart(nint), xend(nint),& 253 & (repsplinein%spLastCoeffs(ii), ii = 1, 6) 254 call checkioerror(iostat, fname, "Error in reading last spline coeffs") 255 repsplinein%cutoff = xend(nint) 256 ! Check on consistenty 257 do jj = 2, nint 258 if (abs(xend(jj-1) - repsplinein%xstart(jj)) > 1e-8_dp) then 259 if (present(isp1) .and. present(isp2)) then 260 write(chdummy, "(A,I2,A,I2,A)") "Repulsive not continuous for species& 261 & pair ", isp1, "-", isp2, "." 262 else 263 write(chdummy, "(A)") "Repulsive not continuous." 264 end if 265 call error(chdummy) 266 end if 267 end do 268 269 end subroutine OldSKData_readsplinerep 270 271 272 !> Reads the RangeSep data from an open file. 273 subroutine readRangeSep(fp, fname, rangeSepSK) 274 275 !> File identifier 276 integer, intent(in) :: fp 277 278 !> File name 279 character(*), intent(in) :: fname 280 281 !> Rangesep data 282 type(TRangeSepSKTag), intent(inout) :: rangeSepSK 283 284 integer :: iostat 285 integer :: nint, ii, jj 286 character(lc) :: chdummy 287 real(dp) :: omega 288 logical :: hasRangeSep 289 real(dp), allocatable :: xend(:) 290 291 !> Seek rangesep part in SK file 292 do 293 read(fp, '(A)', iostat=iostat) chdummy 294 if (iostat /= 0) then 295 hasRangeSep = .false. 296 exit 297 elseif (chdummy == "RangeSep") then 298 hasRangeSep = .true. 299 exit 300 end if 301 end do 302 303 if ( .not. hasRangeSep) then 304 write(chdummy, "(A,A,A)") "RangeSep extension tag not found in file '",& 305 & trim(fname), "'" 306 call error(chdummy) 307 end if 308 309 read(fp, *, iostat=iostat) chdummy, omega 310 call checkioerror(iostat, fname, "Error in reading range-sep method and range-sep parameter") 311 312 if (chdummy /= "LC") then 313 write(chdummy, "(A,A,A)") "Unknown range-separation method in SK file '", trim(fname), "'" 314 call error(chdummy) 315 end if 316 317 if (omega < 0.0_dp) then 318 write(chdummy, "(A)") "Range-separation parameter is negative" 319 call error(chdummy) 320 end if 321 322 rangeSepSK%omega = omega 323 324 end subroutine ReadRangeSep 325 326 327 !> Checks for IO errors and prints message. 328 subroutine checkIOError(iostat, fname, msg) 329 330 !> Flag of the IO operation. 331 integer, intent(in) :: iostat 332 333 !> Name of the file. 334 character(*), intent(in) :: fname 335 336 !> Message to print if IO operation flag is non-zero. 337 character(*), intent(in) :: msg 338 339 if (iostat /= 0) then 340 call error("IO error in file '" // trim(fname) // "': " // trim(msg)) 341 end if 342 343 end subroutine checkIOError 344 345end module dftbp_oldskdata 346