diff --git a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f index ce343f9b..3a68ffc6 100644 --- a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f +++ b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -113,7 +113,6 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] ! --- allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num)) - tmp_2 = 0.d0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & @@ -347,7 +346,7 @@ END_PROVIDER ! --- -BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a_old, (mo_num, mo_num)] BEGIN_DOC ! @@ -366,17 +365,18 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] PROVIDE mo_l_coef mo_r_coef PROVIDE fock_3e_uhf_mo_cs - !print *, ' Providing fock_3e_uhf_mo_a ...' - !call wall_time(ti) + print *, ' Providing fock_3e_uhf_mo_a_old ...' + call wall_time(ti) o = elec_beta_num + 1 call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij) - fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs + PROVIDE fock_3e_uhf_mo_cs + fock_3e_uhf_mo_a_old = fock_3e_uhf_mo_cs !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) & - !$OMP SHARED (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_a) + !$OMP SHARED (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_a_old) allocate(tmp(mo_num,mo_num)) tmp = 0.d0 @@ -460,7 +460,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] !$OMP CRITICAL do a = 1, mo_num do b = 1, mo_num - fock_3e_uhf_mo_a(b,a) += tmp(b,a) + fock_3e_uhf_mo_a_old(b,a) += tmp(b,a) enddo enddo !$OMP END CRITICAL @@ -468,19 +468,21 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] deallocate(tmp) !$OMP END PARALLEL - !call wall_time(tf) - !print *, ' Wall time for fock_3e_uhf_mo_a =', tf - ti + call wall_time(tf) + print *, ' Wall time for fock_3e_uhf_mo_a_old =', tf - ti END_PROVIDER ! --- -BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b_old, (mo_num, mo_num)] BEGIN_DOC + ! ! BETA part of the Fock matrix from three-electron terms ! ! WARNING :: non hermitian if bi-ortho MOS used + ! END_DOC implicit none @@ -491,17 +493,18 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] PROVIDE mo_l_coef mo_r_coef - !print *, ' PROVIDING fock_3e_uhf_mo_b ...' - !call wall_time(ti) + print *, ' PROVIDING fock_3e_uhf_mo_b_old ...' + call wall_time(ti) o = elec_beta_num + 1 call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij) - fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs + PROVIDE fock_3e_uhf_mo_cs_old + fock_3e_uhf_mo_b_old = fock_3e_uhf_mo_cs_old !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) & - !$OMP SHARED (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_b) + !$OMP SHARED (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_b_old) allocate(tmp(mo_num,mo_num)) tmp = 0.d0 @@ -575,7 +578,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] !$OMP CRITICAL do a = 1, mo_num do b = 1, mo_num - fock_3e_uhf_mo_b(b,a) += tmp(b,a) + fock_3e_uhf_mo_b_old(b,a) += tmp(b,a) enddo enddo !$OMP END CRITICAL @@ -583,8 +586,8 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] deallocate(tmp) !$OMP END PARALLEL - !call wall_time(tf) - !print *, ' total Wall time for fock_3e_uhf_mo_b =', tf - ti + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_mo_b_old =', tf - ti END_PROVIDER @@ -760,5 +763,421 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)] END_PROVIDER +! --- + + BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] +&BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] + + BEGIN_DOC + ! + ! Open Shell part of the Fock matrix from three-electron terms + ! + ! WARNING :: non hermitian if bi-ortho MOS used + ! + END_DOC + + implicit none + integer :: a, b, i, j, ipoint, o + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: loc_1, loc_2, loc_3, loc_4 + double precision :: ti, tf + double precision, allocatable :: tmp(:,:) + double precision, allocatable :: Okappa(:), Jkappa(:,:), Obarkappa(:), Jbarkappa(:,:) + double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:) + double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:) + double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:) + + PROVIDE mo_l_coef mo_r_coef + PROVIDE fock_3e_uhf_mo_cs + + print *, ' Providing fock_3e_uhf_mo_a and fock_3e_uhf_mo_b ...' + call wall_time(ti) + + o = elec_beta_num + 1 + call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, I_bij_aij) + + PROVIDE fock_3e_uhf_mo_cs + fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs + fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs + + ! --- + + allocate(Jkappa(n_points_final_grid,3), Okappa(n_points_final_grid)) + allocate(Jbarkappa(n_points_final_grid,3), Obarkappa(n_points_final_grid)) + Jkappa = 0.d0 + Okappa = 0.d0 + Jbarkappa = 0.d0 + Obarkappa = 0.d0 + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, tmp_omp_d1, tmp_omp_d2) & + !$OMP SHARED (n_points_final_grid, elec_beta_num, elec_alpha_num, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP int2_grad1_u12_bimo_t, Okappa, Jkappa, Obarkappa, Jbarkappa) + + allocate(tmp_omp_d2(n_points_final_grid,3), tmp_omp_d1(n_points_final_grid)) + + tmp_omp_d2 = 0.d0 + tmp_omp_d1 = 0.d0 + !$OMP DO + do i = 1, elec_beta_num + do ipoint = 1, n_points_final_grid + tmp_omp_d2(ipoint,1) += int2_grad1_u12_bimo_t(ipoint,1,i,i) + tmp_omp_d2(ipoint,2) += int2_grad1_u12_bimo_t(ipoint,2,i,i) + tmp_omp_d2(ipoint,3) += int2_grad1_u12_bimo_t(ipoint,3,i,i) + tmp_omp_d1(ipoint) += mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i) + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + do ipoint = 1, n_points_final_grid + Jkappa(ipoint,1) += tmp_omp_d2(ipoint,1) + Jkappa(ipoint,2) += tmp_omp_d2(ipoint,2) + Jkappa(ipoint,3) += tmp_omp_d2(ipoint,3) + Okappa(ipoint) += tmp_omp_d1(ipoint) + enddo + !$OMP END CRITICAL + + tmp_omp_d2 = 0.d0 + tmp_omp_d1 = 0.d0 + !$OMP DO + do i = elec_beta_num+1, elec_alpha_num + do ipoint = 1, n_points_final_grid + tmp_omp_d2(ipoint,1) += int2_grad1_u12_bimo_t(ipoint,1,i,i) + tmp_omp_d2(ipoint,2) += int2_grad1_u12_bimo_t(ipoint,2,i,i) + tmp_omp_d2(ipoint,3) += int2_grad1_u12_bimo_t(ipoint,3,i,i) + tmp_omp_d1(ipoint) += mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i) + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + do ipoint = 1, n_points_final_grid + Jbarkappa(ipoint,1) += tmp_omp_d2(ipoint,1) + Jbarkappa(ipoint,2) += tmp_omp_d2(ipoint,2) + Jbarkappa(ipoint,3) += tmp_omp_d2(ipoint,3) + Obarkappa(ipoint) += tmp_omp_d1(ipoint) + enddo + !$OMP END CRITICAL + + deallocate(tmp_omp_d2, tmp_omp_d1) + !$OMP END PARALLEL + + ! --- + + allocate(tmp_1(n_points_final_grid,5)) + + do ipoint = 1, n_points_final_grid + + loc_1 = -2.d0 * Okappa (ipoint) + loc_2 = -2.d0 * Obarkappa(ipoint) + + tmp_1(ipoint,1) = loc_1 * Jbarkappa(ipoint,1) + loc_2 * Jkappa(ipoint,1) + tmp_1(ipoint,2) = loc_1 * Jbarkappa(ipoint,2) + loc_2 * Jkappa(ipoint,2) + tmp_1(ipoint,3) = loc_1 * Jbarkappa(ipoint,3) + loc_2 * Jkappa(ipoint,3) + + tmp_1(ipoint,4) = Obarkappa(ipoint) + tmp_1(ipoint,5) = -loc_1 + enddo + + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, i, j, loc_1, loc_2, tmp_omp_d2) & + !$OMP SHARED (n_points_final_grid, elec_beta_num, elec_alpha_num, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP int2_grad1_u12_bimo_t, tmp_1) + + allocate(tmp_omp_d2(n_points_final_grid,3)) + tmp_omp_d2 = 0.d0 + + !$OMP DO COLLAPSE(2) + do i = 1, elec_beta_num + do j = elec_beta_num+1, elec_alpha_num + do ipoint = 1, n_points_final_grid + + loc_1 = mos_l_in_r_array_transp(ipoint,j) * mos_r_in_r_array_transp(ipoint,i) + loc_2 = mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,j) + + tmp_omp_d2(ipoint,1) += loc_1 * int2_grad1_u12_bimo_t(ipoint,1,i,j) + loc_2 * int2_grad1_u12_bimo_t(ipoint,1,j,i) + tmp_omp_d2(ipoint,2) += loc_1 * int2_grad1_u12_bimo_t(ipoint,2,i,j) + loc_2 * int2_grad1_u12_bimo_t(ipoint,2,j,i) + tmp_omp_d2(ipoint,3) += loc_1 * int2_grad1_u12_bimo_t(ipoint,3,i,j) + loc_2 * int2_grad1_u12_bimo_t(ipoint,3,j,i) + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do ipoint = 1, n_points_final_grid + tmp_1(ipoint,1) += tmp_omp_d2(ipoint,1) + tmp_1(ipoint,2) += tmp_omp_d2(ipoint,2) + tmp_1(ipoint,3) += tmp_omp_d2(ipoint,3) + enddo + !$OMP END CRITICAL + + deallocate(tmp_omp_d2) + !$OMP END PARALLEL + + ! --- + + allocate(tmp_2(n_points_final_grid,5,mo_num,mo_num)) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, a, b) & + !$OMP SHARED (n_points_final_grid, mo_num, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, & + !$OMP tmp_2) + !$OMP DO COLLAPSE(2) + do a = 1, mo_num + do b = 1, mo_num + do ipoint = 1, n_points_final_grid + tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a) + tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a) + tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, a, b, i) & + !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t, & + !$OMP tmp_2) + !$OMP DO COLLAPSE(2) + do a = 1, mo_num + do b = 1, mo_num + + tmp_2(:,4,b,a) = 0.d0 + do i = 1, elec_beta_num + do ipoint = 1, n_points_final_grid + tmp_2(ipoint,4,b,a) += final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) & + + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) & + + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) ) + enddo + enddo + + tmp_2(:,5,b,a) = 0.d0 + do i = elec_beta_num+1, elec_alpha_num + do ipoint = 1, n_points_final_grid + tmp_2(ipoint,5,b,a) += final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) & + + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) & + + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) ) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + ! --- + + call dgemv( 'T', 5*n_points_final_grid, mo_num*mo_num, 1.d0 & + , tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) & + , tmp_1(1,1), 1 & + , 1.d0, fock_3e_uhf_mo_a(1,1), 1) + + deallocate(tmp_1, tmp_2) + + ! --- + + allocate(tmp_3(n_points_final_grid,8,mo_num), tmp_4(n_points_final_grid,8,mo_num)) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, b, loc_1, loc_2) & + !$OMP SHARED (n_points_final_grid, mo_num, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP final_weight_at_r_vector, Jkappa, Jbarkappa, tmp_3, tmp_4) + !$OMP DO + do b = 1, mo_num + tmp_3(:,:,b) = 0.d0 + tmp_4(:,:,b) = 0.d0 + do ipoint = 1, n_points_final_grid + + tmp_3(ipoint,1,b) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,b) + + loc_1 = -2.d0 * mos_r_in_r_array_transp(ipoint,b) + tmp_4(ipoint,1,b) = loc_1 * ( Jbarkappa(ipoint,1) * Jkappa(ipoint,1) & + + Jbarkappa(ipoint,2) * Jkappa(ipoint,2) & + + Jbarkappa(ipoint,3) * Jkappa(ipoint,3) ) + + tmp_4(ipoint,8,b) = mos_r_in_r_array_transp(ipoint,b) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, b, i, loc_1, loc_2, loc_3, loc_4) & + !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP Jkappa, Jbarkappa, tmp_3, tmp_4) + !$OMP DO + do b = 1, mo_num + + do i = 1, elec_beta_num + do ipoint = 1, n_points_final_grid + + loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) + loc_2 = mos_r_in_r_array_transp(ipoint,i) + + tmp_3(ipoint,2,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,1,b,i) + tmp_3(ipoint,3,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,2,b,i) + tmp_3(ipoint,4,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,3,b,i) + tmp_3(ipoint,8,b) += loc_1 * ( Jbarkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,b,i) & + + Jbarkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,b,i) & + + Jbarkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,b,i) ) + + tmp_4(ipoint,5,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,1,i,b) + tmp_4(ipoint,6,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,2,i,b) + tmp_4(ipoint,7,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,3,i,b) + tmp_4(ipoint,1,b) += loc_2 * ( Jbarkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,i,b) & + + Jbarkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,i,b) & + + Jbarkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,i,b) ) + enddo + enddo + + do i = elec_beta_num+1, elec_alpha_num + do ipoint = 1, n_points_final_grid + + loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) + loc_3 = 2.d0 * loc_1 + loc_2 = mos_r_in_r_array_transp(ipoint,i) + loc_4 = 2.d0 * loc_2 + + tmp_3(ipoint,5,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,1,b,i) + tmp_3(ipoint,6,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,2,b,i) + tmp_3(ipoint,7,b) -= loc_1 * int2_grad1_u12_bimo_t(ipoint,3,b,i) + tmp_3(ipoint,8,b) += loc_3 * ( Jkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,b,i) & + + Jkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,b,i) & + + Jkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,b,i) ) + + tmp_4(ipoint,2,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,1,i,b) + tmp_4(ipoint,3,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,2,i,b) + tmp_4(ipoint,4,b) += loc_2 * int2_grad1_u12_bimo_t(ipoint,3,i,b) + tmp_4(ipoint,1,b) += loc_4 * ( Jkappa(ipoint,1) * int2_grad1_u12_bimo_t(ipoint,1,i,b) & + + Jkappa(ipoint,2) * int2_grad1_u12_bimo_t(ipoint,2,i,b) & + + Jkappa(ipoint,3) * int2_grad1_u12_bimo_t(ipoint,3,i,b) ) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint, b, i, j, loc_1, loc_2, loc_3) & + !$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, elec_alpha_num, & + !$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t, & + !$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, & + !$OMP tmp_3, tmp_4) + !$OMP DO + do b = 1, mo_num + do i = 1, elec_beta_num + do j = elec_beta_num+1, elec_alpha_num + do ipoint = 1, n_points_final_grid + + loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,j) + loc_2 = mos_r_in_r_array_transp(ipoint,b) + loc_3 = mos_r_in_r_array_transp(ipoint,i) + + tmp_3(ipoint,8,b) -= loc_1 * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,j) & + + int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,j) & + + int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,j) ) + + tmp_4(ipoint,1,b) -= loc_3 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,b) & + + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,b) & + + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,b) ) + tmp_4(ipoint,1,b) += loc_2 * ( int2_grad1_u12_bimo_t(ipoint,1,i,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) & + + int2_grad1_u12_bimo_t(ipoint,2,i,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) & + + int2_grad1_u12_bimo_t(ipoint,3,i,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) ) + + loc_1 = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) + loc_3 = mos_r_in_r_array_transp(ipoint,j) + + tmp_3(ipoint,8,b) -= loc_1 * ( int2_grad1_u12_bimo_t(ipoint,1,b,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) & + + int2_grad1_u12_bimo_t(ipoint,2,b,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) & + + int2_grad1_u12_bimo_t(ipoint,3,b,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) ) + + tmp_4(ipoint,1,b) -= loc_3 * ( int2_grad1_u12_bimo_t(ipoint,1,j,i) * int2_grad1_u12_bimo_t(ipoint,1,i,b) & + + int2_grad1_u12_bimo_t(ipoint,2,j,i) * int2_grad1_u12_bimo_t(ipoint,2,i,b) & + + int2_grad1_u12_bimo_t(ipoint,3,j,i) * int2_grad1_u12_bimo_t(ipoint,3,i,b) ) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + ! --- + + call dgemm( 'T', 'N', mo_num, mo_num, 8*n_points_final_grid, 1.d0 & + , tmp_3(1,1,1), 8*n_points_final_grid & + , tmp_4(1,1,1), 8*n_points_final_grid & + , 1.d0, fock_3e_uhf_mo_a(1,1), mo_num) + + deallocate(tmp_3, tmp_4) + deallocate(Jkappa, Okappa) + + ! --- + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (a, b, i, j, I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia, tmp) & + !$OMP SHARED (mo_num, o, elec_alpha_num, elec_beta_num, fock_3e_uhf_mo_a) + + allocate(tmp(mo_num,mo_num)) + tmp = 0.d0 + + !$OMP DO + do a = 1, mo_num + do b = 1, mo_num + + do j = o, elec_alpha_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + tmp(b,a) -= 0.5d0 * ( I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - I_bij_iaj & + - I_bij_jia ) + + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do a = 1, mo_num + do b = 1, mo_num + fock_3e_uhf_mo_a(b,a) += tmp(b,a) + enddo + enddo + !$OMP END CRITICAL + + deallocate(tmp) + !$OMP END PARALLEL + + call wall_time(tf) + print *, ' Wall time for fock_3e_uhf_mo_a =', tf - ti + +END_PROVIDER + ! --- diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 0ae515bb..fcca29ac 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -190,30 +190,14 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] if(bi_ortho) then - !allocate(tmp(ao_num,ao_num)) - !tmp = Fock_matrix_tc_ao_alpha - !if(three_body_h_tc) then - ! tmp += fock_3e_uhf_ao_a - !endif - !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1)) - !deallocate(tmp) - PROVIDE mo_l_coef mo_r_coef - !call wall_time(tt0) call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) - !call wall_time(tt1) - !print*, ' 2-e term:', tt1-tt0 if(three_body_h_tc) then - !call wall_time(tt0) - !PROVIDE fock_a_tot_3e_bi_orth - !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth PROVIDE fock_3e_uhf_mo_a Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a - !call wall_time(tt1) - !print*, ' 3-e term:', tt1-tt0 endif else @@ -241,19 +225,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] if(bi_ortho) then - !allocate(tmp(ao_num,ao_num)) - !tmp = Fock_matrix_tc_ao_beta - !if(three_body_h_tc) then - ! tmp += fock_3e_uhf_ao_b - !endif - !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1)) - !deallocate(tmp) - call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) if(three_body_h_tc) then - !PROVIDE fock_b_tot_3e_bi_orth - !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth PROVIDE fock_3e_uhf_mo_b Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b endif diff --git a/src/tc_scf/fock_three_bi_ortho.irp.f b/src/tc_scf/fock_three_bi_ortho.irp.f index 5d2f199c..8475c387 100644 --- a/src/tc_scf/fock_three_bi_ortho.irp.f +++ b/src/tc_scf/fock_three_bi_ortho.irp.f @@ -34,7 +34,7 @@ BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] enddo !call wall_time(t1) - !print*, ' Wall time for fock_a_tot_3e_bi_orth =', t1-t0 + !print*, ' Wall time for fock_a_tot_3e_bi_orth =', t1 - t0 END_PROVIDER diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index d7780497..b925c9df 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -57,6 +57,7 @@ program test_ints ! call test_old_ints call test_fock_3e_uhf_mo_cs() + call test_fock_3e_uhf_mo_a() end @@ -1109,6 +1110,13 @@ subroutine test_fock_3e_uhf_mo_cs() double precision :: I_old, I_new double precision :: diff_tot, diff, thr_ih, norm +! double precision :: t0, t1 +! print*, ' Providing fock_a_tot_3e_bi_orth ...' +! call wall_time(t0) +! PROVIDE fock_a_tot_3e_bi_orth +! call wall_time(t1) +! print*, ' Wall time for fock_a_tot_3e_bi_orth =', t1 - t0 + PROVIDE fock_3e_uhf_mo_cs fock_3e_uhf_mo_cs_old thr_ih = 1d-10 @@ -1123,7 +1131,7 @@ subroutine test_fock_3e_uhf_mo_cs() diff = dabs(I_old - I_new) if(diff .gt. thr_ih) then - print *, ' problem on ', j, i + print *, ' problem in fock_3e_uhf_mo_cs on ', j, i print *, ' old value = ', I_old print *, ' new value = ', I_new stop @@ -1139,4 +1147,83 @@ subroutine test_fock_3e_uhf_mo_cs() return end subroutine test_fock_3e_uhf_mo_cs +! --- + +subroutine test_fock_3e_uhf_mo_a() + + implicit none + integer :: i, j + double precision :: I_old, I_new + double precision :: diff_tot, diff, thr_ih, norm + + PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_a_old + + thr_ih = 1d-10 + norm = 0.d0 + diff_tot = 0.d0 + + do i = 1, mo_num + do j = 1, mo_num + + I_old = fock_3e_uhf_mo_a_old(j,i) + I_new = fock_3e_uhf_mo_a (j,i) + + diff = dabs(I_old - I_new) + if(diff .gt. thr_ih) then + print *, ' problem in fock_3e_uhf_mo_a on ', j, i + print *, ' old value = ', I_old + print *, ' new value = ', I_new + stop + endif + + norm += dabs(I_old) + diff_tot += diff + enddo + enddo + + print *, ' diff tot (%) = ', 100.d0 * diff_tot / norm + + return +end subroutine test_fock_3e_uhf_mo_a + +! --- + +subroutine test_fock_3e_uhf_mo_b() + + implicit none + integer :: i, j + double precision :: I_old, I_new + double precision :: diff_tot, diff, thr_ih, norm + + PROVIDE fock_3e_uhf_mo_b fock_3e_uhf_mo_b_old + + thr_ih = 1d-10 + norm = 0.d0 + diff_tot = 0.d0 + + do i = 1, mo_num + do j = 1, mo_num + + I_old = fock_3e_uhf_mo_b_old(j,i) + I_new = fock_3e_uhf_mo_b (j,i) + + diff = dabs(I_old - I_new) + if(diff .gt. thr_ih) then + print *, ' problem in fock_3e_uhf_mo_b on ', j, i + print *, ' old value = ', I_old + print *, ' new value = ', I_new + stop + endif + + norm += dabs(I_old) + diff_tot += diff + enddo + enddo + + print *, ' diff tot (%) = ', 100.d0 * diff_tot / norm + + return +end subroutine test_fock_3e_uhf_mo_b + +! ---