diff --git a/src/tc_bi_ortho/no_dressing_naive.irp.f b/src/tc_bi_ortho/no_dressing_naive.irp.f new file mode 100644 index 00000000..10c80cec --- /dev/null +++ b/src/tc_bi_ortho/no_dressing_naive.irp.f @@ -0,0 +1,92 @@ + +! --- + +BEGIN_PROVIDER [double precision, no_0_naive] + + implicit none + integer :: ii, jj, kk + integer :: i, j, k + double precision :: sigma_i, sigma_j, sigma_k + double precision :: tmp + double precision :: I_ijk_ijk, I_ijk_kij, I_ijk_jki, I_ijk_jik, I_ijk_kji, I_ijk_ikj + double precision :: t0, t1 + logical, external :: is_same_spin + + print*, " Providing no_0_naive ..." + call wall_time(t0) + + + tmp = 0.d0 + do ii = 1, elec_num + + if(ii .le. elec_beta_num) then + i = ii + sigma_i = -1.d0 + else + i = ii - elec_beta_num + sigma_i = +1.d0 + endif + + do jj = 1, elec_num + + if(jj .le. elec_beta_num) then + j = jj + sigma_j = -1.d0 + else + j = jj - elec_beta_num + sigma_j = +1.d0 + endif + + do kk = 1, elec_num + + if(kk .le. elec_beta_num) then + k = kk + sigma_k = -1.d0 + else + k = kk - elec_beta_num + sigma_k = +1.d0 + endif + + call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k & + , i, sigma_i, j, sigma_j, k, sigma_k & + , I_ijk_ijk) + + call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k & + , k, sigma_k, i, sigma_i, j, sigma_j & + , I_ijk_kij) + + call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k & + , j, sigma_j, k, sigma_k, i, sigma_i & + , I_ijk_jki) + + call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k & + , j, sigma_j, i, sigma_i, k, sigma_k & + , I_ijk_jik) + + call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k & + , k, sigma_k, j, sigma_j, i, sigma_i & + , I_ijk_kji) + + call give_integrals_3_body_bi_ort_spin( i, sigma_i, j, sigma_j, k, sigma_k & + , i, sigma_i, k, sigma_k, j, sigma_j & + , I_ijk_ikj) + + + tmp = tmp + I_ijk_ijk + I_ijk_kij + I_ijk_jki - I_ijk_jik - I_ijk_kji - I_ijk_ikj + !tmp = tmp + I_ijk_ijk + 2.d0 * I_ijk_kij - 3.d0 * I_ijk_jik + enddo + enddo + enddo + + no_0_naive = -1.d0 * (-tmp) / 6.d0 + + call wall_time(t1) + print*, " Wall time for no_0_naive (sec) = ", (t1 - t0)/60.d0 + + print*, " no_0_naive = ", no_0_naive + +END_PROVIDER + +! --- + + diff --git a/src/tc_bi_ortho/no_dressing_v0.irp.f b/src/tc_bi_ortho/no_dressing_v0.irp.f new file mode 100644 index 00000000..9a070dab --- /dev/null +++ b/src/tc_bi_ortho/no_dressing_v0.irp.f @@ -0,0 +1,105 @@ + +! --- + +BEGIN_PROVIDER [double precision, no_0_v0] + + implicit none + integer :: i, j, k + double precision :: tmp + double precision :: I_ijk_ijk, I_ijk_kij, I_ijk_jik, I_ijk_jki, I_ijk_ikj, I_ijk_kji + double precision :: t0, t1 + + call wall_time(t0) + print*, " Providing no_0_v0 ..." + + if(elec_alpha_num .eq. elec_beta_num) then + + tmp = 0.d0 + + do i = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(i, j, k, i, j, k, I_ijk_ijk) + call give_integrals_3_body_bi_ort(i, j, k, k, i, j, I_ijk_kij) + call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik) + + tmp = tmp + 4.d0 * (2.d0 * I_ijk_ijk + I_ijk_kij - 3.d0 * I_ijk_jik) + enddo + enddo + enddo + + no_0_v0 = -1.d0 * (-tmp) / 6.d0 + + else + + tmp = 0.d0 + + do i = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(i, j, k, i, j, k, I_ijk_ijk) + call give_integrals_3_body_bi_ort(i, j, k, k, i, j, I_ijk_kij) + call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik) + + tmp = tmp + 4.d0 * (2.d0 * I_ijk_ijk + I_ijk_kij - 3.d0 * I_ijk_jik) + enddo + enddo + enddo + + do i = elec_beta_num+1, elec_alpha_num + do j = elec_beta_num+1, elec_alpha_num + do k = elec_beta_num+1, elec_alpha_num + + call give_integrals_3_body_bi_ort(i, j, k, i, j, k, I_ijk_ijk) + call give_integrals_3_body_bi_ort(i, j, k, k, i, j, I_ijk_kij) + call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik) + + tmp = tmp + I_ijk_ijk + 2.d0 * I_ijk_kij - 3.d0 * I_ijk_jik + enddo + enddo + enddo + + do i = elec_beta_num+1, elec_alpha_num + do j = 1, elec_beta_num + + do k = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(i, j, k, i, j, k, I_ijk_ijk) + call give_integrals_3_body_bi_ort(i, j, k, j, k, i, I_ijk_jki) + call give_integrals_3_body_bi_ort(i, j, k, i, k, j, I_ijk_ikj) + call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik) + call give_integrals_3_body_bi_ort(i, j, k, k, j, i, I_ijk_kji) + + tmp = tmp + 6.d0 * (2.d0 * I_ijk_ijk + I_ijk_jki - I_ijk_ikj - I_ijk_jik - I_ijk_kji) + enddo + + do k = elec_beta_num+1, elec_alpha_num + + call give_integrals_3_body_bi_ort(i, j, k, i, j, k, I_ijk_ijk) + call give_integrals_3_body_bi_ort(i, j, k, j, k, i, I_ijk_jki) + call give_integrals_3_body_bi_ort(i, j, k, i, k, j, I_ijk_ikj) + call give_integrals_3_body_bi_ort(i, j, k, j, i, k, I_ijk_jik) + call give_integrals_3_body_bi_ort(i, j, k, k, j, i, I_ijk_kji) + + tmp = tmp + 3.d0 * (2.d0 * I_ijk_ijk + 2.d0 * I_ijk_jki - I_ijk_ikj - I_ijk_jik - 2.d0 * I_ijk_kji) + enddo + + enddo + enddo + + no_0_v0 = -1.d0 * (-tmp) / 6.d0 + + endif + + call wall_time(t1) + print*, " Wall time for no_0_v0 (sec) = ", (t1 - t0)/60.d0 + + print*, " no_0_v0 = ", no_0_v0 + +END_PROVIDER + +! --- + + diff --git a/src/tc_bi_ortho/tc_utils.irp.f b/src/tc_bi_ortho/tc_utils.irp.f index 53fe5884..ba5ffff9 100644 --- a/src/tc_bi_ortho/tc_utils.irp.f +++ b/src/tc_bi_ortho/tc_utils.irp.f @@ -66,3 +66,25 @@ end ! --- + +logical function is_same_spin(sigma_1, sigma_2) + + BEGIN_DOC + ! + ! true if sgn(sigma_1) = sgn(sigma_2) + ! + END_DOC + + implicit none + double precision, intent(in) :: sigma_1, sigma_2 + + if((sigma_1 * sigma_2) .gt. 0.d0) then + is_same_spin = .true. + else + is_same_spin = .false. + endif + +end function is_same_spin + +! --- + diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f index d25a1f70..226854ed 100644 --- a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -35,7 +35,9 @@ program tc_bi_ortho !call test_no_aaa() !call test_no() - call test_no_v0() + !call test_no_v0() + + call test_no_0() end @@ -315,7 +317,7 @@ subroutine test_no_v0() print*, ' accu (%) = ', 100.d0*accu/norm return -end subroutine test_no +end subroutine test_no_0 ! --- @@ -498,3 +500,23 @@ end ! --- +subroutine test_no_0() + + implicit none + double precision :: accu, norm + + print*, ' testing test_no_0 ...' + + PROVIDE no_0_naive + PROVIDE no_0_v0 + + accu = dabs(no_0_naive - no_0_v0) + norm = dabs(no_0_naive) + + print*, ' accu (%) = ', 100.d0*accu/norm + + return +end + +! --- +