From 207e52d2202bca52eebc0ba8f394150a9bd98b48 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Wed, 26 Apr 2023 08:52:06 +0200 Subject: [PATCH] fixed bug in TC integrals --- src/bi_ort_ints/one_e_bi_ort.irp.f | 17 +- src/bi_ort_ints/semi_num_ints_mo.irp.f | 86 ++++--- src/bi_ort_ints/three_body_ints_bi_ort.irp.f | 21 +- src/bi_ort_ints/total_twoe_pot.irp.f | 3 + src/bi_ortho_mos/mos_rl.irp.f | 4 + src/non_h_ints_mu/grad_squared.irp.f | 8 +- src/non_h_ints_mu/new_grad_tc.irp.f | 139 ++++++----- src/non_h_ints_mu/tc_integ.irp.f | 246 +++++++++---------- src/non_h_ints_mu/total_tc_int.irp.f | 30 ++- src/tc_scf/diis_tcscf.irp.f | 8 + src/tc_scf/fock_3e_bi_ortho_uhf.irp.f | 14 +- src/tc_scf/fock_tc.irp.f | 14 +- src/tc_scf/fock_tc_mo_tot.irp.f | 218 ++++++++-------- src/tc_scf/fock_three_bi_ortho.irp.f | 8 +- src/tc_scf/tc_scf_energy.irp.f | 1 + 15 files changed, 453 insertions(+), 364 deletions(-) diff --git a/src/bi_ort_ints/one_e_bi_ort.irp.f b/src/bi_ort_ints/one_e_bi_ort.irp.f index 7f89899b..b0455570 100644 --- a/src/bi_ort_ints/one_e_bi_ort.irp.f +++ b/src/bi_ort_ints/one_e_bi_ort.irp.f @@ -11,14 +11,17 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)] provide j1b_type if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then + + print *, ' do things properly !' + stop - do i = 1, ao_num - do j = 1, ao_num - ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) & - + j1b_gauss_hermII (j,i) & - + j1b_gauss_nonherm(j,i) ) - enddo - enddo + !do i = 1, ao_num + ! do j = 1, ao_num + ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) & + ! + j1b_gauss_hermII (j,i) & + ! + j1b_gauss_nonherm(j,i) ) + ! enddo + !enddo endif diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f index 4694a998..0d727785 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -110,27 +110,36 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, print *, ' providing int2_grad1_u12_ao_transp ...' call wall_time(wall0) - if(test_cycle_tc)then - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, ao_num - int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1) - int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2) - int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3) - enddo - enddo - enddo + if(test_cycle_tc) then + + PROVIDE int2_grad1_u12_ao_test + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1) + int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2) + int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3) + enddo + enddo + enddo + else - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, ao_num - int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1) - int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2) - int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3) - enddo - enddo - enddo + + PROVIDE int2_grad1_u12_ao + + do ipoint = 1, n_points_final_grid + do i = 1, ao_num + do j = 1, ao_num + int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1) + int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2) + int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3) + enddo + enddo + enddo + endif + call wall_time(wall1) print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0 @@ -144,9 +153,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, integer :: ipoint double precision :: wall0, wall1 - !print *, ' providing int2_grad1_u12_bimo_transp' + PROVIDE mo_l_coef mo_r_coef + PROVIDE int2_grad1_u12_ao_transp + + !print *, ' providing int2_grad1_u12_bimo_transp' + !call wall_time(wall0) - call wall_time(wall0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint) & @@ -163,25 +175,31 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, !$OMP END DO !$OMP END PARALLEL - call wall_time(wall1) + !call wall_time(wall1) !print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0 END_PROVIDER ! --- -BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3, mo_num, mo_num )] - implicit none - integer :: i, j, ipoint - do ipoint = 1, n_points_final_grid - do i = 1, mo_num - do j = 1, mo_num - int2_grad1_u12_bimo_t(ipoint,1,j,i) = int2_grad1_u12_bimo_transp(j,i,1,ipoint) - int2_grad1_u12_bimo_t(ipoint,2,j,i) = int2_grad1_u12_bimo_transp(j,i,2,ipoint) - int2_grad1_u12_bimo_t(ipoint,3,j,i) = int2_grad1_u12_bimo_transp(j,i,3,ipoint) - enddo - enddo - enddo +BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)] + + implicit none + integer :: i, j, ipoint + + PROVIDE mo_l_coef mo_r_coef + PROVIDE int2_grad1_u12_bimo_transp + + do ipoint = 1, n_points_final_grid + do i = 1, mo_num + do j = 1, mo_num + int2_grad1_u12_bimo_t(ipoint,1,j,i) = int2_grad1_u12_bimo_transp(j,i,1,ipoint) + int2_grad1_u12_bimo_t(ipoint,2,j,i) = int2_grad1_u12_bimo_transp(j,i,2,ipoint) + int2_grad1_u12_bimo_t(ipoint,3,j,i) = int2_grad1_u12_bimo_transp(j,i,3,ipoint) + enddo + enddo + enddo + END_PROVIDER ! --- diff --git a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f index 48fa84f7..e8b56307 100644 --- a/src/bi_ort_ints/three_body_ints_bi_ort.irp.f +++ b/src/bi_ort_ints/three_body_ints_bi_ort.irp.f @@ -81,21 +81,24 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) integer :: ipoint double precision :: weight + PROVIDE mo_l_coef mo_r_coef + PROVIDE int2_grad1_u12_bimo_t + integral = 0.d0 do ipoint = 1, n_points_final_grid weight = final_weight_at_r_vector(ipoint) - integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & - * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) & - + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) & + integral += weight * mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) & + * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,l,j) & + + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,l,j) & + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,l,j) ) - integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & - * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & - + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + integral += weight * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,j) & + * ( int2_grad1_u12_bimo_t(ipoint,1,n,m) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & + + int2_grad1_u12_bimo_t(ipoint,2,n,m) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + int2_grad1_u12_bimo_t(ipoint,3,n,m) * int2_grad1_u12_bimo_t(ipoint,3,k,i) ) - integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & - * ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & - + int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + integral += weight * mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,m) & + * ( int2_grad1_u12_bimo_t(ipoint,1,l,j) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & + + int2_grad1_u12_bimo_t(ipoint,2,l,j) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + int2_grad1_u12_bimo_t(ipoint,3,l,j) * int2_grad1_u12_bimo_t(ipoint,3,k,i) ) enddo diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index bdebe890..721ea0c8 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -20,6 +20,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, a enddo END_PROVIDER + ! --- BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ] @@ -66,6 +67,8 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n else + PROVIDE ao_tc_int_chemist + do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num diff --git a/src/bi_ortho_mos/mos_rl.irp.f b/src/bi_ortho_mos/mos_rl.irp.f index d51999fc..c69309d1 100644 --- a/src/bi_ortho_mos/mos_rl.irp.f +++ b/src/bi_ortho_mos/mos_rl.irp.f @@ -17,6 +17,8 @@ subroutine ao_to_mo_bi_ortho(A_ao, LDA_ao, A_mo, LDA_mo) double precision, intent(out) :: A_mo(LDA_mo,mo_num) double precision, allocatable :: T(:,:) + PROVIDE mo_l_coef mo_r_coef + allocate ( T(ao_num,mo_num) ) !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T @@ -54,6 +56,8 @@ subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao) double precision, intent(out) :: A_ao(LDA_ao,ao_num) double precision, allocatable :: tmp_1(:,:), tmp_2(:,:) + PROVIDE mo_l_coef mo_r_coef + ! ao_overlap x mo_r_coef allocate( tmp_1(ao_num,mo_num) ) call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 & diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index 8b801f9d..ea4bd36c 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -352,7 +352,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao BEGIN_DOC ! - ! tc_grad_square_ao(k,i,l,j) = 1/2 + ! tc_grad_square_ao(k,i,l,j) = -1/2 ! END_DOC @@ -373,6 +373,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao else + PROVIDE int2_grad1_u12_square_ao + allocate(b_mat(n_points_final_grid,ao_num,ao_num)) b_mat = 0.d0 @@ -392,8 +394,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao !$OMP END PARALLEL tc_grad_square_ao = 0.d0 - call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & - , int2_grad1_u12_squared_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & , 0.d0, tc_grad_square_ao, ao_num*ao_num) deallocate(b_mat) diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 064eb9f1..24e7e743 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -1,68 +1,68 @@ ! --- -BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)] - - BEGIN_DOC - ! - ! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1) - ! - ! where r1 = r(ipoint) - ! - ! if J(r1,r2) = u12: - ! - ! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1) - ! = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ] - ! = -int2_grad1_u12_ao(i,j,ipoint,:) - ! - ! if J(r1,r2) = u12 x v1 x v2 - ! - ! int1_grad2_u12_ao(:,i,j,ipoint) = v2 x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ] - ! - \grad_2 v2 x [ \int dr1 u12 v1 \phi_i(r1) \phi_j(r1) ] - ! = -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) - ! + 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) - ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) - ! - ! - END_DOC - - implicit none - integer :: ipoint, i, j - double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 - - PROVIDE j1b_type - - if(j1b_type .eq. 3) then - - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) - - tmp0 = 0.5d0 * v_1b(ipoint) - tmp_x = v_1b_grad(1,ipoint) - tmp_y = v_1b_grad(2,ipoint) - tmp_z = v_1b_grad(3,ipoint) - - do j = 1, ao_num - do i = 1, ao_num - - tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) - tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) - - int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x - int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y - int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z - enddo - enddo - enddo - - else - - int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao - - endif - -END_PROVIDER +!BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)] +! +! BEGIN_DOC +! ! +! ! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1) +! ! +! ! where r1 = r(ipoint) +! ! +! ! if J(r1,r2) = u12: +! ! +! ! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1) +! ! = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ] +! ! = -int2_grad1_u12_ao(i,j,ipoint,:) +! ! +! ! if J(r1,r2) = u12 x v1 x v2 +! ! +! ! int1_grad2_u12_ao(:,i,j,ipoint) = v2 x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ] +! ! - \grad_2 v2 x [ \int dr1 u12 v1 \phi_i(r1) \phi_j(r1) ] +! ! = -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) +! ! + 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) +! ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) +! ! +! ! +! END_DOC +! +! implicit none +! integer :: ipoint, i, j +! double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 +! +! PROVIDE j1b_type +! +! if(j1b_type .eq. 3) then +! +! do ipoint = 1, n_points_final_grid +! x = final_grid_points(1,ipoint) +! y = final_grid_points(2,ipoint) +! z = final_grid_points(3,ipoint) +! +! tmp0 = 0.5d0 * v_1b(ipoint) +! tmp_x = v_1b_grad(1,ipoint) +! tmp_y = v_1b_grad(2,ipoint) +! tmp_z = v_1b_grad(3,ipoint) +! +! do j = 1, ao_num +! do i = 1, ao_num +! +! tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) +! tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) +! +! int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x +! int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y +! int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z +! enddo +! enddo +! enddo +! +! else +! +! int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao +! +! endif +! +!END_PROVIDER ! --- @@ -192,7 +192,10 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ! ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) . \grad_1 | ij > ! - ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ! = -1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 (-1) \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) + ! + ! -1 in \int dr2 ! ! This is obtained by integration by parts. ! @@ -207,8 +210,6 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, print*, ' providing tc_grad_and_lapl_ao ...' call wall_time(time0) - PROVIDE int2_grad1_u12_ao - if(read_tc_integ) then open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao', action="read") @@ -217,6 +218,8 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, else + PROVIDE int2_grad1_u12_ao + allocate(b_mat(n_points_final_grid,ao_num,ao_num,3)) b_mat = 0.d0 @@ -247,10 +250,10 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, do m = 1, 3 call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & , int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & - , 0.d0, tc_grad_and_lapl_ao, ao_num*ao_num) + , 1.d0, tc_grad_and_lapl_ao, ao_num*ao_num) enddo deallocate(b_mat) - + call sum_A_At(tc_grad_and_lapl_ao(1,1,1,1), ao_num*ao_num) ! !$OMP PARALLEL & ! !$OMP DEFAULT (NONE) & diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f index 865522e0..41209a36 100644 --- a/src/non_h_ints_mu/tc_integ.irp.f +++ b/src/non_h_ints_mu/tc_integ.irp.f @@ -1,8 +1,7 @@ ! --- - BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)] -&BEGIN_PROVIDER [double precision, int2_grad1_u12_squared_ao, (ao_num, ao_num, n_points_final_grid)] +BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_DOC ! @@ -23,10 +22,6 @@ ! - 0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) ! - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint) ! - ! - ! - ! int2_grad1_u12_squared_ao = int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 - ! END_DOC implicit none @@ -34,25 +29,22 @@ double precision :: time0, time1 double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 - print*, ' providing int2_grad1_u12_ao & int2_grad1_u12_squared_ao ...' + print*, ' providing int2_grad1_u12_ao ...' + call wall_time(time0) PROVIDE j1b_type if(read_tc_integ) then - call wall_time(time0) open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read") read(11) int2_grad1_u12_ao - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_ao =', time1-time0 endif if(j1b_type .eq. 3) then - ! --- - if(.not.read_tc_integ) then - call wall_time(time0) - PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b + + PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b + int2_grad1_u12_ao = 0.d0 ! TODO OPENMP do ipoint = 1, n_points_final_grid @@ -73,42 +65,14 @@ enddo enddo enddo - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_ao =', time1-time0 + endif - ! --- - - call wall_time(time0) - PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 - int2_grad1_u12_squared_ao = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, ipoint) & - !$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - do j = 1, ao_num - do i = 1, ao_num - int2_grad1_u12_squared_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_squared_ao =', time1-time0 - - ! --- - elseif(j1b_type .ge. 100) then - ! --- - - call wall_time(time0) PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra - PROVIDE grad1_u12_squared_num PROVIDE grad1_u12_num + double precision, allocatable :: tmp(:,:,:) allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) tmp = 0.d0 @@ -130,34 +94,14 @@ if(.not.read_tc_integ) then int2_grad1_u12_ao = 0.d0 do m = 1, 3 - call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, 1.d0 & + !call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 & + ! this work also because of the symmetry in K(1,2) and sign compensation in L(1,2,3) + call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -1.d0 & , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid & , 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num) enddo - !! DEBUG - !PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b - !int2_grad1_u12_ao = 0.d0 - !do ipoint = 1, n_points_final_grid - ! x = final_grid_points(1,ipoint) - ! y = final_grid_points(2,ipoint) - ! z = final_grid_points(3,ipoint) - ! tmp0 = 0.5d0 * v_1b(ipoint) - ! tmp_x = v_1b_grad(1,ipoint) - ! tmp_y = v_1b_grad(2,ipoint) - ! tmp_z = v_1b_grad(3,ipoint) - ! do j = 1, ao_num - ! do i = 1, ao_num - ! tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) - ! tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) - ! int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x - ! int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y - ! int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z - ! enddo - ! enddo - !enddo - - ! these dgemm are equivalen to + !! these dgemm are equivalen to !!$OMP PARALLEL & !!$OMP DEFAULT (NONE) & !!$OMP PRIVATE (j, i, ipoint, jpoint, w) & @@ -169,7 +113,9 @@ ! do j = 1, ao_num ! do i = 1, ao_num ! do jpoint = 1, n_points_extra_final_grid - ! w = tmp(jpoint,i,j) + ! w = -tmp(jpoint,i,j) + ! !w = tmp(jpoint,i,j) this work also because of the symmetry in K(1,2) + ! ! and sign compensation in L(1,2,3) ! int2_grad1_u12_ao(i,j,ipoint,1) += w * grad1_u12_num(jpoint,ipoint,1) ! int2_grad1_u12_ao(i,j,ipoint,2) += w * grad1_u12_num(jpoint,ipoint,2) ! int2_grad1_u12_ao(i,j,ipoint,3) += w * grad1_u12_num(jpoint,ipoint,3) @@ -179,67 +125,15 @@ !enddo !!$OMP END DO !!$OMP END PARALLEL - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_ao =', time1-time0 endif - ! --- - - call wall_time(time0) - int2_grad1_u12_squared_ao = 0.d0 - call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 & - , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid & - , 0.d0, int2_grad1_u12_squared_ao(1,1,1), ao_num*ao_num) - - !! DEBUG - !PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 - !int2_grad1_u12_squared_ao = 0.d0 - !!$OMP PARALLEL & - !!$OMP DEFAULT (NONE) & - !!$OMP PRIVATE (i, j, ipoint) & - !!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) - !!$OMP DO SCHEDULE (static) - !do ipoint = 1, n_points_final_grid - ! do j = 1, ao_num - ! do i = 1, ao_num - ! int2_grad1_u12_squared_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) - ! enddo - ! enddo - !enddo - !!$OMP END DO - !!$OMP END PARALLEL - !call wall_time(time1) - - !! this dgemm is equivalen to - !!$OMP PARALLEL & - !!$OMP DEFAULT (NONE) & - !!$OMP PRIVATE (i, j, ipoint, jpoint, w) & - !!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, & - !!$OMP n_points_extra_final_grid, final_weight_at_r_vector_extra, & - !!$OMP aos_in_r_array_extra_transp, grad1_u12_squared_num, tmp) - !!$OMP DO SCHEDULE (static) - !do ipoint = 1, n_points_final_grid - ! do j = 1, ao_num - ! do i = 1, ao_num - ! do jpoint = 1, n_points_extra_final_grid - ! w = 0.5d0 * tmp(jpoint,i,j) - ! int2_grad1_u12_squared_ao(i,j,ipoint) += w * grad1_u12_squared_num(jpoint,ipoint) - ! enddo - ! enddo - ! enddo - !enddo - !!$OMP END DO - !!$OMP END PARALLEL - call wall_time(time1) - print*, ' wall time for int2_grad1_u12_squared_ao =', time1-time0 - - ! --- - deallocate(tmp) else + print *, ' j1b_type = ', j1b_type, 'not implemented yet' stop + endif if(write_tc_integ.and.mpi_master) then @@ -250,6 +144,112 @@ call ezfio_set_tc_keywords_io_tc_integ('Read') endif + call wall_time(time1) + print*, ' wall time for int2_grad1_u12_ao =', time1-time0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_points_final_grid)] + + BEGIN_DOC + ! + ! int2_grad1_u12_square_ao = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2 + ! + END_DOC + + implicit none + integer :: ipoint, i, j, m, jpoint + double precision :: time0, time1 + double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 + + print*, ' providing int2_grad1_u12_square_ao ...' + call wall_time(time0) + + PROVIDE j1b_type + + if(j1b_type .eq. 3) then + + PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12 + + int2_grad1_u12_square_ao = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, ipoint) & + !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + elseif(j1b_type .ge. 100) then + + PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra + PROVIDE grad1_u12_squared_num + + double precision, allocatable :: tmp(:,:,:) + allocate(tmp(n_points_extra_final_grid,ao_num,ao_num)) + tmp = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, i, jpoint) & + !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp) + !$OMP DO SCHEDULE (static) + do j = 1, ao_num + do i = 1, ao_num + do jpoint = 1, n_points_extra_final_grid + tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + int2_grad1_u12_square_ao = 0.d0 + call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 & + , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid & + , 0.d0, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num) + + !! this dgemm is equivalen to + !!$OMP PARALLEL & + !!$OMP DEFAULT (NONE) & + !!$OMP PRIVATE (i, j, ipoint, jpoint, w) & + !!$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, & + !!$OMP n_points_extra_final_grid, final_weight_at_r_vector_extra, & + !!$OMP aos_in_r_array_extra_transp, grad1_u12_squared_num, tmp) + !!$OMP DO SCHEDULE (static) + !do ipoint = 1, n_points_final_grid + ! do j = 1, ao_num + ! do i = 1, ao_num + ! do jpoint = 1, n_points_extra_final_grid + ! w = -0.5d0 * tmp(jpoint,i,j) + ! int2_grad1_u12_square_ao(i,j,ipoint) += w * grad1_u12_squared_num(jpoint,ipoint) + ! enddo + ! enddo + ! enddo + !enddo + !!$OMP END DO + !!$OMP END PARALLEL + + deallocate(tmp) + + else + + print *, ' j1b_type = ', j1b_type, 'not implemented yet' + stop + + endif + + call wall_time(time1) + print*, ' wall time for int2_grad1_u12_square_ao =', time1-time0 + END_PROVIDER ! --- diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f index 4f8dc74d..a60c99da 100644 --- a/src/non_h_ints_mu/total_tc_int.irp.f +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -48,9 +48,14 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao print *, ' providing ao_tc_int_chemist ...' call wall_time(wall0) - if(test_cycle_tc)then + if(test_cycle_tc) then + ao_tc_int_chemist = ao_tc_int_chemist_test + else + + PROVIDE tc_grad_square_ao tc_grad_and_lapl_ao ao_two_e_coul + do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num @@ -68,27 +73,34 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao END_PROVIDER -BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)] ! --- + +BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)] + implicit none integer :: i, j, k, l double precision :: wall1, wall0 + print *, ' providing ao_tc_int_chemist_no_cycle ...' call wall_time(wall0) - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) -! ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j) - enddo + + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) + !ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j) enddo enddo enddo + enddo + call wall_time(wall1) print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0 END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)] implicit none diff --git a/src/tc_scf/diis_tcscf.irp.f b/src/tc_scf/diis_tcscf.irp.f index 0b08f784..4ec70de5 100644 --- a/src/tc_scf/diis_tcscf.irp.f +++ b/src/tc_scf/diis_tcscf.irp.f @@ -92,17 +92,22 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)] allocate(F(ao_num,ao_num)) if(var_tc) then + do i = 1, ao_num do j = 1, ao_num F(j,i) = Fock_matrix_vartc_ao_tot(j,i) enddo enddo + else + + PROVIDE Fock_matrix_tc_ao_tot do i = 1, ao_num do j = 1, ao_num F(j,i) = Fock_matrix_tc_ao_tot(j,i) enddo enddo + endif allocate(tmp(ao_num,ao_num)) @@ -139,6 +144,9 @@ BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)] implicit none + PROVIDE mo_r_coef mo_l_coef + PROVIDE FQS_SQF_ao + call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) & , FQS_SQF_mo, size(FQS_SQF_mo, 1) ) 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 d8b962d7..a67a3705 100644 --- a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f +++ b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -50,19 +50,23 @@ END_PROVIDER BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] BEGIN_DOC -! ALPHA part of the Fock matrix from three-electron terms -! -! WARNING :: non hermitian if bi-ortho MOS used + ! + ! ALPHA 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, o double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia double precision :: ti, tf PROVIDE mo_l_coef mo_r_coef + PROVIDE fock_3e_uhf_mo_cs !print *, ' PROVIDING fock_3e_uhf_mo_a ...' - call wall_time(ti) + !call wall_time(ti) o = elec_beta_num + 1 @@ -142,7 +146,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] enddo enddo - call wall_time(tf) + !call wall_time(tf) !print *, ' total 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 1d651c4e..52eeb694 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -18,6 +18,8 @@ double precision :: density, density_a, density_b double precision :: t0, t1 + PROVIDE ao_two_e_tc_tot + !print*, ' providing two_e_tc_non_hermit_integral_seq ...' !call wall_time(t0) @@ -80,6 +82,10 @@ END_PROVIDER double precision :: t0, t1 double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + PROVIDE ao_two_e_tc_tot + PROVIDE mo_l_coef mo_r_coef + PROVIDE TCSCF_density_matrix_ao_alpha TCSCF_density_matrix_ao_beta + !print*, ' providing two_e_tc_non_hermit_integral ...' !call wall_time(t0) @@ -142,7 +148,7 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)] implicit none - Fock_matrix_tc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_alpha + Fock_matrix_tc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_tc_non_hermit_integral_alpha END_PROVIDER @@ -181,14 +187,17 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] !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 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) ) if(three_body_h_tc) then !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 endif else + call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) @@ -276,6 +285,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ] implicit none + PROVIDE mo_l_coef mo_r_coef + PROVIDE Fock_matrix_tc_mo_tot + call mo_to_ao_bi_ortho( Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) & , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) ) diff --git a/src/tc_scf/fock_tc_mo_tot.irp.f b/src/tc_scf/fock_tc_mo_tot.irp.f index a03a0624..78f4c9b0 100644 --- a/src/tc_scf/fock_tc_mo_tot.irp.f +++ b/src/tc_scf/fock_tc_mo_tot.irp.f @@ -1,107 +1,120 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ] &BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)] - implicit none - BEGIN_DOC - ! TC-Fock matrix on the MO basis. WARNING !!! NON HERMITIAN !!! - ! For open shells, the ROHF Fock Matrix is :: - ! - ! | F-K | F + K/2 | F | - ! |---------------------------------| - ! | F + K/2 | F | F - K/2 | - ! |---------------------------------| - ! | F | F - K/2 | F + K | - ! - ! - ! F = 1/2 (Fa + Fb) - ! - ! K = Fb - Fa - ! - END_DOC - integer :: i,j,n - if (elec_alpha_num == elec_beta_num) then - Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha - else - do j=1,elec_beta_num - ! F-K - do i=1,elec_beta_num !CC - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + BEGIN_DOC + ! TC-Fock matrix on the MO basis. WARNING !!! NON HERMITIAN !!! + ! For open shells, the ROHF Fock Matrix is :: + ! + ! | F-K | F + K/2 | F | + ! |---------------------------------| + ! | F + K/2 | F | F - K/2 | + ! |---------------------------------| + ! | F | F - K/2 | F + K | + ! + ! + ! F = 1/2 (Fa + Fb) + ! + ! K = Fb - Fa + ! + END_DOC + + implicit none + integer :: i, j, n + + if(elec_alpha_num == elec_beta_num) then + + PROVIDE Fock_matrix_tc_mo_alpha + + Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha + + else + + PROVIDE Fock_matrix_tc_mo_beta Fock_matrix_tc_mo_alpha + + do j = 1, elec_beta_num + ! F-K + do i = 1, elec_beta_num !CC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& - (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - enddo - ! F+K/2 - do i=elec_beta_num+1,elec_alpha_num !CA - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + enddo + ! F+K/2 + do i = elec_beta_num+1, elec_alpha_num !CA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - enddo - ! F - do i=elec_alpha_num+1, mo_num !CV - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) - enddo - enddo + enddo + ! F + do i = elec_alpha_num+1, mo_num !CV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + enddo - do j=elec_beta_num+1,elec_alpha_num - ! F+K/2 - do i=1,elec_beta_num !AC - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + do j = elec_beta_num+1, elec_alpha_num + ! F+K/2 + do i = 1, elec_beta_num !AC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - enddo - ! F - do i=elec_beta_num+1,elec_alpha_num !AA - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) - enddo - ! F-K/2 - do i=elec_alpha_num+1, mo_num !AV - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + enddo + ! F + do i = elec_beta_num+1, elec_alpha_num !AA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + ! F-K/2 + do i = elec_alpha_num+1, mo_num !AV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - enddo - enddo + enddo + enddo - do j=elec_alpha_num+1, mo_num - ! F - do i=1,elec_beta_num !VC - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) - enddo - ! F-K/2 - do i=elec_beta_num+1,elec_alpha_num !VA - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + do j = elec_alpha_num+1, mo_num + ! F + do i = 1, elec_beta_num !VC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + ! F-K/2 + do i = elec_beta_num+1, elec_alpha_num !VA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - enddo - ! F+K - do i=elec_alpha_num+1,mo_num !VV - Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) & + enddo + ! F+K + do i = elec_alpha_num+1, mo_num !VV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) & + (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) - enddo - enddo - if(three_body_h_tc)then + enddo + enddo + + if(three_body_h_tc) then + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + ! C-O do j = 1, elec_beta_num - do i = elec_beta_num+1, elec_alpha_num - Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) - Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) - enddo + do i = elec_beta_num+1, elec_alpha_num + Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo enddo ! C-V do j = 1, elec_beta_num - do i = elec_alpha_num+1, mo_num - Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) - Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) - enddo + do i = elec_alpha_num+1, mo_num + Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo enddo ! O-V do j = elec_beta_num+1, elec_alpha_num - do i = elec_alpha_num+1, mo_num - Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) - Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) - enddo + do i = elec_alpha_num+1, mo_num + Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo enddo - endif + endif - endif + endif - do i = 1, mo_num - Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i) - enddo + do i = 1, mo_num + Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i) + enddo if(frozen_orb_scf)then @@ -116,28 +129,29 @@ enddo endif - if(no_oa_or_av_opt)then - do i = 1, n_act_orb - iorb = list_act(i) - do j = 1, n_inact_orb - jorb = list_inact(j) - Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 - Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 - enddo - do j = 1, n_virt_orb - jorb = list_virt(j) - Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 - Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 - enddo - do j = 1, n_core_orb - jorb = list_core(j) - Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 - Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 - enddo - enddo - endif + if(no_oa_or_av_opt)then + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + if(.not.bi_ortho .and. three_body_h_tc)then - Fock_matrix_tc_mo_tot += fock_3_mat + Fock_matrix_tc_mo_tot += fock_3_mat endif END_PROVIDER diff --git a/src/tc_scf/fock_three_bi_ortho.irp.f b/src/tc_scf/fock_three_bi_ortho.irp.f index 7c776ce5..cca4b5aa 100644 --- a/src/tc_scf/fock_three_bi_ortho.irp.f +++ b/src/tc_scf/fock_three_bi_ortho.irp.f @@ -4,14 +4,16 @@ BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] BEGIN_DOC -! Alpha part of the Fock matrix from three-electron terms -! -! WARNING :: non hermitian if bi-ortho MOS used + ! Alpha part of the Fock matrix from three-electron terms + ! + ! WARNING :: non hermitian if bi-ortho MOS used END_DOC + implicit none integer :: i, a PROVIDE mo_l_coef mo_r_coef + PROVIDE fock_cs_3e_bi_orth fock_a_tmp1_bi_ortho fock_a_tmp2_bi_ortho fock_a_tot_3e_bi_orth = 0.d0 diff --git a/src/tc_scf/tc_scf_energy.irp.f b/src/tc_scf/tc_scf_energy.irp.f index c3de0322..1fb09828 100644 --- a/src/tc_scf/tc_scf_energy.irp.f +++ b/src/tc_scf/tc_scf_energy.irp.f @@ -11,6 +11,7 @@ integer :: i, j PROVIDE mo_l_coef mo_r_coef + PROVIDE two_e_tc_non_hermit_integral_alpha two_e_tc_non_hermit_integral_beta TC_HF_energy = nuclear_repulsion TC_HF_one_e_energy = 0.d0