UKRMOL_module.f90

Go to the documentation of this file.

! Copyright 2019
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-in (UKRmol+ suite).
!
!     UKRmol-in is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     UKRmol-in is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with  UKRmol-in (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
 
module ukrmol_module
 
    use precisn,             only: longint, wp
    use const_gbl,           only: stdout
    use baseintegral_module, only: baseintegral
    use options_module,      only: options
 
    implicit none
 
    public ukrmolintegral
 
    private
 
    type, extends(baseintegral) :: ukrmolintegral
        integer  :: num_csf
        integer  :: num_one_electron_integrals
        integer  :: num_two_electron_integrals
        integer  :: symmetry_type
        integer  :: matrix_io_unit
        integer  :: num_symmetry
        real(wp) :: exotic_mass
        integer  :: hamiltonian_unit
        integer,          allocatable :: num_orbitals_sym(:)
        real(wp),         allocatable :: xnuc(:), ynuc(:), znuc(:), charge(:)
        character(len=8), allocatable, dimension(:) :: cname
    contains
        procedure, public :: initialize_self    => initialize_ukrmol
        procedure, public :: finalize_self      => finalize_ukrmol
        procedure, public :: load_integrals     => load_ukrmol
        procedure, public :: get_integral_ijklm => get_integral_ukrmol
        procedure, public :: write_geometries   => write_geometries_ukrmol
        procedure, public :: destroy_integrals  => destroy_integral_ukrmol
    end type
 
contains
 
    subroutine initialize_ukrmol (this, option)
        class(ukrmolintegral)      :: this
        class(options), intent(in) :: option
 
        write (stdout, "('Using UKRMOL+ integrals')")
 
        this % num_csf = option % num_csfs
        this % num_symmetry = option % num_syms
 
        allocate(this % num_orbitals_sym(this % num_symmetry))
 
        this % num_orbitals_sym(:) = option % num_electronic_orbitals_sym(:)
        this % hamiltonian_unit = option % hamiltonian_unit
        this % exotic_mass = option % exotic_mass
        this % matrix_io_unit = option % hamiltonian_unit
        this % name = option % name
 
    end subroutine initialize_ukrmol
 
 
    subroutine finalize_ukrmol (this)
 
        use ukrmol_interface_gbl, only: destroy_ukrmolp
 
        class(ukrmolintegral) :: this
 
        if (allocated(this % num_orbitals_sym)) deallocate (this % num_orbitals_sym)
        if (allocated(this % xnuc))   deallocate (this % xnuc)
        if (allocated(this % ynuc))   deallocate (this % ynuc)
        if (allocated(this % znuc))   deallocate (this % znuc)
        if (allocated(this % charge)) deallocate (this % charge)
        if (allocated(this % cname))  deallocate (this % cname)
 
        call destroy_ukrmolp
 
    end subroutine finalize_ukrmol
 
 
    subroutine load_ukrmol (this, iounit)
        use ukrmol_interface_gbl, only: read_ukrmolp_ints, get_geom
 
        class(ukrmolintegral) :: this
        integer, intent(in)   :: iounit
        integer :: nalm = 0, lembf, matrix_size = 0,dummy_int
 
        this % num_one_electron_integrals = -1
        this % num_two_electron_integrals = -1
 
        call read_ukrmolp_ints(this % matrix_io_unit, iounit, this % number_of_matrix_records, &
                               this % num_one_electron_integrals, this % num_two_electron_integrals, &
                               this % num_csf, stdout, this % symmetry_type, this % num_symmetry, this % num_orbitals_sym, &
                               this % positron_flag, 1.0_wp, this % name, &
                               nalm, .false.)
 
        !Get geometries
        call get_geom(this % nnuc, this % cname, this % xnuc, this % ynuc, this % znuc, this % charge)
 
        !Read header files locally to store to memory
        rewind(this % matrix_io_unit)
        read (this % matrix_io_unit) matrix_size, this % number_of_matrix_records , dummy_int, matrix_size, &
                                     dummy_int, this % num_symmetries, dummy_int, dummy_int, dummy_int, dummy_int, this % nnuc, &
                                     dummy_int, this % NAME, this % NHE, this % DTNUC
        this % core_energy = this % DTNUC(1)
        close (this % matrix_io_unit)
 
    end subroutine load_ukrmol
 
 
    function get_integral_ukrmol (this, i, j, k, l, m) result(coeff)
        use ukrmol_interface_gbl, only: get_integral
 
        class(ukrmolintegral) :: this
        integer, intent(in)   :: i, j, k, l, m
        integer  :: a, b, c, d
        real(wp) :: coeff
 
        if (i == 0) then
            a = this % orbital_mapping(j)
            b = this % orbital_mapping(l)
            coeff = get_integral(a, b, 0, 0, m)
        else
            a = this % orbital_mapping(i)
            b = this % orbital_mapping(j)
            c = this % orbital_mapping(k)
            d = this % orbital_mapping(l)
            coeff = get_integral(a, b, c, d, m)
        end if
 
    end function get_integral_ukrmol
 
 
    subroutine write_geometries_ukrmol (this, iounit)
 
        class(ukrmolintegral), intent(in) :: this
        integer,               intent(in) :: iounit
        integer :: ido, ifail, i
 
        do ido = 1, this % nnuc
            write (iounit) this % cname(ido), this % xnuc(ido), this % ynuc(ido), this % znuc(ido), this % charge(ido)
        end do
 
        write (stdout, 120) (this % cname(i), this % xnuc(i), this % ynuc(i), this % znuc(i), this % charge(i), i = 1, this % nnuc)
        120 FORMAT(/,' Nuclear data     X         Y         Z       Charge',/,(3x,a8,2x,4f10.6))
 
    end subroutine write_geometries_ukrmol
 
 
    subroutine destroy_integral_ukrmol (this)
        use ukrmol_interface_gbl, only: destroy_ukrmolp
 
        class(ukrmolintegral) :: this
 
        !Where is the destory integral or cleanup function?
        !MADE ONE!
        call destroy_ukrmolp()
 
    end subroutine destroy_integral_ukrmol
 
end module ukrmol_module