Orbital_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 orbital_module
 
    use const_gbl,       only: stdout
    use consts_mpi_ci,   only: symtype_cinfv, symtype_dinfh, symtype_d2h
    use mpi_gbl,         only: mpi_xermsg
    use precisn,         only: wp
    use scatci_routines, only: mkorbs, pmkorbs
    use symbolic_module, only: symbolicelementvector
 
    implicit none
 
    private
 
    public spinorbital, orbitaltable
 
    type :: spinorbital
        integer :: orbital_number
        integer :: gerude
        integer :: m_quanta
        integer :: spin
        integer :: positron_flag
        integer :: electron_idx
    end type spinorbital
 
 
    type :: orbitaltable
        logical :: initialized = .false.    
        integer :: total_num_spin_orbitals
        integer :: total_num_orbitals
        integer :: symmetry_type
        integer :: num_symmetries
        integer :: positron_flag
        integer :: mflg = 0
        integer,           allocatable :: mcorb(:), mcon(:), nsrb(:)
        type(spinorbital), allocatable :: spin_orbitals(:)
        integer,           allocatable :: orbital_map(:)
    contains
        procedure, public    ::    initialize => initialize_table
        procedure, public    ::    construct => construct_table
        procedure, public    ::    get_orbital_number
        procedure, public    ::    get_spin
        procedure, public    ::    get_gerude
        procedure, public    ::    get_electron_number
        procedure, public    ::    compute_electron_index
        procedure, public    ::    check_max_mcon_in_determinants
        procedure, public    ::    evaluate_ijkl_and_coeffs
        procedure, public    ::    get_mcon
        procedure, public    ::    get_minimum_mcon
        procedure, public    ::    get_two_minimum_mcon
        procedure, public    ::    add_positron
        procedure, private   ::    evaluate_case_one
        procedure, private   ::    evaluate_case_two
        procedure, private   ::    evaluate_case_three
        procedure, private   ::    evaluate_case_other
    end type orbitaltable
 
contains
 
    subroutine initialize_table (this, nsrb, norb, symtype, nsym, positron_flag)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: nsrb, nsym, norb, symtype, positron_flag
 
        this % total_num_spin_orbitals = nsrb
        this % total_num_orbitals      = norb
        this % symmetry_type           = symtype
        this % num_symmetries          = nsym
        this % positron_flag           = positron_flag
 
        if (allocated(this % spin_orbitals))    deallocate (this % spin_orbitals)
        if (allocated(this % mcon))             deallocate (this % mcon)
        if (allocated(this % nsrb))             deallocate (this % nsrb)
        if (allocated(this % mcorb))            deallocate (this % mcorb)
        if (allocated(this % orbital_map))      deallocate (this % orbital_map)
 
        allocate(this % spin_orbitals(this % total_num_spin_orbitals))
        allocate(this % mcon(this % total_num_spin_orbitals))
        allocate(this % nsrb(this % total_num_spin_orbitals))
        allocate(this % mcorb(this % total_num_orbitals))
        allocate(this % orbital_map(this % total_num_orbitals))
 
        this % mcorb(:) = 0
 
    end subroutine initialize_table
 
 
    subroutine construct_table (this, num_orbital_target_sym, num_orbital_target_sym_dinf, &
                                num_orbitals, num_elec_orbitals, num_orbitals_congen)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: num_orbital_target_sym(this % num_symmetries),        &
                                              num_orbital_target_sym_dinf(this % num_symmetries),   &
                                              num_orbitals(this % num_symmetries),                  &
                                              num_elec_orbitals(this % num_symmetries),             &
                                              num_orbitals_congen(this % num_symmetries)
 
        integer :: MN(this % total_num_spin_orbitals), MG(this % total_num_spin_orbitals), &
                   MM(this % total_num_spin_orbitals), MS(this % total_num_spin_orbitals), MPOS(this % total_num_spin_orbitals)
        integer :: dummy_iposit = 0, dummy_lusme = 0, ido
 
        if (this % positron_flag /= 0 .and. this % symmetry_type == symtype_d2h) then
            call pmkorbs(num_orbitals, num_elec_orbitals, num_orbitals_congen, this % num_symmetries,    &
                         mn, mg, mm, ms, this % mcon, this % mcorb, this % total_num_orbitals,           &
                         this % total_num_spin_orbitals, this % orbital_map, mpos, this % positron_flag, &
                         this % symmetry_type, dummy_lusme, stdout)
        else
            call mkorbs(this % num_symmetries, mn, mg, mm, ms, this % total_num_orbitals,      &
                        this % total_num_spin_orbitals, this % orbital_map, mpos, this % mcon, &
                        this % mcorb, this % positron_flag, num_orbital_target_sym,            &
                        num_orbital_target_sym_dinf, this % symmetry_type, dummy_lusme, stdout)
        end if
 
        do ido = 1, this % total_num_spin_orbitals
            this % spin_orbitals(ido) % orbital_number = mn(ido)
            this % spin_orbitals(ido) % gerude         = mg(ido)
            this % spin_orbitals(ido) % m_quanta       = mm(ido)
            this % spin_orbitals(ido) % spin           = ms(ido)
            this % spin_orbitals(ido) % positron_flag  = mpos(ido)
            this % spin_orbitals(ido) % electron_idx   = 0
        end do
 
    end subroutine construct_table
 
 
    subroutine compute_electron_index (this, num_electrons, reference_determinants)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: num_electrons
        integer,             intent(in)    :: reference_determinants(num_electrons)
 
        integer :: ido
 
        do ido = 1, num_electrons
            this % spin_orbitals(reference_determinants(ido)) % electron_idx = ido
        end do
 
    end subroutine compute_electron_index
 
 
    integer function check_max_mcon_in_determinants (this, n, determinants)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: determinants(:), n
 
        integer :: i
 
        check_max_mcon_in_determinants = 0
        do i = 1, n
            check_max_mcon_in_determinants = max(check_max_mcon_in_determinants, this % mcon(determinants(i)))
        end do
 
    end function check_max_mcon_in_determinants
 
 
    subroutine evaluate_ijkl_and_coeffs (this, dtrs, coeff, symmetry_type, symbol, flag)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: dtrs(4), symmetry_type, flag
        real(wp),            intent(in)    :: coeff
        class(symbolicelementvector)       :: symbol
 
        type(spinorbital) :: P,R            !Determinants A
        type(spinorbital) :: Q,S            !Determinants B
        integer  :: KA, KB, lpositron = 0
        integer  :: NSFA, NSFB, MLA, MLB    !I'm not 100% what these do alpha/beta?
        real(wp) :: SIGN
 
        if (this % positron_flag /= 0) lpositron = this % add_positron(dtrs, 1, 4)
 
        !Get our spin orbitals
        p = this % spin_orbitals(dtrs(1))
        q = this % spin_orbitals(dtrs(2))
        r = this % spin_orbitals(dtrs(3))
        s = this % spin_orbitals(dtrs(4))
 
        !If there are opposing spins
        if(p % spin + q % spin == 1 .or. r % spin + s % spin == 1 .or. p % positron_flag /= q % positron_flag) then
            nsfa = 0
        else
            nsfa = 1
            mla = 0
            if (symmetry_type < symtype_d2h) then
                ka = max(dtrs(1), dtrs(2))
                kb = max(dtrs(3), dtrs(4))
                if (ka < kb) then
                    if (r % m_quanta * s % m_quanta < 0) mla = 1
                else
                    if(p % m_quanta * q % m_quanta < 0) mla = 1
                end if
            end if
        end if
 
        !If there are opposing spins or if there is a single positron
        if (p % spin + s % spin == 1 .or. q % spin + r % spin == 1 .or. r % positron_flag /= s % positron_flag) then
            nsfb = 0
        else
            nsfb = 1
            mlb = 0
            if (symmetry_type < symtype_d2h) then
                ka = max(dtrs(1), dtrs(4))
                kb = max(dtrs(3), dtrs(2))
                if (ka < kb) then
                    if (q % m_quanta * r % m_quanta < 0) mlb = 1
                else
                    if (p % m_quanta * s % m_quanta < 0) mlb = 1
                end if
            end if
        end if
 
        !No exchange so ignore
        if (nsfa == 0 .and. nsfb == 0) return
 
        sign = 1_wp
        if (lpositron == 1) sign = -1_wp
 
        if (flag == 1) then
            if (nsfa /= 0) then
                call this % evaluate_case_other(p % orbital_number, q % orbital_number, r % orbital_number, &
                                                s % orbital_number, mla, 0, sign * coeff, symbol)
            end if
            if (nsfb /= 0) then
                call this % evaluate_case_other(p % orbital_number, s % orbital_number, r % orbital_number, &
                                                q % orbital_number, mlb, 1, sign * coeff, symbol)
            end if
        else if (p % orbital_number == q % orbital_number .and. r % orbital_number == s % orbital_number) then
 
            call this % evaluate_case_one(p % orbital_number, r % orbital_number, nsfa, nsfb, mla, mlb, sign * coeff, symbol)
 
        else if (p % orbital_number == q % orbital_number .or. r % orbital_number == s % orbital_number) then
            if (nsfa /= 0) then
                call this % evaluate_case_other(p % orbital_number, q % orbital_number, r % orbital_number, &
                                                s % orbital_number, mla, 0, sign * coeff, symbol)
            end if
            if (nsfb /= 0) then
                call this % evaluate_case_other(p % orbital_number, s % orbital_number, r % orbital_number, &
                                                q % orbital_number, mlb, 1, sign * coeff, symbol)
            end if
        else if (p % orbital_number == r % orbital_number .and. q % orbital_number == s % orbital_number) then
 
            call this % evaluate_case_three(p % orbital_number, r % orbital_number, nsfa, nsfb, mla, mlb, sign * coeff, symbol)
 
        else if (p % orbital_number == r % orbital_number) then
            if (nsfa /= 0) then
                call this % evaluate_case_other(p % orbital_number, q % orbital_number, r % orbital_number, &
                                                s % orbital_number, mla, 0, sign * coeff, symbol)
            end if
            if (nsfb /= 0) then
                call this % evaluate_case_other(p % orbital_number, s % orbital_number, r % orbital_number, &
                                                q % orbital_number, mlb, 1, sign * coeff, symbol)
            end if
        else if(p % orbital_number == s % orbital_number .and. q % orbital_number == r % orbital_number) then
 
            call this % evaluate_case_two(p % orbital_number, r % orbital_number, nsfa, nsfb, mla, mlb, sign * coeff, symbol)
 
        else
            if (nsfa /= 0) then
                call this % evaluate_case_other(p % orbital_number, q % orbital_number, r % orbital_number, &
                                                s % orbital_number, mla, 0, sign * coeff, symbol)
            end if
            if (nsfb /= 0) then
                call this % evaluate_case_other(p % orbital_number, s % orbital_number, r % orbital_number, &
                                                q % orbital_number, mlb, 1, sign * coeff, symbol)
            end if
        end if
 
    end subroutine evaluate_ijkl_and_coeffs
 
 
    subroutine evaluate_case_one (this, p, r, nsfa, nsfb, mla, mlb, coeff, symbol)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: p, r, NSFA, NSFB, MLA, MLB
        real(wp),            intent(in)    :: coeff
        class(symbolicelementvector)       :: symbol
 
        integer :: I, J
 
        if (p < r) then
            i = r
            j = p
        else
            i = p
            j = r
        end if
 
        if (nsfa /= 0) then
            if (mla == 0) then
                call symbol % insert_ijklm_symbol(i, i, j, j, 0, coeff)
            else
                call symbol % insert_ijklm_symbol(i, i, j, j, 1, coeff)
            end if
        end if
 
        if (nsfb /= 0) then
            if (mlb == 0) then
                call symbol % insert_ijklm_symbol(i, j, i, j, 0, -coeff)
            else
                call symbol % insert_ijklm_symbol(i, j, i, j, 1, -coeff)
            end if
        end if
 
    end subroutine evaluate_case_one
 
 
    subroutine evaluate_case_two (this, p, q, nsfa, nsfb, mla, mlb, coeff, symbol)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: p, q, NSFA, NSFB, MLA, MLB
        real(wp),            intent(in)    :: coeff
        class(symbolicelementvector)       :: symbol
 
        integer :: I, J
 
        if (p < q) then
            i = q
            j = p
        else
            i = p
            j = q
        end if
 
        if (nsfa /= 0) then
            if (mla == 0) then
                call symbol % insert_ijklm_symbol(i, j, i, j, 1, coeff)
            else
                call symbol % insert_ijklm_symbol(i, j, i, j, 0, coeff)
            end if
        end if
 
        if (nsfb /= 0) then
            if (mlb == 0) then
                call symbol % insert_ijklm_symbol(i, i, j, j, 1, -coeff)
            else
                call symbol % insert_ijklm_symbol(i, i, j, j, 0, -coeff)
            end if
        end if
 
    end subroutine evaluate_case_two
 
 
    subroutine evaluate_case_three (this, p, q, nsfa, nsfb, mla, mlb, coeff, symbol)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: p, q, NSFA, NSFB, MLA, MLB
        real(wp),            intent(in)    :: coeff
        class(symbolicelementvector)       :: symbol
 
        integer :: I, J
 
        if (p < q) then
            i = q
            j = p
        else
            i = p
            j = q
        end if
 
        if (nsfa /= 0) then
            if (mla == 0) then
                call symbol % insert_ijklm_symbol(i, j, i, j, 1, coeff)
            else
                call symbol % insert_ijklm_symbol(i, j, i, j, 0, coeff)
            end if
        end if
 
        if (nsfb /= 0) then
            if (mlb == 0) then
                call symbol % insert_ijklm_symbol(i, j, i, j, 1, -coeff)
            else
                call symbol % insert_ijklm_symbol(i, j, i, j, 0, -coeff)
            end if
        end if
 
    end subroutine evaluate_case_three
 
 
    subroutine evaluate_case_other (this, p, q, r, s, mla, ms, coeff, symbol)
        class(orbitaltable), intent(inout) :: this
        integer,             intent(in)    :: p, q, r, s, mla, ms
        real(wp),            intent(in)    :: coeff
        class(symbolicelementvector)       :: symbol
 
        real(wp) :: cfd
        integer  :: I, J, K, L
 
        if (p < q) then
            i = q
            j = p
        else
            i = p
            j = q
        end if
        if (r < s) then
            k = s
            l = r
        else
            k = r
            l = s
        end if
 
        cfd = coeff
        if (ms /= 0) cfd = -cfd
 
        if (i > k .or. (i == k .and. j > l)) then
            call symbol % insert_ijklm_symbol(i, j, k, l, mla, cfd)
        else
            call symbol % insert_ijklm_symbol(k, l, i, j, mla, cfd)
        end if
 
    end subroutine evaluate_case_other
 
 
    integer function get_orbital_number (this, spin_orbital)
 
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: spin_orbital
 
        if (spin_orbital > this % total_num_spin_orbitals .or. spin_orbital < 1) then
            write (stdout, "('Orbital selected: ',2i12)") spin_orbital, this % total_num_spin_orbitals
            flush(stdout)
            call mpi_xermsg('Orbital_module', 'get_orbital_number', 'Selected orbital not within range', 1, 1)
        end if
 
        get_orbital_number = this % spin_orbitals(spin_orbital) % orbital_number
 
    end function get_orbital_number
 
 
    integer function get_spin (this, spin_orbital)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: spin_orbital
 
        if (spin_orbital > this % total_num_spin_orbitals .or. spin_orbital < 1) then
            write (stdout, "('Orbital selected: ',2i12)") spin_orbital, this % total_num_spin_orbitals
            flush(stdout)
            call mpi_xermsg('Orbital_module', 'get_spin', 'Selected orbital not within range', 1, 1)
        end if
 
        get_spin = this % spin_orbitals(spin_orbital) % spin
 
    end function get_spin
 
 
    integer function get_gerude (this, spin_orbital)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: spin_orbital
 
        if (spin_orbital > this % total_num_spin_orbitals .or. spin_orbital < 1) then
            write (stdout, "('Orbital selected: ',2i12)") spin_orbital, this % total_num_spin_orbitals
            flush(stdout)
            call mpi_xermsg('Orbital_module', 'get_gerude', 'Selected orbital not within range', 1, 1)
        end if
 
        get_gerude = this % spin_orbitals(spin_orbital) % gerude
 
    end function get_gerude
 
 
    integer function get_electron_number (this, spin_orbital)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: spin_orbital
 
        if (spin_orbital > this % total_num_spin_orbitals .or. spin_orbital < 1) then
            write (stdout, "('Orbital selected: ',2i12)") spin_orbital, this % total_num_spin_orbitals
            flush(stdout)
            call mpi_xermsg('Orbital_module', 'get_electron_number', 'Selected orbital not within range', 1, 1)
        end if
 
        get_electron_number = this % spin_orbitals(spin_orbital) % electron_idx
 
    end function get_electron_number
 
 
    integer function get_minimum_mcon (this, determinants)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: determinants(4)
 
        !get_minimum_mcon = min(this%get_two_minimum_mcon(determinants(1),determinants(2)),  &
        !           &    this%get_two_minimum_mcon(determinants(3),determinants(4)))
 
        get_minimum_mcon = min(this % mcon(determinants(1)), &
                               this % mcon(determinants(2)), &
                               this % mcon(determinants(3)), &
                               this % mcon(determinants(4)))
 
    end function get_minimum_mcon
 
 
    integer function get_mcon (this, spin_orbital)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: spin_orbital
 
        if (spin_orbital > this % total_num_spin_orbitals .or. spin_orbital < 1) then
            write (stdout, "('Orbital selected: ',2i12)") spin_orbital, this % total_num_spin_orbitals
            flush(stdout)
            call mpi_xermsg('Orbital_module', 'get_mcon', 'Selected orbital not within range', 1, 1)
        end if
 
        get_mcon = this % mcon(spin_orbital)
 
    end function get_mcon
 
 
    integer function get_two_minimum_mcon (this, determinant_one, determinant_two)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: determinant_one, determinant_two
 
        get_two_minimum_mcon = min(this % mcon(determinant_one), this % mcon(determinant_two))
 
    end function get_two_minimum_mcon
 
 
    integer function add_positron (this, determinants, ia, ib)
        class(orbitaltable), intent(in) :: this
        integer,             intent(in) :: determinants(4), ia, ib
        integer :: icount, ido, ip
 
        add_positron = 0
        icount = 0
 
        do ido = ia, ib
            ip = determinants(ido)
            if (this % spin_orbitals(ip) % positron_flag /= 0) icount = icount + 1
        end do
 
        if (icount == 0) return
        if (icount == 2) then
            add_positron = 1
            return
        end if
 
        write (stdout, "(' ICOUNT ',i4,' NDTC = ',4i4)") icount, determinants(:)
        flush(stdout)
        call mpi_xermsg('Orbital_module', 'add_positron', 'Invalid count', 1, 1)
 
    end function add_positron
 
end module orbital_module