diff --git a/src/tc_scf/diago_bi_ort_tcfock.irp.f b/src/tc_scf/diago_bi_ort_tcfock.irp.f index 726169d9..02545315 100644 --- a/src/tc_scf/diago_bi_ort_tcfock.irp.f +++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f @@ -38,33 +38,9 @@ , fock_tc_leigvec_mo, fock_tc_reigvec_mo & , n_real_tc, eigval_right_tmp ) - !if(max_ov_tc_scf)then - ! call non_hrmt_fock_mat( mo_num, F_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag & - ! , fock_tc_leigvec_mo, fock_tc_reigvec_mo & - ! , n_real_tc, eigval_right_tmp ) - !else - ! call non_hrmt_diag_split_degen_bi_orthog( mo_num, F_tmp & - ! , fock_tc_leigvec_mo, fock_tc_reigvec_mo & - ! , n_real_tc, eigval_right_tmp ) - !endif - deallocate(F_tmp) - -! if(n_real_tc .ne. mo_num)then -! print*,'n_real_tc ne mo_num ! ',n_real_tc -! stop -! endif - eigval_fock_tc_mo = eigval_right_tmp -! print*,'Eigenvalues of Fock_matrix_tc_mo_tot' -! do i = 1, elec_alpha_num -! print*, i, eigval_fock_tc_mo(i) -! enddo -! do i = elec_alpha_num+1, mo_num -! print*, i, eigval_fock_tc_mo(i) - level_shift_tcscf -! enddo -! deallocate( eigval_right_tmp ) ! L.T x R call dgemm( "T", "N", mo_num, mo_num, mo_num, 1.d0 & 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 fccfd837..d8b962d7 100644 --- a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f +++ b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -49,6 +49,11 @@ 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 + 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 @@ -145,6 +150,11 @@ END_PROVIDER ! --- BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (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 integer :: a, b, i, j, o diff --git a/src/tc_scf/fock_for_right.irp.f b/src/tc_scf/fock_hermit.irp.f similarity index 100% rename from src/tc_scf/fock_for_right.irp.f rename to src/tc_scf/fock_hermit.irp.f diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 6796666d..e21938de 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -6,10 +6,11 @@ BEGIN_DOC ! - ! two_e_tc_non_hermit_integral_seq_alpha(k,i) = + ! two_e_tc_non_hermit_integral_seq_alpha(k,i) = ON THE AO BASIS ! - ! where F^tc is the two-body part of the TC Fock matrix and k,i are AO basis functions + ! where F^tc is the TWO-BODY part of the TC Fock matrix and k,i are AO basis functions ! + ! works in SEQUENTIAL END_DOC implicit none @@ -17,8 +18,6 @@ double precision :: density, density_a, density_b double precision :: t0, t1 - !print*, ' providing two_e_tc_non_hermit_integral_seq ...' - !call wall_time(t0) two_e_tc_non_hermit_integral_seq_alpha = 0.d0 two_e_tc_non_hermit_integral_seq_beta = 0.d0 @@ -32,24 +31,6 @@ 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_seq_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_seq_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_seq_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_seq_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) - !! 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_seq_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j) ! rho(l,j) * < k l| T | i j> @@ -64,8 +45,6 @@ enddo enddo - !call wall_time(t1) - !print*, ' wall time for two_e_tc_non_hermit_integral_seq after = ', t1 - t0 END_PROVIDER @@ -76,9 +55,9 @@ END_PROVIDER BEGIN_DOC ! - ! two_e_tc_non_hermit_integral_alpha(k,i) = + ! two_e_tc_non_hermit_integral_alpha(k,i) = ON THE AO BASIS ! - ! where F^tc is the two-body part of the TC Fock matrix and k,i are AO basis functions + ! where F^tc is the TWO-BODY part of the TC Fock matrix and k,i are AO basis functions ! END_DOC @@ -88,8 +67,6 @@ END_PROVIDER double precision :: t0, t1 double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) - !print*, ' providing two_e_tc_non_hermit_integral ...' - !call wall_time(t0) two_e_tc_non_hermit_integral_alpha = 0.d0 two_e_tc_non_hermit_integral_beta = 0.d0 @@ -135,8 +112,6 @@ END_PROVIDER deallocate(tmp_a, tmp_b) !$OMP END PARALLEL - !call wall_time(t1) - !print*, ' wall time for two_e_tc_non_hermit_integral after = ', t1 - t0 END_PROVIDER @@ -181,14 +156,6 @@ 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) - 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 @@ -217,14 +184,6 @@ 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 diff --git a/src/tc_scf/fock_tc_mo_tot.irp.f b/src/tc_scf/fock_tc_mo_tot.irp.f index 2f33cd17..a03a0624 100644 --- a/src/tc_scf/fock_tc_mo_tot.irp.f +++ b/src/tc_scf/fock_tc_mo_tot.irp.f @@ -3,7 +3,7 @@ &BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)] implicit none BEGIN_DOC - ! Fock matrix on the MO basis. + ! TC-Fock matrix on the MO basis. WARNING !!! NON HERMITIAN !!! ! For open shells, the ROHF Fock Matrix is :: ! ! | F-K | F + K/2 | F | diff --git a/src/tc_scf/fock_three_bi_ortho.irp.f b/src/tc_scf/fock_three_bi_ortho.irp.f index 279670b8..7c776ce5 100644 --- a/src/tc_scf/fock_three_bi_ortho.irp.f +++ b/src/tc_scf/fock_three_bi_ortho.irp.f @@ -1,178 +1,296 @@ -BEGIN_PROVIDER [ double precision, fock_a_abb_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC -! fock_a_abb_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,beta,beta contribution - END_DOC - fock_a_abb_3e_bi_orth_old = 0.d0 - integer :: i,a,j,k - double precision :: direct_int, exch_23_int - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = j+1, elec_beta_num - ! see contrib_3e_soo - 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 - fock_a_abb_3e_bi_orth_old(a,i) += direct_int - exch_23_int + +! --- + +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 + END_DOC + implicit none + integer :: i, a + + PROVIDE mo_l_coef mo_r_coef + + 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 - enddo - fock_a_abb_3e_bi_orth_old = - fock_a_abb_3e_bi_orth_old + END_PROVIDER -BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC -! fock_a_aba_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,beta contribution - END_DOC - fock_a_aba_3e_bi_orth_old = 0.d0 - integer :: i,a,j,k - double precision :: direct_int, exch_13_int - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_alpha_num ! a - do k = 1, elec_beta_num ! b - ! a b a a b a - 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 - fock_a_aba_3e_bi_orth_old(a,i) += direct_int - exch_13_int +! --- + +BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (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 + integer :: i, a + + PROVIDE mo_l_coef mo_r_coef + + 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 - enddo - fock_a_aba_3e_bi_orth_old = - fock_a_aba_3e_bi_orth_old + END_PROVIDER -BEGIN_PROVIDER [ double precision, fock_a_aaa_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC -! fock_a_aaa_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution - END_DOC - fock_a_aaa_3e_bi_orth_old = 0.d0 - integer :: i,a,j,k - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_alpha_num - do k = j+1, elec_alpha_num - ! positive terms :: cycle contrib - 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 > - fock_a_aaa_3e_bi_orth_old(a,i) += 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 - fock_a_aaa_3e_bi_orth_old(a,i) += - exch_13_int - exch_23_int - exch_12_int +! --- + +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 + + PROVIDE mo_l_coef mo_r_coef + + 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_a_aaa_3e_bi_orth_old = - fock_a_aaa_3e_bi_orth_old -END_PROVIDER - -BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC - ! fock_a_tot_3e_bi_orth_old = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions - END_DOC - fock_a_tot_3e_bi_orth_old = fock_a_abb_3e_bi_orth_old + fock_a_aba_3e_bi_orth_old + fock_a_aaa_3e_bi_orth_old + + fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth END_PROVIDER -BEGIN_PROVIDER [ double precision, fock_b_baa_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC -! fock_b_baa_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,alpha contribution - END_DOC - fock_b_baa_3e_bi_orth_old = 0.d0 - integer :: i,a,j,k - double precision :: direct_int, exch_23_int - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_alpha_num - do k = j+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, i, j, k, exch_23_int)! < a k j | i j k > : E_23 - fock_b_baa_3e_bi_orth_old(a,i) += direct_int - exch_23_int +! --- + +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 + + PROVIDE mo_l_coef mo_r_coef + + 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 - enddo - fock_b_baa_3e_bi_orth_old = - fock_b_baa_3e_bi_orth_old + + fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho + END_PROVIDER -BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC -! fock_b_bab_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,beta contribution - END_DOC - fock_b_bab_3e_bi_orth_old = 0.d0 - integer :: i,a,j,k - double precision :: direct_int, exch_13_int - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = 1, elec_alpha_num - ! b a b b a b - 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 - fock_b_bab_3e_bi_orth_old(a,i) += direct_int - exch_13_int +! --- + +BEGIN_PROVIDER [double precision, fock_a_tmp2_bi_ortho, (mo_num, mo_num)] + + implicit none + integer :: i, a, j, k + double precision :: contrib_sss + + PROVIDE mo_l_coef mo_r_coef + + 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 - fock_b_bab_3e_bi_orth_old = - fock_b_bab_3e_bi_orth_old + END_PROVIDER -BEGIN_PROVIDER [ double precision, fock_b_bbb_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC -! fock_b_bbb_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution - END_DOC - fock_b_bbb_3e_bi_orth_old = 0.d0 - integer :: i,a,j,k - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - do i = 1, mo_num - do a = 1, mo_num - - do j = 1, elec_beta_num - do k = j+1, elec_beta_num - ! positive terms :: cycle contrib - 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 > - fock_b_bbb_3e_bi_orth_old(a,i) += 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 - fock_b_bbb_3e_bi_orth_old(a,i) += - exch_13_int - exch_23_int - exch_12_int +! --- + +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 + + PROVIDE mo_l_coef mo_r_coef + + 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 - enddo - fock_b_bbb_3e_bi_orth_old = - fock_b_bbb_3e_bi_orth_old -END_PROVIDER -BEGIN_PROVIDER [ double precision, fock_b_tot_3e_bi_orth_old, (mo_num, mo_num)] - implicit none - BEGIN_DOC - ! fock_b_tot_3e_bi_orth_old = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions - END_DOC - fock_b_tot_3e_bi_orth_old = fock_b_bbb_3e_bi_orth_old + fock_b_bab_3e_bi_orth_old + fock_b_baa_3e_bi_orth_old + 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 + + PROVIDE mo_l_coef mo_r_coef + + 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 + +END_PROVIDER + +! --- + +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 + + PROVIDE mo_l_coef mo_r_coef + + 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) + + 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 + + PROVIDE mo_l_coef mo_r_coef + + 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) + + BEGIN_DOC + ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + PROVIDE mo_l_coef mo_r_coef + + 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/fock_three_bi_ortho_new_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f deleted file mode 100644 index f73171a3..00000000 --- a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f +++ /dev/null @@ -1,286 +0,0 @@ - -! --- - -BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] - - implicit none - integer :: i, a - - PROVIDE mo_l_coef mo_r_coef - - 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 - -END_PROVIDER - -! --- - -BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)] - - implicit none - integer :: i, a - - PROVIDE mo_l_coef mo_r_coef - - 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 - -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 - - PROVIDE mo_l_coef mo_r_coef - - 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 - - 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 - - PROVIDE mo_l_coef mo_r_coef - - 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 - - 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 - - PROVIDE mo_l_coef mo_r_coef - - 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 - -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 - - PROVIDE mo_l_coef mo_r_coef - - 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 - - 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 - - PROVIDE mo_l_coef mo_r_coef - - 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 - -END_PROVIDER - -! --- - -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 - - PROVIDE mo_l_coef mo_r_coef - - 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) - - 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 - - PROVIDE mo_l_coef mo_r_coef - - 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) - - BEGIN_DOC - ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k - END_DOC - - PROVIDE mo_l_coef mo_r_coef - - 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/fock_three.irp.f b/src/tc_scf/fock_three_hermit.irp.f similarity index 68% rename from src/tc_scf/fock_three.irp.f rename to src/tc_scf/fock_three_hermit.irp.f index 424eeffd..5c48970b 100644 --- a/src/tc_scf/fock_three.irp.f +++ b/src/tc_scf/fock_three_hermit.irp.f @@ -227,3 +227,144 @@ BEGIN_PROVIDER [ double precision, fock_3_mat_b_op_sh, (mo_num, mo_num)] enddo END_PROVIDER + + +BEGIN_PROVIDER [ double precision, fock_3_w_kk_sum, (n_points_final_grid,3)] + implicit none + integer :: mm, ipoint,k + double precision :: w_kk + fock_3_w_kk_sum = 0.d0 + do k = 1, elec_beta_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + w_kk = x_W_ij_erf_rk(ipoint,mm,k,k) + fock_3_w_kk_sum(ipoint,mm) += w_kk + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_w_ki_mos_k, (n_points_final_grid,3,mo_num)] + implicit none + integer :: mm, ipoint,k,i + double precision :: w_ki, mo_k + fock_3_w_ki_mos_k = 0.d0 + do i = 1, mo_num + do k = 1, elec_beta_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + w_ki = x_W_ij_erf_rk(ipoint,mm,k,i) + mo_k = mos_in_r_array(k,ipoint) + fock_3_w_ki_mos_k(ipoint,mm,i) += w_ki * mo_k + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_w_kl_w_kl, (n_points_final_grid,3)] + implicit none + integer :: k,j,ipoint,mm + double precision :: w_kj + fock_3_w_kl_w_kl = 0.d0 + do j = 1, elec_beta_num + do k = 1, elec_beta_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + w_kj = x_W_ij_erf_rk(ipoint,mm,k,j) + fock_3_w_kl_w_kl(ipoint,mm) += w_kj * w_kj + enddo + enddo + enddo + enddo + + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_rho_beta, (n_points_final_grid)] + implicit none + integer :: ipoint,k + fock_3_rho_beta = 0.d0 + do ipoint = 1, n_points_final_grid + do k = 1, elec_beta_num + fock_3_rho_beta(ipoint) += mos_in_r_array(k,ipoint) * mos_in_r_array(k,ipoint) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_w_kl_mo_k_mo_l, (n_points_final_grid,3)] + implicit none + integer :: ipoint,k,l,mm + double precision :: mos_k, mos_l, w_kl + fock_3_w_kl_mo_k_mo_l = 0.d0 + do k = 1, elec_beta_num + do l = 1, elec_beta_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + mos_k = mos_in_r_array_transp(ipoint,k) + mos_l = mos_in_r_array_transp(ipoint,l) + w_kl = x_W_ij_erf_rk(ipoint,mm,l,k) + fock_3_w_kl_mo_k_mo_l(ipoint,mm) += w_kl * mos_k * mos_l + enddo + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_w_ki_wk_a, (n_points_final_grid,3,mo_num, mo_num)] + implicit none + integer :: ipoint,i,a,k,mm + double precision :: w_ki,w_ka + fock_3_w_ki_wk_a = 0.d0 + do i = 1, mo_num + do a = 1, mo_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + do k = 1, elec_beta_num + w_ki = x_W_ij_erf_rk(ipoint,mm,k,i) + w_ka = x_W_ij_erf_rk(ipoint,mm,k,a) + fock_3_w_ki_wk_a(ipoint,mm,a,i) += w_ki * w_ka + enddo + enddo + enddo + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_trace_w_tilde, (n_points_final_grid,3)] + implicit none + integer :: ipoint,k,mm + fock_3_trace_w_tilde = 0.d0 + do k = 1, elec_beta_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + fock_3_trace_w_tilde(ipoint,mm) += fock_3_w_ki_wk_a(ipoint,mm,k,k) + enddo + enddo + enddo + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_3_w_kl_wla_phi_k, (n_points_final_grid,3,mo_num)] + implicit none + integer :: ipoint,a,k,mm,l + double precision :: w_kl,w_la, mo_k + fock_3_w_kl_wla_phi_k = 0.d0 + do a = 1, mo_num + do k = 1, elec_beta_num + do l = 1, elec_beta_num + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + w_kl = x_W_ij_erf_rk(ipoint,mm,l,k) + w_la = x_W_ij_erf_rk(ipoint,mm,l,a) + mo_k = mos_in_r_array_transp(ipoint,k) + fock_3_w_kl_wla_phi_k(ipoint,mm,a) += w_kl * w_la * mo_k + enddo + enddo + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/tc_scf/fock_three_utils.irp.f b/src/tc_scf/fock_three_utils.irp.f deleted file mode 100644 index 5aec1d9e..00000000 --- a/src/tc_scf/fock_three_utils.irp.f +++ /dev/null @@ -1,140 +0,0 @@ - -BEGIN_PROVIDER [ double precision, fock_3_w_kk_sum, (n_points_final_grid,3)] - implicit none - integer :: mm, ipoint,k - double precision :: w_kk - fock_3_w_kk_sum = 0.d0 - do k = 1, elec_beta_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - w_kk = x_W_ij_erf_rk(ipoint,mm,k,k) - fock_3_w_kk_sum(ipoint,mm) += w_kk - enddo - enddo - enddo -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_w_ki_mos_k, (n_points_final_grid,3,mo_num)] - implicit none - integer :: mm, ipoint,k,i - double precision :: w_ki, mo_k - fock_3_w_ki_mos_k = 0.d0 - do i = 1, mo_num - do k = 1, elec_beta_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - w_ki = x_W_ij_erf_rk(ipoint,mm,k,i) - mo_k = mos_in_r_array(k,ipoint) - fock_3_w_ki_mos_k(ipoint,mm,i) += w_ki * mo_k - enddo - enddo - enddo - enddo - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_w_kl_w_kl, (n_points_final_grid,3)] - implicit none - integer :: k,j,ipoint,mm - double precision :: w_kj - fock_3_w_kl_w_kl = 0.d0 - do j = 1, elec_beta_num - do k = 1, elec_beta_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - w_kj = x_W_ij_erf_rk(ipoint,mm,k,j) - fock_3_w_kl_w_kl(ipoint,mm) += w_kj * w_kj - enddo - enddo - enddo - enddo - - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_rho_beta, (n_points_final_grid)] - implicit none - integer :: ipoint,k - fock_3_rho_beta = 0.d0 - do ipoint = 1, n_points_final_grid - do k = 1, elec_beta_num - fock_3_rho_beta(ipoint) += mos_in_r_array(k,ipoint) * mos_in_r_array(k,ipoint) - enddo - enddo -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_w_kl_mo_k_mo_l, (n_points_final_grid,3)] - implicit none - integer :: ipoint,k,l,mm - double precision :: mos_k, mos_l, w_kl - fock_3_w_kl_mo_k_mo_l = 0.d0 - do k = 1, elec_beta_num - do l = 1, elec_beta_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - mos_k = mos_in_r_array_transp(ipoint,k) - mos_l = mos_in_r_array_transp(ipoint,l) - w_kl = x_W_ij_erf_rk(ipoint,mm,l,k) - fock_3_w_kl_mo_k_mo_l(ipoint,mm) += w_kl * mos_k * mos_l - enddo - enddo - enddo - enddo - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_w_ki_wk_a, (n_points_final_grid,3,mo_num, mo_num)] - implicit none - integer :: ipoint,i,a,k,mm - double precision :: w_ki,w_ka - fock_3_w_ki_wk_a = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - do k = 1, elec_beta_num - w_ki = x_W_ij_erf_rk(ipoint,mm,k,i) - w_ka = x_W_ij_erf_rk(ipoint,mm,k,a) - fock_3_w_ki_wk_a(ipoint,mm,a,i) += w_ki * w_ka - enddo - enddo - enddo - enddo - enddo -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_trace_w_tilde, (n_points_final_grid,3)] - implicit none - integer :: ipoint,k,mm - fock_3_trace_w_tilde = 0.d0 - do k = 1, elec_beta_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - fock_3_trace_w_tilde(ipoint,mm) += fock_3_w_ki_wk_a(ipoint,mm,k,k) - enddo - enddo - enddo - -END_PROVIDER - -BEGIN_PROVIDER [ double precision, fock_3_w_kl_wla_phi_k, (n_points_final_grid,3,mo_num)] - implicit none - integer :: ipoint,a,k,mm,l - double precision :: w_kl,w_la, mo_k - fock_3_w_kl_wla_phi_k = 0.d0 - do a = 1, mo_num - do k = 1, elec_beta_num - do l = 1, elec_beta_num - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - w_kl = x_W_ij_erf_rk(ipoint,mm,l,k) - w_la = x_W_ij_erf_rk(ipoint,mm,l,a) - mo_k = mos_in_r_array_transp(ipoint,k) - fock_3_w_kl_wla_phi_k(ipoint,mm,a) += w_kl * w_la * mo_k - enddo - enddo - enddo - enddo - enddo -END_PROVIDER - diff --git a/src/tc_scf/minimize_tc_angles.irp.f b/src/tc_scf/minimize_tc_angles.irp.f index cb729eb2..1363e62b 100644 --- a/src/tc_scf/minimize_tc_angles.irp.f +++ b/src/tc_scf/minimize_tc_angles.irp.f @@ -1,10 +1,11 @@ program print_angles implicit none + BEGIN_DOC +! program that minimizes the angle between left- and right-orbitals when degeneracies are found in the TC-Fock matrix + END_DOC 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 = 14 ! small grid for quick debug touch my_n_pt_r_grid my_n_pt_a_grid ! call sort_by_tc_fock call minimize_tc_orb_angles diff --git a/src/tc_scf/print_angle_tc_orb.irp.f b/src/tc_scf/print_angle_tc_orb.irp.f deleted file mode 100644 index 09260395..00000000 --- a/src/tc_scf/print_angle_tc_orb.irp.f +++ /dev/null @@ -1,9 +0,0 @@ -program print_angles - implicit none - 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 = 14 ! small grid for quick debug - call print_angles_tc -end diff --git a/src/tc_scf/rotate_tcscf_orbitals.irp.f b/src/tc_scf/rotate_tcscf_orbitals.irp.f index fc4a7935..31999c18 100644 --- a/src/tc_scf/rotate_tcscf_orbitals.irp.f +++ b/src/tc_scf/rotate_tcscf_orbitals.irp.f @@ -4,7 +4,7 @@ program rotate_tcscf_orbitals BEGIN_DOC - ! TODO : Put the documentation of the program here + ! TODO : Rotate the bi-orthonormal orbitals in order to minimize left-right angles when degenerate END_DOC implicit none diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f index 596ae500..3c12118f 100644 --- a/src/tc_scf/routines_rotates.irp.f +++ b/src/tc_scf/routines_rotates.irp.f @@ -1,7 +1,54 @@ + +! --- + +subroutine LTxSxR(n, m, L, S, R, C) + + implicit none + integer, intent(in) :: n, m + double precision, intent(in) :: L(n,m), S(n,n), R(n,m) + double precision, intent(out) :: C(m,m) + integer :: i, j + double precision :: accu_d, accu_nd + double precision, allocatable :: tmp(:,:) + + ! L.T x S x R + allocate(tmp(m,n)) + call dgemm( 'T', 'N', m, n, n, 1.d0 & + , L, size(L, 1), S, size(S, 1) & + , 0.d0, tmp, size(tmp, 1) ) + call dgemm( 'N', 'N', m, m, n, 1.d0 & + , tmp, size(tmp, 1), R, size(R, 1) & + , 0.d0, C, size(C, 1) ) + deallocate(tmp) + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(j.eq.i) then + accu_d += dabs(C(j,i)) + else + accu_nd += C(j,i) * C(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) + + print*, ' accu_d = ', accu_d + print*, ' accu_nd = ', accu_nd + +end subroutine LTxR + +! --- + + ! --- subroutine minimize_tc_orb_angles() + BEGIN_DOC + ! routine that minimizes the angle between left- and right-orbitals when degeneracies are found + END_DOC implicit none logical :: good_angles diff --git a/src/tc_scf/tc_scf_dm.irp.f b/src/tc_scf/tc_scf_dm.irp.f index 90719f47..07da8a58 100644 --- a/src/tc_scf/tc_scf_dm.irp.f +++ b/src/tc_scf/tc_scf_dm.irp.f @@ -2,6 +2,9 @@ BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num) ] + BEGIN_DOC + ! TC-SCF transition density matrix on the AO basis for BETA electrons + END_DOC implicit none if(bi_ortho) then @@ -16,6 +19,9 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_num) ] + BEGIN_DOC + ! TC-SCF transition density matrix on the AO basis for ALPHA electrons + END_DOC implicit none if(bi_ortho) then @@ -31,6 +37,9 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_tot, (ao_num, ao_num) ] implicit none + BEGIN_DOC + ! TC-SCF transition density matrix on the AO basis for ALPHA+BETA electrons + END_DOC TCSCF_density_matrix_ao_tot = TCSCF_density_matrix_ao_beta + TCSCF_density_matrix_ao_alpha END_PROVIDER diff --git a/src/tc_scf/tc_scf_utils.irp.f b/src/tc_scf/tc_scf_utils.irp.f deleted file mode 100644 index dde477c4..00000000 --- a/src/tc_scf/tc_scf_utils.irp.f +++ /dev/null @@ -1,43 +0,0 @@ - -! --- - -subroutine LTxSxR(n, m, L, S, R, C) - - implicit none - integer, intent(in) :: n, m - double precision, intent(in) :: L(n,m), S(n,n), R(n,m) - double precision, intent(out) :: C(m,m) - integer :: i, j - double precision :: accu_d, accu_nd - double precision, allocatable :: tmp(:,:) - - ! L.T x S x R - allocate(tmp(m,n)) - call dgemm( 'T', 'N', m, n, n, 1.d0 & - , L, size(L, 1), S, size(S, 1) & - , 0.d0, tmp, size(tmp, 1) ) - call dgemm( 'N', 'N', m, m, n, 1.d0 & - , tmp, size(tmp, 1), R, size(R, 1) & - , 0.d0, C, size(C, 1) ) - deallocate(tmp) - - accu_d = 0.d0 - accu_nd = 0.d0 - do i = 1, m - do j = 1, m - if(j.eq.i) then - accu_d += dabs(C(j,i)) - else - accu_nd += C(j,i) * C(j,i) - endif - enddo - enddo - accu_nd = dsqrt(accu_nd) - - print*, ' accu_d = ', accu_d - print*, ' accu_nd = ', accu_nd - -end subroutine LTxR - -! --- - diff --git a/src/tc_scf/test_Ne.sh b/src/tc_scf/test_Ne.sh deleted file mode 100755 index a6422931..00000000 --- a/src/tc_scf/test_Ne.sh +++ /dev/null @@ -1,13 +0,0 @@ -QP_ROOT=/home/eginer/new_qp2/qp2 -source ${QP_ROOT}/quantum_package.rc - echo Ne > Ne.xyz - echo $QP_ROOT - qp create_ezfio -b cc-pcvdz Ne.xyz -o Ne_tc_scf - qp run scf - qp set tc_keywords bi_ortho True - qp set ao_two_e_erf_ints mu_erf 0.87 - qp set tc_keywords j1b_pen [1.5] - qp set tc_keywords j1b_type 3 - qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out - grep "TC energy =" Ne.ezfio.tc_scf.out | tail -1 - eref=-128.552134 diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f deleted file mode 100644 index a14c4126..00000000 --- a/src/tc_scf/test_int.irp.f +++ /dev/null @@ -1,1003 +0,0 @@ -program test_ints - - BEGIN_DOC -! TODO : Put the documentation of the program here - END_DOC - - implicit none - - print *, ' starting test_ints ...' - - my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 -! my_n_pt_r_grid = 15 ! 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 - - my_extra_grid_becke = .True. - my_n_pt_r_extra_grid = 30 - my_n_pt_a_extra_grid = 50 ! small extra_grid for quick debug - touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid - -!! OK -!call routine_int2_u_grad1u_j1b2 -!! OK -!call routine_v_ij_erf_rk_cst_mu_j1b -!! OK -! call routine_x_v_ij_erf_rk_cst_mu_j1b -!! OK -! call routine_v_ij_u_cst_mu_j1b - -!! OK -!call routine_int2_u2_j1b2 - -!! OK -!call routine_int2_u_grad1u_x_j1b2 - -!! OK -! call routine_int2_grad1u2_grad2u2_j1b2 -! call routine_int2_u_grad1u_j1b2 -! call test_total_grad_lapl -! call test_total_grad_square -! call test_ao_tc_int_chemist -! call test_grid_points_ao -! call test_tc_scf - !call test_int_gauss - - !call test_fock_3e_uhf_ao() - !call test_fock_3e_uhf_mo() - - !call test_tc_grad_and_lapl_ao() - !call test_tc_grad_square_ao() - - call test_two_e_tc_non_hermit_integral() - -end - -! --- - -subroutine test_tc_scf - implicit none - integer :: i -! provide int2_u_grad1u_x_j1b2_test - provide x_v_ij_erf_rk_cst_mu_j1b_test -! provide x_v_ij_erf_rk_cst_mu_j1b_test -! print*,'TC_HF_energy = ',TC_HF_energy -! print*,'grad_non_hermit = ',grad_non_hermit -end - -subroutine test_ao_tc_int_chemist - implicit none - provide ao_tc_int_chemist -! provide ao_tc_int_chemist_test -! provide tc_grad_square_ao_test -! provide tc_grad_and_lapl_ao_test -end - -! --- - -subroutine routine_test_j1b - implicit none - integer :: i,icount,j - icount = 0 - do i = 1, List_all_comb_b3_size - if(dabs(List_all_comb_b3_coef(i)).gt.1.d-10)then - print*,'' - print*,List_all_comb_b3_expo(i),List_all_comb_b3_coef(i) - print*,List_all_comb_b3_cent(1:3,i) - print*,'' - icount += 1 - endif - - enddo - print*,'List_all_comb_b3_coef,icount = ',List_all_comb_b3_size,icount - do i = 1, ao_num - do j = 1, ao_num - do icount = 1, List_comb_thr_b3_size(j,i) - print*,'',j,i - print*,List_comb_thr_b3_expo(icount,j,i),List_comb_thr_b3_coef(icount,j,i) - print*,List_comb_thr_b3_cent(1:3,icount,j,i) - print*,'' - enddo -! enddo - enddo - enddo - print*,'max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size,List_all_comb_b3_size - -end - -subroutine routine_int2_u_grad1u_j1b2 - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - array(j,i,l,k) += int2_u_grad1u_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += int2_u_grad1u_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - -subroutine routine_v_ij_erf_rk_cst_mu_j1b - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) -! print*,'ao_overlap_abs = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) -! enddo -! print*,'center = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) -! enddo -! print*,'sigma = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) -! enddo - - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - array(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - - -subroutine routine_x_v_ij_erf_rk_cst_mu_j1b - implicit none - integer :: i,j,ipoint,k,l,m - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) -! print*,'ao_overlap_abs = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) -! enddo -! print*,'center = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) -! enddo -! print*,'sigma = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) -! enddo - - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - do m = 1, 3 - array(j,i,l,k) += x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += x_v_ij_erf_rk_cst_mu_j1b (j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - - - -subroutine routine_v_ij_u_cst_mu_j1b_test - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) -! print*,'ao_overlap_abs = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) -! enddo -! print*,'center = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) -! enddo -! print*,'sigma = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) -! enddo - - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - -subroutine routine_int2_grad1u2_grad2u2_j1b2 - implicit none - integer :: i,j,ipoint,k,l - integer :: ii , jj - double precision :: weight,accu_relat, accu_abs, contrib - 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 - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight -! !array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight -! array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight -! !array(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight -! array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight -! if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then -! if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then -! print*,j,i,ipoint -! print*,int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) , int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)) -! print*,int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) , int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) - int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint)) -! stop -! endif -! endif - enddo - enddo - enddo - enddo - enddo - double precision :: e_ref, e_new - accu_relat = 0.d0 - accu_abs = 0.d0 - e_ref = 0.d0 - e_new = 0.d0 - do ii = 1, elec_alpha_num - do jj = ii, elec_alpha_num - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - e_ref += mo_coef(j,ii) * mo_coef(i,ii) * array_ref(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj) - e_new += mo_coef(j,ii) * mo_coef(i,ii) * array(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj) - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib -! if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then -! accu_relat += contrib/dabs(array_ref(j,i,l,k)) -! endif - enddo - enddo - enddo - enddo - - enddo - enddo - print*,'e_ref = ',e_ref - print*,'e_new = ',e_new -! print*,'accu_abs = ',accu_abs/dble(ao_num)**4 -! print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - -subroutine routine_int2_u2_j1b2 - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) -! print*,'ao_overlap_abs = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) -! enddo -! print*,'center = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) -! enddo -! print*,'sigma = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) -! enddo - - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - array(j,i,l,k) += int2_u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += int2_u2_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - - -subroutine routine_int2_u_grad1u_x_j1b2 - implicit none - integer :: i,j,ipoint,k,l,m - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) -! print*,'ao_overlap_abs = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) -! enddo -! print*,'center = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) -! enddo -! print*,'sigma = ' -! do i = 1, ao_num -! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) -! enddo - - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - do m = 1, 3 - array(j,i,l,k) += int2_u_grad1u_x_j1b2_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += int2_u_grad1u_x_j1b2 (j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - - -end - -subroutine routine_v_ij_u_cst_mu_j1b - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) - - allocate(array(ao_num, ao_num, ao_num, ao_num)) - array = 0.d0 - allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) - array_ref = 0.d0 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight - enddo - enddo - enddo - enddo - enddo - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) - accu_abs += contrib - if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(array_ref(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - -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 fock_3e_uhf_ao_b - - ! --- - - 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) - - ! --- - - 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 - print *, ' ' - - 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 - -! --- - -subroutine test_total_grad_lapl - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(tc_grad_and_lapl_ao_test(j,i,l,k) - tc_grad_and_lapl_ao(j,i,l,k)) - accu_abs += contrib - if(dabs(tc_grad_and_lapl_ao(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(tc_grad_and_lapl_ao(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - -end - -subroutine test_total_grad_square - implicit none - integer :: i,j,ipoint,k,l - double precision :: weight,accu_relat, accu_abs, contrib - accu_relat = 0.d0 - accu_abs = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do j = 1, ao_num - contrib = dabs(tc_grad_square_ao_test(j,i,l,k) - tc_grad_square_ao(j,i,l,k)) - accu_abs += contrib - if(dabs(tc_grad_square_ao(j,i,l,k)).gt.1.d-10)then - accu_relat += contrib/dabs(tc_grad_square_ao(j,i,l,k)) - endif - enddo - enddo - enddo - enddo - print*,'accu_abs = ',accu_abs/dble(ao_num)**4 - print*,'accu_relat = ',accu_relat/dble(ao_num)**4 - - -end - -subroutine test_grid_points_ao - implicit none - integer :: i,j,ipoint,icount,icount_good, icount_bad,icount_full - double precision :: thr - thr = 1.d-10 -! print*,'max_n_pts_grid_ao_prod = ',max_n_pts_grid_ao_prod -! print*,'n_pts_grid_ao_prod' - do i = 1, ao_num - do j = i, ao_num - icount = 0 - icount_good = 0 - icount_bad = 0 - icount_full = 0 - do ipoint = 1, n_points_final_grid -! if(dabs(int2_u_grad1u_x_j1b2_test(j,i,ipoint,1)) & -! + dabs(int2_u_grad1u_x_j1b2_test(j,i,ipoint,2)) & -! + dabs(int2_u_grad1u_x_j1b2_test(j,i,ipoint,3)) ) -! if(dabs(int2_u2_j1b2_test(j,i,ipoint)).gt.thr)then -! icount += 1 -! endif - if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then - icount_full += 1 - endif - if(dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)).gt.thr)then - icount += 1 - if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then - icount_good += 1 - else - print*,j,i,ipoint - print*,dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint))/dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)) - icount_bad += 1 - endif - endif -! if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr)then -! endif - enddo - print*,'' - print*,j,i - print*,icount,icount_full, icount_bad!,n_pts_grid_ao_prod(j,i) - print*,dble(icount)/dble(n_points_final_grid),dble(icount_full)/dble(n_points_final_grid) -! dble(n_pts_grid_ao_prod(j,i))/dble(n_points_final_grid) -! if(icount.gt.n_pts_grid_ao_prod(j,i))then -! print*,'pb !!' -! endif - enddo - enddo -end - -subroutine test_int_gauss - implicit none - integer :: i,j - print*,'center' - do i = 1, ao_num - do j = i, ao_num - print*,j,i - print*,ao_prod_sigma(j,i),ao_overlap_abs_grid(j,i) - print*,ao_prod_center(1:3,j,i) - enddo - enddo - print*,'' - double precision :: weight, r(3),integral_1,pi,center(3),f_r,alpha,distance,integral_2 - center = 0.d0 - pi = dacos(-1.d0) - integral_1 = 0.d0 - integral_2 = 0.d0 - alpha = 0.75d0 - do i = 1, n_points_final_grid - ! you get x, y and z of the ith grid point - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - weight = final_weight_at_r_vector(i) - distance = dsqrt( (r(1) - center(1))**2 + (r(2) - center(2))**2 + (r(3) - center(3))**2 ) - f_r = dexp(-alpha * distance*distance) - ! you add the contribution of the grid point to the integral - integral_1 += f_r * weight - integral_2 += f_r * distance * weight - enddo - print*,'integral_1 =',integral_1 - print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5 - print*,'integral_2 =',integral_2 - print*,'(pi/alpha)**1.5 =',2.d0*pi / (alpha)**2 - - -end - -! --- - -subroutine test_tc_grad_and_lapl_ao() - - implicit none - integer :: i, j, k, l - double precision :: diff_tot, diff, thr_ih, norm - - thr_ih = 1d-10 - - PROVIDE tc_grad_and_lapl_ao tc_grad_and_lapl_ao_loop - - norm = 0.d0 - diff_tot = 0.d0 - do i = 1, ao_num - do j = 1, ao_num - do k = 1, ao_num - do l = 1, ao_num - - diff = dabs(tc_grad_and_lapl_ao_loop(l,k,j,i) - tc_grad_and_lapl_ao(l,k,j,i)) - if(diff .gt. thr_ih) then - print *, ' difference on ', l, k, j, i - print *, ' loops : ', tc_grad_and_lapl_ao_loop(l,k,j,i) - print *, ' lapack: ', tc_grad_and_lapl_ao (l,k,j,i) - !stop - endif - - norm += dabs(tc_grad_and_lapl_ao_loop(l,k,j,i)) - diff_tot += diff - enddo - enddo - enddo - enddo - - print *, ' diff tot = ', diff_tot / norm - print *, ' norm = ', norm - print *, ' ' - - return - -end - -! --- - -subroutine test_tc_grad_square_ao() - - implicit none - integer :: i, j, k, l - double precision :: diff_tot, diff, thr_ih, norm - - thr_ih = 1d-10 - - PROVIDE tc_grad_square_ao tc_grad_square_ao_loop - - norm = 0.d0 - diff_tot = 0.d0 - do i = 1, ao_num - do j = 1, ao_num - do k = 1, ao_num - do l = 1, ao_num - - diff = dabs(tc_grad_square_ao_loop(l,k,j,i) - tc_grad_square_ao(l,k,j,i)) - if(diff .gt. thr_ih) then - print *, ' difference on ', l, k, j, i - print *, ' loops : ', tc_grad_square_ao_loop(l,k,j,i) - print *, ' lapack: ', tc_grad_square_ao (l,k,j,i) - !stop - endif - - norm += dabs(tc_grad_square_ao_loop(l,k,j,i)) - diff_tot += diff - enddo - enddo - enddo - enddo - - print *, ' diff tot = ', diff_tot / norm - print *, ' norm = ', norm - print *, ' ' - - return - -end - -! --- - -subroutine test_two_e_tc_non_hermit_integral() - - implicit none - integer :: i, j - double precision :: diff_tot, diff, thr_ih, norm - - thr_ih = 1d-10 - - PROVIDE two_e_tc_non_hermit_integral_beta two_e_tc_non_hermit_integral_alpha - PROVIDE two_e_tc_non_hermit_integral_seq_beta two_e_tc_non_hermit_integral_seq_alpha - - ! --- - - norm = 0.d0 - diff_tot = 0.d0 - do i = 1, ao_num - do j = 1, ao_num - - diff = dabs(two_e_tc_non_hermit_integral_seq_alpha(j,i) - two_e_tc_non_hermit_integral_alpha(j,i)) - if(diff .gt. thr_ih) then - print *, ' difference on ', j, i - print *, ' seq : ', two_e_tc_non_hermit_integral_seq_alpha(j,i) - print *, ' // : ', two_e_tc_non_hermit_integral_alpha (j,i) - !stop - endif - - norm += dabs(two_e_tc_non_hermit_integral_seq_alpha(j,i)) - diff_tot += diff - enddo - enddo - - print *, ' diff tot a = ', diff_tot / norm - print *, ' norm a = ', norm - print *, ' ' - - ! --- - - norm = 0.d0 - diff_tot = 0.d0 - do i = 1, ao_num - do j = 1, ao_num - - diff = dabs(two_e_tc_non_hermit_integral_seq_beta(j,i) - two_e_tc_non_hermit_integral_beta(j,i)) - if(diff .gt. thr_ih) then - print *, ' difference on ', j, i - print *, ' seq : ', two_e_tc_non_hermit_integral_seq_beta(j,i) - print *, ' // : ', two_e_tc_non_hermit_integral_beta (j,i) - !stop - endif - - norm += dabs(two_e_tc_non_hermit_integral_seq_beta(j,i)) - diff_tot += diff - enddo - enddo - - print *, ' diff tot b = ', diff_tot / norm - print *, ' norm b = ', norm - print *, ' ' - - ! --- - - return - -end - -! --- - ->>>>>>> 92a4e33f8a21717cab0c0e4f8412ed6903afb04a