ALCHEMY_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/>.
 
!> \brief   ALCHEMY integral module
!> \authors A Al-Refaie
!> \date    2017
!>
!> \note 16/01/2019 - Jakub Benda: Unifom coding style and expanded documentation.
!>
module alchemy_module
 
    use precisn,              only: longint, wp
    use const_gbl,            only: stdout
    use consts_mpi_ci,        only: nidx
    use baseintegral_module,  only: baseintegral
    use memorymanager_module, only: master_memory
    use options_module,       only: options
    use global_utils,         only: indfunc
    use scatci_routines,      only: rdint, rdints
    use utility_module,       only: pack_ints, unpack_ints
 
    implicit none
 
    public alchemyintegral
 
    private
 
    type, extends(baseintegral) :: alchemyintegral
 
        !Our integrals
        real(wp), pointer :: one_electron_integral(:)
        real(wp), pointer :: two_electron_integral(:)
        integer           :: num_one_electron_integrals
        integer           :: num_two_electron_integrals
        integer           :: integral_ordering
        integer           :: use_scf = 0
 
        !> How many IJKL pairs we have
        integer              :: num_unique_pairs
        !> The list of unique labels
        integer(longint), allocatable :: pair_labels(:,:)
        !> The number of labels per symmetry
        integer, allocatable :: num_orbitals_sym(:)
 
        integer              :: max_number_pair_sets
        integer, allocatable :: num_two_electron_blocks(:)
        integer, allocatable :: num_one_electron_blocks(:)
        integer, allocatable :: one_electron_pointer(:)
        integer, allocatable :: two_electron_pointer(:)
 
        integer :: num_pq, num_rs
        integer :: num_pair_idx
 
        !>the pair id
        integer, allocatable :: pair_idx(:)
        integer, allocatable :: orbital_idx(:)
        integer, allocatable :: symmetry_idx(:)
 
    contains
 
        procedure, public  :: initialize_self    => initialize_alchemy
        procedure, public  :: finalize_self      => finalize_alchemy
        procedure, public  :: load_integrals     => load_integrals_alchemy
        procedure, public  :: get_integral_ijklm => get_integral_alchemy
        procedure, public  :: destroy_integrals  => destroy_integral_alchemy
        procedure, public  :: write_geometries   => write_geometries_alchemy
        procedure, private :: count_num_pairs
        procedure, private :: generate_pairs
        procedure, private :: generate_pointer_table
        procedure, private :: generate_pair_index
        procedure, private :: generate_orbital_index
       !procedure, private :: read_one_electron_integrals
       !procedure, private :: read_two_electron_integrals
        procedure, private :: get_one_electron_index
        procedure, private :: get_two_electron_index
 
    end type alchemyintegral
 
contains
 
    !> \brief   ?
    !> \authors A Al-Refaie
    !> \date    2017
    !>
    integer function count_num_pairs (this)
        class(alchemyintegral) :: this
        integer :: pair_number, mdo, ndo, ido, jdo, nn(5), msym, begin, pass, istart, num_sym
 
        num_sym = this % num_symmetries
        msym = num_sym + num_sym - 1
        count_num_pairs = 0
 
        if (this % use_SCF == 0) then
            !Standard Alchemy list
            do mdo = 1, msym
                begin = mdo / 2 + 1
                nn(1) = mdo - 1
                do ndo = begin, num_sym
                    nn(2) = ndo - 1
                    nn(3) = abs(ndo - mdo)
                    do ido = begin, ndo
                        nn(4) = ido - 1
                        nn(5) = abs(ido - mdo)
                        count_num_pairs = count_num_pairs + 1
                    end do
                end do
            end do
        else
            do pass = 1, 2
                do mdo = pass, msym
                    begin = mdo / 2 + 1
                    nn(1) = mdo - 1
                    do ndo = mdo / 2 + pass, num_sym
                        nn(2) = ndo - 1
                        nn(3) = abs(ndo - mdo)
                        istart = ndo + 1 - pass
                        if (mdo /= 1 .and. pass == 1) begin = istart
                        do ido = begin, istart
                            nn(4) = ido - 1
                            nn(5) = abs(ido - mdo)
                            count_num_pairs = count_num_pairs + 1
                        end do
                    end do
                end do
                msym = msym - 2
            end do
        end if
 
    end function count_num_pairs
 
 
    !> \brief   ?
    !> \authors A Al-Refaie
    !> \date    2017
    !>
    subroutine generate_pairs (this)
        class(alchemyintegral) :: this
        integer :: pair_number, npqrs, nn(5), mdo, ndo, ido, jdo, I, IBGN, IF, IPASS, J, K, M, MSYM, N, NBGN, begin, pass, num_sym
        integer(longint) :: packed_label(NIDX)
        integer,allocatable :: lpqrs(:,:)
 
        num_sym = this % num_symmetries
        msym = num_sym + num_sym - 1
 
        allocate(lpqrs(5, this % num_unique_pairs))
 
        pair_number = 0
 
        if (this % use_SCF == 0) then
            !Standard Alchemy list
            do mdo = 1, msym
                begin = mdo / 2 + 1
                nn(1) = mdo - 1
                do ndo = begin, num_sym
                    nn(2) = ndo - 1
                    nn(3) = abs(ndo - mdo)
                    do ido = begin, ndo
                        nn(4) = ido - 1
                        nn(5) = abs(ido - mdo)
                        pair_number = pair_number + 1
                        do jdo = 1, 5
                            lpqrs(jdo,pair_number) = nn(jdo)
                        end do
                    end do
                end do
            end do
        else
            do ipass = 1, 2
                do m = ipass, msym
                    ibgn = m / 2 + 1
                    nn(1) = m - 1
                    do n = m / 2 + ipass, num_sym
                        nn(2) = n - 1
                        nn(3) = abs(n - m)
                        IF = n + 1 - ipass
                        if (m /= 1 .and. ipass == 1) ibgn = IF
                        do i = ibgn, IF
                            nn(4) = i - 1
                            nn(5) = abs(i - m)
                            pair_number = pair_number + 1
                            do j = 1, 5
                                lpqrs(j,pair_number) = nn(j)
                            end do
                        end do
                     end do
                 end do
                 msym = msym - 2
            end do
        end if
 
        if (pair_number /= this % num_unique_pairs) stop "Error in pair count"
 
        do ido = 1, this % num_unique_pairs
            call pack_ints(lpqrs(2,ido), lpqrs(3,ido), lpqrs(4,ido), lpqrs(5,ido), 0, 0, 0, 0, packed_label)
            if (lpqrs(2,ido) + lpqrs(3,ido) /= lpqrs(1,ido)) packed_label(1) = -packed_label(1)
            this % pair_labels(:,ido) = packed_label
        end do
 
        write (stdout, "(' NMPRS=',I5//' LPQRS'/(2X,5(2X,I3)))") this % num_unique_pairs, &
            ((lpqrs(ido,jdo), ido = 1, 5), jdo = 1, this % num_unique_pairs)
 
        deallocate(lpqrs)
 
    end subroutine generate_pairs
 
 
    subroutine generate_pointer_table (this)
        class(alchemyintegral) :: this
        integer                :: label(8), num_orbs(4), size_two_pointer
        integer(longint)       :: packed_label(NIDX)
        integer                :: ido, jdo, nam1, nam2, num_pq, num_rs, num_pr
        integer                :: block_number, mvl, md, mm, m, mpa, mpb, mp
 
        !Lets allocate our one elctron pointers
 
        size_two_pointer = this % num_symmetries * (this % num_symmetries + 1) / 2
        do ido = 2, this % num_symmetries
            size_two_pointer = size_two_pointer + (this % num_symmetries - ido + 2) * (this % num_symmetries - ido + 1)
        end do
 
        allocate(this % num_two_electron_blocks(2 * this % num_symmetries + 1))
 
        this % num_two_electron_blocks(1) = 0
        this % num_two_electron_blocks(2) = indfunc(this % num_symmetries + 1, 0)
 
        mvl = 2 * this% num_symmetries - 1
 
        do ido = 2, mvl
            md = (ido - 2) / 2
            mm = this % num_symmetries - md - 1
            this% num_two_electron_blocks(ido + 1) = this % num_two_electron_blocks(ido) + indfunc(mm + 1, 0)
        end do
 
        allocate(this % one_electron_pointer(this % num_symmetries + 1))
 
        this % one_electron_pointer(1) = 1
 
        do ido = 1, this % num_symmetries
            this % one_electron_pointer(ido + 1) = this % one_electron_pointer(ido) + indfunc(this % num_orbitals_sym(ido) + 1, 0)
        end do
 
        allocate(this % two_electron_pointer(this % num_two_electron_blocks(mvl + 1)))
 
        do ido =1, this % num_two_electron_blocks(mvl + 1)
            this % two_electron_pointer(ido) = 0
        end do
 
        this % num_one_electron_integrals = this % one_electron_pointer(this % num_symmetries + 1) - 1
        this % num_two_electron_integrals = 1
 
        write (stdout, "(//,5X,'Integral Storage Table follows : ',/)")
 
        do ido = 1, this % num_unique_pairs
            !> Get our lsabels
            packed_label(:) = this % pair_labels(:,ido)
 
            if (packed_label(1) > 0) then
                md = 0
            else
                md = 1
                packed_label(1) = -packed_label(1)
            end if
 
            call unpack_ints(packed_label, label)
 
            do jdo = 1, 4
                num_orbs(jdo) = this % num_orbitals_sym(label(jdo) + 1)
            end do
 
            if (md == 0) then
                m = label(1) + label(2)
            else
                m = abs(label(1) - label(2))
            end if
 
            if (m == 0) then
                md = -1
            else
                md = (m - 1) / 2
            end if
            mpa = label(1) - md
            mpb = label(3) - md
            if (mpa < mpb) then
                !MP=(MPB*(MPB-1))/2+MPA+this%num_two_electron_blocks(M+1)
                mp = (mpb * (mpb - 1)) / 2 + mpa + this % num_two_electron_blocks(m + 1)
            else
                !MP=(MPA*(MPA-1))/2+MPB+this%num_two_electron_blocks(M+1)
                mp = (mpa * (mpa - 1)) / 2 + mpb + this % num_two_electron_blocks(m + 1)
            end if
 
            this % two_electron_pointer(mp) = this % num_two_electron_integrals
            !write(stdout,"('ELEC ',3i8)") mp,this%num_two_electron_integrals,this%two_electron_pointer(mp)
            if (label(1) == label(2)) then
                num_pq = indfunc(num_orbs(1) + 1, 0)
            else
                num_pq = num_orbs(1) * num_orbs(2)
            end if
 
            if (label(3) == label(4)) then
                num_rs = indfunc(num_orbs(3) + 1, 0)
            else
                num_rs = num_orbs(3) * num_orbs(4)
            end if
            this % max_number_pair_sets = max(this % max_number_pair_sets, num_pq, num_rs)
            if (label(1) == label(3) .and. label(2) == label(4)) then
                this % num_two_electron_integrals = this % num_two_electron_integrals + indfunc(num_pq + 1, 0)
            else
                this % num_two_electron_integrals = this % num_two_electron_integrals + num_pq * num_rs
            end if
        end do
        !Move to the end of the block
        this % num_two_electron_integrals = this % num_two_electron_integrals - 1
 
        write (stdout, "('No of 1 electron integrals = ',i8)") this % num_one_electron_integrals
        write (stdout, "('No of 2 electron integrals = ',i8)") this % num_two_electron_integrals
 
    end subroutine generate_pointer_table
 
 
    subroutine generate_orbital_index (this)
        class(alchemyintegral) :: this
        integer :: total_num_orbitals, orbital, iorb, isym, max_orbitals, ifail
 
        total_num_orbitals = sum(this % num_orbitals_sym(:))
 
        allocate(this % orbital_idx(total_num_orbitals), &
                 this % symmetry_idx(total_num_orbitals), stat = ifail)
 
        call master_memory % track_memory(storage_size(this%orbital_idx)/8,  size(this % orbital_idx),  ifail, 'SWEDEN::pair_idx')
        call master_memory % track_memory(storage_size(this%symmetry_idx)/8, size(this % symmetry_idx), ifail, &
                                          'SWEDEN::symmetry_idx')
 
        this % orbital_idx(:) = 0
        this % symmetry_idx(:) = 0
 
        max_orbitals = 0
        orbital = 0
 
        do isym = 1,this % num_symmetries
            max_orbitals = max(max_orbitals, this % num_orbitals_sym(isym)**2)
            do iorb = 1, this % num_orbitals_sym(isym)
                orbital = orbital + 1
                this % orbital_idx(orbital) = iorb
                this % symmetry_idx(orbital) = isym - 1
            end do
        end do
 
        this % num_pair_idx = max((this% num_symmetries**2 + this % num_symmetries) / 2 + 1, &
                                  this % max_number_pair_sets, &
                                  300, &
                                  max_orbitals)
 
    end subroutine generate_orbital_index
 
 
    subroutine generate_pair_index (this)
        class(alchemyintegral) :: this
        integer                :: ido, idx, ifail
 
        allocate(this % pair_idx(0:this % num_pair_idx), stat = ifail)
 
        call master_memory % track_memory(storage_size(this % pair_idx)/8, size(this % pair_idx), ifail, 'SWEDEN::pair_idx')
 
        idx = 0
        this % pair_idx(0) = 0
        do ido = 1, this % num_pair_idx
            this % pair_idx(ido) = idx
            idx = idx + ido
        end do
 
    end subroutine generate_pair_index
 
 
    subroutine initialize_alchemy (this, option)
        class(alchemyintegral)     :: this
        class(options), intent(in) :: option
        integer                    :: ido, jdo
 
        write (stdout, "('Using ALCHEMY integral')")
 
        this % num_unique_pairs = 0
        this % num_symmetries = option % num_syms
        this % integral_ordering = option % integral_ordering
        this % use_SCF = option % use_SCF
 
        !Lets begin counting
 
        allocate(this % num_orbitals_sym(this % num_symmetries))
 
        this % num_orbitals_sym(:) = option % num_electronic_orbitals_sym(:)
        this % num_unique_pairs = this % count_num_pairs()
 
        write (stdout, *) 'Number of pairs =',this % num_unique_pairs
 
        allocate(this % pair_labels(2, this % num_unique_pairs))
 
        write (stdout,*) 'ALCHEMY Integral -- Generating pairs'
 
        !>Generate indexes
        call this % generate_pairs
        call this % generate_orbital_index
        call this % generate_pair_index
 
        write (stdout, "('1',/,10x,'D2h Two Electron Integral Box ','Information')")
        write (stdout, "(    /,10x,'No. of Boxes = ',i4,/)") this % num_unique_pairs
        write (stdout, "(      10x,'Box Descriptions: ',/)")
        write (stdout, *) 'ALCHEMY Integral -- Generating block pointers'
 
        call this % generate_pointer_table
 
    end subroutine initialize_alchemy
 
 
    subroutine finalize_alchemy (this)
        class(alchemyintegral) :: this
    end subroutine finalize_alchemy
 
 
    !> This is just a copy from scatci_routines with only the relevant ALCHEMY parts.
    !> This will be replaced when the prototype is completed with a 'caching' system.
    subroutine load_integrals_alchemy (this, iounit)
        use params,          only: ctrans1, cpoly
        use global_utils,    only: intape
        use scatci_routines, only: search, chn2e
 
        class(alchemyintegral)    :: this
        integer, intent(in)       :: iounit
        !     ALCHEMY header arrays
        character(len=4), dimension(4)  :: LABEL
        character(len=4), dimension(33) :: NAM1
        integer, dimension(8)  :: NAO, NCORE
        integer, dimension(20) :: NHE, NOB
        real(wp), allocatable           :: xtempe(:), xtempp(:)
 
        integer  :: I, NSYM, IONEIN, NT, ND, LTRI, NNUC, N, IA, IAT, NALM
        integer  :: ifail, k, l, ind, IMAX, IMIN, NEND, IMAXP
        real(wp) :: EN, SIGN, x1e
        real(wp) :: vsmall = 1.e-10_wp
        real(kind=wp), dimension(20) :: charg, xnuc
 
        write (stdout, *) ' Loading ALCHEMY Integrals into core'
 
        !allocate the integral space, this will be replaced by a caching system
        !One electron integral
        allocate(this % one_electron_integral(2 * this % num_one_electron_integrals), stat = ifail)
        if (ifail /= 0) stop "could not allocate 1 electron integral"
 
        !Two electron integral
        allocate(this % two_electron_integral(this % num_two_electron_integrals), stat = ifail)
        if (ifail /= 0) stop "could not allocate 2 electron integral"
 
        read (iounit) (nam1(i), i = 1, 33), nsym, nt, nnuc, nd, ltri, (nob(i), i = 1, nsym), (nd, i = 1, nt), (nd, i = 1, nt), &
                      (charg(i), i = 1, nnuc), (xnuc(i), i = 1, nnuc)
 
        write (stdout, "(/' Transformed integrals read:',/5X,30A4)") (nam1(i), i = 1, 30)
 
        en = 0.0_wp
        do n = 2, nnuc
            if (abs(charg(n)) < vsmall) cycle
            do i = 1, n - 1
               if (abs(charg(i)) < vsmall) cycle
               en = en + charg(n) * charg(i) / abs(xnuc(n) - xnuc(i))
            end do
        end do
 
        this % core_energy = en
 
        write (stdout, "(' NSYM =',I5,3X,'NOB  =',20I5)") nsym, (nob(i), i = 1, nsym)
        write (stdout, "(/,10x,'Nuclear repulsion energy = ',f15.7)") en
        write (stdout, "(' NNUC =',I5,3X,'CHARG=',10F10.0/21X,10F10.0)") nnuc, (charg(i), i = 1, nnuc)
        write (stdout, "(15X,'XNUC =',10F10.5/21X,10F10.5)") (xnuc(i), i = 1, nnuc)
 
        ia = 1
        call rdints(iounit, nsym, nob, ltri, this % num_one_electron_integrals + 1, ia, this % one_electron_integral, nalm)
        ! print*,this%one_electron_integral(1)
        if (nalm /= 0) stop "UNABLE TO READ ALCHEMY INTEGRALS"
 
        ia = 1
        call rdints(iounit, nsym, nob, ltri, this % num_one_electron_integrals + 1, ia, this % one_electron_integral, nalm)
        ! print*,this%one_electron_integral(1)
        if (nalm /= 0) stop "UNABLE TO READ ALCHEMY INTEGRALS"
 
        iat = ia - 1
        ia = 1
        allocate(xtempp(this % num_one_electron_integrals))
        call rdints(iounit, nsym, nob, ltri, this % num_one_electron_integrals + 1, ia, xtempp, nalm)
        ! print*,xtempp(1)
 
        if (nalm == 0) then
            if (iat /= this % num_one_electron_integrals) stop "INCORRECT NUMBER OF INTEGRALS"
            sign = 1.0
            if(abs(this % positron_flag) == 1) sign = -1.0
            do i = 1, iat
                x1e = this % one_electron_integral(i) + xtempp(i)
                this % one_electron_integral(i) = x1e
            end do
            !this%one_electron_integral(1:) = this%one_electron_integral(:) + xtempp(:)
        end if
 
        !print *,this%one_electron_integral(1)
 
        imax = 0
        imin = 0
        nend = -1
        imaxp = 0
        ia = 1
        call rdint(imin, imax, nend, imaxp, iounit, ia, this % num_two_electron_integrals + 1, this % two_electron_integral)
        if (nend /= 1) stop "UNABLE TO READ ALCHEMY INTEGRALS"
 
        deallocate(xtempp)
 
        !endif
 
        1600 FORMAT(' Integrals read successfully:',       /, &
                    ' 1-electron integrals, NINT1e =',i10, /, &
                    ' 2-electron integrals, NINT2e =',i10)
        write (stdout, 1600) this % num_one_electron_integrals, this % num_two_electron_integrals
 
        this % nhe(:) = nhe(:)
        this % nnuc = nnuc
        this % dtnuc(:) = 0
        this % dtnuc(1) = en
        this % nhe(:) = nob(:)
        this % dtnuc(22:21+nnuc) = xnuc(1:nnuc)
        this % dtnuc(2:1+nnuc) = charg(1:nnuc)
 
    end subroutine load_integrals_alchemy
 
 
    function get_integral_alchemy (this, i, j, k, l, m) result(coeff)
        class(alchemyintegral) :: this
        integer, intent(in)    :: i,j,k,l,m
        real(wp) :: coeff
        integer  :: ia, ib, integral_idx
 
        coeff = 0.0
 
        !One electron case
        if (i == 0) then
            integral_idx = this % get_one_electron_index(j, l, m)
            coeff = this % one_electron_integral(integral_idx)
            !print *,coeff
        else
            integral_idx = this % get_two_electron_index(i, j, k, l, m)
            coeff = this % two_electron_integral(integral_idx)
        end if
 
        !WRITE(stdout,*)i, j, k, l, m
        !write(stdout,*)coeff,integral_idx
        !write(stdout,*) i,j,k,l,m,coeff,integral_idx
 
    end function get_integral_alchemy
 
 
    integer function get_one_electron_index (this, i, j, pos)
        class(alchemyintegral), intent(in) :: this
        integer,                intent(in) :: i, j, pos
        integer :: m, p, q, ii, jj
 
        ii = this % orbital_mapping(i)
        jj = this % orbital_mapping(j)
 
        p = this % orbital_idx(ii)
        q = this % orbital_idx(jj)
 
        m = this % symmetry_idx(ii)
 
        if (p < q) then
            get_one_electron_index = this % pair_idx(q) + p + this % one_electron_pointer(m + 1) - 1
        else
            get_one_electron_index = this % pair_idx(p) + q + this % one_electron_pointer(m + 1) - 1
        end if
 
        if (pos /= 0) get_one_electron_index = get_one_electron_index + this % num_one_electron_integrals
        !write(stdout,"(10i8)") 0,p,0,q,0,get_one_electron_index,m,this%one_electron_pointer(m+1),this%pair_idx(p)
 
    end function get_one_electron_index
 
 
    integer function get_two_electron_index (this, i, j, k, l, m)
        class(alchemyintegral), intent(in) :: this
        integer,                intent(in) :: i, j, k, l, m
        integer :: symmetry, mpp(4), nbp(4), npp(4), mapped(4), npq(4), ido, block_pointer, ipqrs, mv, md, mpq, mrs, mpr
 
        mapped(1) = this % orbital_mapping(i)
        mapped(2) = this % orbital_mapping(j)
        mapped(3) = this % orbital_mapping(k)
        mapped(4) = this % orbital_mapping(l)
 
        do ido = 1, 4
            npp(ido) = this % orbital_idx(mapped(ido))
            symmetry = this % symmetry_idx(mapped(ido))
            mpp(ido) = symmetry
            nbp(ido) = this % num_orbitals_sym(symmetry + 1)
        end do
 
        if (m /= 0) then
            mv = mpp(1) + mpp(2)
        else
            mv = abs(mpp(1) - mpp(2))
        end if
 
        if (mv /= 0) then
            md = (mv - 1) / 2
        else
            md = -1
        end if
 
        mpq = mpp(1) - md
        mrs = mpp(3) - md
        mpr = this % pair_idx(mpq) + mrs + this % num_two_electron_blocks(mv + 1)
 
        block_pointer = this % two_electron_pointer(mpr) - 1
 
        do ido = 1, 4, 2
            if (mpp(ido) /= mpp(ido + 1)) then
                npq(ido)     = nbp(ido) * (npp(ido + 1) - 1) + npp(ido)
                npq(ido + 1) = nbp(ido) *  nbp(ido + 1)
            else
                if (npp(ido) < npp(ido + 1)) then
                    npq(ido) = this % pair_idx(npp(ido + 1)) + npp(ido)
                else
                    npq(ido) = this % pair_idx(npp(ido))     + npp(ido + 1)
                end if
                npq(ido + 1) = this % pair_idx(nbp(ido) + 1)
            end if
        end do
 
        if (mpp(1) == mpp(3) .and. mpp(2) == mpp(4)) then
            if (npq(1) < npq(3)) then
                ipqrs = this % pair_idx(npq(3)) + npq(1)
            else
                ipqrs = this % pair_idx(npq(1)) + npq(3)
            end if
        else
            if (this % integral_ordering == 0) then
                ipqrs = npq(2) * (npq(3) - 1) + npq(1)
            else
                ipqrs = npq(4) * (npq(1) - 1) + npq(3)
            end if
        end if
 
        get_two_electron_index = ipqrs + block_pointer
        !print *,'ALC ',block_pointer,     IPQRS
        !write(stdout,"(6i8)") i,j,k,l,m,get_two_electron_index
 
    end function get_two_electron_index
 
 
    subroutine write_geometries_alchemy (this, iounit)
        use params,          only : cpoly
        use scatci_routines, only : search
 
        class(alchemyintegral), intent(in) :: this
        integer,                intent(in) :: iounit
        integer :: ido, ifail
 
        !REWIND iounit
        !CALL SEARCH(iounit,cpoly,ifail)
        !READ(iounit)nnuc
        ! ALLOCATE(xnuc(nnuc),ynuc(nnuc),znuc(nnuc),charge(nnuc),CNAME(nnuc))
        !DO i=1, nnuc
        !         READ(iounit)cname(i), ii, xnuc(i), ynuc(i), znuc(i), charge(i)
        !END DO
 
    end subroutine write_geometries_alchemy
 
 
    subroutine destroy_integral_alchemy (this)
        class(alchemyintegral) :: this
 
        if (associated(this % one_electron_integral)) deallocate(this % one_electron_integral)
        if (associated(this % two_electron_integral)) deallocate(this % two_electron_integral)
 
        if (allocated(this % pair_labels)) deallocate(this % pair_labels)
 
        !> The number of labels per symmetry
        if (allocated(this % num_orbitals_sym)) deallocate(this % num_orbitals_sym)
 
        if (allocated(this % one_electron_pointer)) deallocate(this % one_electron_pointer)
        if (allocated(this % two_electron_pointer)) deallocate(this % two_electron_pointer)
 
        if (allocated(this % pair_idx)) deallocate(this % pair_idx)
        if (allocated(this % orbital_idx)) deallocate(this % orbital_idx)
        if (allocated(this % symmetry_idx)) deallocate(this % symmetry_idx)
 
        if (allocated(this % num_two_electron_blocks)) deallocate(this % num_two_electron_blocks)
        if (allocated(this % num_one_electron_blocks)) deallocate(this % num_one_electron_blocks)
 
    end subroutine destroy_integral_alchemy
 
end module alchemy_module