From 10b461f5a21306299639c0c0f790847e004620ed Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 02:10:45 +0100 Subject: [PATCH] tc_scf: combined --- src/tc_scf/diago_vartcfock.irp.f | 96 +++ src/tc_scf/diis_tcscf.irp.f | 40 +- src/tc_scf/fock_tc.irp.f | 48 +- src/tc_scf/fock_vartc.irp.f | 287 ++++++++ src/tc_scf/rh_tcscf_diis.irp.f | 72 +- src/tc_scf/rh_vartcscf_simple.irp.f | 89 +++ src/tc_scf/tc_scf.irp.f | 1 + src/tc_scf/tc_scf_energy.irp.f | 29 + src/tc_scf/test_int.irp.f | 1059 +++++++++++++++++++++++++++ 9 files changed, 1659 insertions(+), 62 deletions(-) create mode 100644 src/tc_scf/diago_vartcfock.irp.f create mode 100644 src/tc_scf/fock_vartc.irp.f create mode 100644 src/tc_scf/rh_vartcscf_simple.irp.f create mode 100644 src/tc_scf/test_int.irp.f diff --git a/src/tc_scf/diago_vartcfock.irp.f b/src/tc_scf/diago_vartcfock.irp.f new file mode 100644 index 00000000..0c881dcb --- /dev/null +++ b/src/tc_scf/diago_vartcfock.irp.f @@ -0,0 +1,96 @@ + +! --- + +BEGIN_PROVIDER [ double precision, fock_vartc_eigvec_mo, (mo_num, mo_num)] + + implicit none + + integer :: i, j + integer :: liwork, lwork, n, info + integer, allocatable :: iwork(:) + double precision, allocatable :: work(:), F(:,:), F_save(:,:) + double precision, allocatable :: diag(:) + + PROVIDE mo_r_coef + PROVIDE Fock_matrix_vartc_mo_tot + + allocate( F(mo_num,mo_num), F_save(mo_num,mo_num) ) + allocate (diag(mo_num) ) + + do j = 1, mo_num + do i = 1, mo_num + F(i,j) = Fock_matrix_vartc_mo_tot(i,j) + enddo + enddo + + ! Insert level shift here + do i = elec_beta_num+1, elec_alpha_num + F(i,i) += 0.5d0 * level_shift_tcscf + enddo + do i = elec_alpha_num+1, mo_num + F(i,i) += level_shift_tcscf + enddo + + n = mo_num + lwork = 1+6*n + 2*n*n + liwork = 3 + 5*n + + allocate(work(lwork)) + allocate(iwork(liwork) ) + + lwork = -1 + liwork = -1 + + F_save = F + call dsyevd('V', 'U', mo_num, F, size(F, 1), diag, work, lwork, iwork, liwork, info) + + if (info /= 0) then + print *, irp_here//' DSYEVD failed : ', info + stop 1 + endif + lwork = int(work(1)) + liwork = iwork(1) + deallocate(iwork) + deallocate(work) + + allocate(work(lwork)) + allocate(iwork(liwork) ) + call dsyevd('V', 'U', mo_num, F, size(F, 1), diag, work, lwork, iwork, liwork, info) + deallocate(iwork) + + if (info /= 0) then + F = F_save + call dsyev('V', 'L', mo_num, F, size(F, 1), diag, work, lwork, info) + + if (info /= 0) then + print *, irp_here//' DSYEV failed : ', info + stop 1 + endif + endif + + do i = 1, mo_num + do j = 1, mo_num + fock_vartc_eigvec_mo(j,i) = F(j,i) + enddo + enddo + + deallocate(work, F, F_save, diag) + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, fock_vartc_eigvec_ao, (ao_num, mo_num)] + + implicit none + + PROVIDE mo_r_coef + + call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & + , mo_r_coef, size(mo_r_coef, 1), fock_vartc_eigvec_mo, size(fock_vartc_eigvec_mo, 1) & + , 0.d0, fock_vartc_eigvec_ao, size(fock_vartc_eigvec_ao, 1)) + +END_PROVIDER + +! --- + diff --git a/src/tc_scf/diis_tcscf.irp.f b/src/tc_scf/diis_tcscf.irp.f index ff1077f5..0b08f784 100644 --- a/src/tc_scf/diis_tcscf.irp.f +++ b/src/tc_scf/diis_tcscf.irp.f @@ -1,17 +1,3 @@ -! --- - -BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero_TCSCF ] - - implicit none - - if(threshold_DIIS_TCSCF == 0.d0) then - threshold_DIIS_nonzero_TCSCF = dsqrt(thresh_tcscf) - else - threshold_DIIS_nonzero_TCSCF = threshold_DIIS_TCSCF - endif - ASSERT(threshold_DIIS_nonzero_TCSCF >= 0.d0) - -END_PROVIDER ! --- @@ -100,13 +86,30 @@ END_PROVIDER BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)] implicit none + integer :: i, j double precision, allocatable :: tmp(:,:) + double precision, allocatable :: F(:,:) + + allocate(F(ao_num,ao_num)) + if(var_tc) then + do i = 1, ao_num + do j = 1, ao_num + F(j,i) = Fock_matrix_vartc_ao_tot(j,i) + enddo + enddo + else + do i = 1, ao_num + do j = 1, ao_num + F(j,i) = Fock_matrix_tc_ao_tot(j,i) + enddo + enddo + endif allocate(tmp(ao_num,ao_num)) ! F x Q - call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & - , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), Q_matrix, size(Q_matrix, 1) & + call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & + , F, size(F, 1), Q_matrix, size(Q_matrix, 1) & , 0.d0, tmp, size(tmp, 1) ) ! F x Q x S @@ -121,11 +124,12 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)] , 0.d0, tmp, size(tmp, 1) ) ! F x Q x S - S x Q x F - call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0 & - , tmp, size(tmp, 1), Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & + call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0 & + , tmp, size(tmp, 1), F, size(F, 1) & , 1.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) ) deallocate(tmp) + deallocate(F) END_PROVIDER diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index e21938de..1d651c4e 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -18,6 +18,8 @@ double precision :: density, density_a, density_b double precision :: t0, t1 + !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 @@ -31,6 +33,15 @@ 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_seq_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j) ! rho(l,j) * < k l| T | i j> @@ -45,6 +56,8 @@ enddo enddo + !call wall_time(t1) + !print*, ' wall time for two_e_tc_non_hermit_integral_seq after = ', t1 - t0 END_PROVIDER @@ -67,6 +80,8 @@ 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 @@ -112,6 +127,8 @@ 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 @@ -156,6 +173,14 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] if(bi_ortho) then + !allocate(tmp(ao_num,ao_num)) + !tmp = Fock_matrix_tc_ao_alpha + !if(three_body_h_tc) then + ! tmp += fock_3e_uhf_ao_a + !endif + !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1)) + !deallocate(tmp) + 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 @@ -184,6 +209,14 @@ 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 @@ -216,10 +249,6 @@ END_PROVIDER do k = elec_beta_num+1, elec_alpha_num grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i))) grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k))) - !grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i)) - !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k)) - !grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i) - !grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k) enddo enddo @@ -227,10 +256,6 @@ END_PROVIDER do k = elec_alpha_num+1, mo_num grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i))) grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k))) - !grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i)) - !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k)) - grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i) - grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k) enddo enddo @@ -238,15 +263,10 @@ END_PROVIDER do k = elec_alpha_num+1, mo_num grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i))) grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k))) - !grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i)) - !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k)) - grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i) - grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k) enddo enddo - !grad_non_hermit = dsqrt(grad_non_hermit_left) + dsqrt(grad_non_hermit_right) - grad_non_hermit = grad_non_hermit_left + grad_non_hermit_right + grad_non_hermit = max(grad_non_hermit_left, grad_non_hermit_right) END_PROVIDER diff --git a/src/tc_scf/fock_vartc.irp.f b/src/tc_scf/fock_vartc.irp.f new file mode 100644 index 00000000..03899b07 --- /dev/null +++ b/src/tc_scf/fock_vartc.irp.f @@ -0,0 +1,287 @@ + +! --- + + BEGIN_PROVIDER [ double precision, two_e_vartc_integral_alpha, (ao_num, ao_num)] +&BEGIN_PROVIDER [ double precision, two_e_vartc_integral_beta , (ao_num, ao_num)] + + implicit none + integer :: i, j, k, l + double precision :: density, density_a, density_b, I_coul, I_kjli + double precision :: t0, t1 + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + + two_e_vartc_integral_alpha = 0.d0 + two_e_vartc_integral_beta = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli) & + !$OMP SHARED (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_vartc_tot, & + !$OMP two_e_vartc_integral_alpha, two_e_vartc_integral_beta) + + allocate(tmp_a(ao_num,ao_num), tmp_b(ao_num,ao_num)) + tmp_a = 0.d0 + tmp_b = 0.d0 + + !$OMP DO + do j = 1, ao_num + do l = 1, ao_num + density_a = TCSCF_density_matrix_ao_alpha(l,j) + density_b = TCSCF_density_matrix_ao_beta (l,j) + density = density_a + density_b + do i = 1, ao_num + do k = 1, ao_num + + I_coul = density * ao_two_e_vartc_tot(k,i,l,j) + I_kjli = ao_two_e_vartc_tot(k,j,l,i) + + tmp_a(k,i) += I_coul - density_a * I_kjli + tmp_b(k,i) += I_coul - density_b * I_kjli + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do i = 1, ao_num + do j = 1, ao_num + two_e_vartc_integral_alpha(j,i) += tmp_a(j,i) + two_e_vartc_integral_beta (j,i) += tmp_b(j,i) + enddo + enddo + !$OMP END CRITICAL + + deallocate(tmp_a, tmp_b) + !$OMP END PARALLEL + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_alpha, (ao_num, ao_num)] + + implicit none + + Fock_matrix_vartc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_vartc_integral_alpha + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_beta, (ao_num, ao_num)] + + implicit none + + Fock_matrix_vartc_ao_beta = ao_one_e_integrals_tc_tot + two_e_vartc_integral_beta + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_alpha, (mo_num, mo_num) ] + + implicit none + + call ao_to_mo_bi_ortho( Fock_matrix_vartc_ao_alpha, size(Fock_matrix_vartc_ao_alpha, 1) & + , Fock_matrix_vartc_mo_alpha, size(Fock_matrix_vartc_mo_alpha, 1) ) + if(three_body_h_tc) then + Fock_matrix_vartc_mo_alpha += fock_3e_uhf_mo_a + endif + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_beta, (mo_num,mo_num) ] + + implicit none + + call ao_to_mo_bi_ortho( Fock_matrix_vartc_ao_beta, size(Fock_matrix_vartc_ao_beta, 1) & + , Fock_matrix_vartc_mo_beta, size(Fock_matrix_vartc_mo_beta, 1) ) + if(three_body_h_tc) then + Fock_matrix_vartc_mo_beta += fock_3e_uhf_mo_b + endif + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, grad_vartc] + + implicit none + integer :: i, k + double precision :: grad_left, grad_right + + grad_left = 0.d0 + grad_right = 0.d0 + + do i = 1, elec_beta_num ! doc --> SOMO + do k = elec_beta_num+1, elec_alpha_num + grad_left = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i))) + grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k))) + enddo + enddo + + do i = 1, elec_beta_num ! doc --> virt + do k = elec_alpha_num+1, mo_num + grad_left = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i))) + grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k))) + enddo + enddo + + do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt + do k = elec_alpha_num+1, mo_num + grad_left = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i))) + grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k))) + enddo + enddo + + grad_vartc = grad_left + grad_right + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_tot, (ao_num, ao_num) ] + + implicit none + + call mo_to_ao_bi_ortho( Fock_matrix_vartc_mo_tot, size(Fock_matrix_vartc_mo_tot, 1) & + , Fock_matrix_vartc_ao_tot, size(Fock_matrix_vartc_ao_tot, 1) ) + +END_PROVIDER + +! --- + + BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_tot, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_diag_mo_tot, (mo_num)] + + implicit none + integer :: i, j, n + + if(elec_alpha_num == elec_beta_num) then + Fock_matrix_vartc_mo_tot = Fock_matrix_vartc_mo_alpha + else + + do j = 1, elec_beta_num + ! F-K + do i = 1, elec_beta_num !CC + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))& + - (Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j)) + enddo + ! F+K/2 + do i = elec_beta_num+1, elec_alpha_num !CA + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j)) + enddo + ! F + do i = elec_alpha_num+1, mo_num !CV + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j)) + enddo + enddo + + do j = elec_beta_num+1, elec_alpha_num + ! F+K/2 + do i = 1, elec_beta_num !AC + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j)) + enddo + ! F + do i = elec_beta_num+1, elec_alpha_num !AA + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j)) + enddo + ! F-K/2 + do i = elec_alpha_num+1, mo_num !AV + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j)) + enddo + enddo + + do j = elec_alpha_num+1, mo_num + ! F + do i = 1, elec_beta_num !VC + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j)) + enddo + ! F-K/2 + do i = elec_beta_num+1, elec_alpha_num !VA + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j)) + enddo + ! F+K + do i = elec_alpha_num+1, mo_num !VV + Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j)) & + + (Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j)) + enddo + enddo + if(three_body_h_tc)then + ! C-O + do j = 1, elec_beta_num + do i = elec_beta_num+1, elec_alpha_num + Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo + enddo + ! C-V + do j = 1, elec_beta_num + do i = elec_alpha_num+1, mo_num + Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo + enddo + ! O-V + do j = elec_beta_num+1, elec_alpha_num + do i = elec_alpha_num+1, mo_num + Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo + enddo + endif + + endif + + do i = 1, mo_num + Fock_matrix_vartc_diag_mo_tot(i) = Fock_matrix_vartc_mo_tot(i,i) + enddo + + if(frozen_orb_scf)then + integer :: iorb, jorb + do i = 1, n_core_orb + iorb = list_core(i) + do j = 1, n_act_orb + jorb = list_act(j) + Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + + if(no_oa_or_av_opt)then + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + + !call check_sym(Fock_matrix_vartc_mo_tot, mo_num) + !do i = 1, mo_num + ! write(*,'(100(F15.8, I4))') Fock_matrix_vartc_mo_tot(i,:) + !enddo + +END_PROVIDER + +! --- + diff --git a/src/tc_scf/rh_tcscf_diis.irp.f b/src/tc_scf/rh_tcscf_diis.irp.f index 306c78b3..645742c8 100644 --- a/src/tc_scf/rh_tcscf_diis.irp.f +++ b/src/tc_scf/rh_tcscf_diis.irp.f @@ -21,7 +21,7 @@ subroutine rh_tcscf_diis() dim_DIIS = 0 g_delta_th = 1d0 er_delta_th = 1d0 - rate_th = 100.d0 !0.01d0 !0.2d0 + rate_th = 0.1d0 allocate(mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num)) mo_l_coef_save = 0.d0 @@ -38,17 +38,25 @@ subroutine rh_tcscf_diis() PROVIDE level_shift_TCSCF PROVIDE mo_l_coef mo_r_coef - write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & - '====', '================', '================', '================', '================', '================' & - , '================', '================', '================', '====', '========' + !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + ! '====', '================', '================', '================', '================', '================' & + ! , '================', '================', '================', '====', '========' + !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + ! ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' & + ! , ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' + !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + ! '====', '================', '================', '================', '================', '================' & + ! , '================', '================', '================', '====', '========' - write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '====', '========' + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' & - , ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' - - write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + , ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & '====', '================', '================', '================', '================', '================' & - , '================', '================', '================', '====', '========' + , '================', '================', '====', '========' ! first iteration (HF orbitals) @@ -61,23 +69,26 @@ subroutine rh_tcscf_diis() if(three_body_h_tc) then etc_3e = diag_three_elem_hf endif - tc_grad = grad_non_hermit + !tc_grad = grad_non_hermit er_DIIS = maxval(abs(FQS_SQF_mo)) e_delta = dabs(etc_tot - e_save) e_save = etc_tot - g_save = tc_grad + !g_save = tc_grad er_save = er_DIIS call wall_time(t1) - write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & - it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + !write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + ! it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 ! --- PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot - do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. threshold_DIIS_nonzero_TCSCF)) + !do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. dsqrt(thresh_tcscf))) + do while(er_DIIS .gt. dsqrt(thresh_tcscf)) call wall_time(t0) @@ -118,12 +129,10 @@ subroutine rh_tcscf_diis() ! --- - g_delta = grad_non_hermit - g_save + !g_delta = grad_non_hermit - g_save er_delta = maxval(abs(FQS_SQF_mo)) - er_save - !if((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1)) then - if((g_delta > rate_th * g_delta_th) .and. (it > 1)) then - !if((g_delta > 0.d0) .and. (it > 1)) then + if((er_delta > rate_th * er_save) .and. (it > 1)) then Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_dim_DIIS) call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & @@ -140,15 +149,16 @@ subroutine rh_tcscf_diis() ! --- - g_delta = grad_non_hermit - g_save + !g_delta = grad_non_hermit - g_save er_delta = maxval(abs(FQS_SQF_mo)) - er_save mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) - !do while((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1)) - do while((g_delta > rate_th * g_delta_th) .and. (it > 1)) - print *, ' big or bad step : ', g_delta, rate_th * g_delta_th + do while((er_delta > rate_th * er_save) .and. (it > 1)) + print *, ' big or bad step ' + !print *, g_delta , rate_th * g_save + print *, er_delta, rate_th * er_save mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) @@ -165,7 +175,7 @@ subroutine rh_tcscf_diis() !call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) TOUCH mo_l_coef mo_r_coef - g_delta = grad_non_hermit - g_save + !g_delta = grad_non_hermit - g_save er_delta = maxval(abs(FQS_SQF_mo)) - er_save if(level_shift_TCSCF - level_shift_save > 40.d0) then @@ -189,25 +199,27 @@ subroutine rh_tcscf_diis() if(three_body_h_tc) then etc_3e = diag_three_elem_hf endif - tc_grad = grad_non_hermit + !tc_grad = grad_non_hermit er_DIIS = maxval(abs(FQS_SQF_mo)) e_delta = dabs(etc_tot - e_save) - g_delta = tc_grad - g_save + !g_delta = tc_grad - g_save er_delta = er_DIIS - er_save e_save = etc_tot - g_save = tc_grad + !g_save = tc_grad level_shift_save = level_shift_TCSCF er_save = er_DIIS - g_delta_th = dabs(tc_grad) ! g_delta) + !g_delta_th = dabs(tc_grad) ! g_delta) er_delta_th = dabs(er_DIIS) !er_delta) call wall_time(t1) - write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & - it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + !write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + ! it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 - if(g_delta .lt. 0.d0) then + if(er_delta .lt. 0.d0) then call ezfio_set_tc_scf_bitc_energy(etc_tot) call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) diff --git a/src/tc_scf/rh_vartcscf_simple.irp.f b/src/tc_scf/rh_vartcscf_simple.irp.f new file mode 100644 index 00000000..ecb0709e --- /dev/null +++ b/src/tc_scf/rh_vartcscf_simple.irp.f @@ -0,0 +1,89 @@ +! --- + +subroutine rh_vartcscf_simple() + + implicit none + integer :: i, j, it, dim_DIIS + double precision :: t0, t1 + double precision :: e_save, e_delta, rho_delta + double precision :: etc_tot, etc_1e, etc_2e, etc_3e + double precision :: er_DIIS + + + it = 0 + e_save = 0.d0 + dim_DIIS = 0 + + ! --- + + PROVIDE level_shift_tcscf + PROVIDE mo_r_coef + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '====', '========' + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' & + , ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '====', '========' + + + ! first iteration (HF orbitals) + call wall_time(t0) + + etc_tot = VARTC_HF_energy + etc_1e = VARTC_HF_one_e_energy + etc_2e = VARTC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + etc_3e = diag_three_elem_hf + endif + er_DIIS = maxval(abs(FQS_SQF_mo)) + e_delta = dabs(etc_tot - e_save) + e_save = etc_tot + + call wall_time(t1) + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + + do while(er_DIIS .gt. dsqrt(thresh_tcscf)) + call wall_time(t0) + + it += 1 + if(it > n_it_tcscf_max) then + print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max + stop + endif + + mo_r_coef = fock_vartc_eigvec_ao + mo_l_coef = mo_r_coef + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + TOUCH mo_l_coef mo_r_coef + + etc_tot = VARTC_HF_energy + etc_1e = VARTC_HF_one_e_energy + etc_2e = VARTC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + etc_3e = diag_three_elem_hf + endif + er_DIIS = maxval(abs(FQS_SQF_mo)) + e_delta = dabs(etc_tot - e_save) + e_save = etc_tot + + call ezfio_set_tc_scf_bitc_energy(etc_tot) + + call wall_time(t1) + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + enddo + + print *, ' VAR-TCSCF Simple converged !' + +end + +! --- + diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index deaf8d82..187750ff 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -73,3 +73,4 @@ subroutine create_guess() end subroutine create_guess ! --- + diff --git a/src/tc_scf/tc_scf_energy.irp.f b/src/tc_scf/tc_scf_energy.irp.f index 611b8b4c..c3de0322 100644 --- a/src/tc_scf/tc_scf_energy.irp.f +++ b/src/tc_scf/tc_scf_energy.irp.f @@ -30,5 +30,34 @@ END_PROVIDER +! --- + + BEGIN_PROVIDER [ double precision, VARTC_HF_energy] +&BEGIN_PROVIDER [ double precision, VARTC_HF_one_e_energy] +&BEGIN_PROVIDER [ double precision, VARTC_HF_two_e_energy] + + implicit none + integer :: i, j + + PROVIDE mo_r_coef + + VARTC_HF_energy = nuclear_repulsion + VARTC_HF_one_e_energy = 0.d0 + VARTC_HF_two_e_energy = 0.d0 + + do j = 1, ao_num + do i = 1, ao_num + VARTC_HF_two_e_energy += 0.5d0 * ( two_e_vartc_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) & + + two_e_vartc_integral_beta (i,j) * TCSCF_density_matrix_ao_beta (i,j) ) + VARTC_HF_one_e_energy += ao_one_e_integrals_tc_tot(i,j) & + * (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (i,j) ) + enddo + enddo + + VARTC_HF_energy += VARTC_HF_one_e_energy + VARTC_HF_two_e_energy + VARTC_HF_energy += diag_three_elem_hf + +END_PROVIDER + ! --- diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f new file mode 100644 index 00000000..0866cdaf --- /dev/null +++ b/src/tc_scf/test_int.irp.f @@ -0,0 +1,1059 @@ +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() + + call test_tc_grad_square_ao_test() + + PROVIDE TC_HF_energy VARTC_HF_energy + print *, ' TC_HF_energy = ', TC_HF_energy + print *, ' VARTC_HF_energy = ', VARTC_HF_energy + +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 +! do i = 1, ng_fit_jast +! print*,expo_gauss_1_erf_x_2(i),coef_gauss_1_erf_x_2(i) +! enddo +! provide tc_grad_square_ao_test +! provide tc_grad_and_lapl_ao_test +! provide int2_u_grad1u_x_j1b2_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 + +! --- + +subroutine test_tc_grad_square_ao_test() + + implicit none + integer :: i, j, k, l + double precision :: diff_tot, diff, thr_ih, norm + + print*, ' test_tc_grad_square_ao_test ' + + thr_ih = 1d-7 + + PROVIDE tc_grad_square_ao_test tc_grad_square_ao_test_ref + + 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_test(l,k,j,i) - tc_grad_square_ao_test_ref(l,k,j,i)) + if(diff .gt. thr_ih) then + print *, ' difference on ', l, k, j, i + print *, ' new : ', tc_grad_square_ao_test (l,k,j,i) + print *, ' ref : ', tc_grad_square_ao_test_ref(l,k,j,i) + !stop + endif + + norm += dabs(tc_grad_square_ao_test_ref(l,k,j,i)) + diff_tot += diff + enddo + enddo + enddo + enddo + + print *, ' diff tot = ', diff_tot / norm + print *, ' norm = ', norm + print *, ' ' + + return +end + +! --- + +