basis_library.F90 Source File

module basis_library

  use, intrinsic :: iso_fortran_env, only: real64
  use elements, only: MAX_ELEMENT_Z, get_element_id
  use strings, only: to_upper
  use constants, only: ANGULAR_LABEL, NUM_CART_BF
  use io_constants, only: IW

  implicit none

  private

  type, public :: atom_basis_t
    integer :: nshells = 0
    integer :: nbfs = 0
    integer :: nprims = 0
    integer, allocatable :: ang(:)
    integer, allocatable :: ncontract(:)
    real(real64), allocatable :: ex(:)
    real(real64), allocatable :: cc(:)
  contains
    procedure, private :: reserve => atom_basis_reserve
  end type

  type, public :: basis_library_t
    type(atom_basis_t) :: atoms(MAX_ELEMENT_Z)
  contains
    procedure :: from_file => read_basis_library
    procedure :: calc_req_storage
    procedure :: echo => basis_library_echo
  end type

contains

!-------------------------------------------------------------------------------

  subroutine calc_req_storage(this, z, nshell, ngauss, nbasis)
    class(basis_library_t), target, intent(in) :: this
    integer, intent(in) :: z(:)
    integer, intent(out) :: nshell, ngauss, nbasis

    nshell = sum(this%atoms(z)%nshells)
    nbasis = sum(this%atoms(z)%nbfs)
    ngauss = sum(this%atoms(z)%nprims)

  end subroutine

!-------------------------------------------------------------------------------

  subroutine read_basis_library(this, fname)
    class(basis_library_t), target, intent(inout) :: this
    character(*), intent(in) :: fname
    integer :: iunit
    integer :: stat
    character(len=1024) :: msg

    open (newunit=iunit, file=trim(adjustl(fname)), status='old', action='read', &
          iostat=stat, iomsg=msg)

    if (stat /= 0) then
      write (IW, *) 'Can''t open basis library file ''', trim(adjustl(fname)), ''''
      write (IW, *) 'The error is:'
      write (IW, *) trim(msg)
      flush (IW)
      call abort
    end if

    call read_basis_library_file(this, iunit)

    close (iunit)

  end subroutine

!-------------------------------------------------------------------------------

  subroutine read_basis_library_file(this, iunit)
    class(basis_library_t), target, intent(inout) :: this
    integer, intent(in) :: iunit
    type(atom_basis_t), pointer :: atom
    integer :: z
    integer :: ig
    integer :: ln
    integer :: state
    character(1) :: label
    integer :: ng
    character(len=1024) :: error
    integer :: iost

    integer :: id

    character(1024) :: line

    integer :: l
    integer, parameter :: &
      READ_ATOM = 0, &
      READ_SHELL = 1, &
      READ_GAUSS = 2

    state = READ_ATOM
    ig = 0
    ln = 0

    do
      read (unit=iunit, fmt='(A)', end=100) line
      ln = ln+1
!          Check if line is commented out
      if (skip_string(line)) cycle

      select case (state)
      case (READ_ATOM)
        z = get_element_id(trim(adjustl(line)))

        if (z > 0) then
          state = READ_SHELL
          atom => this%atoms(z)
        end if

      case (READ_SHELL)
        if (len_trim(line) == 0) then
          state = READ_ATOM
          ig = 0
          cycle
        end if

        read (line, *, err=200, iostat=iost, iomsg=error) label, ng

        atom%nshells = atom%nshells+1
        call atom%reserve(nshell=atom%nshells, ngauss=ig+ng)

        l = scan(ANGULAR_LABEL, label)
        atom%ang(atom%nshells) = l-1
        atom%nbfs = atom%nbfs+NUM_CART_BF(l-1)

        atom%ncontract(atom%nshells) = ng
        atom%nprims = atom%nprims+ng

        if (atom%ang(atom%nshells) < 0) then
          write (IW, *) "Unknown basis function type: '", label, "'"
          goto 200
        end if

        state = READ_GAUSS

      case (READ_GAUSS)
        ig = ig+1
        read (line, fmt=*, err=200, iostat=iost, iomsg=error) id, atom%ex(ig), atom%cc(ig)

        ng = ng-1
        if (ng == 0) state = READ_SHELL

      end select
    end do

100 continue
    return

200 continue
    write (IW, '("Error while reading basis set file at line ",I0,":")') ln
    write (IW, '(A)') trim(line)
    call abort

  end subroutine

!-------------------------------------------------------------------------------

  subroutine basis_library_echo(this)
    use elements, only: MAX_ELEMENT_Z, ELEMENTS_LONG_NAME
    class(basis_library_t), intent(in) :: this
    integer :: ia, sh, ig, g

    do ia = 1, MAX_ELEMENT_Z
      associate (at => this%atoms(ia))
        if (at%nshells == 0) cycle
        write (*, *) trim(ELEMENTS_LONG_NAME(ia)), at%nshells, at%nprims, at%nbfs
        ig = 0
        do sh = 1, at%nshells
          write (*, '(A1,i10)') ANGULAR_LABEL(at%ang(sh):at%ang(sh)), at%ncontract(sh)
          do g = 1, at%ncontract(sh)
            ig = ig+1
            write (*, '(i4,2ES25.15)') g, at%ex(ig), at%cc(ig)
          end do
        end do
      end associate
    end do
  end subroutine

!-------------------------------------------------------------------------------

  subroutine atom_basis_reserve(this, nshell, ngauss)
    class(atom_basis_t), intent(inout) :: this
    integer, optional, intent(in) :: nshell, ngauss

    if (present(nshell)) then
      call reserve_array_int(this%ang, nshell)
      call reserve_array_int(this%ncontract, nshell)
    end if

    if (present(ngauss)) then
      call reserve_array_real(this%ex, ngauss)
      call reserve_array_real(this%cc, ngauss)
    end if
  end subroutine

  subroutine reserve_array_int(a, n)
    integer, intent(in) :: n
    integer, allocatable, intent(inout) :: a(:)
    integer, allocatable :: tmp(:)
    integer :: lda
    lda = 0
    if (allocated(a)) lda = ubound(a, 1)
    if (n <= lda) return
    allocate (tmp(max(n, lda+16)))
    if (allocated(a)) tmp(:) = a(:)
    call move_alloc(from=tmp, to=a)
  end subroutine

  subroutine reserve_array_real(a, n)
    integer, intent(in) :: n
    real(real64), allocatable, intent(inout) :: a(:)
    real(real64), allocatable :: tmp(:)
    integer :: lda
    lda = 0
    if (allocated(a)) lda = ubound(a, 1)
    if (n <= lda) return
    allocate (tmp(max(n, lda+16)))
    if (allocated(a)) tmp(:lda) = a(:)
    call move_alloc(from=tmp, to=a)
  end subroutine

  logical function skip_string(str) result(res)
    character(*) :: str
    character(*), parameter :: TO_SKIP = '!#&$/\'
    integer :: pos
    res = .false.
    pos = verify(str, ' ')
    if (pos == 0) return
    res = scan(TO_SKIP, str(pos:pos)) /= 0
  end function

end module