Symbolic_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 symbolic_module
 
    use const_gbl,              only: stdout
    use consts_mpi_ci,          only: default_integral_threshold
    use containers,             only: bstree
    use integer_packing,        only: pack8ints, unpack8ints
    use iso_c_binding,          only: c_loc, c_ptr, c_f_pointer
    use mpi_gbl,                only: mpi_reduceall_max, mpi_mod_rotate_arrays_around_ring, mpi_xermsg
    use precisn,                only: longint, wp
    use memorymanager_module,   only: master_memory
    use parallelization_module, only: grid => process_grid
 
    implicit none
 
    public symbolicelementvector
 
    private
 
    type, extends(bstree) :: symbolicelementvector
        private
        integer(longint), allocatable :: electron_integral(:,:)
        real(wp),         allocatable :: coeffs(:)
        integer                       :: n = 0
        logical                       :: constructed = .false.   
        real(wp)                      :: threshold  = 0.0
        integer                       :: max_capacity = 0
        integer                       :: expand_size = 10
    contains
        ! type-bound procedures used by the bstree interface
        procedure, public  :: compare => bstree_compare
 
        ! own type-bound procedures
        procedure, public  :: construct
        procedure, public  :: insert_ijklm_symbol
        procedure, public  :: insert_symbol
        procedure, public  :: remove_symbol_at
        procedure, public  :: is_empty
        procedure, public  :: clear
        procedure, public  :: check_same_integral
        procedure, public  :: get_integral_label
        procedure, public  :: synchronize_symbols
        procedure, public  :: modify_coeff
        procedure, public  :: get_coefficient
        procedure, public  :: get_coeff_and_integral
        procedure, public  :: get_size
        procedure, public  :: add_symbols
        procedure, public  :: print => print_symbols
        procedure, private :: expand_array
        procedure, private :: check_bounds
        procedure, public  :: destroy
        procedure, private :: check_constructed
    end type  symbolicelementvector
 
contains
 
    integer function bstree_compare (this, i, j, data) result (verdict)
 
        class(symbolicelementvector), intent(in) :: this
        integer,                      intent(in) :: i, j
        type(c_ptr),        optional, intent(in) :: data
 
        integer(longint)          :: ii(2), jj(2)
        integer(longint), pointer :: kk(:)
 
        ! get the default index set to compare with (if provided)
        if (present(data)) call c_f_pointer(data, kk, (/ 2 /))
 
        ! obtain the index sets to compare
        if (i <= 0) then; ii = kk; else; ii = this % electron_integral(:,i); end if
        if (j <= 0) then; jj = kk; else; jj = this % electron_integral(:,j); end if
 
        ! compare the index sets using lexicographical comparison
        if (ii(1) < jj(1)) then
            verdict = -1
        else if (ii(1) > jj(1)) then
            verdict = +1
        else if (ii(2) < jj(2)) then
            verdict = -1
        else if (ii(2) > jj(2)) then
            verdict = +1
        else
            verdict = 0
        end if
 
    end function bstree_compare
 
 
    subroutine check_constructed (this)
        class(symbolicelementvector) :: this
        if (.not. this % constructed) then
            write (stdout, "('Vector::constructed - Vector is not constructed')")
            stop "Vector::constructed - Vector is not constructed"
        end if
    end subroutine check_constructed
 
 
    subroutine construct (this, threshold, initial_size)
        class(symbolicelementvector)   :: this
        real(wp), optional, intent(in) :: threshold
        integer,  optional, intent(in) :: initial_size
        integer :: err
 
        if (present(threshold)) then
            this % threshold = threshold
        else
            this % threshold = default_integral_threshold
        end if
 
        this % expand_size = 100
        this % max_capacity = this % expand_size
 
        !Allocate the vectors
        allocate(this % electron_integral(2, this % max_capacity), this % coeffs(this % max_capacity), stat = err)
        call master_memory % track_memory(kind(this % electron_integral), size(this % electron_integral), &
                                          err, 'SYMBOLIC::ELECTRONINT')
        call master_memory % track_memory(kind(this % coeffs), size(this % coeffs), err, 'SYMBOLIC::ELECTRONCOEFF')
        if (err /= 0) then
            write (stdout, "('SymbolicVector::construct- arrays not allocated')")
            stop "SymbolicVector arrays not allocated"
        end if
 
        !Do an intial clear
        call this % clear
 
        !This has been succesfully constructed
        this % constructed = .true.
    end subroutine construct
 
 
    logical function check_same_integral (this, integral, idx) result (found)
 
        class(symbolicelementvector), intent(in)  :: this
        integer(longint), target,     intent(in)  :: integral(2)
        integer,                      intent(out) :: idx
 
        call this % check_constructed
 
        idx   = this % locate(-1, c_loc(integral))
        found = idx > 0
 
    end function check_same_integral
 
 
    subroutine insert_ijklm_symbol (this, i, j, k, l, m, coeff, check_same_)
        class(symbolicelementvector)  :: this
        integer,           intent(in) :: i, j, k, l, m
        real(wp),          intent(in) :: coeff
        logical, optional, intent(in) :: check_same_
        integer(longint) :: integral_label(2)
        logical          :: check_same
 
        if (present(check_same_)) then
            check_same = check_same_
        else
            check_same = .true.
        end if
 
        call pack8ints(i, j, k, l, m, 0, 0, 0, integral_label)
        call this % insert_symbol(integral_label, coeff, check_same)
    end subroutine insert_ijklm_symbol
 
 
    subroutine insert_symbol (this, integral_label, coeff, check_same_)
        class(symbolicelementvector)  :: this
        integer(longint),  intent(in) :: integral_label(2)
        real(wp),          intent(in) :: coeff
        logical, optional, intent(in) :: check_same_
        integer :: idx = 0, check_idx
        logical :: check_same
 
        if (present(check_same_)) then
            check_same = check_same_
        else
            check_same = .true.
        end if
 
        ! Filter small coefficients
        if (abs(coeff) < this % threshold) return
 
        if (.not.check_same) then
 
            this % n = this % n + 1                                     ! update number of elements
            idx = this % n                                              ! the idx is the last element
 
            if (this % n > this % max_capacity) then                    ! if we've exceeded capacity then expand
                call this % expand_array()
            end if
 
            this % electron_integral(:,idx) = integral_label(:)         ! place into the array
            this % coeffs(idx) = 0.0_wp                                 ! zero the coefficient since we can guareentee the compiler will do so
            call this % insert(idx)                                     ! insert idx into the binary tree
 
        !Check if we have the same integral, if we do get the index
        else if (.not. this % check_same_integral(integral_label, idx)) then
 
            this % n = this % n + 1                                     ! update number of elements
            idx = this % n                                              ! the idx is the last element
 
            if (this % n > this % max_capacity) then                    ! if we've exceeded capacity then expand
                call this % expand_array()
            end if
 
            this % electron_integral(:,idx) = integral_label(:)         ! place into the array
            this % coeffs(idx) = 0.0_wp                                 ! zero the coefficient since we can guareentee the compiler will do so
            call this % insert(idx)                                     ! insert index into the binary tree
        end if
 
        !Now we have an index, lets add in the coeffcient
        this % coeffs(idx) = this % coeffs(idx) + coeff
        !write(111,"('COEFF :',2es16.8)") coeff,this%coeffs(idx)
 
    end subroutine insert_symbol
 
 
    subroutine expand_array(this)
        class(symbolicelementvector)  :: this
        integer(longint), allocatable :: temp_integral(:,:)
        real(wp),         allocatable :: temp_coeffs(:)
        integer :: temp_capacity,ido,err
 
        !Of course lets check if we are constructed
        call this % check_constructed
 
        !Now the tempcapacity is the current max capacity
        temp_capacity = this % max_capacity
 
        allocate(temp_integral(2, this % max_capacity), temp_coeffs(this % max_capacity))
 
        !Copy over the elements into the temporary array
        do ido = 1, temp_capacity
            temp_integral(1,ido) = this % electron_integral(1,ido)
            temp_integral(2,ido) = this % electron_integral(2,ido)
        end do
 
        temp_coeffs(1:temp_capacity) = this % coeffs(1:temp_capacity)
 
        call master_memory % free_memory(kind(this % electron_integral),size(this % electron_integral))
        call master_memory % free_memory(kind(this % coeffs), size(this % coeffs))
        !deallocate our old values
        deallocate(this % electron_integral)
        deallocate(this % coeffs)
 
        !Increase the max capacity
        this % max_capacity = this % max_capacity + this % expand_size
        !Allocate our new array with the new max capacity
        allocate(this % electron_integral(2, this % max_capacity), this % coeffs(this % max_capacity), stat = err)
        call master_memory % track_memory(kind(this % electron_integral), size(this % electron_integral), &
                                          err, 'SYMBOLIC::EXP::ELECTRONINT')
        call master_memory % track_memory(kind(this % coeffs), size(this % coeffs), err, 'SYMBOLIC::EXP::ELECTRONCOEFF')
 
        !Zero the elements
        this % electron_integral(:,:) = -1
        this % coeffs(:) = 0.0_wp
 
        do ido = 1, temp_capacity
            !Copy our old values back
            this % electron_integral(1,ido) = temp_integral(1,ido)
            this % electron_integral(2,ido) = temp_integral(2,ido)
        end do
        this % coeffs(1:temp_capacity) = temp_coeffs(1:temp_capacity)
        deallocate(temp_integral, temp_coeffs)
 
    end subroutine expand_array
 
 
    subroutine add_symbols (this, rhs, alpha_)
        class(symbolicelementvector)             :: this
        class(symbolicelementvector), intent(in) :: rhs
        real(wp), optional,           intent(in) :: alpha_
        real(wp) :: alpha
        integer  :: ido
 
        if (present(alpha_)) then
            alpha = alpha_
        else
            alpha = 1.0_wp
        end if
 
        !Loop through each rhs symbolic and add it to me scaled by alpha
        do ido = 1, rhs % n
            call this % insert_symbol(rhs % electron_integral(1:2,ido), rhs % coeffs(ido) * alpha)
        end do
    end subroutine add_symbols
 
 
    subroutine synchronize_symbols (this)
        class(symbolicelementvector) :: this
        real(wp),         allocatable         :: coeffs(:)
        integer(longint), allocatable, target :: labels(:,:)
        integer(longint), pointer             :: label_ptr(:)
        integer(longint)                      :: my_num_symbols, largest_num_symbols, procs_num_of_symbols
        integer                               :: ido, proc_id, ierr
 
        if (grid % gprocs <= 1) then
            return
        end if
 
        my_num_symbols = this % get_size()
 
        call mpi_reduceall_max(my_num_symbols, largest_num_symbols, grid % gcomm)
 
        if (largest_num_symbols == 0) return
 
        call master_memory % track_memory(kind(labels), int(largest_num_symbols * 2), 0, 'SYMBOLIC::MPISYNCH::LABELS')
        call master_memory % track_memory(kind(coeffs), int(largest_num_symbols), 0, 'SYMBOLIC::MPISYNCH::COEFFS')
        allocate(labels(2,largest_num_symbols), coeffs(largest_num_symbols), stat = ierr)
 
        label_ptr(1:largest_num_symbols*2) => labels(:,:)
        labels = 0
        coeffs = 0
 
        if (my_num_symbols > 0) then
            labels(1,1:my_num_symbols) = this % electron_integral(1,1:my_num_symbols)
            labels(2,1:my_num_symbols) = this % electron_integral(2,1:my_num_symbols)
            coeffs(1:my_num_symbols) = this % coeffs(1:my_num_symbols)
        end if
 
        procs_num_of_symbols = my_num_symbols
 
        do proc_id = 1, grid % gprocs - 1
            call mpi_mod_rotate_arrays_around_ring(procs_num_of_symbols, label_ptr, coeffs, largest_num_symbols, grid % gcomm)
            do ido = 1, procs_num_of_symbols
                call this % insert_symbol(labels(1:2,ido), coeffs(ido))
            end do
        end do
 
        call master_memory % free_memory(kind(labels), size(labels))
        call master_memory % free_memory(kind(coeffs), size(coeffs))
        deallocate(labels, coeffs)
 
    end subroutine synchronize_symbols
 
 
    logical function check_bounds (this, i)
        class(symbolicelementvector) :: this
        integer, intent(in) :: i
 
        if (i <= 0 .or. i > this % n) then
            write (stdout, "('SymbolicVector::check_bounds - Out of Bounds access')")
            stop "SymbolicVector::check_bounds - Out of Bounds access"
            check_bounds = .false.
        else
            check_bounds = .true.
        end if
 
    end function check_bounds
 
 
    subroutine remove_symbol_at (this, idx)
        class(symbolicelementvector) :: this
        integer, intent(in)          :: idx
        integer                      :: ido
 
        if (this % check_bounds(idx)) then
            do ido = idx + 1, this % n
                this % electron_integral(1,ido-1) = this % electron_integral(1,ido)
                this % electron_integral(2,ido-1) = this % electron_integral(2,ido)
                this % coeffs(ido-1)              = this % coeffs(ido)
            end do
            this % electron_integral(:, this % n) = 0
            this % coeffs(this % n) = 0
            this % n = this % n - 1
        end if
 
    end subroutine remove_symbol_at
 
 
    logical function is_empty (this)
        class(symbolicelementvector) :: this
 
        is_empty = (this % n == 0)
 
    end function is_empty
 
 
    subroutine clear (this)
        class(symbolicelementvector) :: this
 
        this % n = 0
        call this % bstree % destroy
 
    end subroutine clear
 
 
    subroutine modify_coeff (this, idx, coeff)
        class(symbolicelementvector) :: this
        integer, intent(in) :: idx
        real(wp) :: coeff
 
        this % coeffs(idx) = this % coeffs(idx) + coeff
 
    end subroutine modify_coeff
 
 
    function get_integral_label (this, idx)
        class(symbolicelementvector)   :: this
        integer, intent(in)            :: idx
        integer(longint), dimension(2) :: get_integral_label
 
        if (this % check_bounds(idx)) get_integral_label(:) = this % electron_integral(:,idx)
 
    end function get_integral_label
 
 
    real(wp) function get_coefficient (this, idx)
        class(symbolicelementvector) :: this
        integer, intent(in)          :: idx
 
        if (this % check_bounds(idx)) get_coefficient = this % coeffs(idx)
 
    end function get_coefficient
 
 
    subroutine get_coeff_and_integral (this, idx, coeff, label)
        class(symbolicelementvector), intent(in)  :: this
        integer,                      intent(in)  :: idx
        real(wp),                     intent(out) :: coeff
        integer(longint),             intent(out) :: label(2)
 
        if (this % check_bounds(idx)) then
            label(1:2) = this % electron_integral(1:2,idx)
            coeff = this % coeffs(idx)
        end if
 
    end subroutine get_coeff_and_integral
 
 
    integer function get_size(this)
        class(symbolicelementvector) :: this
        get_size = this % n
    end function get_size
 
 
    subroutine destroy (this)
        class(symbolicelementvector) :: this
 
        if (allocated(this % electron_integral)) then
            call master_memory % free_memory(kind(this % electron_integral), size(this % electron_integral))
            deallocate(this % electron_integral)
        end if
        if (allocated(this % coeffs)) then
            call master_memory % free_memory(kind(this % coeffs), size(this % coeffs))
            deallocate(this % coeffs)
        end if
 
        call this % bstree % destroy
 
        this % constructed = .false.
 
    end subroutine destroy
 
 
    subroutine print_symbols (this)
        class(symbolicelementvector) :: this
        integer  :: labels(8), ido
        real(wp) :: coeff
 
        if (.not. this % is_empty()) write (stdout, "('Outputting symbolic elements....')")
        do ido = 1, this % n
            call unpack8ints(this % electron_integral(1,ido), labels)
            write (stdout, "(5i4,' -- ',es14.3)") labels(1:5),this % coeffs(ido)
        end do
 
    end subroutine print_symbols
 
end module symbolic_module