multidip_tests.f90

Go to the documentation of this file.

! Copyright 2020
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-out (UKRmol+ suite).
!
!     UKRmol-out 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-out 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-out (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
!
module multidip_tests

    use precisn_gbl, only: wp

    implicit none

contains

    subroutine run_tests

        print '(/,A)', 'Running all tests'

        call test_permutations
        call test_coul
        call test_1p_eint
        call test_1p_cint
        call test_2p_cint
        call test_2p_gint

    end subroutine run_tests


    subroutine test_permutations

        use multidip_special, only: next_permutation

        integer :: i, order(4)

        order = -1

        print '(/,A,/)', 'Testing permutations of 4 elements...'

        i = 0
        do while (next_permutation(order))
            i = i + 1
            print '(A,I2,A,*(1x,I0))', '  ', i, ':', order
        end do

    end subroutine test_permutations


    subroutine test_coul

        use multidip_params,  only: ntermsasy
        use multidip_special, only: coul_gsl, coul_ukrmol

        real(wp) :: r = 100, f, fp, g1, g1p, g2, g2p, g3, g3p, k, a, b
        real(wp) :: Ek(4) = [ -1.000, -0.100, -0.010, -0.001 ]
        integer  :: ls(3) = [ 0, 1, 2 ], ipw, ie, n, l

        print '(/,A)', 'Testing Whittaker functions'
        print '(/,2x,A,3x,A,7x,A,20x,A,19x,A,17x,A,16x,A,20x,A)', 'l', 'Ek', 'GSL W', 'GSL W''', 'UKRmol W', 'UKRmol W''', &
            'asy W', 'asy W'''

        do ie = 1, size(ek)
            k = sqrt(-2*ek(ie))
            do ipw = 1, size(ls)
                l = ls(ipw)
                call coul_gsl(l, ek(ie), r, f, fp, g1, g1p)
                call coul_ukrmol(l, ek(ie), r, f, fp, g2, g2p)
                f = 1; g3 = 0
                do n = 0, ntermsasy
                    if (n > 0) f = f * (l + 1 - 1/k + n - 1) * (-l - 1/k + n - 1) / n / (-2*k*r)
                    g3 = g3 + f
                    g3p = g3p + (1/k - n)/r * f
                end do
                g3p = exp(-k*r) * (2*k*r)**(1/k) * (-k*g3 + g3p)
                g3 = exp(-k*r) * (2*k*r)**(1/k) * g3
                print '(I3,F8.3,6E25.15)', l, ek(ie), g1, g1p, g2, g2p, g3, g3p
            end do
        end do

    end subroutine test_coul


    subroutine test_1p_eint

        use multidip_integ, only: nested_exp_integ

        real(wp) :: a = 150, c = 0
        real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
        real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]

        integer  :: ie, m, s2, ms(1), ss(2)

        complex(wp) :: ks(2), integ
        complex(wp) :: reference(2,3,2) = reshape([(-4.772929e-01_wp, +8.885647e+01_wp), (-8.424650e-04_wp, +5.923864e-01_wp),  &
                                                   (+9.979744e-06_wp, +3.949186e-03_wp), (+1.307372e+00_wp, +3.240230e+02_wp),  &
                                                   (+3.965361e-02_wp, +2.158463e+00_wp), (+4.697322e-04_wp, +1.437268e-02_wp),  &
                                                   (-2.211202e-01_wp, +7.029786e+01_wp), (-9.899674e-06_wp, +4.686525e-01_wp),  &
                                                   (+9.695004e-06_wp, +3.124290e-03_wp), (-1.752969e+00_wp, +1.992839e+02_wp),  &
                                                   (+8.061387e-05_wp, +1.328557e+00_wp), (+7.894722e-05_wp, +8.855647e-03_wp)], &
                                                   [2, 3, 2])

        print '(/,A)', 'Testing one-photon exponential integrals'
        print '(/,2x,A,1x,A,1x,A,2x,A,12x,A,13x,A,2x,A,3x,A,4x,A)', 'm', 's2', 's1', 'k2', 'k1', &
            're calculated', 'im calculated', 're reference', 'im reference'

        do ie = 1, size(k1)
            do m = 0, 2
                do s2 = 1, 2

                    ms = [ 1 - m ]
                    ss = [ (-1)**(s2+1), +1 ]
                    ks = [ k1(ie), k2(ie) ]

                    integ = nested_exp_integ(a, c, 1, ms, ss, ks)

                    print '(SP,3I3,SS,2F14.10,SP,4E16.7)', ms(1), ss(2), ss(1), real(ks(2),wp), real(ks(1),wp), &
                          real(integ, wp), aimag(integ), real(reference(s2,m+1,ie), wp), aimag(reference(s2,m+1,ie))
                end do
            end do
        end do

    end subroutine test_1p_eint


    subroutine test_1p_cint

        use multidip_integ, only: nested_coul_integ

        real(wp) :: a = 150, c = 0
        real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
        real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]

        integer  :: ie, ipw, s2, ms(1), ss(2), ls(2)
        integer  :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
        integer  :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]

        complex(wp) :: ks(2), integ
        complex(wp) :: reference(2,8,2) = reshape([(-3.454030e+13_wp, +1.290212e+14_wp), (+4.518674e+14_wp, +1.770164e+14_wp),  &
                                                   (-1.054015e+14_wp, -8.206518e+13_wp), (-3.390989e+14_wp, +3.473927e+14_wp),  &
                                                   (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp),  &
                                                   (-1.054015e+14_wp, -8.206518e+13_wp), (-3.390989e+14_wp, +3.473927e+14_wp),  &
                                                   (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp),  &
                                                   (+1.307948e+14_wp, +3.038993e+13_wp), (+1.111367e+14_wp, -4.753884e+14_wp),  &
                                                   (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp),  &
                                                   (+1.307948e+14_wp, +3.038993e+13_wp), (+1.111367e+14_wp, -4.753884e+14_wp),  &
                                                   (-6.382646e+01_wp, +2.799838e+01_wp), (+7.071896e+01_wp, +1.844497e+02_wp),  &
                                                   (+1.584632e+01_wp, -6.787748e+01_wp), (-1.900621e+02_wp, -5.395706e+01_wp),  &
                                                   (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp),  &
                                                   (+1.584632e+01_wp, -6.787748e+01_wp), (-1.900621e+02_wp, -5.395706e+01_wp),  &
                                                   (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp),  &
                                                   (-6.443654e+01_wp, -2.668020e+01_wp), (-1.150300e+02_wp, -1.606738e+02_wp),  &
                                                   (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp),  &
                                                   (-6.443654e+01_wp, -2.668020e+01_wp), (-1.150300e+02_wp, -1.606738e+02_wp)], &
                                                   [2, 8, 2])

        print '(/,A)', 'Testing one-photon Coulomb integrals'
        print '(/,1x,A,1x,A,1x,A,2x,A,12x,A,13x,A,2x,A,3x,A,4x,A)', 's2', 'l2', 'l1', 'k2', 'k1', &
            're calculated', 'im calculated', 're reference', 'im reference'

        do ie = 1, size(k1)
            do ipw = 1, size(l1)
                do s2 = 1, 2

                    ms = [ 1 ]
                    ss = [ (-1)**(s2+1), +1 ]
                    ls = [ l1(ipw), l2(ipw) ]
                    ks = [ k1(ie), k2(ie) ]

                    integ = nested_coul_integ(a, c, 1, ms, ss, ls, ks)

                    print '(SP,I3,SS,2I3,2F14.10,SP,4E16.7)', (-1)**(s2+1), l2(ipw), l1(ipw), k2(ie), k1(ie), &
                          real(integ, wp), aimag(integ), real(reference(s2,ipw,ie), wp), aimag(reference(s2,ipw,ie))
                end do
            end do
        end do

    end subroutine test_1p_cint


    subroutine test_2p_cint

        use multidip_integ, only: nested_coul_integ

        real(wp) :: a = 150, c = 0

        real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
        real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
        real(wp) :: k3(2) = [ 1.511902774652_wp, 1.917220644579_wp ]

        integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
        integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
        integer :: l3(8) = [ 1, 1, 1, 3, 3, 3, 3, 3 ]

        integer :: ie, ipw, s4, ms(2), ss(4), ls(4)

        complex(wp) :: ks(4), integ1, integ2, integ
        complex(wp) :: reference(2,8,2) = reshape([(-1.465883e+16_wp, -4.950701e+15_wp), (-3.842765e+16_wp,  8.620168e+16_wp),  &
                                                   (-1.466146e+16_wp, -4.951551e+15_wp), (-3.843793e+16_wp,  8.623263e+16_wp),  &
                                                   (-1.443664e+16_wp, -5.780723e+15_wp), (-3.378946e+16_wp,  8.865409e+16_wp),  &
                                                   ( 1.520251e+16_wp, -2.902864e+15_wp), (-8.964058e+15_wp, -9.400575e+16_wp),  &
                                                   ( 1.541476e+16_wp, -2.070261e+15_wp), (-1.420349e+16_wp, -9.382680e+16_wp),  &
                                                   ( 1.542117e+16_wp, -2.071208e+15_wp), (-1.422236e+16_wp, -9.390101e+16_wp),  &
                                                   ( 1.541476e+16_wp, -2.070261e+15_wp), (-1.420349e+16_wp, -9.382680e+16_wp),  &
                                                   ( 1.542117e+16_wp, -2.071208e+15_wp), (-1.422236e+16_wp, -9.390101e+16_wp),  &
                                                   (-4.322904e+03_wp, -4.890496e+03_wp), (-2.301503e+04_wp,  1.829823e+04_wp),  &
                                                   (-4.323345e+03_wp, -4.890988e+03_wp), (-2.301886e+04_wp,  1.830140e+04_wp),  &
                                                   ( 6.409175e+03_wp, -1.254518e+03_wp), ( 2.742317e+04_wp,  1.058896e+04_wp),  &
                                                   ( 5.914821e+03_wp,  2.763232e+03_wp), ( 1.382079e+04_wp, -2.596191e+04_wp),  &
                                                   (-5.378829e+03_wp,  3.705000e+03_wp), (-2.937519e+04_wp,  1.216512e+03_wp),  &
                                                   (-5.380092e+03_wp,  3.705885e+03_wp), (-2.938678e+04_wp,  1.217206e+03_wp),  &
                                                   (-5.378829e+03_wp,  3.705000e+03_wp), (-2.937519e+04_wp,  1.216512e+03_wp),  &
                                                   (-5.380092e+03_wp,  3.705885e+03_wp), (-2.938678e+04_wp,  1.217206e+03_wp)], &
                                                   [2, 8, 2])

        print '(/,A)', 'Testing two-photon Coulomb integrals'
        print '(/,4(1x,A),2x,A,12x,A,12x,A,12x,A,3x,A,3x,A,4x,A)', 's4', 'l3', 'l2', 'l1', 'k3', 'k2', 'k1', &
            're calculated', 'im calculated', 're reference', 'im reference'

        do ie = 1, size(k1)
            do ipw = 1, size(l1)
                do s4 = 1, 2

                    ms = [ 1, 1 ]

                    ss = [ (-1)**(s4+1), +1, -1, +1 ]
                    ls = [ l1(ipw), l2(ipw), l2(ipw), l3(ipw) ]
                    ks = [ k1(ie), k2(ie), k2(ie), k3(ie) ]

                    integ1 = nested_coul_integ(a, c, 2, ms, ss, ls, ks)

                    ss = [ +1, +1, -1, (-1)**(s4+1) ]
                    ls = [ l3(ipw), l2(ipw), l2(ipw), l1(ipw) ]
                    ks = [ k3(ie), k2(ie), k2(ie), k1(ie) ]

                    integ2 = nested_coul_integ(a, c, 2, ms, ss, ls, ks)

                    integ = integ1 + integ2

                    print '(SP,I3,SS,3I3,3F14.10,SP,4E16.7)', (-1)**(s4+1), l3(ie), l2(ipw), l1(ipw), k3(ie), k2(ie), k1(ie), &
                          real(integ, wp), aimag(integ), real(reference(s4,ipw,ie), wp), aimag(reference(s4,ipw,ie))
                end do
            end do
        end do

    end subroutine test_2p_cint


    subroutine test_2p_gint

        use multidip_integ, only: nested_cgreen_integ

        real(wp) :: a = 150, c = 0

        real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
        real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
        real(wp) :: k3(2) = [ 1.511902774652_wp, 1.917220644579_wp ]

        integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
        integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
        integer :: l3(8) = [ 1, 1, 1, 3, 3, 3, 3, 3 ]

        integer :: ie, ipw, s4, ms(2), ls(3)

        complex(wp) :: ks(3), integ1, integ2
        complex(wp) :: reference(2,2,8,2) = reshape([(-8.310927e+15_wp,  7.523384e+15_wp), (-5.909915e+16_wp,  6.931934e+16_wp),  &
                                                     ( 5.679523e+16_wp,  4.672527e+16_wp), ( 8.225897e+24_wp,  1.848188e+24_wp),  &
                                                     (-6.506509e+15_wp,  2.066669e+16_wp), (-1.618385e+16_wp,  4.224605e+15_wp),  &
                                                     ( 1.050096e+17_wp,  4.554189e+16_wp), (-4.312479e+24_wp, -7.240598e+24_wp),  &
                                                     (-7.673705e+15_wp,  2.037601e+16_wp), (-1.647265e+16_wp,  3.342092e+15_wp),  &
                                                     ( 1.078882e+17_wp,  3.987695e+16_wp), (-4.728274e+24_wp, -7.026109e+24_wp),  &
                                                     (-4.503787e+15_wp, -2.119630e+16_wp), ( 1.619053e+16_wp,  4.287343e+15_wp),  &
                                                     (-1.138495e+17_wp,  1.201762e+16_wp), ( 2.556995e+23_wp,  8.432496e+24_wp),  &
                                                     (-3.344463e+15_wp, -2.151751e+16_wp), ( 1.600768e+16_wp,  5.199000e+15_wp),  &
                                                     (-1.135684e+17_wp,  1.836697e+16_wp), ( 7.238276e+23_wp,  8.446832e+24_wp),  &
                                                     (-8.612323e+15_wp, -1.723199e+16_wp), ( 3.535435e+16_wp,  4.008848e+16_wp),  &
                                                     (-9.807095e+16_wp,  3.755844e+16_wp), ( 7.275593e+24_wp,  4.341652e+24_wp),  &
                                                     (-3.344463e+15_wp, -2.151751e+16_wp), ( 1.600768e+16_wp,  5.199000e+15_wp),  &
                                                     (-1.135684e+17_wp,  1.836697e+16_wp), ( 7.238276e+23_wp,  8.446832e+24_wp),  &
                                                     (-8.612323e+15_wp, -1.723199e+16_wp), ( 3.535435e+16_wp,  4.008848e+16_wp),  &
                                                     (-9.807095e+16_wp,  3.755844e+16_wp), ( 7.275593e+24_wp,  4.341652e+24_wp),  &
                                                     (-1.247516e+03_wp,  3.118370e+03_wp), (-2.763140e+04_wp,  1.784564e+04_wp),  &
                                                     ( 1.332880e+04_wp,  9.920923e+03_wp), ( 6.345060e+03_wp,  4.578023e+04_wp),  &
                                                     (-3.896753e+03_wp,  9.080720e+02_wp), (-6.662204e+03_wp,  3.014237e+04_wp),  &
                                                     ( 6.717303e+03_wp,  1.712697e+04_wp), ( 4.397613e+04_wp, -1.193480e+04_wp),  &
                                                     ( 1.189054e+03_wp, -3.822275e+03_wp), (-2.263210e+04_wp, -2.099869e+04_wp),  &
                                                     ( 1.136951e+04_wp, -1.445428e+04_wp), (-3.253996e+04_wp, -3.189015e+04_wp),  &
                                                     ( 3.212357e+03_wp, -2.385969e+03_wp), ( 1.813020e+04_wp, -2.499891e+04_wp),  &
                                                     (-1.298806e+04_wp, -1.303291e+04_wp), (-4.510261e+04_wp, -6.582287e+03_wp),  &
                                                     ( 4.327121e+02_wp,  3.979850e+03_wp), ( 1.239242e+04_wp,  2.828936e+04_wp),  &
                                                     (-4.667920e+03_wp,  1.779027e+04_wp), ( 1.713954e+04_wp,  4.222912e+04_wp),  &
                                                     ( 1.147174e+03_wp,  5.355515e+03_wp), ( 1.429665e+04_wp,  1.753491e+04_wp),  &
                                                     (-5.010378e+03_wp,  2.217364e+04_wp), ( 3.781810e+04_wp,  1.923663e+04_wp),  &
                                                     ( 4.327121e+02_wp,  3.979850e+03_wp), ( 1.239242e+04_wp,  2.828936e+04_wp),  &
                                                     (-4.667920e+03_wp,  1.779027e+04_wp), ( 1.713954e+04_wp,  4.222912e+04_wp),  &
                                                     ( 1.147174e+03_wp,  5.355515e+03_wp), ( 1.429665e+04_wp,  1.753491e+04_wp),  &
                                                     (-5.010378e+03_wp,  2.217364e+04_wp), ( 3.781810e+04_wp,  1.923663e+04_wp)], &
                                                     [2,2,8,2])

        print '(/,A)', 'Testing two-photon Green integrals'
        print '(/,5(1x,A),2x,A,12x,A,12x,A,12x,A,3x,A,3x,A,4x,A)', 's4', 's1', 'l3', 'l2', 'l1', 'k3', 'k2', 'k1', &
            're calculated', 'im calculated', 're reference', 'im reference'

        do ie = 1, size(k1)
            do ipw = 1, size(l1)
                do s4 = 1, 2

                    ms = [ 1, 1 ]
                    ls = [ l1(ipw), l2(ipw), l3(ipw) ]
                    ks = [ k1(ie), k2(ie), k3(ie) ]

                    integ1 = nested_cgreen_integ(a, c, 2, (-1)**(s4+1), +1, ms, ls, ks)
                    integ2 = nested_cgreen_integ(a, c, 2, (-1)**(s4+1), -1, ms, ls, ks)

                    print '(SP,2I3,SS,3I3,3F14.10,SP,8E16.7)', (-1)**(s4+1), +1, l3(ipw), l2(ipw), l1(ipw), k3(ie), k2(ie), k1(ie),&
                          real(integ1, wp), aimag(integ1), real(reference(1,s4,ipw,ie), wp), aimag(reference(1,s4,ipw,ie))
                    print '(SP,2I3,SS,3I3,3F14.10,SP,8E16.7)', (-1)**(s4+1), -1, l3(ipw), l2(ipw), l1(ipw), k3(ie), k2(ie), k1(ie),&
                          real(integ2, wp), aimag(integ2), real(reference(2,s4,ipw,ie), wp), aimag(reference(2,s4,ipw,ie))
                end do
            end do
        end do
    end subroutine test_2p_gint

end module multidip_tests