Hamiltonian_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 hamiltonian_module
 
    use integer_packing,        only: unpack8ints
    use precisn,                only: longint, wp
    use baseintegral_module,    only: baseintegral
    use csf_module,             only: csfobject, csforbital
    use options_module,         only: options
    use orbital_module,         only: orbitaltable, spinorbital
    use parallelization_module, only: grid => process_grid
    use symbolic_module,        only: symbolicelementvector
 
    implicit none
 
    public basehamiltonian
 
    type, abstract :: basehamiltonian
        class(orbitaltable), pointer :: orbitals
        class(options),      pointer :: options
        class(baseintegral), pointer :: integral
        class(csfobject),    pointer :: csfs(:)
 
        integer  :: nflg = 0
        integer  :: diagonal_flag
        integer  :: positron_flag
        integer  :: phase_flag
        logical  :: constructed = .false.   
        logical  :: initialized = .false.   
        integer  :: job_id = 0
        integer  :: number_of_integrals = 0
        real(wp) :: element_one = 0.0
 
        type(symbolicelementvector) :: reference_symbol
    contains
        !Constructor
        procedure, public  :: construct => construct_base_hamiltonian
        procedure(generic_build), deferred :: build_hamiltonian
 
        !-----------------------Private procedures-----------------!
        procedure, public  :: slater_rules
       !procedure, public  :: obey_slater_rules
        procedure, public  :: evaluate_integrals
        procedure, public  :: evaluate_integrals_singular
       !procedure, public  :: store_first_element
 
        !-------Slater rule functions------------------!
        procedure, private :: handle_two_pair
        procedure, private :: handle_one_pair
        procedure, private :: handle_same
        procedure, public  :: my_job
    end type basehamiltonian
 
    abstract interface
 
        subroutine generic_build(this,matrix_elements)
            use basematrix_module, only: basematrix
            import :: basehamiltonian
            class(basehamiltonian)           :: this
            class(basematrix), intent(inout) :: matrix_elements
        end subroutine generic_build
    end interface
 
contains
 
    subroutine construct_base_hamiltonian (this, option, csfs, orbitals, integral)
        class(basehamiltonian),      intent(inout) :: this
        class(options),      target, intent(in)    :: option
        class(csfobject),    target, intent(in)    :: csfs(:)
        class(orbitaltable), target, intent(in)    :: orbitals
        class(baseintegral), target, intent(in)    :: integral
 
        !Assign our member pointers for each class
        this % options  => option
        this % csfs     => csfs
        this % orbitals => orbitals
        this % integral => integral
        this % phase_flag = this % options % phase_correction_flag
        this % positron_flag = this % options % positron_flag
 
        this % job_id = 0
        this % diagonal_flag = min(this % options % matrix_eval, 2)
        this % constructed = .true.
 
    end subroutine construct_base_hamiltonian
 
 
    subroutine slater_rules (this, csf_a, csf_b, sym_elements, flag, job_)
        class(basehamiltonian),       intent(in)    :: this
        class(csfobject),             intent(in)    :: csf_a, csf_b
        class(symbolicelementvector), intent(inout) :: sym_elements
        integer,           intent(in)  :: flag
        logical, optional, intent(in)  :: job_
        logical  :: job
        integer  :: diag, num_differences, dtr_diff(4), dtr_diff_fast(4), dtr_idx_a, dtr_idx_b
        real(wp) :: coeff
 
        !Check for job parameter
        if (present(job_)) then
            job = job_
        else
            job = .true.
        end if
 
        !Clear the symbolic elements ahead of time
        call sym_elements % clear()
 
        !Check if we should be doing this
        if (job) then
            !If not then return
            if (.not. this % my_job()) return
        end if
 
        !Check if we are not the same
        if (csf_a % id /= csf_b % id) then
            !If we are not then quickly check if we need to do a slater calculation
            if (csf_a % check_slater(csf_b) > 4) return
            !if(this%obey_slater_rules(csf_a,csf_b) > 4) return
        end if
 
        !Loop through the first set of determinants
        do dtr_idx_a = 1, csf_a % num_orbitals
            !Loop through the second set of determinants
            do dtr_idx_b = 1, csf_b % num_orbitals
 
                !Use the fast slater rules method to check the number of differences and get the differing determinants
                call csf_a % orbital(dtr_idx_a) % compare_excitations_fast(csf_b % orbital(dtr_idx_b), &
                                                                           this % options % num_electrons, &
                                                                           coeff, &
                                                                           num_differences, &
                                                                           dtr_diff)
 
                !Depending on the number of differences, we handle it appropriately
                select case (num_differences)
                    case(3:)
                        cycle
                    case(2)
                        call this % handle_two_pair(dtr_diff, coeff, sym_elements, flag)
                    case(1)
                        call this % handle_one_pair(dtr_diff, coeff, sym_elements, csf_a % orbital(dtr_idx_a), flag)
                    case(0)
                        call this % handle_same(dtr_diff, coeff, sym_elements, csf_a % orbital(dtr_idx_a), flag)
                end select
 
            end do
        end do
 
    end subroutine slater_rules
 
 
    subroutine handle_two_pair (this, dtrs, coeff, sym_element, flag)
        class(basehamiltonian)       ::    this
        class(symbolicelementvector) ::    sym_element
        integer,  intent(in)         ::    dtrs(4)
        real(wp), intent(in)         ::    coeff
        integer,  intent(in)         :: flag
 
        if (this % diagonal_flag> 1 .and. this % orbitals % get_minimum_mcon(dtrs) > 0) return
 
        !Evaluate for symbolic elements
        call this % orbitals % evaluate_IJKL_and_coeffs(dtrs, coeff, 0, sym_element, flag)
 
    end subroutine handle_two_pair
 
 
    subroutine handle_one_pair (this, dtrs, coeff, sym_element, csf, flag)
        class(basehamiltonian)           :: this
        class(symbolicelementvector)     :: sym_element
        class(csforbital), intent(in)    :: csf
        integer,           intent(inout) :: dtrs(4)
        real(wp),          intent(in)    :: coeff
        integer,           intent(in)    :: flag
 
        type(spinorbital) :: so1, so2
 
        integer :: ref_dtrs(this % options % num_electrons)
        integer :: num_electrons, N, M, ido, ia1, ia2, lnadd
        integer :: nchk !Whether its continuum maybe?
 
        !Reset the positron flag
        lnadd = 0
 
        !Get number of electrons
        num_electrons = this % options % num_electrons
 
        !Copt the reference determinant to a seperate array
        !ref_dtrs(:) = this%options%reference_dtrs(:)
 
        !Loop through the reference and replace the spin orbitals with ones from the CSF
        !do ido=1,csf%num_replaced_dtrs
        !    N = csf%replacing(ido)
        !    M = this%orbitals%get_electron_number(N)
        !    ref_dtrs(M) = csf%replaced(ido)
        !enddo
 
        call csf % get_determinants(ref_dtrs, num_electrons)
 
        !Not sure what this does exactly
        if (this % diagonal_flag > 1 .and. this % orbitals % get_two_minimum_mcon(dtrs(3), dtrs(4)) > 0) then
            nchk = 1
        else
            nchk = 0
        end if
 
        !Loop through all the spin orbitals
        do ido = 1, num_electrons
            !If we have the same spin orbitals then skip
            if (ref_dtrs(ido) == dtrs(3)) cycle
            !Again not quite sure what this does
            if (nchk == 1 .and. this % orbitals % get_mcon(ref_dtrs(ido)) > 0) cycle
            !Place the reference into our determinant array
            dtrs(1) = ref_dtrs(ido)
            dtrs(2) = ref_dtrs(ido)
            !Now evaluate for symbols
            call this % orbitals % evaluate_IJKL_and_coeffs(dtrs, coeff, 0, sym_element, flag)
        end do
 
        !If this run contains positrons, then add them in if necessary
        if (this % positron_flag /= 0) lnadd = this % orbitals % add_positron(dtrs, 3, 4)
        if (nchk == 1) return
 
        !Get our original slater replacements
        so1 = this % orbitals % spin_orbitals(dtrs(3))
        so2 = this % orbitals % spin_orbitals(dtrs(4))
 
        !Check if they posses the same spin and m quanta
        if (so1 % m_quanta /= so2 % m_quanta .or. so1 % spin /= so2 % spin) return
        if (so1 % orbital_number < so2 % orbital_number) then
            ia1 = so2 % orbital_number
            ia2 = so1 % orbital_number
        else
            ia1 = so1 % orbital_number
            ia2 = so2 % orbital_number
        end if
 
        !if they do then add a 0I0J + P symbol
        call sym_element % insert_ijklm_symbol(0, ia1, 0, ia2, lnadd, coeff)
 
    end subroutine handle_one_pair
 
 
    subroutine handle_same (this, dtrs, coeff, sym_element, csf, flag)
        class(basehamiltonian)           :: this
        class(symbolicelementvector)     :: sym_element
        class(csforbital), intent(in)    :: csf
        integer,           intent(inout) :: dtrs(4)
        real(wp),          intent(in)    :: coeff
        integer,           intent(in)    :: flag
 
        type(spinorbital) :: so1, so2
 
        integer :: ref_dtrs(this % options % num_electrons)
        integer :: num_electrons, N, M, ido, jdo, ia1, ia2, lnadd
        integer :: nchk   !Whether its continuum maybe?
 
        !Reset the positron flag
        lnadd = 0
 
        !Get number of electrons
        num_electrons = this % options % num_electrons
 
        !Copt the reference determinant to a seperate array
        !ref_dtrs(:) = this%options%reference_dtrs(:)
 
        !Loop through the reference and replace the spin orbitals with ones from the CSF
        !do ido=1,csf%num_replaced_dtrs
        !    N = csf%replacing(ido)
        !    M = this%orbitals%get_electron_number(N)
        !    ref_dtrs(M) = csf%replaced(ido)
        !enddo
 
        call csf % get_determinants(ref_dtrs, num_electrons)
 
        do ido = 2, num_electrons
            dtrs(1) = ref_dtrs(ido)
            dtrs(2) = dtrs(1)
            do jdo = 1, ido - 1
                dtrs(3) = ref_dtrs(jdo)
                dtrs(4) = dtrs(3)
                if (this % diagonal_flag > 1 .and. this % orbitals % get_two_minimum_mcon(dtrs(1), dtrs(3)) > 0) cycle
                !Now evaluate for symbols
                call this % orbitals % evaluate_IJKL_and_coeffs(dtrs, coeff, 0, sym_element, flag)
            end do
        end do
 
        if (this % NFLG /= 0) this % NFLG = 2
        do ido = 1, num_electrons
            n = ref_dtrs(ido)
            if (this % diagonal_flag > 1 .and. this % orbitals%get_mcon(n) /= 0) cycle
 
            m = this % orbitals % get_orbital_number(n)
            if (this % positron_flag /= 0) lnadd = this % orbitals % add_positron(dtrs, 1, 2)
            call sym_element % insert_ijklm_symbol(0, m, 0, m, lnadd, coeff)
        end do
 
    end subroutine handle_same
 
 
    real(wp) function evaluate_integrals_singular (this, label, coeff, dummy_orb)
        class(basehamiltonian), intent(inout) :: this
        integer,                intent(in)    :: dummy_orb
        real(wp),               intent(in)    :: coeff
        integer(longint),       intent(in)    :: label(2)
 
        integer :: ido, jdo, num_elements, lwd(8)
 
        evaluate_integrals_singular = 0.0
 
        if (coeff == 0.0_wp) return
 
        !Search for the dummy orbital
        if (this % phase_flag > 0) then
            call unpack8ints(label, lwd)
            do jdo = 1, 4
                if (lwd(jdo) == this % phase_flag) return
            end do
        else if (this % phase_flag == 0) then
            call unpack8ints(label, lwd)
            do jdo = 1, 4
                if (lwd(jdo) <= 0) cycle
                if (this % orbitals % mcorb(lwd(jdo)) == 0) return
            end do
        end if
 
        evaluate_integrals_singular = coeff * this % integral % get_integralf(label)
 
        !$OMP ATOMIC
        this % number_of_integrals = this % number_of_integrals + 1
 
    end function evaluate_integrals_singular
 
 
    real(wp) function evaluate_integrals (this, symbolic_elements, dummy_orb)
        class(basehamiltonian),       intent(inout) :: this
        class(symbolicelementvector), intent(in)    :: symbolic_elements
        integer,                      intent(in)    :: dummy_orb
 
        integer          :: ido, num_elements
        integer(longint) :: label(2)
        real(wp)         :: coeff
 
        num_elements = symbolic_elements % get_size()
 
        evaluate_integrals = 0.0_wp
        if (num_elements == 0) return
        do ido = 1, num_elements
            call symbolic_elements % get_coeff_and_integral(ido, coeff, label)
            evaluate_integrals = evaluate_integrals + this % evaluate_integrals_singular(label, coeff, dummy_orb)
        end do
 
    end function evaluate_integrals
 
 
    logical function my_job (this)
        class(basehamiltonian) :: this
 
        my_job = .true.
 
        if (grid % gprocs == 1 .or. grid % grank == -1) return
 
        if (this % job_id /= grid % grank) then
            my_job = .false.
        end if
 
        this % job_id = this % job_id + 1
        if (this % job_id == grid % gprocs) then
            this % job_id = 0
        end if
 
    end function my_job
 
 
    integer function fast_slat (rep1, rep2, fin1, fin2, num_rep1, num_rep2, spin_orbs, tot_num_spin_orbs)
        integer, intent(in)    :: num_rep1, rep1(num_rep1), fin1(num_rep1)
        integer, intent(in)    :: num_rep2, rep2(num_rep2), fin2(num_rep2), tot_num_spin_orbs
        integer, intent(inout) :: spin_orbs(tot_num_spin_orbs)
 
        integer :: i_orbital
        integer :: i_dtrs
        integer :: i_spin_orbital
        integer :: orb_add
        integer :: orb_remove
 
        do i_dtrs = 1, num_rep1
            orb_remove = (rep1(i_dtrs) + 1) / 2
            orb_add = (fin1(i_dtrs) + 1) / 2
            spin_orbs(orb_remove) = spin_orbs(orb_remove) - 1
            spin_orbs(orb_add) = spin_orbs(orb_add) + 1
        end do
 
        do i_dtrs = 1, num_rep2
            orb_remove = (rep2(i_dtrs) + 1) / 2
            orb_add = (fin2(i_dtrs) + 1) / 2
            spin_orbs(orb_remove) = spin_orbs(orb_remove) + 1
            spin_orbs(orb_add) = spin_orbs(orb_add) - 1
        end do
 
        fast_slat = 0
 
        do i_dtrs = 1, num_rep1
            orb_remove = (rep1(i_dtrs) + 1) / 2
            orb_add = (fin1(i_dtrs) + 1) / 2
            fast_slat = fast_slat + spin_orbs(orb_remove)
            spin_orbs(orb_remove) = 0
            fast_slat = fast_slat + spin_orbs(orb_add)
            spin_orbs(orb_add) = 0
        end do
 
        do i_dtrs = 1, num_rep2
            orb_remove = (rep2(i_dtrs) + 1) / 2
            orb_add = (fin2(i_dtrs) + 1) / 2
            fast_slat = fast_slat + spin_orbs(orb_remove)
            spin_orbs(orb_remove) = 0
            fast_slat = fast_slat + spin_orbs(orb_add)
            spin_orbs(orb_add) = 0
        end do
 
    end function fast_slat
 
end module hamiltonian_module