From 8669b46de1a561fc1760554944c503ff228d8119 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Thu, 17 Oct 2024 00:25:02 +0200 Subject: [PATCH 1/5] improv cpx erf calculation --- src/hartree_fock/deb_ao_2e_int.irp.f | 102 +++++++++++++++++++- src/utils/cpx_boys.irp.f | 138 +++++++++++++++++---------- src/utils/cpx_erf.irp.f | 107 +++++++++++++++------ 3 files changed, 266 insertions(+), 81 deletions(-) diff --git a/src/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f index 29608603..94dc0bda 100644 --- a/src/hartree_fock/deb_ao_2e_int.irp.f +++ b/src/hartree_fock/deb_ao_2e_int.irp.f @@ -3,10 +3,11 @@ program deb_ao_2e_int implicit none - call check_ao_two_e_integral_cgtos() + !call check_ao_two_e_integral_cgtos() !call check_crint1() !call check_crint2() !call check_crint3() + call check_crint4() end @@ -344,3 +345,102 @@ end ! --- +subroutine check_crint4() + + implicit none + + integer :: i_test, n_test + integer :: i, seed_size, clock_time + double precision :: xr(1), x, shift + double precision :: yr(1), y + double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im + double precision :: delta_ref + double precision :: t1, t2, t_int1, t_int2 + complex*16 :: rho + complex*16 :: int1, int2, int3 + integer, allocatable :: seed(:) + + + + call random_seed(size=seed_size) + allocate(seed(seed_size)) + call system_clock(count=clock_time) + seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /) + !seed = 123456789 + call random_seed(put=seed) + + + t_int1 = 0.d0 + t_int2 = 0.d0 + + n_test = 100 + shift = 15.d0 + + acc_re = 0.d0 + nrm_re = 0.d0 + acc_im = 0.d0 + nrm_im = 0.d0 + do i_test = 1, n_test + + call random_number(xr) + call random_number(yr) + + x = 1.d0 * (2.d0 * shift * xr(1) - shift) + y = 1.d0 * (2.d0 * shift * yr(1) - shift) + + rho = x + (0.d0, 1.d0) * y + + call wall_time(t1) + call zboysfun00_1(rho, int1) + call wall_time(t2) + t_int1 = t_int1 + t2 - t1 + + call wall_time(t1) + call zboysfun00_2(rho, int2) + call wall_time(t2) + t_int2 = t_int2 + t2 - t1 + + dif_re = dabs(real(int1) - real(int2)) + dif_im = dabs(aimag(int1) - aimag(int2)) + if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then + print*, ' important error found: ' + print*, " rho = ", x, y + print*, real(int1), real(int2), dif_re + print*, aimag(int1), aimag(int2), dif_im + call crint_quad_12(0, rho, 10000000, int3) + if(zabs(int1 - int3) .lt. zabs(int2 - int3)) then + print*, ' implementation 2 seems to be wrong' + else + print*, ' implementation 1 seems to be wrong' + print*, ' quad 10000000:', real(int3), aimag(int3) + call crint_quad_12(0, rho, 100000000, int3) + print*, ' quad 100000000:', real(int3), aimag(int3) + endif + !print*, ' quad:', real(int3), aimag(int3) + !stop + endif + + if((real(int1) /= real(int1)) .or. (aimag(int1) /= aimag(int1)) .or. & + (real(int2) /= real(int2)) .or. (aimag(int2) /= aimag(int2)) ) then + cycle + else + acc_re += dif_re + acc_im += dif_im + nrm_re += dabs(real(int1)) + nrm_im += dabs(aimag(int1)) + endif + enddo + + print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re + 1d-15) + print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im + 1d-15) + + print*, "zerf_1 wall time (sec) = ", t_int1 + print*, "zerf_2 wall time (sec) = ", t_int2 + + + deallocate(seed) + +end + +! --- + diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f index 687433e4..acac1f6f 100644 --- a/src/utils/cpx_boys.irp.f +++ b/src/utils/cpx_boys.irp.f @@ -10,44 +10,100 @@ complex*16 function crint_1(n, rho) complex*16, intent(in) :: rho integer :: i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im - double precision :: n_tmp + double precision :: rho_mod + double precision :: tmp complex*16 :: sq_rho, rho_inv, rho_exp - complex*16 :: crint_smallz, cpx_erf_1 + complex*16 :: crint_smallz - rho_re = real (rho) - rho_im = aimag(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) + rho_mod = zabs(rho) + + if(rho_mod < 3.5d0) then + + if(rho_mod .lt. 0.35d0) then + + select case(n) + case(0) + crint_1 = (((((((((1.3122532963802805073d-08 * rho & + - 1.450385222315046877d-07) * rho & + + 1.458916900093370682d-06) * rho & + - 0.132275132275132275d-04) * rho & + + 0.106837606837606838d-03) * rho & + - 0.757575757575757576d-03) * rho & + + 0.462962962962962963d-02) * rho & + - 0.238095238095238095d-01) * rho & + + 0.10000000000000000000d0) * rho & + - 0.33333333333333333333d0) * rho & + + 1.0d0 + case(1) + crint_1 = (((((((((1.198144314086343d-08 * rho & + - 1.312253296380281d-07) * rho & + + 1.305346700083542d-06) * rho & + - 1.167133520074696d-05) * rho & + + 9.259259259259259d-05) * rho & + - 6.410256410256410d-04) * rho & + + 3.787878787878788d-03) * rho & + - 1.851851851851852d-02) * rho & + + 7.142857142857142d-02) * rho & + - 2.000000000000000d-01) * rho & + + 3.333333333333333d-01 + case(2) + crint_1 = (((((((((1.102292768959436d-08 * rho & + - 1.198144314086343d-07) * rho & + + 1.181027966742252d-06) * rho & + - 1.044277360066834d-05) * rho & + + 8.169934640522875d-05) * rho & + - 5.555555555555556d-04) * rho & + + 3.205128205128205d-03) * rho & + - 1.515151515151515d-02) * rho & + + 5.555555555555555d-02) * rho & + - 1.428571428571428d-01) * rho & + + 2.000000000000000d-01 + case(3) + crint_1 = (((((((((1.020641452740218d-08 * rho & + - 1.102292768959436d-07) * rho & + + 1.078329882677709d-06) * rho & + - 9.448223733938020d-06) * rho & + + 7.309941520467836d-05) * rho & + - 4.901960784313725d-04) * rho & + + 2.777777777777778d-03) * rho & + - 1.282051282051282d-02) * rho & + + 4.545454545454546d-02) * rho & + - 1.111111111111111d-01) * rho & + + 1.428571428571428d-01 + case default + tmp = dble(n + n + 1) + crint_1 = (((((((((2.755731922398589d-07 * rho / (tmp + 20.d0) & + - 2.755731922398589d-06 / (tmp + 18.d0)) * rho & + + 2.480158730158730d-05 / (tmp + 16.d0)) * rho & + - 1.984126984126984d-04 / (tmp + 14.d0)) * rho & + + 1.388888888888889d-03 / (tmp + 12.d0)) * rho & + - 8.333333333333333d-03 / (tmp + 10.d0)) * rho & + + 4.166666666666666d-02 / (tmp + 8.d0)) * rho & + - 1.666666666666667d-01 / (tmp + 6.d0)) * rho & + + 5.000000000000000d-01 / (tmp + 4.d0)) * rho & + - 1.000000000000000d+00 / (tmp + 2.d0)) * rho & + + 1.0d0 / tmp + end select - if(rho_mod < 10.d0) then - ! small z - if(rho_mod .lt. 1.d-15) then - crint_1 = 1.d0 / dble(n + n + 1) else + crint_1 = crint_smallz(n, rho) + endif else - ! large z if(rho_mod .gt. 40.d0) then - n_tmp = dble(n) + 0.5d0 - crint_1 = 0.5d0 * gamma(n_tmp) / (rho**n_tmp) + crint_1 = 0.5d0 * gamma(dble(n) + 0.5d0) / (rho**n * zsqrt(rho)) else - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - rho_exp = 0.5d0 * zexp(-rho) rho_inv = (1.d0, 0.d0) / rho - crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + call zboysfun00_1(rho, crint_1) mmax = n if(mmax .gt. 0) then do i = 0, mmax-1 @@ -56,9 +112,9 @@ complex*16 function crint_1(n, rho) endif endif - endif + return end ! --- @@ -72,27 +128,15 @@ complex*16 function crint_sum_1(n_pt_out, rho, d1) complex*16, intent(in) :: rho, d1(0:n_pt_out) integer :: n, i, mmax - double precision :: rho_mod, rho_re, rho_im - double precision :: sq_rho_re, sq_rho_im + double precision :: rho_mod complex*16 :: sq_rho, F0 complex*16 :: rho_tmp, rho_inv, rho_exp complex*16, allocatable :: Fm(:) - complex*16 :: crint_smallz, cpx_erf_1 + complex*16 :: crint_smallz - rho_re = real (rho) - rho_im = aimag(rho) - rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) - -! ! debug -! double precision :: d1_real(0:n_pt_out) -! double precision :: rint_sum -! do i = 0, n_pt_out -! d1_real(i) = real(d1(i)) -! enddo -! crint_sum_1 = rint_sum(n_pt_out, rho_re, d1_real) -! return + rho_mod = zabs(rho) if(rho_mod < 10.d0) then ! small z @@ -132,12 +176,7 @@ complex*16 function crint_sum_1(n_pt_out, rho, d1) else - ! get \sqrt(rho) - sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) - sq_rho_im = 0.5d0 * rho_im / sq_rho_re - sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im - - F0 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + call zboysfun00_1(rho, F0) crint_sum_1 = F0 * d1(0) rho_exp = 0.5d0 * zexp(-rho) @@ -219,15 +258,16 @@ complex*16 function crint_2(n, rho) integer, intent(in) :: n complex*16, intent(in) :: rho - double precision :: tmp - complex*16 :: rho2 + double precision :: tmp, abs_rho complex*16 :: vals(0:n) complex*16, external :: crint_smallz - if(abs(rho) < 3.5d0) then + abs_rho = zabs(rho) - if(abs(rho) .lt. 0.35d0) then + if(abs_rho < 3.5d0) then + + if(abs_rho .lt. 0.35d0) then select case(n) case(0) @@ -343,7 +383,7 @@ subroutine zboysfun(n_max, x, vals) integer :: n complex*16 :: yy, x_inv - call zboysfun00(x, vals(0)) + call zboysfun00_2(x, vals(0)) yy = 0.5d0 * zexp(-x) x_inv = (1.d0, 0.d0) / x @@ -434,7 +474,7 @@ subroutine crint_2_vec(n_max, rho, vals) integer :: n double precision :: tmp, abs_rho - complex*16 :: rho2, rho3, erho + complex*16 :: erho abs_rho = abs(rho) @@ -855,8 +895,6 @@ subroutine crint_quad_12_vec(n_max, rho, vals) complex*16, intent(out) :: vals(0:n_max) integer :: n - double precision :: tmp, abs_rho - complex*16 :: rho2, rho3, erho do n = 0, n_max call crint_quad_12(n, rho, 10000000, vals(n)) diff --git a/src/utils/cpx_erf.irp.f b/src/utils/cpx_erf.irp.f index 914a5157..6afac49e 100644 --- a/src/utils/cpx_erf.irp.f +++ b/src/utils/cpx_erf.irp.f @@ -1,6 +1,37 @@ ! --- +subroutine zboysfun00_1(rho, F0) + + implicit none + + include 'constants.include.F' + + complex*16, intent(in) :: rho + complex*16, intent(out) :: F0 + + double precision :: rho_re, rho_im, rho_mod + double precision :: sq_rho_re, sq_rho_im + complex*16 :: sq_rho + + complex*16, external :: cpx_erf_1 + + + rho_re = real(rho) + rho_im = aimag(rho) + rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) + + sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod) + sq_rho_im = 0.5d0 * rho_im / sq_rho_re + sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im + + F0 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho + + return +end + +! --- + complex*16 function cpx_erf_1(x, y) BEGIN_DOC @@ -43,6 +74,7 @@ complex*16 function cpx_erf_1(x, y) cpx_erf_1 = conjg(erf_tot) endif + return end ! --- @@ -54,22 +86,42 @@ complex*16 function erf_E(x, yabs) double precision, intent(in) :: x, yabs - if((dabs(x).gt.6.d0) .or. (x==0.d0)) then + double precision :: xy + + + if(dabs(x) .gt. 6.d0) then erf_E = (0.d0, 0.d0) return endif - if(dabs(x) .lt. 1.d-7) then + xy = x * yabs - erf_E = -inv_pi * (0.d0, 1.d0) * yabs + if(dabs(xy) .lt. 0.1d0) then + + erf_E = ((((((((((((-(0.d0, 9.3968347936601903d-8) * xy + (6.5777843555621328d-7, 0.d0) & + ) * xy + (0.d0, 4.2755598311153869d-6) & + ) * xy - (2.5653358986692322d-5, 0.d0) & + ) * xy - (0.d0, 0.00014109347442680775d0) & + ) * xy + (0.00070546737213403885d0, 0.d0) & + ) * xy + (0.d0, 0.0031746031746031746d0) & + ) * xy - (0.012698412698412698d0, 0.d0) & + ) * xy - (0.d0, 0.044444444444444446d0) & + ) * xy + (0.13333333333333333d0, 0.d0) & + ) * xy + (0.d0, 0.33333333333333331d0) & + ) * xy - (0.66666666666666663d0, 0.d0) & + ) * xy - (0.d0, 1.d0) & + ) * xy + (1.d0, 0.d0) + + erf_E = erf_E * (0.d0, 1.d0) * yabs * inv_pi * dexp(-x*x) else erf_E = 0.5d0 * inv_pi * dexp(-x*x) & - * ((1.d0, 0.d0) - zexp(-(2.d0, 0.d0) * x * yabs)) / x + * ((1.d0, 0.d0) - zexp(-(0.d0, 2.d0) * xy)) / x endif + return end ! --- @@ -84,31 +136,30 @@ double precision function erf_F(x, yabs) integer, parameter :: Nmax = 13 integer :: i - double precision :: tmp1, tmp2, x2, ct + double precision :: tmp1, tmp2, x2 if(dabs(x) .gt. 5.8d0) then - erf_F = 0.d0 + return + endif - else - x2 = x * x - ct = x * inv_pi + x2 = x * x - erf_F = 0.d0 - do i = 1, Nmax + erf_F = 0.d0 + do i = 1, Nmax - tmp1 = 0.25d0 * dble(i) * dble(i) + x2 - tmp2 = dexp(-tmp1) / tmp1 - erf_F = erf_F + tmp2 + tmp1 = 0.25d0 * dble(i*i) + x2 + tmp2 = dexp(-tmp1) / tmp1 + erf_F = erf_F + tmp2 - if(dabs(tmp2) .lt. 1d-15) exit - enddo - erf_F = ct * erf_F + if(tmp2 .lt. 1d-15) exit + enddo - endif + erf_F = x * inv_pi * erf_F + return end ! --- @@ -149,6 +200,7 @@ complex*16 function erf_G(x, yabs) enddo + return end ! --- @@ -163,19 +215,14 @@ complex*16 function erf_H(x, yabs) integer, parameter :: Nmax = 13 integer :: i - double precision :: tmp0, tmp1, tmp_mod, x2, ct, idble + double precision :: tmp0, tmp1, x2, idble complex*16 :: tmp2 - if(x .eq. 0.d0) then - erf_H = (0.d0, 0.d0) - return - endif - - if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then + if((dabs(x) .lt. 107d0) .and. (yabs .lt. 6.1d0)) then + !if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then x2 = x * x - ct = 0.5d0 * inv_pi erf_H = 0.d0 do i = 1, Nmax @@ -186,10 +233,10 @@ complex*16 function erf_H(x, yabs) tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1 erf_H = erf_H + tmp2 - tmp_mod = dsqrt(real(tmp2)*real(tmp2) + aimag(tmp2)*aimag(tmp2)) - if(tmp_mod .lt. 1d-15) exit + if(zabs(tmp2) .lt. 1d-15) exit enddo - erf_H = ct * erf_H + + erf_H = 0.5d0 * inv_pi * erf_H else @@ -201,7 +248,7 @@ end ! --- -subroutine zboysfun00(z, val) +subroutine zboysfun00_2(z, val) BEGIN_DOC ! From ff63afa7978aa2ee39a975c927ea9d66d1806822 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Thu, 17 Oct 2024 19:14:22 +0200 Subject: [PATCH 2/5] working on maney center problem in cgtos --- src/ao_basis/EZFIO.cfg | 6 +- src/ao_basis/README.rst | 2 +- src/ao_one_e_ints/aos_cgtos.irp.f | 79 +-- .../one_e_coul_integrals_cgtos.irp.f | 70 +- .../one_e_kin_integrals_cgtos.irp.f | 80 +-- .../two_e_coul_integrals_cgtos.irp.f | 605 +++++++++--------- src/hartree_fock/deb_ao_2e_int.irp.f | 195 +++++- src/utils/cgtos_one_e.irp.f | 44 +- src/utils/cgtos_utils.irp.f | 72 ++- src/utils/cpx_boys.irp.f | 188 +++--- 10 files changed, 773 insertions(+), 568 deletions(-) diff --git a/src/ao_basis/EZFIO.cfg b/src/ao_basis/EZFIO.cfg index ad36ca6b..bd716383 100644 --- a/src/ao_basis/EZFIO.cfg +++ b/src/ao_basis/EZFIO.cfg @@ -73,7 +73,7 @@ doc: If true, use cgtos for AO integrals interface: ezfio default: False -[ao_expo_im_cgtos] +[ao_expo_im] type: double precision doc: imag part for Exponents for each primitive of each cGTOs |AO| size: (ao_basis.ao_num,ao_basis.ao_prim_num_max) @@ -82,12 +82,12 @@ interface: ezfio, provider [ao_expo_pw] type: double precision doc: plane wave part for each primitive GTOs |AO| -size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max) +size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max) interface: ezfio, provider [ao_expo_phase] type: double precision doc: phase shift for each primitive GTOs |AO| -size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max) +size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max) interface: ezfio, provider diff --git a/src/ao_basis/README.rst b/src/ao_basis/README.rst index 6b9e6c07..0b7e5c3e 100644 --- a/src/ao_basis/README.rst +++ b/src/ao_basis/README.rst @@ -40,6 +40,6 @@ Complex Gaussian-Type Orbitals (cGTOs) are also supported: \chi_i(\mathbf{r}) = x^a y^b z^c \sum_k c_{ki} \left( e^{-\alpha_{ki} \mathbf{r}^2 - \imath \mathbf{k}_{ki} \cdot \mathbf{r} - \imath \phi_{ki}} + \text{C.C.} \right) where: - - :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im_cgtos``), + - :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im``), - :math:`\mathbf{k} = (k_x, k_y, k_z) \in \mathbb{R}^3` (specified by ``ao_expo_pw``), - :math:`\phi = \phi_x + \phi_y + \phi_z \in \mathbb{R}` (specified by ``ao_expo_phase``). diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index 2792d938..b5dd5263 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -21,21 +21,18 @@ END_PROVIDER &BEGIN_PROVIDER [complex*16, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)] implicit none + integer :: i, j, m do j = 1, ao_num do i = 1, ao_prim_num_max + ao_expo_cgtos_ord_transp(i,j) = ao_expo_cgtos_ord(j,i) - do m = 1, 3 + + do m = 1, 4 ao_expo_pw_ord_transp(m,i,j) = ao_expo_pw_ord(m,j,i) ao_expo_phase_ord_transp(m,i,j) = ao_expo_phase_ord(m,j,i) enddo - ao_expo_pw_ord_transp(4,i,j) = ao_expo_pw_ord_transp(1,i,j) * ao_expo_pw_ord_transp(1,i,j) & - + ao_expo_pw_ord_transp(2,i,j) * ao_expo_pw_ord_transp(2,i,j) & - + ao_expo_pw_ord_transp(3,i,j) * ao_expo_pw_ord_transp(3,i,j) - ao_expo_phase_ord_transp(4,i,j) = ao_expo_phase_ord_transp(1,j,i) & - + ao_expo_phase_ord_transp(2,j,i) & - + ao_expo_phase_ord_transp(3,j,i) enddo enddo @@ -50,16 +47,12 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] integer :: i, j, ii, m, powA(3), nz double precision :: norm double precision :: kA2, phiA - complex*16 :: expo, expo_inv, C_A(3) + complex*16 :: expo, expo_inv, C_Ae(3), C_Ap(3) complex*16 :: overlap_x, overlap_y, overlap_z complex*16 :: integ1, integ2, C1, C2 nz = 100 - C_A(1) = (0.d0, 0.d0) - C_A(2) = (0.d0, 0.d0) - C_A(3) = (0.d0, 0.d0) - ao_coef_norm_cgtos = 0.d0 do i = 1, ao_num @@ -74,19 +67,22 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] do j = 1, ao_prim_num(i) - expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im_cgtos(i,j) + expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im(i,j) expo_inv = (1.d0, 0.d0) / expo do m = 1, 3 - C_A(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j) + C_Ap(m) = nucl_coord(ii,m) + C_Ae(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j) enddo - phiA = ao_expo_phase(4,i,j) - KA2 = ao_expo_pw(4,i,j) + phiA = ao_expo_phase(1,i,j) + ao_expo_phase(2,i,j) + ao_expo_phase(3,i,j) + KA2 = ao_expo_pw(1,i,j) * ao_expo_pw(1,i,j) & + + ao_expo_pw(2,i,j) * ao_expo_pw(2,i,j) & + + ao_expo_pw(3,i,j) * ao_expo_pw(3,i,j) C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2) C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2) - call overlap_cgaussian_xyz(C_A, C_A, expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz) - call overlap_cgaussian_xyz(conjg(C_A), C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz) + call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz) + call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz) norm = 2.d0 * real(C1 * integ1 + C2 * integ2) @@ -126,12 +122,17 @@ END_PROVIDER iorder(j) = j d(j,1) = ao_expo(i,j) d(j,2) = ao_coef_norm_cgtos(i,j) - d(j,3) = ao_expo_im_cgtos(i,j) + d(j,3) = ao_expo_im(i,j) - do m = 1, 4 + do m = 1, 3 d(j,3+m) = ao_expo_pw(m,i,j) + enddo + d(j,7) = d(j,4) * d(j,4) + d(j,5) * d(j,5) + d(j,6) * d(j,6) + + do m = 1, 3 d(j,7+m) = ao_expo_phase(m,i,j) enddo + d(j,11) = d(j,8) + d(j,9) + d(j,10) enddo call dsort(d(1,1), iorder, ao_prim_num(i)) @@ -165,8 +166,8 @@ END_PROVIDER double precision :: c, overlap, overlap_x, overlap_y, overlap_z double precision :: KA2(3), phiA(3) double precision :: KB2(3), phiB(3) - complex*16 :: alpha, alpha_inv, A_center(3) - complex*16 :: beta, beta_inv, B_center(3) + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) complex*16 :: C1(1:4), C2(1:4) complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1 complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2 @@ -178,18 +179,18 @@ END_PROVIDER dim1 = 100 - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & - !$OMP alpha, alpha_inv, A_center, power_A, KA2, phiA, & - !$OMP beta, beta_inv, B_center, power_B, KB2, phiB, & - !$OMP overlap_x , overlap_y , overlap_z , overlap, & - !$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, & - !$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) & - !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & - !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & - !$OMP ao_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z, & + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & + !$OMP alpha, alpha_inv, Ae_center, Ap_center, power_A, KA2, phiA, & + !$OMP beta, beta_inv, Be_center, Bp_center, power_B, KB2, phiB, & + !$OMP overlap_x , overlap_y , overlap_z , overlap, & + !$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, & + !$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) & + !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & + !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & + !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & + !$OMP ao_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z, & !$OMP ao_overlap_cgtos) do j = 1, ao_num @@ -212,7 +213,8 @@ END_PROVIDER alpha_inv = (1.d0, 0.d0) / alpha do m = 1, 3 phiA(m) = ao_expo_phase_ord_transp(m,n,j) - A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) + Ap_center(m) = nucl_coord(jj,m) + Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j) enddo @@ -222,7 +224,8 @@ END_PROVIDER beta_inv = (1.d0, 0.d0) / beta do m = 1, 3 phiB(m) = ao_expo_phase_ord_transp(m,l,i) - B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) + Bp_center(m) = nucl_coord(ii,m) + Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i) enddo @@ -238,10 +241,10 @@ END_PROVIDER C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3))) C2(4) = C2(1) * C2(2) * C2(3) - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) - call overlap_cgaussian_xyz(conjg(A_center), B_center, conjg(alpha), beta, power_A, power_B, & + call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, conjg(Ap_center), Bp_center, & overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2) diff --git a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f index 568d4d8f..fed5f29d 100644 --- a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f @@ -17,25 +17,25 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] double precision :: c, Z, C_center(3) double precision :: phiA, KA2 double precision :: phiB, KB2 - complex*16 :: alpha, alpha_inv, A_center(3) - complex*16 :: beta, beta_inv, B_center(3) + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) complex*16 :: C1, C2, I1, I2 complex*16 :: NAI_pol_mult_cgtos ao_integrals_n_e_cgtos = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, & - !$OMP alpha, alpha_inv, A_center, phiA, KA2, power_A, C1, I1, & - !$OMP beta, beta_inv, B_center, phiB, KB2, power_B, C2, I2) & - !$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, & - !$OMP ao_power, nucl_num, nucl_charge, n_pt_max_integrals, & - !$OMP ao_expo_cgtos_ord_transp, ao_coef_cgtos_norm_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & - !$OMP ao_integrals_n_e_cgtos) - !$OMP DO SCHEDULE (dynamic) + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, C1, C2, I1, I2, & + !$OMP alpha, alpha_inv, Ae_center, Ap_center, phiA, KA2, power_A, & + !$OMP beta, beta_inv, Be_center, Bp_center, phiB, KB2, power_B) & + !$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, & + !$OMP ao_power, nucl_num, nucl_charge, n_pt_max_integrals, & + !$OMP ao_expo_cgtos_ord_transp, ao_coef_cgtos_norm_ord_transp, & + !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & + !$OMP ao_integrals_n_e_cgtos) + !$OMP DO SCHEDULE (dynamic) do j = 1, ao_num @@ -51,9 +51,9 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] alpha = ao_expo_cgtos_ord_transp(n,j) alpha_inv = (1.d0, 0.d0) / alpha - do m = 1, 3 - A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) + Ap_center(m) = nucl_coord(jj,m) + Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) enddo phiA = ao_expo_phase_ord_transp(4,n,j) KA2 = ao_expo_pw_ord_transp(4,n,j) @@ -62,9 +62,9 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] beta = ao_expo_cgtos_ord_transp(l,i) beta_inv = (1.d0, 0.d0) / beta - do m = 1, 3 - B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) + Bp_center(m) = nucl_coord(ii,m) + Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) enddo phiB = ao_expo_phase_ord_transp(4,l,i) KB2 = ao_expo_pw_ord_transp(4,l,i) @@ -79,9 +79,11 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] C_center(1:3) = nucl_coord(k,1:3) - I1 = NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_max_integrals) + I1 = NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, & + Ap_center, Bp_center, C_center, n_pt_max_integrals) - I2 = NAI_pol_mult_cgtos(conjg(A_center), B_center, power_A, power_B, conjg(alpha), beta, C_center, n_pt_max_integrals) + I2 = NAI_pol_mult_cgtos(conjg(Ae_center), Be_center, power_A, power_B, conjg(alpha), beta, & + conjg(Ap_center), Bp_center, C_center, n_pt_max_integrals) c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2) enddo @@ -93,14 +95,15 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)] enddo enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER ! --- -complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in) +complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, & + Ap_center, Bp_center, C_center, n_pt_in) BEGIN_DOC ! @@ -115,7 +118,8 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp integer, intent(in) :: n_pt_in, power_A(3), power_B(3) double precision, intent(in) :: C_center(3) - complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3) + complex*16, intent(in) :: alpha, Ae_center(3), Ap_center(3) + complex*16, intent(in) :: beta, Be_center(3), Bp_center(3) integer :: i, n_pt, n_pt_out double precision :: dist_AB, dist_AC @@ -124,20 +128,19 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp complex*16 :: d(0:n_pt_in) complex*16, external :: V_n_e_cgtos - complex*16, external :: crint_2 - complex*16, external :: crint_sum_2 + complex*16, external :: crint_sum dist_AB = 0.d0 dist_AC = 0.d0 do i = 1, 3 - dist_AB += abs(A_center(i) - B_center(i)) - dist_AC += abs(A_center(i) - C_center(i) * (1.d0, 0.d0)) + dist_AB += abs(Ae_center(i) - Be_center(i)) + dist_AC += abs(Ae_center(i) - C_center(i) * (1.d0, 0.d0)) enddo - if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13)) then + if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13) .or. use_pw) then continue @@ -157,8 +160,8 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp dist = (0.d0, 0.d0) dist_integral = (0.d0, 0.d0) do i = 1, 3 - P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv - dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i)) + P_center(i) = (alpha * Ae_center(i) + beta * Be_center(i)) * p_inv + dist += (Ae_center(i) - Be_center(i)) * (Ae_center(i) - Be_center(i)) dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) enddo @@ -175,12 +178,13 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp n_pt = 2 * ((power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3))) if(n_pt == 0) then - NAI_pol_mult_cgtos = coeff * crint_2(0, const) + !NAI_pol_mult_cgtos = coeff * crint_1(0, const) + NAI_pol_mult_cgtos = coeff * crint_sum(0, const, (1.d0, 0.d0)) return endif d(0:n_pt_in) = (0.d0, 0.d0) - call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & + call give_cpolynomial_mult_center_one_e(Ap_center, Bp_center, alpha, beta, & power_A, power_B, C_center, n_pt_in, d, n_pt_out) if(n_pt_out < 0) then @@ -188,7 +192,7 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp return endif - NAI_pol_mult_cgtos = coeff * crint_sum_2(n_pt_out, const, d) + NAI_pol_mult_cgtos = coeff * crint_sum(n_pt_out, const, d) return end diff --git a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f index 7f1f62cb..b55c37f8 100644 --- a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f @@ -10,8 +10,8 @@ double precision :: c, deriv_tmp double precision :: KA2, phiA double precision :: KB2, phiB - complex*16 :: alpha, alpha_inv, A_center(3), C1 - complex*16 :: beta, beta_inv, B_center(3), C2 + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3), C1 + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3), C2 complex*16 :: overlap_x, overlap_y, overlap_z, overlap complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1 complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2 @@ -24,30 +24,32 @@ ! -- Dummy call to provide everything - A_center(:) = (0.0d0, 0.d0) - B_center(:) = (1.0d0, 0.d0) - alpha = (1.0d0, 0.d0) - beta = (0.1d0, 0.d0) - power_A = 1 - power_B = 0 - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + Ae_center(:) = (0.0d0, 0.d0) + Be_center(:) = (1.0d0, 0.d0) + Ap_center(:) = (0.0d0, 0.d0) + Bp_center(:) = (1.0d0, 0.d0) + alpha = (1.0d0, 0.d0) + beta = (0.1d0, 0.d0) + power_A = 1 + power_B = 0 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) ! --- - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & - !$OMP A_center, power_A, alpha, alpha_inv, KA2, phiA, & - !$OMP B_center, power_B, beta, beta_inv, KB2, phiB, & - !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & - !$OMP overlap_x, overlap_y, overlap_z, overlap, & - !$OMP overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2, & - !$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, & - !$OMP overlap_y0_2, overlap_z0_2) & - !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & - !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & + !$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, & + !$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, & + !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & + !$OMP overlap_x, overlap_y, overlap_z, overlap, & + !$OMP overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2, & + !$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, & + !$OMP overlap_y0_2, overlap_z0_2) & + !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & + !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & + !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & !$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) do j = 1, ao_num @@ -73,7 +75,8 @@ alpha = ao_expo_cgtos_ord_transp(n,j) alpha_inv = (1.d0, 0.d0) / alpha do m = 1, 3 - A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) + Ap_center(m) = nucl_coord(jj,m) + Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_ord_transp(m,n,j) enddo phiA = ao_expo_phase_ord_transp(4,n,j) KA2 = ao_expo_pw_ord_transp(4,n,j) @@ -83,7 +86,8 @@ beta = ao_expo_cgtos_ord_transp(l,i) beta_inv = (1.d0, 0.d0) / beta do m = 1, 3 - B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) + Bp_center(m) = nucl_coord(ii,m) + Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_ord_transp(m,l,i) enddo phiB = ao_expo_phase_ord_transp(4,l,i) KB2 = ao_expo_pw_ord_transp(4,l,i) @@ -93,20 +97,20 @@ C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2)) - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) ! --- power_A(1) = power_A(1) - 2 if(power_A(1) > -1) then - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_m2_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_m2_2, overlap_y, overlap_z, overlap, dim1) else overlap_m2_1 = (0.d0, 0.d0) @@ -114,10 +118,10 @@ endif power_A(1) = power_A(1) + 4 - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_p2_1, overlap_y, overlap_z, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_p2_2, overlap_y, overlap_z, overlap, dim1) power_A(1) = power_A(1) - 2 @@ -138,10 +142,10 @@ power_A(2) = power_A(2) - 2 if(power_A(2) > -1) then - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x, overlap_m2_1, overlap_y, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_x, overlap_m2_2, overlap_y, overlap, dim1) else overlap_m2_1 = (0.d0, 0.d0) @@ -149,10 +153,10 @@ endif power_A(2) = power_A(2) + 4 - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x, overlap_p2_1, overlap_y, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_x, overlap_p2_2, overlap_y, overlap, dim1) power_A(2) = power_A(2) - 2 @@ -173,10 +177,10 @@ power_A(3) = power_A(3) - 2 if(power_A(3) > -1) then - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x, overlap_y, overlap_m2_1, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_x, overlap_y, overlap_m2_2, overlap, dim1) else overlap_m2_1 = (0.d0, 0.d0) @@ -184,10 +188,10 @@ endif power_A(3) = power_A(3) + 4 - call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & overlap_x, overlap_y, overlap_p2_1, overlap, dim1) - call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & overlap_x, overlap_y, overlap_p2_2, overlap, dim1) power_A(3) = power_A(3) - 2 diff --git a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f index 3ac6f1a1..9414e851 100644 --- a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f +++ b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f @@ -21,10 +21,10 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) double precision :: KJ2, phiJ double precision :: KK2, phiK double precision :: KL2, phiL - complex*16 :: expo1, expo1_inv, I_center(3) - complex*16 :: expo2, expo2_inv, J_center(3) - complex*16 :: expo3, expo3_inv, K_center(3) - complex*16 :: expo4, expo4_inv, L_center(3) + complex*16 :: expo1, expo1_inv, Ie_center(3), Ip_center(3) + complex*16 :: expo2, expo2_inv, Je_center(3), Jp_center(3) + complex*16 :: expo3, expo3_inv, Ke_center(3), Kp_center(3) + complex*16 :: expo4, expo4_inv, Le_center(3), Lp_center(3) complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv @@ -44,258 +44,238 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l) ao_two_e_integral_cgtos = ao_2e_cgtos_schwartz_accel(i, j, k, l) + return + endif + + dim1 = n_pt_max_integrals + + ii = ao_nucl(i) + jj = ao_nucl(j) + kk = ao_nucl(k) + ll = ao_nucl(l) + + do m = 1, 3 + I_power(m) = ao_power(i,m) + J_power(m) = ao_power(j,m) + K_power(m) = ao_power(k,m) + L_power(m) = ao_power(l,m) + enddo + + + ao_two_e_integral_cgtos = 0.d0 + + if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then + + do p = 1, ao_prim_num(i) + + coef1 = ao_coef_cgtos_norm_ord_transp(p,i) + expo1 = ao_expo_cgtos_ord_transp(p,i) + expo1_inv = (1.d0, 0.d0) / expo1 + do m = 1, 3 + Ip_center(m) = nucl_coord(ii,m) + Ie_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) + enddo + phiI = ao_expo_phase_ord_transp(4,p,i) + KI2 = ao_expo_pw_ord_transp(4,p,i) + + do q = 1, ao_prim_num(j) + + coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) + expo2 = ao_expo_cgtos_ord_transp(q,j) + expo2_inv = (1.d0, 0.d0) / expo2 + do m = 1, 3 + Jp_center(m) = nucl_coord(jj,m) + Je_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) + enddo + phiJ = ao_expo_phase_ord_transp(4,q,j) + KJ2 = ao_expo_pw_ord_transp(4,q,j) + + call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & + expo1, expo2, I_power, J_power, Ie_center, Je_center, Ip_center, Jp_center, dim1) + + p1_inv = (1.d0, 0.d0) / pp1 + + call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & + conjg(expo1), expo2, I_power, J_power, conjg(Ie_center), Je_center, conjg(Ip_center), Jp_center, dim1) + + p2_inv = (1.d0, 0.d0) / pp2 + + do r = 1, ao_prim_num(k) + + coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) + expo3 = ao_expo_cgtos_ord_transp(r,k) + expo3_inv = (1.d0, 0.d0) / expo3 + do m = 1, 3 + Kp_center(m) = nucl_coord(kk,m) + Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) + enddo + phiK = ao_expo_phase_ord_transp(4,r,k) + KK2 = ao_expo_pw_ord_transp(4,r,k) + + do s = 1, ao_prim_num(l) + + coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) + expo4 = ao_expo_cgtos_ord_transp(s,l) + expo4_inv = (1.d0, 0.d0) / expo4 + do m = 1, 3 + Lp_center(m) = nucl_coord(ll,m) + Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) + enddo + phiL = ao_expo_phase_ord_transp(4,s,l) + KL2 = ao_expo_pw_ord_transp(4,s,l) + + call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & + expo3, expo4, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1) + + q1_inv = (1.d0, 0.d0) / qq1 + + call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & + conjg(expo3), expo4, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1) + + q2_inv = (1.d0, 0.d0) / qq2 + + C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & + - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) + C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & + - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) + C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & + - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) + C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & + - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) + C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) & + - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) + C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) & + - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) + C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) & + - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) + C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) & + - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) + + int1 = general_primitive_integral_cgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + + int2 = general_primitive_integral_cgtos(dim1, & + P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + + int3 = general_primitive_integral_cgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + + int4 = general_primitive_integral_cgtos(dim1, & + P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + + int5 = general_primitive_integral_cgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + + int6 = general_primitive_integral_cgtos(dim1, & + conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + + int7 = general_primitive_integral_cgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) + + int8 = general_primitive_integral_cgtos(dim1, & + conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & + Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) + + int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 + + ao_two_e_integral_cgtos = ao_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) + enddo ! s + enddo ! r + enddo ! q + enddo ! p + else - dim1 = n_pt_max_integrals + do p = 1, ao_prim_num(i) - ii = ao_nucl(i) - jj = ao_nucl(j) - kk = ao_nucl(k) - ll = ao_nucl(l) + coef1 = ao_coef_cgtos_norm_ord_transp(p,i) + expo1 = ao_expo_cgtos_ord_transp(p,i) + phiI = ao_expo_phase_ord_transp(4,p,i) - do m = 1, 3 - I_power(m) = ao_power(i,m) - J_power(m) = ao_power(j,m) - K_power(m) = ao_power(k,m) - L_power(m) = ao_power(l,m) - enddo + do q = 1, ao_prim_num(j) + coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) + expo2 = ao_expo_cgtos_ord_transp(q,j) + phiJ = ao_expo_phase_ord_transp(4,q,j) - ao_two_e_integral_cgtos = 0.d0 + do r = 1, ao_prim_num(k) - if(ii /= jj .or. kk /= ll .or. jj /= kk) then + coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) + expo3 = ao_expo_cgtos_ord_transp(r,k) + phiK = ao_expo_phase_ord_transp(4,r,k) - do p = 1, ao_prim_num(i) + do s = 1, ao_prim_num(l) - coef1 = ao_coef_cgtos_norm_ord_transp(p,i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) - enddo - phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) + coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) + expo4 = ao_expo_cgtos_ord_transp(s,l) + phiL = ao_expo_phase_ord_transp(4,s,l) - do q = 1, ao_prim_num(j) + C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL)) + C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL)) + C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL)) + C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL)) + C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL)) + C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL)) + C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL)) + C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL)) - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) - enddo - phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) + int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & - expo1, expo2, I_power, J_power, I_center, J_center, dim1) - p1_inv = (1.d0, 0.d0) / pp1 + int2 = ERI_cgtos(expo1, expo2, conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & - conjg(expo1), expo2, I_power, J_power, conjg(I_center), J_center, dim1) - p2_inv = (1.d0, 0.d0) / pp2 + int3 = ERI_cgtos(conjg(expo1), expo2, expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - do r = 1, ao_prim_num(k) + int4 = ERI_cgtos(conjg(expo1), expo2, conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - expo3_inv = (1.d0, 0.d0) / expo3 - do m = 1, 3 - K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) - enddo - phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) + int5 = ERI_cgtos(expo1, conjg(expo2), expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - do s = 1, ao_prim_num(l) + int6 = ERI_cgtos(expo1, conjg(expo2), conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - expo4_inv = (1.d0, 0.d0) / expo4 - do m = 1, 3 - L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) - enddo - phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) + int7 = ERI_cgtos(conjg(expo1), conjg(expo2), expo3, expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & - expo3, expo4, K_power, L_power, K_center, L_center, dim1) - q1_inv = (1.d0, 0.d0) / qq1 + int8 = ERI_cgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, & + I_power(1), J_power(1), K_power(1), L_power(1), & + I_power(2), J_power(2), K_power(2), L_power(2), & + I_power(3), J_power(3), K_power(3), L_power(3)) - call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & - conjg(expo3), expo4, K_power, L_power, conjg(K_center), L_center, dim1) - q2_inv = (1.d0, 0.d0) / qq2 + int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 & + + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) + ao_two_e_integral_cgtos = ao_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) + enddo ! s + enddo ! r + enddo ! q + enddo ! p - int1 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int2 = general_primitive_integral_cgtos(dim1, & - P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int3 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int4 = general_primitive_integral_cgtos(dim1, & - P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int5 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int6 = general_primitive_integral_cgtos(dim1, & - conjg(P2_new), conjg(P2_center), conjg(fact_p2), conjg(pp2), conjg(p2_inv), iorder_p2, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int7 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1) - - int8 = general_primitive_integral_cgtos(dim1, & - conjg(P1_new), conjg(P1_center), conjg(fact_p1), conjg(pp1), conjg(p1_inv), iorder_p1, & - Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2) - - int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - - ao_two_e_integral_cgtos = ao_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - else - - do p = 1, ao_prim_num(i) - - coef1 = ao_coef_cgtos_norm_ord_transp(p,i) - expo1 = ao_expo_cgtos_ord_transp(p,i) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) - enddo - phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) - - do q = 1, ao_prim_num(j) - - coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) - expo2 = ao_expo_cgtos_ord_transp(q,j) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) - enddo - phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) - - do r = 1, ao_prim_num(k) - - coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) - expo3 = ao_expo_cgtos_ord_transp(r,k) - expo3_inv = (1.d0, 0.d0) / expo3 - do m = 1, 3 - K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) - enddo - phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) - - do s = 1, ao_prim_num(l) - - coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) - expo4 = ao_expo_cgtos_ord_transp(s,l) - expo4_inv = (1.d0, 0.d0) / expo4 - do m = 1, 3 - L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) - enddo - phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) - - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - - int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int2 = ERI_cgtos(expo1, expo2, conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int3 = ERI_cgtos(conjg(expo1), expo2, expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int4 = ERI_cgtos(conjg(expo1), expo2, conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int5 = ERI_cgtos(expo1, conjg(expo2), expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int6 = ERI_cgtos(expo1, conjg(expo2), conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int7 = ERI_cgtos(conjg(expo1), conjg(expo2), expo3, expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int8 = ERI_cgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, & - I_power(1), J_power(1), K_power(1), L_power(1), & - I_power(2), J_power(2), K_power(2), L_power(2), & - I_power(3), J_power(3), K_power(3), L_power(3)) - - int_tot = C1 * int1 + C2 * int2 + C3 * int3 + C4 * int4 & - + C5 * int5 + C6 * int6 + C7 * int7 + C8 * int8 - - ao_two_e_integral_cgtos = ao_two_e_integral_cgtos + coef4 * 2.d0 * real(int_tot) - enddo ! s - enddo ! r - enddo ! q - enddo ! p - - endif ! same centers - endif ! do schwartz + endif end @@ -321,10 +301,10 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) double precision :: KJ2, phiJ double precision :: KK2, phiK double precision :: KL2, phiL - complex*16 :: expo1, expo1_inv, I_center(3) - complex*16 :: expo2, expo2_inv, J_center(3) - complex*16 :: expo3, expo3_inv, K_center(3) - complex*16 :: expo4, expo4_inv, L_center(3) + complex*16 :: expo1, expo1_inv, Ie_center(3), Ip_center(3) + complex*16 :: expo2, expo2_inv, Je_center(3), Jp_center(3) + complex*16 :: expo3, expo3_inv, Ke_center(3), Kp_center(3) + complex*16 :: expo4, expo4_inv, Le_center(3), Lp_center(3) complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv @@ -362,7 +342,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) - if(ii /= jj .or. kk /= ll .or. jj /= kk) then + if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then schwartz_kl(0,0) = 0.d0 do r = 1, ao_prim_num(k) @@ -371,7 +351,8 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) expo1 = ao_expo_cgtos_ord_transp(r,k) expo1_inv = (1.d0, 0.d0) / expo1 do m = 1, 3 - K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) + Kp_center(m) = nucl_coord(kk,m) + Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) KK2 = ao_expo_pw_ord_transp(4,r,k) @@ -383,17 +364,20 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) expo2 = ao_expo_cgtos_ord_transp(s,l) expo2_inv = (1.d0, 0.d0) / expo2 do m = 1, 3 - L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) + Lp_center(m) = nucl_coord(ll,m) + Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) KL2 = ao_expo_pw_ord_transp(4,s,l) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & - expo1, expo2, K_power, L_power, K_center, L_center, dim1) + expo1, expo2, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1) + p1_inv = (1.d0, 0.d0) / pp1 call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & - conjg(expo1), expo2, K_power, L_power, conjg(K_center), L_center, dim1) + conjg(expo1), expo2, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1) + p2_inv = (1.d0, 0.d0) / pp2 C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2)) @@ -457,7 +441,8 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) expo1 = ao_expo_cgtos_ord_transp(p,i) expo1_inv = (1.d0, 0.d0) / expo1 do m = 1, 3 - I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) + Ip_center(m) = nucl_coord(ii,m) + Ie_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) enddo phiI = ao_expo_phase_ord_transp(4,p,i) KI2 = ao_expo_pw_ord_transp(4,p,i) @@ -468,17 +453,19 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) expo2 = ao_expo_cgtos_ord_transp(q,j) expo2_inv = (1.d0, 0.d0) / expo2 do m = 1, 3 - J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) + Jp_center(m) = nucl_coord(jj,m) + Je_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) enddo phiJ = ao_expo_phase_ord_transp(4,q,j) KJ2 = ao_expo_pw_ord_transp(4,q,j) call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & - expo1, expo2, I_power, J_power, I_center, J_center, dim1) + expo1, expo2, I_power, J_power, Ie_center, Je_center, Ip_center, Jp_center, dim1) + p1_inv = (1.d0, 0.d0) / pp1 call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & - conjg(expo1), expo2, I_power, J_power, conjg(I_center), J_center, dim1) + conjg(expo1), expo2, I_power, J_power, conjg(Ie_center), Je_center, conjg(Ip_center), Jp_center, dim1) p2_inv = (1.d0, 0.d0) / pp2 C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2)) @@ -538,7 +525,8 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) expo3 = ao_expo_cgtos_ord_transp(r,k) expo3_inv = (1.d0, 0.d0) / expo3 do m = 1, 3 - K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) + Kp_center(m) = nucl_coord(kk,m) + Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) enddo phiK = ao_expo_phase_ord_transp(4,r,k) KK2 = ao_expo_pw_ord_transp(4,r,k) @@ -550,17 +538,20 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) expo4 = ao_expo_cgtos_ord_transp(s,l) expo4_inv = (1.d0, 0.d0) / expo4 do m = 1, 3 - L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) + Lp_center(m) = nucl_coord(ll,m) + Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) enddo phiL = ao_expo_phase_ord_transp(4,s,l) KL2 = ao_expo_pw_ord_transp(4,s,l) call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & - expo3, expo4, K_power, L_power, K_center, L_center, dim1) + expo3, expo4, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1) + q1_inv = (1.d0, 0.d0) / qq1 call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & - conjg(expo3), expo4, K_power, L_power, conjg(K_center), L_center, dim1) + conjg(expo3), expo4, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1) + q2_inv = (1.d0, 0.d0) / qq2 C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & @@ -627,31 +618,21 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) coef1 = ao_coef_cgtos_norm_ord_transp(r,k) * ao_coef_cgtos_norm_ord_transp(r,k) expo1 = ao_expo_cgtos_ord_transp(r,k) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) - enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) schwartz_kl(0,r) = 0.d0 do s = 1, ao_prim_num(l) coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(s,l) * ao_coef_cgtos_norm_ord_transp(s,l) expo2 = ao_expo_cgtos_ord_transp(s,l) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) - enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) - C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2)) - C2 = zexp(-(0.d0, 2.d0) * phiL - 0.5d0 * (real(expo1_inv) * KK2 + expo2_inv * KL2)) + C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL)) + C2 = zexp(-(0.d0, 2.d0) * phiL) !C3 = C2 - C4 = zexp((0.d0, 2.d0) * (phiK - phiL) - 0.5d0 * (conjg(expo1_inv) * KK2 + expo2_inv * KL2)) - C5 = zexp(-(0.d0, 2.d0) * phiK - 0.5d0 * (expo1_inv * KK2 + real(expo2_inv) * KL2)) - C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) + C4 = zexp((0.d0, 2.d0) * (phiK - phiL)) + C5 = zexp(-(0.d0, 2.d0) * phiK) + C6 = (1.d0, 0.d0) !C7 = C6 !C8 = conjg(C5) @@ -711,30 +692,20 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) coef1 = ao_coef_cgtos_norm_ord_transp(p,i) expo1 = ao_expo_cgtos_ord_transp(p,i) - expo1_inv = (1.d0, 0.d0) / expo1 - do m = 1, 3 - I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) - enddo phiI = ao_expo_phase_ord_transp(4,p,i) - KI2 = ao_expo_pw_ord_transp(4,p,i) do q = 1, ao_prim_num(j) coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) expo2 = ao_expo_cgtos_ord_transp(q,j) - expo2_inv = (1.d0, 0.d0) / expo2 - do m = 1, 3 - J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) - enddo phiJ = ao_expo_phase_ord_transp(4,q,j) - KJ2 = ao_expo_pw_ord_transp(4,q,j) - C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2)) - C2 = zexp(-(0.d0, 2.d0) * phiJ - 0.5d0 * (real(expo1_inv) * KI2 + expo2_inv * KJ2)) + C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ)) + C2 = zexp(-(0.d0, 2.d0) * phiJ) !C3 = C2 - C4 = zexp((0.d0, 2.d0) * (phiI - phiJ) - 0.5d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2)) - C5 = zexp(-(0.d0, 2.d0) * phiI - 0.5d0 * (expo1_inv * KI2 + real(expo2_inv) * KJ2)) - C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) + C4 = zexp((0.d0, 2.d0) * (phiI - phiJ)) + C5 = zexp(-(0.d0, 2.d0) * phiI) + C6 = (1.d0, 0.d0) !C7 = C6 !C8 = conjg(C5) @@ -791,41 +762,23 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l) coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) expo3 = ao_expo_cgtos_ord_transp(r,k) - expo3_inv = (1.d0, 0.d0) / expo3 - do m = 1, 3 - K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) - enddo phiK = ao_expo_phase_ord_transp(4,r,k) - KK2 = ao_expo_pw_ord_transp(4,r,k) do s = 1, ao_prim_num(l) if(schwartz_kl(s,r)*schwartz_ij < thr) cycle coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) expo4 = ao_expo_cgtos_ord_transp(s,l) - expo4_inv = (1.d0, 0.d0) / expo4 - do m = 1, 3 - L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) - enddo phiL = ao_expo_phase_ord_transp(4,s,l) - KL2 = ao_expo_pw_ord_transp(4,s,l) - C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) & - - 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) - C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) - C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) & - - 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) + C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL)) + C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL)) + C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL)) + C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL)) + C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL)) + C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL)) + C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL)) + C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL)) int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & I_power(1), J_power(1), K_power(1), L_power(1), & @@ -937,7 +890,7 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ complex*16 :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim) complex*16 :: d1(0:max_dim), d_poly(0:max_dim) - complex*16 :: crint_sum_2 + complex*16 :: crint_sum !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol @@ -1076,7 +1029,7 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_ !DIR$ FORCEINLINE call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) - accu = crint_sum_2(n_pt_out, const, d1) + accu = crint_sum(n_pt_out, const, d1) general_primitive_integral_cgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq @@ -1108,25 +1061,25 @@ complex*16 function ERI_cgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_y, ERI_cgtos = (0.d0, 0.d0) - ASSERT (REAL(alpha) >= 0.d0) - ASSERT (REAL(beta ) >= 0.d0) - ASSERT (REAL(delta) >= 0.d0) - ASSERT (REAL(gama ) >= 0.d0) + ASSERT (real(alpha) >= 0.d0) + ASSERT (real(beta ) >= 0.d0) + ASSERT (real(delta) >= 0.d0) + ASSERT (real(gama ) >= 0.d0) nx = a_x + b_x + c_x + d_x - if(iand(nx,1) == 1) then + if(iand(nx, 1) == 1) then ERI_cgtos = (0.d0, 0.d0) return endif ny = a_y + b_y + c_y + d_y - if(iand(ny,1) == 1) then + if(iand(ny, 1) == 1) then ERI_cgtos = (0.d0, 0.d0) return endif nz = a_z + b_z + c_z + d_z - if(iand(nz,1) == 1) then + if(iand(nz, 1) == 1) then ERI_cgtos = (0.d0, 0.d0) return endif @@ -1727,3 +1680,21 @@ end ! --- +BEGIN_PROVIDER [logical, use_pw] + + implicit none + + logical :: exist + + use_pw = .false. + + call ezfio_has_ao_basis_ao_expo_pw(exist) + if(exist) then + PROVIDE ao_expo_pw + if(maxval(dabs(ao_expo_pw(4,:,:))) .gt. 1d-15) use_pw = .true. + endif + +END_PROVIDER + +! --- + diff --git a/src/hartree_fock/deb_ao_2e_int.irp.f b/src/hartree_fock/deb_ao_2e_int.irp.f index 94dc0bda..99b5b36a 100644 --- a/src/hartree_fock/deb_ao_2e_int.irp.f +++ b/src/hartree_fock/deb_ao_2e_int.irp.f @@ -7,7 +7,9 @@ program deb_ao_2e_int !call check_crint1() !call check_crint2() !call check_crint3() - call check_crint4() + !call check_crint4() + call check_crint5() + !call check_crint6() end @@ -354,7 +356,6 @@ subroutine check_crint4() double precision :: xr(1), x, shift double precision :: yr(1), y double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im - double precision :: delta_ref double precision :: t1, t2, t_int1, t_int2 complex*16 :: rho complex*16 :: int1, int2, int3 @@ -444,3 +445,193 @@ end ! --- +subroutine check_crint5() + + implicit none + + integer :: i_test, n_test + integer :: i, seed_size, clock_time + integer :: n + double precision :: xr(1), yr(1), nr(1), x, shift, y + double precision :: dif1_re, dif1_im, acc1_re, acc1_im + double precision :: dif2_re, dif2_im, acc2_re, acc2_im + double precision :: nrm_re, nrm_im + double precision :: t1, t2, t_int1, t_int2 + complex*16 :: rho + complex*16 :: int1, int2, int_ref + integer, allocatable :: seed(:) + + complex*16, external :: crint_1, crint_2 + + + + call random_seed(size=seed_size) + allocate(seed(seed_size)) + call system_clock(count=clock_time) + seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /) + !seed = 123456789 + call random_seed(put=seed) + + + t_int1 = 0.d0 + t_int2 = 0.d0 + + n_test = 100 + + acc1_re = 0.d0 + acc1_im = 0.d0 + acc2_re = 0.d0 + acc2_im = 0.d0 + nrm_re = 0.d0 + nrm_im = 0.d0 + do i_test = 1, n_test + + call random_number(xr) + call random_number(yr) + call random_number(nr) + + x = 1.d+1 * (30.d0 * xr(1) - 15.d0) + y = 1.d+1 * (30.d0 * yr(1) - 15.d0) + n = int(16.d0 * nr(1)) + + rho = x + (0.d0, 1.d0) * y + + call wall_time(t1) + int1 = crint_1(n, rho) + call wall_time(t2) + t_int1 = t_int1 + t2 - t1 + + call wall_time(t1) + int2 = crint_2(n, rho) + call wall_time(t2) + t_int2 = t_int2 + t2 - t1 + + call crint_quad_12(n, rho, 10000000, int_ref) + + dif1_re = dabs(real(int1) - real(int_ref)) + dif1_im = dabs(aimag(int1) - aimag(int_ref)) + + dif2_re = dabs(real(int2) - real(int_ref)) + dif2_im = dabs(aimag(int2) - aimag(int_ref)) + + if((dif2_re .gt. 1d-7) .or. (dif2_im .gt. 1d-7)) then + print*, ' important error found: ' + print*, " n, rho = ", n, x, y + print*, real(int1), real(int2), real(int_ref) + print*, aimag(int1), aimag(int2), aimag(int_ref) + !stop + endif + + acc1_re += dif1_re + acc1_im += dif1_im + + acc2_re += dif2_re + acc2_im += dif2_im + + nrm_re += dabs(real(int_ref)) + nrm_im += dabs(aimag(int_ref)) + enddo + + print*, "accuracy on boys_1 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc1_im / (nrm_im + 1d-15) + print*, "accuracy on boys_2 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc2_im / (nrm_im + 1d-15) + + print*, "boys_1 wall time (sec) = ", t_int1 + print*, "boys_2 wall time (sec) = ", t_int2 + + + deallocate(seed) + +end + +! --- + +subroutine check_crint6() + + implicit none + + integer :: i_test, n_test + integer :: i, seed_size, clock_time + integer :: n + double precision :: xr(1), yr(1), nr(1), x, shift, y + double precision :: dif_re, dif_im, acc_re, acc_im + double precision :: nrm_re, nrm_im + double precision :: t1, t2, t_int1, t_int2 + complex*16 :: rho + complex*16 :: int1, int2, int3 + integer, allocatable :: seed(:) + + complex*16, external :: crint_1, crint_2 + + + + call random_seed(size=seed_size) + allocate(seed(seed_size)) + call system_clock(count=clock_time) + seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /) + !seed = 123456789 + call random_seed(put=seed) + + + t_int1 = 0.d0 + t_int2 = 0.d0 + + n_test = 100 + + acc_re = 0.d0 + acc_im = 0.d0 + nrm_re = 0.d0 + nrm_im = 0.d0 + do i_test = 1, n_test + + call random_number(xr) + call random_number(yr) + call random_number(nr) + + x = 1.d0 * (30.d0 * xr(1) - 15.d0) + y = 1.d0 * (30.d0 * yr(1) - 15.d0) + n = int(16.d0 * nr(1)) + + rho = x + (0.d0, 1.d0) * y + + call wall_time(t1) + int1 = crint_1(n, rho) + call wall_time(t2) + t_int1 = t_int1 + t2 - t1 + + call wall_time(t1) + int2 = crint_2(n, rho) + call wall_time(t2) + t_int2 = t_int2 + t2 - t1 + + dif_re = dabs(real(int1) - real(int2)) + dif_im = dabs(aimag(int1) - aimag(int2)) + + if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then + print*, ' important error found: ' + print*, " n, rho = ", n, x, y + print*, real(int1), real(int2), dif_re + print*, aimag(int1), aimag(int2), dif_im + call crint_quad_12(n, rho, 100000000, int3) + print*, ' quad 100000000:', real(int3), aimag(int3) + !print*, ' quad 100000000:', dabs(real(int1) - real(int3)), dabs(aimag(int1) - aimag(int3)) + !stop + endif + + acc_re += dif_re + acc_im += dif_im + nrm_re += dabs(real(int1)) + nrm_im += dabs(aimag(int1)) + enddo + + print*, "diff (%):", 100.d0 * acc_re / (nrm_re + 1d-15), 100.d0 * acc_im / (nrm_im + 1d-15) + + print*, "boys_1 wall time (sec) = ", t_int1 + print*, "boys_2 wall time (sec) = ", t_int2 + + + deallocate(seed) + +end + +! --- + diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index 80b82eaf..a8773ce0 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -1,11 +1,12 @@ ! --- -complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A, power_B, dim) +complex*16 function overlap_cgaussian_x(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, dim) BEGIN_DOC ! - ! \int_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha (x-A_x)^2) exp(- beta(x-B_X)^2) dx + ! \int_{-infty}^{+infty} (x - Ap_x)^ax (x - Bp_x)^bx exp(-alpha (x - Ae_x)^2) exp(-beta (x - Be_X)^2) dx + ! ! with complex arguments ! END_DOC @@ -14,20 +15,19 @@ complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A include 'constants.include.F' integer, intent(in) :: dim, power_A, power_B - complex*16, intent(in) :: A_center, B_center, alpha, beta + complex*16, intent(in) :: Ae_center, alpha, Ap_center + complex*16, intent(in) :: Be_center, beta, Bp_center integer :: i, iorder_p - double precision :: fact_p_mod complex*16 :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p complex*16 :: Fc_integral - call give_explicit_cpoly_and_cgaussian_x( P_new, P_center, p, fact_p, iorder_p & - , alpha, beta, power_A, power_B, A_center, B_center, dim) + call give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_p, iorder_p, & + alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim) - fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p)) - if(fact_p_mod .lt. 1.d-14) then + if(zabs(fact_p) .lt. 1.d-14) then overlap_cgaussian_x = (0.d0, 0.d0) return endif @@ -37,22 +37,24 @@ complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A overlap_cgaussian_x = (0.d0, 0.d0) do i = 0, iorder_p - overlap_cgaussian_x += P_new(i) * Fc_integral(i, inv_sq_p) + overlap_cgaussian_x = overlap_cgaussian_x + P_new(i) * Fc_integral(i, inv_sq_p) enddo - overlap_cgaussian_x *= fact_p + overlap_cgaussian_x = overlap_cgaussian_x * fact_p end ! --- -subroutine overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & - overlap_x, overlap_y, overlap_z, overlap, dim ) +subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & + overlap_x, overlap_y, overlap_z, overlap, dim) BEGIN_DOC ! - ! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx + ! S_x = \int (x - Ap_x)^{a_x} exp(-\alpha (x - Ae_x)^2) (x - Bp_x)^{b_x} exp(-beta (x - Be_x)^2) dx + ! ! S = S_x S_y S_z + ! ! for complex arguments ! END_DOC @@ -61,20 +63,20 @@ subroutine overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power include 'constants.include.F' integer, intent(in) :: dim, power_A(3), power_B(3) - complex*16, intent(in) :: A_center(3), B_center(3), alpha, beta + complex*16, intent(in) :: Ae_center(3), alpha, Ap_center(3) + complex*16, intent(in) :: Be_center(3), beta, Bp_center(3) complex*16, intent(out) :: overlap_x, overlap_y, overlap_z, overlap integer :: i, nmax, iorder_p(3) - double precision :: fact_p_mod complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p complex*16 :: F_integral_tab(0:max_dim) complex*16 :: Fc_integral - call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim) + call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, & + alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim) - fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p)) - if(fact_p_mod .lt. 1.d-14) then + if(zabs(fact_p) .lt. 1.d-14) then overlap_x = (1.d-10, 0.d0) overlap_y = (1.d-10, 0.d0) overlap_z = (1.d-10, 0.d0) @@ -96,19 +98,19 @@ subroutine overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power do i = 1, iorder_p(1) overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) enddo - call cgaussian_product_x(alpha, A_center(1), beta, B_center(1), fact_p, p, P_center(1)) + call cgaussian_product_x(alpha, Ap_center(1), beta, Bp_center(1), fact_p, p, P_center(1)) overlap_x *= fact_p do i = 1, iorder_p(2) overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) enddo - call cgaussian_product_x(alpha, A_center(2), beta, B_center(2), fact_p, p, P_center(2)) + call cgaussian_product_x(alpha, Ap_center(2), beta, Bp_center(2), fact_p, p, P_center(2)) overlap_y *= fact_p do i = 1, iorder_p(3) overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) enddo - call cgaussian_product_x(alpha, A_center(3), beta, B_center(3), fact_p, p, P_center(3)) + call cgaussian_product_x(alpha, Ap_center(3), beta, Bp_center(3), fact_p, p, P_center(3)) overlap_z *= fact_p overlap = overlap_x * overlap_y * overlap_z diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index 76e96107..40c8a838 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -1,13 +1,19 @@ ! --- -subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim) +subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorder, & + alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim) BEGIN_DOC + ! ! Transform the product of - ! (x-x_A)^a (x-x_B)^b exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) + ! + ! (x - x_Ap)^a (x - x_Bp)^b exp(-alpha (r - Ae)^2) exp(-beta (r - Be)^2) + ! ! into - ! fact_k \sum_{i=0}^{iorder} (x-x_P)^i exp(-p(r-P)^2) + ! + ! fact_k \sum_{i=0}^{iorder} (x - x_P)^i exp(-p (r - P)^2) + ! END_DOC implicit none @@ -15,13 +21,13 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde integer, intent(in) :: dim integer, intent(in) :: a, b - complex*16, intent(in) :: alpha, beta, A_center, B_center + complex*16, intent(in) :: alpha, Ae_center, Ap_center + complex*16, intent(in) :: beta, Be_center, Bp_center integer, intent(out) :: iorder complex*16, intent(out) :: p, P_center, fact_k complex*16, intent(out) :: P_new(0:max_dim) integer :: n_new, i, j - double precision :: tmp_mod complex*16 :: P_a(0:max_dim), P_b(0:max_dim) complex*16 :: p_inv, ab, d_AB, tmp @@ -35,13 +41,12 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde ! new center p_inv = (1.d0, 0.d0) / p ab = alpha * beta - P_center = (alpha * A_center + beta * B_center) * p_inv + P_center = (alpha * Ae_center + beta * Be_center) * p_inv ! get the factor - d_AB = (A_center - B_center) * (A_center - B_center) - tmp = ab * p_inv * d_AB - tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) - if(tmp_mod .lt. 50.d0) then + d_AB = (Ae_center - Be_center) * (Ae_center - Be_center) + tmp = ab * p_inv * d_AB + if(zabs(tmp) .lt. 50.d0) then fact_k = zexp(-tmp) else fact_k = (0.d0, 0.d0) @@ -49,7 +54,7 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde ! Recenter the polynomials P_a and P_b on P_center !DIR$ FORCEINLINE - call recentered_cpoly2(P_a(0), A_center, P_center, a, P_b(0), B_center, P_center, b) + call recentered_cpoly2(P_a(0), Ap_center, P_center, a, P_b(0), Bp_center, P_center, b) n_new = 0 !DIR$ FORCEINLINE @@ -60,31 +65,38 @@ end ! --- -subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim) +subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, & + alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim) BEGIN_DOC + ! ! Transforms the product of - ! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) + ! + ! (x - x_Ap)^a(1) (x - x_Bp)^b(1) exp(-alpha (x - x_Ae)^2) exp(-beta (x - x_Be)^2) x + ! (y - y_Ap)^a(2) (y - y_Bp)^b(2) exp(-alpha (y - y_Ae)^2) exp(-beta (y - y_Be)^2) x + ! (z - z_Ap)^a(3) (z - z_Bp)^b(3) exp(-alpha (z - z_Ae)^2) exp(-beta (z - z_Be)^2) + ! ! into - ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) - ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) - ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) + ! fact_k * [sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x] exp (-p (x-P_center(1))^2) + ! * [sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y] exp (-p (y-P_center(2))^2) + ! * [sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z] exp (-p (z-P_center(3))^2) ! ! WARNING ::: IF fact_k is too smal then: ! returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef + ! END_DOC implicit none include 'constants.include.F' integer, intent(in) :: dim, a(3), b(3) - complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3) + complex*16, intent(in) :: alpha, Ap_center(3), Ae_center(3) + complex*16, intent(in) :: beta, Bp_center(3), Be_center(3) integer, intent(out) :: iorder(3) complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3) integer :: n_new, i, j - double precision :: tmp_mod complex*16 :: P_a(0:max_dim,3), P_b(0:max_dim,3) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b @@ -97,12 +109,11 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, P_new(0,3) = (0.d0, 0.d0) !DIR$ FORCEINLINE - call cgaussian_product(alpha, A_center, beta, B_center, fact_k, p, P_center) + call cgaussian_product(alpha, Ae_center, beta, Be_center, fact_k, p, P_center) ! IF fact_k is too smal then: returns a "s" function centered in zero ! with an inifinite exponent and a zero polynom coef - tmp_mod = dsqrt(real(fact_k)*real(fact_k) + aimag(fact_k)*aimag(fact_k)) - if(tmp_mod < 1d-14) then + if(zabs(fact_k) < 1d-14) then iorder = 0 p = (1.d+14, 0.d0) fact_k = (0.d0 , 0.d0) @@ -112,7 +123,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, endif !DIR$ FORCEINLINE - call recentered_cpoly2(P_a(0,1), A_center(1), P_center(1), a(1), P_b(0,1), B_center(1), P_center(1), b(1)) + call recentered_cpoly2(P_a(0,1), Ap_center(1), P_center(1), a(1), P_b(0,1), Bp_center(1), P_center(1), b(1)) iorder(1) = a(1) + b(1) do i = 0, iorder(1) P_new(i,1) = 0.d0 @@ -122,7 +133,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, call multiply_cpoly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new) !DIR$ FORCEINLINE - call recentered_cpoly2(P_a(0,2), A_center(2), P_center(2), a(2), P_b(0,2), B_center(2), P_center(2), b(2)) + call recentered_cpoly2(P_a(0,2), Ap_center(2), P_center(2), a(2), P_b(0,2), Bp_center(2), P_center(2), b(2)) iorder(2) = a(2) + b(2) do i = 0, iorder(2) P_new(i,2) = 0.d0 @@ -132,7 +143,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, call multiply_cpoly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new) !DIR$ FORCEINLINE - call recentered_cpoly2(P_a(0,3), A_center(3), P_center(3), a(3), P_b(0,3), B_center(3), P_center(3), b(3)) + call recentered_cpoly2(P_a(0,3), Ap_center(3), P_center(3), a(3), P_b(0,3), Bp_center(3), P_center(3), b(3)) iorder(3) = a(3) + b(3) do i = 0, iorder(3) P_new(i,3) = 0.d0 @@ -199,15 +210,15 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) END_DOC implicit none + complex*16, intent(in) :: a, b, xa, xb complex*16, intent(out) :: p, k, xp - double precision :: tmp_mod complex*16 :: p_inv complex*16 :: xab, ab - ASSERT (REAL(a) > 0.) - ASSERT (REAL(b) > 0.) + ASSERT (real(a) > 0.) + ASSERT (real(b) > 0.) ! new center p = a + b @@ -217,9 +228,8 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) p_inv = (1.d0, 0.d0) / p ab = a * b * p_inv - k = ab * xab*xab - tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k)) - if(tmp_mod > 40.d0) then + k = ab * xab*xab + if(zabs(k) > 40.d0) then k = (0.d0, 0.d0) xp = (0.d0, 0.d0) return @@ -337,7 +347,7 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) enddo tmp = d(nd) - tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) + tmp_mod = dsqrt(real(tmp)*real(tmp) + aimag(tmp)*aimag(tmp)) do while(tmp_mod .lt. 1.d-15) nd -= 1 if(nd < 0) exit diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f index acac1f6f..844a86a4 100644 --- a/src/utils/cpx_boys.irp.f +++ b/src/utils/cpx_boys.irp.f @@ -12,7 +12,7 @@ complex*16 function crint_1(n, rho) integer :: i, mmax double precision :: rho_mod double precision :: tmp - complex*16 :: sq_rho, rho_inv, rho_exp + complex*16 :: rho_inv, rho_exp complex*16 :: crint_smallz @@ -94,24 +94,14 @@ complex*16 function crint_1(n, rho) else - if(rho_mod .gt. 40.d0) then + rho_exp = 0.5d0 * zexp(-rho) + rho_inv = (1.d0, 0.d0) / rho - crint_1 = 0.5d0 * gamma(dble(n) + 0.5d0) / (rho**n * zsqrt(rho)) + call zboysfun00_1(rho, crint_1) + do i = 1, n + crint_1 = ((dble(i) - 0.5d0) * crint_1 - rho_exp) * rho_inv + enddo - else - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - call zboysfun00_1(rho, crint_1) - mmax = n - if(mmax .gt. 0) then - do i = 0, mmax-1 - crint_1 = ((dble(i) + 0.5d0) * crint_1 - rho_exp) * rho_inv - enddo - endif - - endif endif return @@ -119,91 +109,120 @@ end ! --- -complex*16 function crint_sum_1(n_pt_out, rho, d1) +subroutine crint_1_vec(n_max, rho, vals) implicit none include 'constants.include.F' - integer, intent(in) :: n_pt_out - complex*16, intent(in) :: rho, d1(0:n_pt_out) + integer, intent(in) :: n_max + complex*16, intent(in) :: rho + complex*16, intent(out) :: vals(0:n_max) - integer :: n, i, mmax - double precision :: rho_mod - complex*16 :: sq_rho, F0 - complex*16 :: rho_tmp, rho_inv, rho_exp - complex*16, allocatable :: Fm(:) - - complex*16 :: crint_smallz + integer :: n + double precision :: rho_mod + double precision :: tmp + complex*16 :: rho_inv, rho_exp + complex*16 :: crint_smallz rho_mod = zabs(rho) - if(rho_mod < 10.d0) then - ! small z + if(rho_mod < 3.5d0) then - if(rho_mod .lt. 1.d-15) then + if(rho_mod .lt. 0.35d0) then - crint_sum_1 = d1(0) - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum_1 = crint_sum_1 + d1(i) / dble(n+n+1) - enddo + vals(0) = (((((((((1.3122532963802805073d-08 * rho & + - 1.450385222315046877d-07) * rho & + + 1.458916900093370682d-06) * rho & + - 0.132275132275132275d-04) * rho & + + 0.106837606837606838d-03) * rho & + - 0.757575757575757576d-03) * rho & + + 0.462962962962962963d-02) * rho & + - 0.238095238095238095d-01) * rho & + + 0.10000000000000000000d0) * rho & + - 0.33333333333333333333d0) * rho & + + 1.0d0 + + if(n > 0) then + + vals(1) = (((((((((1.198144314086343d-08 * rho & + - 1.312253296380281d-07) * rho & + + 1.305346700083542d-06) * rho & + - 1.167133520074696d-05) * rho & + + 9.259259259259259d-05) * rho & + - 6.410256410256410d-04) * rho & + + 3.787878787878788d-03) * rho & + - 1.851851851851852d-02) * rho & + + 7.142857142857142d-02) * rho & + - 2.000000000000000d-01) * rho & + + 3.333333333333333d-01 + + if(n > 1) then + + vals(2) = (((((((((1.102292768959436d-08 * rho & + - 1.198144314086343d-07) * rho & + + 1.181027966742252d-06) * rho & + - 1.044277360066834d-05) * rho & + + 8.169934640522875d-05) * rho & + - 5.555555555555556d-04) * rho & + + 3.205128205128205d-03) * rho & + - 1.515151515151515d-02) * rho & + + 5.555555555555555d-02) * rho & + - 1.428571428571428d-01) * rho & + + 2.000000000000000d-01 + + if(n > 2) then + + vals(3) = (((((((((1.020641452740218d-08 * rho & + - 1.102292768959436d-07) * rho & + + 1.078329882677709d-06) * rho & + - 9.448223733938020d-06) * rho & + + 7.309941520467836d-05) * rho & + - 4.901960784313725d-04) * rho & + + 2.777777777777778d-03) * rho & + - 1.282051282051282d-02) * rho & + + 4.545454545454546d-02) * rho & + - 1.111111111111111d-01) * rho & + + 1.428571428571428d-01 + + do n = 4, n_max + tmp = dble(n + n + 1) + vals(n) = (((((((((2.755731922398589d-07 * rho / (tmp + 20.d0) & + - 2.755731922398589d-06 / (tmp + 18.d0)) * rho & + + 2.480158730158730d-05 / (tmp + 16.d0)) * rho & + - 1.984126984126984d-04 / (tmp + 14.d0)) * rho & + + 1.388888888888889d-03 / (tmp + 12.d0)) * rho & + - 8.333333333333333d-03 / (tmp + 10.d0)) * rho & + + 4.166666666666666d-02 / (tmp + 8.d0)) * rho & + - 1.666666666666667d-01 / (tmp + 6.d0)) * rho & + + 5.000000000000000d-01 / (tmp + 4.d0)) * rho & + - 1.000000000000000d+00 / (tmp + 2.d0)) * rho & + + 1.0d0 / tmp + enddo + + endif ! n > 2 + endif ! n > 1 + endif ! n > 0 else - crint_sum_1 = d1(0) * crint_smallz(0, rho) - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum_1 = crint_sum_1 + d1(i) * crint_smallz(n, rho) - enddo + call crint_smallz_vec(n_max, rho, vals) endif else - ! large z - if(rho_mod .gt. 40.d0) then + rho_exp = 0.5d0 * zexp(-rho) + rho_inv = (1.d0, 0.d0) / rho - rho_inv = (1.d0, 0.d0) / rho - rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv) + call zboysfun00_1(rho, vals(0)) + do n = 1, n_max + vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - rho_exp) * rho_inv + enddo - crint_sum_1 = rho_tmp * d1(0) - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv - crint_sum_1 = crint_sum_1 + rho_tmp * d1(i) - enddo - - else - - call zboysfun00_1(rho, F0) - crint_sum_1 = F0 * d1(0) - - rho_exp = 0.5d0 * zexp(-rho) - rho_inv = (1.d0, 0.d0) / rho - - mmax = shiftr(n_pt_out, 1) - if(mmax .gt. 0) then - - allocate(Fm(mmax)) - Fm(1:mmax) = (0.d0, 0.d0) - - do n = 0, mmax-1 - F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv - Fm(n+1) = F0 - enddo - - do i = 2, n_pt_out, 2 - n = shiftr(i, 1) - crint_sum_1 = crint_sum_1 + Fm(n) * d1(i) - enddo - - deallocate(Fm) - endif - - endif ! rho_mod - endif ! rho_mod + endif + return end ! --- @@ -433,7 +452,7 @@ end ! --- -complex*16 function crint_sum_2(n_pt_out, rho, d1) +complex*16 function crint_sum(n_pt_out, rho, d1) implicit none @@ -448,13 +467,14 @@ complex*16 function crint_sum_2(n_pt_out, rho, d1) n_max = shiftr(n_pt_out, 1) allocate(vals(0:n_max)) + call crint_1_vec(n_max, rho, vals) !call crint_2_vec(n_max, rho, vals) ! FOR DEBUG - call crint_quad_12_vec(n_max, rho, vals) + !call crint_quad_12_vec(n_max, rho, vals) - crint_sum_2 = d1(0) * vals(0) + crint_sum = d1(0) * vals(0) do i = 2, n_pt_out, 2 - crint_sum_2 += d1(i) * vals(shiftr(i, 1)) + crint_sum += d1(i) * vals(shiftr(i, 1)) enddo deallocate(vals) From 4da6a2aea07465e56dcea04a0a1d5b362447e788 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 18 Oct 2024 09:50:39 +0200 Subject: [PATCH 3/5] fixed bug in cGTOS overlaps --- src/ao_basis/aos.irp.f | 4 -- src/ao_one_e_ints/aos_cgtos.irp.f | 13 ++++-- src/utils/cgtos_one_e.irp.f | 66 ++++++++++++++++++++----------- src/utils/cgtos_utils.irp.f | 38 +++++++----------- 4 files changed, 66 insertions(+), 55 deletions(-) diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index 1cbd3976..09604487 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -75,10 +75,6 @@ END_PROVIDER enddo endif - powA(1) = ao_power(i,1) - powA(2) = ao_power(i,2) - powA(3) = ao_power(i,3) - ! Normalization of the contracted basis functions if (ao_normalized) then norm = 0.d0 diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index b5dd5263..f89dfd09 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -4,6 +4,7 @@ BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (ao_prim_num_max, ao_num)] implicit none + integer :: i, j do j = 1, ao_num @@ -62,9 +63,9 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] powA(2) = ao_power(i,2) powA(3) = ao_power(i,3) - ! Normalization of the primitives if(primitives_normalized) then + ! Normalization of the primitives do j = 1, ao_prim_num(i) expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im(i,j) @@ -81,11 +82,15 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2) C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2) - call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz) - call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz) + call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, & + C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz) + + call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, & + conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz) norm = 2.d0 * real(C1 * integ1 + C2 * integ2) + !ao_coef_norm_cgtos(i,j) = 1.d0 / dsqrt(norm) ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm) enddo @@ -95,7 +100,7 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)] ao_coef_norm_cgtos(i,j) = ao_coef(i,j) enddo - endif + endif ! primitives_normalized enddo diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index a8773ce0..dffb4e47 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -46,12 +46,13 @@ end ! --- -subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x, overlap_y, overlap_z, overlap, dim) +subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_z, overlap, dim) BEGIN_DOC ! - ! S_x = \int (x - Ap_x)^{a_x} exp(-\alpha (x - Ae_x)^2) (x - Bp_x)^{b_x} exp(-beta (x - Be_x)^2) dx + ! S_x = \int (x - Ap_x)^{a_x} exp(-\alpha (x - Ae_x)^2) + ! (x - Bp_x)^{b_x} exp(-\beta (x - Be_x)^2) dx ! ! S = S_x S_y S_z ! @@ -70,11 +71,12 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow integer :: i, nmax, iorder_p(3) complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p complex*16 :: F_integral_tab(0:max_dim) + complex*16 :: ab, arg - complex*16 :: Fc_integral + complex*16, external :: Fc_integral call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, & - alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim) + alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim) if(zabs(fact_p) .lt. 1.d-14) then overlap_x = (1.d-10, 0.d0) @@ -85,36 +87,52 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow endif nmax = maxval(iorder_p) - inv_sq_p = (1.d0, 0.d0) / zsqrt(p) do i = 0, nmax F_integral_tab(i) = Fc_integral(i, inv_sq_p) enddo - overlap_x = P_new(0,1) * F_integral_tab(0) - overlap_y = P_new(0,2) * F_integral_tab(0) - overlap_z = P_new(0,3) * F_integral_tab(0) + ab = alpha * beta * inv_sq_p - do i = 1, iorder_p(1) - overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) - enddo - call cgaussian_product_x(alpha, Ap_center(1), beta, Bp_center(1), fact_p, p, P_center(1)) - overlap_x *= fact_p + arg = ab * (Ae_center(1) - Be_center(1)) & + * (Ae_center(1) - Be_center(1)) + if(real(arg) > 40.d0) then + overlap_x = (0.d0, 0.d0) + else + overlap_x = P_new(0,1) * F_integral_tab(0) + do i = 1, iorder_p(1) + overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) + enddo + overlap_x = overlap_x * zexp(-arg) + endif - do i = 1, iorder_p(2) - overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) - enddo - call cgaussian_product_x(alpha, Ap_center(2), beta, Bp_center(2), fact_p, p, P_center(2)) - overlap_y *= fact_p + arg = ab * (Ae_center(2) - Be_center(2)) & + * (Ae_center(2) - Be_center(2)) + if(real(arg) > 40.d0) then + overlap_y = (0.d0, 0.d0) + else + overlap_y = P_new(0,2) * F_integral_tab(0) + do i = 1, iorder_p(2) + overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) + enddo + overlap_y = overlap_y * zexp(-arg) + endif - do i = 1, iorder_p(3) - overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) - enddo - call cgaussian_product_x(alpha, Ap_center(3), beta, Bp_center(3), fact_p, p, P_center(3)) - overlap_z *= fact_p + arg = ab * (Ae_center(3) - Be_center(3)) & + * (Ae_center(3) - Be_center(3)) + if(real(arg) > 40.d0) then + overlap_z = (0.d0, 0.d0) + else + overlap_z = P_new(0,3) * F_integral_tab(0) + do i = 1, iorder_p(3) + overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) + enddo + overlap_z = overlap_z * zexp(-arg) + endif overlap = overlap_x * overlap_y * overlap_z + return end ! --- diff --git a/src/utils/cgtos_utils.irp.f b/src/utils/cgtos_utils.irp.f index 40c8a838..075d9bbc 100644 --- a/src/utils/cgtos_utils.irp.f +++ b/src/utils/cgtos_utils.irp.f @@ -66,7 +66,7 @@ end ! --- subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, & - alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim) + alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim) BEGIN_DOC ! @@ -91,8 +91,8 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, include 'constants.include.F' integer, intent(in) :: dim, a(3), b(3) - complex*16, intent(in) :: alpha, Ap_center(3), Ae_center(3) - complex*16, intent(in) :: beta, Bp_center(3), Be_center(3) + complex*16, intent(in) :: alpha, Ae_center(3), Ap_center(3) + complex*16, intent(in) :: beta, Be_center(3), Bp_center(3) integer, intent(out) :: iorder(3) complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3) @@ -167,13 +167,12 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) complex*16, intent(in) :: a, b, xa(3), xb(3) complex*16, intent(out) :: p, k, xp(3) - double precision :: tmp_mod complex*16 :: p_inv, xab(3), ab !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab - ASSERT (REAL(a) > 0.) - ASSERT (REAL(b) > 0.) + ASSERT (real(a) > 0.) + ASSERT (real(b) > 0.) ! new exponent p = a + b @@ -185,9 +184,8 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp) p_inv = (1.d0, 0.d0) / p ab = a * b * p_inv - k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3)) - tmp_mod = dsqrt(real(k)*real(k) + aimag(k)*aimag(k)) - if(tmp_mod .gt. 40.d0) then + k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3)) + if(real(k) .gt. 40.d0) then k = (0.d0, 0.d0) xp(1:3) = (0.d0, 0.d0) return @@ -228,8 +226,8 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp) p_inv = (1.d0, 0.d0) / p ab = a * b * p_inv - k = ab * xab*xab - if(zabs(k) > 40.d0) then + k = ab * xab * xab + if(real(k) > 40.d0) then k = (0.d0, 0.d0) xp = (0.d0, 0.d0) return @@ -300,7 +298,6 @@ subroutine add_cpoly(b, nb, c, nc, d, nd) complex*16, intent(out) :: d(0:nb+nc) integer :: ib - double precision :: tmp_mod complex*16 :: tmp nd = nb + nc @@ -308,12 +305,10 @@ subroutine add_cpoly(b, nb, c, nc, d, nd) d(ib) = d(ib) + c(ib) + b(ib) enddo - tmp = d(nd) - tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) - do while( (tmp_mod .lt. 1.d-15) .and. (nd >= 0) ) + tmp = d(nd) + do while( (zabs(tmp) .lt. 1.d-15) .and. (nd >= 0) ) nd -= 1 - tmp = d(nd) - tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) + tmp = d(nd) if(nd < 0) exit enddo @@ -336,7 +331,6 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) complex*16, intent(inout) :: d(0:max(nb, nd)) integer :: ib - double precision :: tmp_mod complex*16 :: tmp nd = max(nd, nb) @@ -346,13 +340,11 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd) d(ib) = d(ib) + cst * b(ib) enddo - tmp = d(nd) - tmp_mod = dsqrt(real(tmp)*real(tmp) + aimag(tmp)*aimag(tmp)) - do while(tmp_mod .lt. 1.d-15) + tmp = d(nd) + do while(zabs(tmp) .lt. 1.d-15) nd -= 1 if(nd < 0) exit - tmp = d(nd) - tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) + tmp = d(nd) enddo endif From 45481ac08e3b8377161650ff57b7972108f0a7d3 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 18 Oct 2024 13:42:51 +0200 Subject: [PATCH 4/5] fixed pw-centering-related bug in cGTOs kinetic integrals --- src/ao_one_e_ints/aos_cgtos.irp.f | 8 +- .../one_e_kin_integrals_cgtos.irp.f | 232 +++++++++++------- src/utils/cgtos_one_e.irp.f | 8 +- 3 files changed, 157 insertions(+), 91 deletions(-) diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index f89dfd09..f926c846 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -246,11 +246,11 @@ END_PROVIDER C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3))) C2(4) = C2(1) * C2(2) * C2(3) - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) - call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, conjg(Ap_center), Bp_center, & - overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) + call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, & + conjg(Ap_center), Bp_center, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2) overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2) diff --git a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f index b55c37f8..f2557fc3 100644 --- a/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f @@ -10,13 +10,15 @@ double precision :: c, deriv_tmp double precision :: KA2, phiA double precision :: KB2, phiB + double precision :: aa complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3), C1 complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3), C2 + complex*16 :: xa complex*16 :: overlap_x, overlap_y, overlap_z, overlap complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1 complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2 - complex*16 :: overlap_m2_1, overlap_p2_1 - complex*16 :: overlap_m2_2, overlap_p2_2 + complex*16 :: overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1 + complex*16 :: overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2 complex*16 :: deriv_tmp_1, deriv_tmp_2 @@ -32,26 +34,27 @@ beta = (0.1d0, 0.d0) power_A = 1 power_B = 0 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) ! --- - !$OMP PARALLEL DO SCHEDULE(GUIDED) & - !$OMP DEFAULT(NONE) & - !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & - !$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, & - !$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, & - !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & - !$OMP overlap_x, overlap_y, overlap_z, overlap, & - !$OMP overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2, & - !$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, & - !$OMP overlap_y0_2, overlap_z0_2) & - !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & - !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & - !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & - !$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) - + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, aa, xa, C1, C2, & + !$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, & + !$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, & + !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & + !$OMP overlap_x, overlap_y, overlap_z, overlap, & + !$OMP overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1, & + !$OMP overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2, & + !$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, & + !$OMP overlap_y0_2, overlap_z0_2) & + !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & + !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & + !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & + !$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) + !$OMP DO SCHEDULE(GUIDED) do j = 1, ao_num jj = ao_nucl(j) @@ -97,42 +100,62 @@ C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2)) - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) ! --- - power_A(1) = power_A(1) - 2 + power_A(1) = power_A(1) - 1 if(power_A(1) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_m2_1, overlap_y, overlap_z, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_m2_2, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_m1_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_m1_2, overlap_y, overlap_z, overlap, dim1) + power_A(1) = power_A(1) - 1 + if(power_A(1) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_m2_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_m2_2, overlap_y, overlap_z, overlap, dim1) + else + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(1) = power_A(1) + 1 else + overlap_m1_1 = (0.d0, 0.d0) + overlap_m1_2 = (0.d0, 0.d0) overlap_m2_1 = (0.d0, 0.d0) overlap_m2_2 = (0.d0, 0.d0) endif + power_A(1) = power_A(1) + 1 - power_A(1) = power_A(1) + 4 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_p2_1, overlap_y, overlap_z, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_p2_2, overlap_y, overlap_z, overlap, dim1) - + power_A(1) = power_A(1) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_p1_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_p1_2, overlap_y, overlap_z, overlap, dim1) + power_A(1) = power_A(1) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_p2_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_p2_2, overlap_y, overlap_z, overlap, dim1) power_A(1) = power_A(1) - 2 - deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_1 & - + dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_1 & - + 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1 + aa = dble(power_A(1)) + xa = Ap_center(1) - Ae_center(1) - deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_2 & - + dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_2 & - + 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2 + deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_1 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) + deriv_tmp_1 = deriv_tmp_1 * overlap_y0_1 * overlap_z0_1 + + deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_2 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) + deriv_tmp_2 = deriv_tmp_2 * overlap_y0_2 * overlap_z0_2 deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) @@ -140,34 +163,55 @@ ! --- - power_A(2) = power_A(2) - 2 + power_A(2) = power_A(2) - 1 if(power_A(2) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x, overlap_m2_1, overlap_y, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_x, overlap_m2_2, overlap_y, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_m1_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_m1_2, overlap_z, overlap, dim1) + power_A(2) = power_A(2) - 1 + if(power_A(2) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_m2_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_m2_2, overlap_z, overlap, dim1) + else + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(2) = power_A(2) + 1 else + overlap_m1_1 = (0.d0, 0.d0) + overlap_m1_2 = (0.d0, 0.d0) overlap_m2_1 = (0.d0, 0.d0) overlap_m2_2 = (0.d0, 0.d0) endif + power_A(2) = power_A(2) + 1 - power_A(2) = power_A(2) + 4 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x, overlap_p2_1, overlap_y, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_x, overlap_p2_2, overlap_y, overlap, dim1) - + power_A(2) = power_A(2) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_p1_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_p1_2, overlap_z, overlap, dim1) + power_A(2) = power_A(2) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_p2_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_p2_2, overlap_z, overlap, dim1) power_A(2) = power_A(2) - 2 - deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_1 & - + dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_1 & - + 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1 + aa = dble(power_A(2)) + xa = Ap_center(2) - Ae_center(2) - deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_2 & - + dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_2 & - + 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2 + deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_1 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) + deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_z0_1 + + deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_2 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) + deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_z0_2 deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) @@ -175,34 +219,55 @@ ! --- - power_A(3) = power_A(3) - 2 + power_A(3) = power_A(3) - 1 if(power_A(3) > -1) then - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x, overlap_y, overlap_m2_1, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_x, overlap_y, overlap_m2_2, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_m1_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m1_2, overlap, dim1) + power_A(3) = power_A(3) - 1 + if(power_A(3) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_m2_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m2_2, overlap, dim1) + else + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(3) = power_A(3) + 1 else + overlap_m1_1 = (0.d0, 0.d0) + overlap_m1_2 = (0.d0, 0.d0) overlap_m2_1 = (0.d0, 0.d0) overlap_m2_2 = (0.d0, 0.d0) endif + power_A(3) = power_A(3) + 1 - power_A(3) = power_A(3) + 4 - call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, & - overlap_x, overlap_y, overlap_p2_1, overlap, dim1) - - call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, Ap_center, conjg(Bp_center), & - overlap_x, overlap_y, overlap_p2_2, overlap, dim1) - + power_A(3) = power_A(3) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_p1_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p1_2, overlap, dim1) + power_A(3) = power_A(3) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_p2_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p2_2, overlap, dim1) power_A(3) = power_A(3) - 2 - - deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_1 & - + dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_1 & - + 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1 - deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_2 & - + dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_2 & - + 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2 + aa = dble(power_A(3)) + xa = Ap_center(3) - Ae_center(3) + + deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_1 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) + deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_y0_1 + + deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_2 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) + deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_y0_2 deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) @@ -214,7 +279,8 @@ enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + !$OMP END PARALLEL END_PROVIDER diff --git a/src/utils/cgtos_one_e.irp.f b/src/utils/cgtos_one_e.irp.f index dffb4e47..83616969 100644 --- a/src/utils/cgtos_one_e.irp.f +++ b/src/utils/cgtos_one_e.irp.f @@ -79,10 +79,10 @@ subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, pow alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim) if(zabs(fact_p) .lt. 1.d-14) then - overlap_x = (1.d-10, 0.d0) - overlap_y = (1.d-10, 0.d0) - overlap_z = (1.d-10, 0.d0) - overlap = (1.d-10, 0.d0) + overlap_x = (0.d0, 0.d0) + overlap_y = (0.d0, 0.d0) + overlap_z = (0.d0, 0.d0) + overlap = (0.d0, 0.d0) return endif From f6728533f988fefea323acdbded36a90d5ecaa62 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 18 Oct 2024 19:04:49 +0200 Subject: [PATCH 5/5] fixed index-4 bug in use_pw --- src/ao_one_e_ints/aos_cgtos.irp.f | 4 ++-- .../one_e_coul_integrals_cgtos.irp.f | 4 ++-- .../two_e_coul_integrals_cgtos.irp.f | 4 ++-- src/utils/cpx_boys.irp.f | 24 +++++++++---------- 4 files changed, 18 insertions(+), 18 deletions(-) diff --git a/src/ao_one_e_ints/aos_cgtos.irp.f b/src/ao_one_e_ints/aos_cgtos.irp.f index f926c846..62c312f8 100644 --- a/src/ao_one_e_ints/aos_cgtos.irp.f +++ b/src/ao_one_e_ints/aos_cgtos.irp.f @@ -18,8 +18,8 @@ END_PROVIDER ! --- BEGIN_PROVIDER [complex*16, ao_expo_cgtos_ord_transp, (ao_prim_num_max, ao_num)] -&BEGIN_PROVIDER [complex*16, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)] -&BEGIN_PROVIDER [complex*16, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)] +&BEGIN_PROVIDER [double precision, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)] +&BEGIN_PROVIDER [double precision, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)] implicit none diff --git a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f index fed5f29d..9294b139 100644 --- a/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f +++ b/src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f @@ -129,6 +129,7 @@ complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, a complex*16, external :: V_n_e_cgtos complex*16, external :: crint_sum + complex*16, external :: crint_1 @@ -178,8 +179,7 @@ complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, a n_pt = 2 * ((power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3))) if(n_pt == 0) then - !NAI_pol_mult_cgtos = coeff * crint_1(0, const) - NAI_pol_mult_cgtos = coeff * crint_sum(0, const, (1.d0, 0.d0)) + NAI_pol_mult_cgtos = coeff * crint_1(0, const) return endif diff --git a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f index 9414e851..3cdb6b73 100644 --- a/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f +++ b/src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f @@ -1690,8 +1690,8 @@ BEGIN_PROVIDER [logical, use_pw] call ezfio_has_ao_basis_ao_expo_pw(exist) if(exist) then - PROVIDE ao_expo_pw - if(maxval(dabs(ao_expo_pw(4,:,:))) .gt. 1d-15) use_pw = .true. + PROVIDE ao_expo_pw_ord_transp + if(maxval(dabs(ao_expo_pw_ord_transp(4,:,:))) .gt. 1d-15) use_pw = .true. endif END_PROVIDER diff --git a/src/utils/cpx_boys.irp.f b/src/utils/cpx_boys.irp.f index 844a86a4..6ddc15db 100644 --- a/src/utils/cpx_boys.irp.f +++ b/src/utils/cpx_boys.irp.f @@ -143,7 +143,7 @@ subroutine crint_1_vec(n_max, rho, vals) - 0.33333333333333333333d0) * rho & + 1.0d0 - if(n > 0) then + if(n_max > 0) then vals(1) = (((((((((1.198144314086343d-08 * rho & - 1.312253296380281d-07) * rho & @@ -157,7 +157,7 @@ subroutine crint_1_vec(n_max, rho, vals) - 2.000000000000000d-01) * rho & + 3.333333333333333d-01 - if(n > 1) then + if(n_max > 1) then vals(2) = (((((((((1.102292768959436d-08 * rho & - 1.198144314086343d-07) * rho & @@ -171,7 +171,7 @@ subroutine crint_1_vec(n_max, rho, vals) - 1.428571428571428d-01) * rho & + 2.000000000000000d-01 - if(n > 2) then + if(n_max > 2) then vals(3) = (((((((((1.020641452740218d-08 * rho & - 1.102292768959436d-07) * rho & @@ -200,9 +200,9 @@ subroutine crint_1_vec(n_max, rho, vals) + 1.0d0 / tmp enddo - endif ! n > 2 - endif ! n > 1 - endif ! n > 0 + endif ! n_max > 2 + endif ! n_max > 1 + endif ! n_max > 0 else @@ -515,7 +515,7 @@ subroutine crint_2_vec(n_max, rho, vals) - 0.33333333333333333333d0) * rho & + 1.0d0 - if(n > 0) then + if(n_max > 0) then vals(1) = (((((((((1.198144314086343d-08 * rho & - 1.312253296380281d-07) * rho & @@ -529,7 +529,7 @@ subroutine crint_2_vec(n_max, rho, vals) - 2.000000000000000d-01) * rho & + 3.333333333333333d-01 - if(n > 1) then + if(n_max > 1) then vals(2) = (((((((((1.102292768959436d-08 * rho & - 1.198144314086343d-07) * rho & @@ -543,7 +543,7 @@ subroutine crint_2_vec(n_max, rho, vals) - 1.428571428571428d-01) * rho & + 2.000000000000000d-01 - if(n > 2) then + if(n_max > 2) then vals(3) = (((((((((1.020641452740218d-08 * rho & - 1.102292768959436d-07) * rho & @@ -572,9 +572,9 @@ subroutine crint_2_vec(n_max, rho, vals) + 1.0d0 / tmp enddo - endif ! n > 2 - endif ! n > 1 - endif ! n > 0 + endif ! n_max > 2 + endif ! n_max > 1 + endif ! n_max > 0 else