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