From 76d502bd353feb37541a282d0583e665c2850677 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Wed, 21 Dec 2022 00:40:24 +0100 Subject: [PATCH] added UHF Fock matrices --- src/bi_ort_ints/semi_num_ints_mo.irp.f | 21 + src/bi_ort_ints/three_body_ints_bi_ort.irp.f | 43 ++- src/hartree_fock/fock_matrix_hf.irp.f | 25 +- src/non_h_ints_mu/new_grad_tc.irp.f | 106 +++++- src/tc_keywords/EZFIO.cfg | 2 +- src/tc_scf/fock_3e_bi_ortho_uhf.irp.f | 377 ++++++++++++++++++ src/tc_scf/fock_tc.irp.f | 44 ++- src/tc_scf/fock_three_bi_ortho_new_new.irp.f | 380 +++++++++++-------- src/tc_scf/tc_scf.irp.f | 69 +++- src/tc_scf/test_int.irp.f | 185 ++++++++- 10 files changed, 1058 insertions(+), 194 deletions(-) create mode 100644 src/tc_scf/fock_3e_bi_ortho_uhf.irp.f 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 4762c25e..89098676 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -170,6 +170,27 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,3 enddo END_PROVIDER +! --- + +BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_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_ao_t(ipoint,1,j,i) = int2_grad1_u12_ao(1,j,i,ipoint) + int2_grad1_u12_ao_t(ipoint,2,j,i) = int2_grad1_u12_ao(2,j,i,ipoint) + int2_grad1_u12_ao_t(ipoint,3,j,i) = int2_grad1_u12_ao(3,j,i,ipoint) + enddo + enddo + enddo + +END_PROVIDER + +! --- + BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo, (3, mo_num, mo_num, n_points_final_grid)] BEGIN_DOC 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 c1c27f06..48fa84f7 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 @@ -15,7 +15,7 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n character*(128) :: name_file three_body_ints_bi_ort = 0.d0 - print*,'Providing the three_body_ints_bi_ort ...' + print *, ' Providing the three_body_ints_bi_ort ...' call wall_time(wall0) name_file = 'six_index_tensor' @@ -71,7 +71,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS ! END_DOC @@ -104,12 +104,11 @@ end subroutine give_integrals_3_body_bi_ort ! --- - subroutine give_integrals_3_body_bi_ort_old(n, l, k, m, j, i, integral) BEGIN_DOC ! - ! < n l k | -L | m j i > with a BI-ORTHONORMAL ORBITALS + ! < n l k | -L | m j i > with a BI-ORTHONORMAL MOLECULAR ORBITALS ! END_DOC @@ -170,3 +169,39 @@ end subroutine give_integrals_3_body_bi_ort_old ! --- +subroutine give_integrals_3_body_bi_ort_ao(n, l, k, m, j, i, integral) + + BEGIN_DOC + ! + ! < n l k | -L | m j i > with a BI-ORTHONORMAL ATOMIC ORBITALS + ! + END_DOC + + implicit none + integer, intent(in) :: n, l, k, m, j, i + double precision, intent(out) :: integral + integer :: ipoint + double precision :: weight + + integral = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + + integral += weight * aos_in_r_array_transp(ipoint,k) * aos_in_r_array_transp(ipoint,i) & + * ( int2_grad1_u12_ao_t(ipoint,1,n,m) * int2_grad1_u12_ao_t(ipoint,1,l,j) & + + int2_grad1_u12_ao_t(ipoint,2,n,m) * int2_grad1_u12_ao_t(ipoint,2,l,j) & + + int2_grad1_u12_ao_t(ipoint,3,n,m) * int2_grad1_u12_ao_t(ipoint,3,l,j) ) + integral += weight * aos_in_r_array_transp(ipoint,l) * aos_in_r_array_transp(ipoint,j) & + * ( int2_grad1_u12_ao_t(ipoint,1,n,m) * int2_grad1_u12_ao_t(ipoint,1,k,i) & + + int2_grad1_u12_ao_t(ipoint,2,n,m) * int2_grad1_u12_ao_t(ipoint,2,k,i) & + + int2_grad1_u12_ao_t(ipoint,3,n,m) * int2_grad1_u12_ao_t(ipoint,3,k,i) ) + integral += weight * aos_in_r_array_transp(ipoint,n) * aos_in_r_array_transp(ipoint,m) & + * ( int2_grad1_u12_ao_t(ipoint,1,l,j) * int2_grad1_u12_ao_t(ipoint,1,k,i) & + + int2_grad1_u12_ao_t(ipoint,2,l,j) * int2_grad1_u12_ao_t(ipoint,2,k,i) & + + int2_grad1_u12_ao_t(ipoint,3,l,j) * int2_grad1_u12_ao_t(ipoint,3,k,i) ) + + enddo + +end subroutine give_integrals_3_body_bi_ort_ao + +! --- diff --git a/src/hartree_fock/fock_matrix_hf.irp.f b/src/hartree_fock/fock_matrix_hf.irp.f index d7d8fa7d..cb698fbb 100644 --- a/src/hartree_fock/fock_matrix_hf.irp.f +++ b/src/hartree_fock/fock_matrix_hf.irp.f @@ -1,12 +1,27 @@ +! --- BEGIN_PROVIDER [ double precision, ao_two_e_integral_alpha, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta , (ao_num, ao_num) ] - use map_module - implicit none +&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta , (ao_num, ao_num) ] + BEGIN_DOC - ! Alpha and Beta Fock matrices in AO basis set + ! + ! 2-e part of alpha and beta Fock matrices (F^{a} & F^{b}) in AO basis set + ! + ! F^{a} = h + G^{a} + ! F^{b} = h + G^{b} + ! + ! where : + ! F^{a} = J^{a} + J^{b} - K^{a} ==> G_{ij}^{a} = \sum_{k,l} P_{kl} (kl|ij) - P_{kl}^{a} (ki|lj) + ! F^{b} = J^{a} + J^{b} - K^{b} ==> G_{ij}^{b} = \sum_{k,l} P_{kl} (kl|ij) - P_{kl}^{b} (ki|lj) + ! + ! and P_{kl} = P_{kl}^{a} + P_{kl}^{b} + ! END_DOC + use map_module + + implicit none + integer :: i,j,k,l,k1,r,s integer :: i0,j0,k0,l0 integer*8 :: p,q @@ -153,6 +168,8 @@ END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num, ao_num) ] &BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num, ao_num) ] implicit none 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 d34e629c..484e3850 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -82,11 +82,77 @@ 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_tmp_j1b(1,i,j,ipoint) - tmp2 * tmp_x + int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(2,i,j,ipoint) - tmp2 * tmp_y + int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_tmp_j1b(3,i,j,ipoint) - tmp2 * tmp_z + enddo + enddo + enddo + + else + + int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao + + endif + +END_PROVIDER + +! --- + BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, ao_num)] BEGIN_DOC ! - ! tc_grad_and_lapl_ao(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij > + ! 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) ! @@ -98,11 +164,14 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, integer :: ipoint, i, j, k, l double precision :: weight1, contrib_x, contrib_y, contrib_z, tmp_x, tmp_y, tmp_z double precision :: ao_k_r, ao_i_r, ao_i_dx, ao_i_dy, ao_i_dz + double precision :: ao_j_r, ao_l_r, ao_l_dx, ao_l_dy, ao_l_dz double precision, allocatable :: ac_mat(:,:,:,:) allocate(ac_mat(ao_num,ao_num,ao_num,ao_num)) ac_mat = 0.d0 + ! --- + do ipoint = 1, n_points_final_grid weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) @@ -132,12 +201,47 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, enddo enddo enddo + + ! --- + + !do ipoint = 1, n_points_final_grid + ! weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) + + ! do l = 1, ao_num + ! ao_l_r = weight1 * aos_in_r_array_transp (ipoint,l) + ! ao_l_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,1) + ! ao_l_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,2) + ! ao_l_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,3) + + ! do j = 1, ao_num + ! ao_j_r = aos_in_r_array_transp(ipoint,j) + + ! tmp_x = ao_j_r * ao_l_dx - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,1) + ! tmp_y = ao_j_r * ao_l_dy - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,2) + ! tmp_z = ao_j_r * ao_l_dz - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,3) + + ! do i = 1, ao_num + ! do k = 1, ao_num + + ! contrib_x = int2_grad1_u12_ao(1,k,i,ipoint) * tmp_x + ! contrib_y = int2_grad1_u12_ao(2,k,i,ipoint) * tmp_y + ! contrib_z = int2_grad1_u12_ao(3,k,i,ipoint) * tmp_z + + ! ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z + ! enddo + ! enddo + ! enddo + ! enddo + !enddo + + ! --- do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num do k = 1, ao_num tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + !tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) enddo enddo enddo diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index a206dfa9..45af9723 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -158,7 +158,7 @@ default: 0. type: character*(32) doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] interface: ezfio,provider,ocaml -default: DIIS +default: Simple [im_thresh_tcscf] type: Threshold diff --git a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f new file mode 100644 index 00000000..607140f9 --- /dev/null +++ b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -0,0 +1,377 @@ + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + !print *, ' PROVIDING fock_3e_uhf_mo_cs ...' + call wall_time(ti) + + fock_3e_uhf_mo_cs = 0.d0 + + do a = 1, mo_num + do b = 1, mo_num + + do j = 1, elec_beta_num + do i = 1, elec_beta_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) + + fock_3e_uhf_mo_cs(b,a) -= 0.5d0 * ( 4.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - 2.d0 * I_bij_aji & + - 2.d0 * I_bij_iaj & + - 2.d0 * I_bij_jia ) + + enddo + enddo + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] + + 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 + + !print *, ' PROVIDING fock_3e_uhf_mo_a ...' + call wall_time(ti) + + o = elec_beta_num + 1 + + fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs + + do a = 1, mo_num + do b = 1, mo_num + + ! --- + + do j = o, elec_alpha_num + do i = 1, elec_beta_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) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - I_bij_iaj & + - 2.d0 * I_bij_jia ) + + enddo + enddo + + ! --- + + do j = 1, elec_beta_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) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - 2.d0 * I_bij_iaj & + - I_bij_jia ) + + enddo + enddo + + ! --- + + 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) + + fock_3e_uhf_mo_a(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 + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] + + 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 + + !print *, ' PROVIDING fock_3e_uhf_mo_b ...' + call wall_time(ti) + + o = elec_beta_num + 1 + + fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs + + do a = 1, mo_num + do b = 1, mo_num + + ! --- + + do j = o, elec_alpha_num + do i = 1, elec_beta_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) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + - I_bij_aji & + - I_bij_iaj ) + + enddo + enddo + + ! --- + + do j = 1, elec_beta_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) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + - I_bij_aji & + - I_bij_jia ) + + enddo + enddo + + ! --- + + 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) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( I_bij_aij & + - I_bij_aji ) + + enddo + enddo + + ! --- + + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_b =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)] + + BEGIN_DOC + ! + ! Equations (B6) and (B7) + ! + ! g <--> gamma + ! d <--> delta + ! e <--> eta + ! k <--> kappa + ! + END_DOC + + implicit none + integer :: g, d, e, k, mu, nu + double precision :: dm_ge_a, dm_ge_b, dm_ge + double precision :: dm_dk_a, dm_dk_b, dm_dk + double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu + double precision :: ti, tf + + print *, ' PROVIDING fock_3e_uhf_ao_a ...' + call wall_time(ti) + + fock_3e_uhf_ao_a = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, & + !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & + !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a) + !$OMP DO + do g = 1, ao_num + do e = 1, ao_num + dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) + dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) + dm_ge = dm_ge_a + dm_ge_b + + do d = 1, ao_num + do k = 1, ao_num + dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) + dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) + dm_dk = dm_dk_a + dm_dk_b + + do mu = 1, ao_num + do nu = 1, ao_num + + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) + + fock_3e_uhf_ao_a(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & + + dm_ge_a * dm_dk_a * i_mugd_eknu & + + dm_ge_a * dm_dk_a * i_mugd_knue & + - dm_ge * dm_dk_a * i_mugd_kenu & + - dm_ge_a * dm_dk * i_mugd_enuk & + - dm_ge_a * dm_dk_a * i_mugd_nuke & + - dm_ge_b * dm_dk_b * i_mugd_nuke ) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)] + + BEGIN_DOC + ! + ! Equations (B6) and (B7) + ! + ! g <--> gamma + ! d <--> delta + ! e <--> eta + ! k <--> kappa + ! + END_DOC + + implicit none + integer :: g, d, e, k, mu, nu + double precision :: dm_ge_a, dm_ge_b, dm_ge + double precision :: dm_dk_a, dm_dk_b, dm_dk + double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu + double precision :: ti, tf + + print *, ' PROVIDING fock_3e_uhf_ao_b ...' + call wall_time(ti) + + fock_3e_uhf_ao_b = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, & + !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & + !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b) + !$OMP DO + do g = 1, ao_num + do e = 1, ao_num + dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) + dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) + dm_ge = dm_ge_a + dm_ge_b + + do d = 1, ao_num + do k = 1, ao_num + dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) + dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) + dm_dk = dm_dk_a + dm_dk_b + + do mu = 1, ao_num + do nu = 1, ao_num + + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) + + fock_3e_uhf_ao_b(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & + + dm_ge_b * dm_dk_b * i_mugd_eknu & + + dm_ge_b * dm_dk_b * i_mugd_knue & + - dm_ge * dm_dk_b * i_mugd_kenu & + - dm_ge_b * dm_dk * i_mugd_enuk & + - dm_ge_b * dm_dk_b * i_mugd_nuke & + - dm_ge_a * dm_dk_a * i_mugd_nuke ) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti + +END_PROVIDER + +! --- + diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index c3642a7e..5981791c 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -31,13 +31,22 @@ density_b = TCSCF_density_matrix_ao_beta (l,j) density = density_a + density_b + !! rho(l,j) * < k l| T | i j> + !two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(l,j,k,i) + !! rho(l,j) * < k l| T | i j> + !two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(l,j,k,i) + !! rho_a(l,j) * < l k| T | i j> + !two_e_tc_non_hermit_integral_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i) + !! rho_b(l,j) * < l k| T | i j> + !two_e_tc_non_hermit_integral_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i) + ! rho(l,j) * < k l| T | i j> - two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(l,j,k,i) + two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j) ! rho(l,j) * < k l| T | i j> - two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(l,j,k,i) - ! rho_a(l,j) * < l k| T | i j> + two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(k,i,l,j) + ! rho_a(l,j) * < k l| T | j i> two_e_tc_non_hermit_integral_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i) - ! rho_b(l,j) * < l k| T | i j> + ! rho_b(l,j) * < k l| T | j i> two_e_tc_non_hermit_integral_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i) enddo @@ -84,13 +93,23 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] END_DOC implicit none + double precision, allocatable :: tmp(:,:) 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) + 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 + !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth + Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a endif else @@ -110,14 +129,23 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] END_DOC implicit none + double precision, allocatable :: tmp(:,:) if(bi_ortho) then - 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) ) + !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 - Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b endif else diff --git a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f index b0345957..859c06ca 100644 --- a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f +++ b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f @@ -1,202 +1,266 @@ +! --- + BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo,contrib - fock_a_tot_3e_bi_orth = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - fock_a_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth(a,i) - fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp1_bi_ortho(a,i) - fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp2_bi_ortho(a,i) + + implicit none + integer :: i, a + + fock_a_tot_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + fock_a_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth (a,i) + fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp1_bi_ortho(a,i) + fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp2_bi_ortho(a,i) + enddo enddo - enddo + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo,contrib - fock_b_tot_3e_bi_orth = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - fock_b_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth(a,i) - fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp2_bi_ortho(a,i) - fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp1_bi_ortho(a,i) + + implicit none + integer :: i, a + + fock_b_tot_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + fock_b_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth (a,i) + fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp2_bi_ortho(a,i) + fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp1_bi_ortho(a,i) + enddo enddo - enddo + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_cs_3e_bi_orth, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo, contrib - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - double precision :: new - fock_cs_3e_bi_orth = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = 1, elec_beta_num -! call contrib_3e_sss(a,i,j,k,contrib_sss) -! call contrib_3e_soo(a,i,j,k,contrib_soo) -! call contrib_3e_sos(a,i,j,k,contrib_sos) -! contrib = 0.5d0 * (contrib_sss + contrib_soo) + contrib_sos - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > - call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > - ! negative terms :: exchange contrib - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 - new = 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) & - -1.5d0 * exch_13_int - exch_23_int - fock_cs_3e_bi_orth(a,i) += new + implicit none + integer :: i, a, j, k + double precision :: contrib_sss, contrib_sos, contrib_soo, contrib + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + double precision :: new + + fock_cs_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = 1, elec_beta_num + + !!call contrib_3e_sss(a,i,j,k,contrib_sss) + !!call contrib_3e_soo(a,i,j,k,contrib_soo) + !!call contrib_3e_sos(a,i,j,k,contrib_sos) + !!contrib = 0.5d0 * (contrib_sss + contrib_soo) + contrib_sos + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + + new = 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) -1.5d0 * exch_13_int - exch_23_int + + fock_cs_3e_bi_orth(a,i) += new + enddo + enddo enddo - enddo - enddo - enddo - fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth + + fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth END_PROVIDER +! --- BEGIN_PROVIDER [double precision, fock_a_tmp1_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss, contrib_sos, contrib_soo, contrib - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - double precision :: new - fock_a_tmp1_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - - do j = elec_beta_num + 1, elec_alpha_num - do k = 1, elec_beta_num - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > - call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 - fock_a_tmp1_bi_ortho(a,i) += 1.5d0 * (direct_int - exch_13_int) & - + 0.5d0 * (c_3_int + c_minus_3_int - exch_23_int - exch_12_int) - enddo - enddo + implicit none + integer :: i, a, j, k + double precision :: contrib_sss, contrib_sos, contrib_soo, contrib + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + double precision :: new + + fock_a_tmp1_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + + do j = elec_beta_num + 1, elec_alpha_num + do k = 1, elec_beta_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + + fock_a_tmp1_bi_ortho(a,i) += 1.5d0 * (direct_int - exch_13_int) + 0.5d0 * (c_3_int + c_minus_3_int - exch_23_int - exch_12_int) + enddo + enddo + enddo enddo - enddo - fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho + + fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_a_tmp2_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_sss - fock_a_tmp2_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - do j = 1, elec_alpha_num - do k = elec_beta_num+1, elec_alpha_num - call contrib_3e_sss(a,i,j,k,contrib_sss) - fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss + + implicit none + integer :: i, a, j, k + double precision :: contrib_sss + + fock_a_tmp2_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = 1, elec_alpha_num + do k = elec_beta_num+1, elec_alpha_num + call contrib_3e_sss(a, i, j, k, contrib_sss) + + fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss + enddo + enddo enddo - enddo enddo - enddo + END_PROVIDER - - - +! --- BEGIN_PROVIDER [double precision, fock_b_tmp1_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int - double precision :: new - fock_b_tmp1_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = elec_beta_num+1, elec_alpha_num - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - fock_b_tmp1_bi_ortho(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int - enddo - enddo + implicit none + integer :: i, a, j, k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int + double precision :: new + + fock_b_tmp1_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = 1, elec_beta_num + do k = elec_beta_num+1, elec_alpha_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + + fock_b_tmp1_bi_ortho(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int + enddo + enddo + enddo enddo - enddo - fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho + + fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho + END_PROVIDER +! --- + BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)] - implicit none - integer :: i,a,j,k - double precision :: contrib_soo - fock_b_tmp2_bi_ortho = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - do j = elec_beta_num + 1, elec_alpha_num - do k = 1, elec_alpha_num - call contrib_3e_soo(a,i,j,k,contrib_soo) - fock_b_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_soo + + implicit none + integer :: i, a, j, k + double precision :: contrib_soo + + fock_b_tmp2_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = elec_beta_num + 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_soo(a, i, j, k, contrib_soo) + + fock_b_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_soo + enddo + enddo enddo - enddo enddo - enddo + END_PROVIDER -subroutine contrib_3e_sss(a,i,j,k,integral) - implicit none - integer, intent(in) :: a,i,j,k - BEGIN_DOC - ! returns the pure same spin contribution to F(a,i) from two orbitals j,k - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > - call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > - integral = direct_int + c_3_int + c_minus_3_int - ! negative terms :: exchange contrib - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 - integral += - exch_13_int - exch_23_int - exch_12_int - integral = -integral +! --- + +subroutine contrib_3e_sss(a, i, j, k, integral) + + BEGIN_DOC + ! returns the pure same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + integral = direct_int + c_3_int + c_minus_3_int + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + integral += - exch_13_int - exch_23_int - exch_12_int + + integral = -integral + end +! --- + subroutine contrib_3e_soo(a,i,j,k,integral) - implicit none - integer, intent(in) :: a,i,j,k - BEGIN_DOC - ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_23_int - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 - integral = direct_int - exch_23_int - integral = -integral + + BEGIN_DOC + ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_23_int + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + integral = direct_int - exch_23_int + + integral = -integral + end -subroutine contrib_3e_sos(a,i,j,k,integral) - implicit none - integer, intent(in) :: a,i,j,k - BEGIN_DOC - ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_13_int - call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > - call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 - integral = direct_int - exch_13_int - integral = -integral +! --- + +subroutine contrib_3e_sos(a, i, j, k, integral) + + BEGIN_DOC + ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + integral = direct_int - exch_13_int + + integral = -integral + end + +! --- + diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index c84c837f..5090cc97 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -145,8 +145,10 @@ subroutine simple_tcscf() endif e_delta = dabs(TC_HF_energy - e_save) - print *, ' delta E = ', e_delta - print *, ' gradient = ', grad_non_hermit + print *, ' delta E = ', e_delta + print *, ' gradient = ', grad_non_hermit + print *, ' max TC DIIS error = ', maxval(abs(FQS_SQF_mo)) + !print *, ' gradient= ', grad_non_hermit_right !rho_new = TCSCF_bi_ort_dm_ao @@ -168,6 +170,8 @@ subroutine simple_tcscf() TOUCH mo_l_coef mo_r_coef call ezfio_set_tc_scf_bitc_energy(TC_HF_energy) + !call test_fock_3e_uhf_mo() + print *, ' ***' print *, '' @@ -202,3 +206,64 @@ end subroutine simple_tcscf ! --- +subroutine test_fock_3e_uhf_mo() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + + thr_ih = 1d-12 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_b + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_a(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_a (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' norm_a = ', norm + print *, ' ' + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_b(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_b (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' norm_b = ', norm + print *, ' ' + + ! --- + +end subroutine test_fock_3e_uhf_mo() + diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f index 6961d2f0..e12e629b 100644 --- a/src/tc_scf/test_int.irp.f +++ b/src/tc_scf/test_int.irp.f @@ -9,22 +9,29 @@ program test_ints print *, 'starting ...' my_grid_becke = .True. -! my_n_pt_r_grid = 30 -! my_n_pt_a_grid = 50 - my_n_pt_r_grid = 10 ! small grid for quick debug - my_n_pt_a_grid = 26 ! small grid for quick debug + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + ! my_n_pt_r_grid = 10 ! small grid for quick debug + ! my_n_pt_a_grid = 26 ! small grid for quick debug touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid - call routine_int2_u_grad1u_j1b2 - call routine_v_ij_erf_rk_cst_mu_j1b - call routine_x_v_ij_erf_rk_cst_mu_tmp_j1b - call routine_v_ij_u_cst_mu_j1b + !call routine_int2_u_grad1u_j1b2 + !call routine_v_ij_erf_rk_cst_mu_j1b + !call routine_x_v_ij_erf_rk_cst_mu_tmp_j1b + !call routine_v_ij_u_cst_mu_j1b ! ! call routine_test_j1b -! call routine_int2_grad1u2_grad2u2_j1b2 + !call routine_int2_grad1u2_grad2u2_j1b2 + + + !call test_fock_3e_uhf_ao() + call test_fock_3e_uhf_mo() + end +! --- + subroutine routine_test_j1b implicit none integer :: i,icount,j @@ -286,13 +293,13 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2 double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) double precision, allocatable :: ints(:,:,:) allocate(ints(ao_num, ao_num, n_points_final_grid)) - do ipoint = 1, n_points_final_grid - do i = 1, ao_num - do j = 1, ao_num - read(33,*)ints(j,i,ipoint) - enddo - enddo - enddo +! do ipoint = 1, n_points_final_grid +! do i = 1, ao_num +! do j = 1, ao_num +! read(33,*)ints(j,i,ipoint) +! enddo +! enddo +! enddo allocate(array(ao_num, ao_num, ao_num, ao_num)) array = 0.d0 @@ -344,3 +351,149 @@ subroutine routine_int2_grad1u2_grad2u2_j1b2 end + +! --- + +subroutine test_fock_3e_uhf_ao() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + double precision, allocatable :: fock_3e_uhf_ao_a_mo(:,:), fock_3e_uhf_ao_b_mo(:,:) + + thr_ih = 1d-7 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + + ! --- + + PROVIDE fock_3e_uhf_ao_a + + allocate(fock_3e_uhf_ao_a_mo(mo_num,mo_num)) + call ao_to_mo_bi_ortho( fock_3e_uhf_ao_a , size(fock_3e_uhf_ao_a , 1) & + , fock_3e_uhf_ao_a_mo, size(fock_3e_uhf_ao_a_mo, 1) ) + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_ao_a_mo(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_ao_a_mo (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' ' + + deallocate(fock_3e_uhf_ao_a_mo) + + ! --- + + PROVIDE fock_3e_uhf_ao_b + + allocate(fock_3e_uhf_ao_b_mo(mo_num,mo_num)) + call ao_to_mo_bi_ortho( fock_3e_uhf_ao_b , size(fock_3e_uhf_ao_b , 1) & + , fock_3e_uhf_ao_b_mo, size(fock_3e_uhf_ao_b_mo, 1) ) + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_ao_b_mo(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_ao_b_mo (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + + deallocate(fock_3e_uhf_ao_b_mo) + + ! --- + +end subroutine test_fock_3e_uhf_ao() + +! --- + +subroutine test_fock_3e_uhf_mo() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + + thr_ih = 1d-12 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_b + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_a(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_a (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' norm_a = ', norm + print *, ' ' + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_b(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_b (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' norm_b = ', norm + print *, ' ' + + ! --- + +end subroutine test_fock_3e_uhf_mo() + +! --- + + + + +