From 10b461f5a21306299639c0c0f790847e004620ed Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 02:10:45 +0100 Subject: [PATCH 1/6] 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 + +! --- + + From 5b58b062d9fb7173135c6882115e0664578f8377 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 02:19:15 +0100 Subject: [PATCH 2/6] non_h_ints_mu: combined --- src/non_h_ints_mu/grad_squared.irp.f | 139 +++++++------ src/non_h_ints_mu/grad_squared_manu.irp.f | 174 ++++++++++++++++- src/non_h_ints_mu/new_grad_tc.irp.f | 222 +++++++++++++-------- src/non_h_ints_mu/new_grad_tc_manu.irp.f | 227 +++++++++++++--------- src/non_h_ints_mu/total_tc_int.irp.f | 58 +++++- 5 files changed, 579 insertions(+), 241 deletions(-) diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f index ff3d11f3..7925fa7c 100644 --- a/src/non_h_ints_mu/grad_squared.irp.f +++ b/src/non_h_ints_mu/grad_squared.irp.f @@ -299,7 +299,6 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_ END_PROVIDER -! --- ! --- BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ] @@ -364,70 +363,100 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao integer :: ipoint, i, j, k, l double precision :: weight1, ao_ik_r, ao_i_r double precision :: time0, time1 - double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:), tmp(:,:,:) + double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:) print*, ' providing tc_grad_square_ao ...' call wall_time(time0) - allocate(ac_mat(ao_num,ao_num,ao_num,ao_num), b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid)) + if(read_tc_integ) then - b_mat = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, k, ipoint) & - !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) - !$OMP DO SCHEDULE (static) - do i = 1, ao_num - do k = 1, ao_num - do ipoint = 1, n_points_final_grid - b_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - tmp = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (j, l, ipoint) & - !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) - !$OMP DO SCHEDULE (static) - do ipoint = 1, n_points_final_grid - do j = 1, ao_num - do l = 1, ao_num - tmp(l,j,ipoint) = u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) + 0.5d0 * grad12_j12(l,j,ipoint) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - - ac_mat = 0.d0 - call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & - , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & - , 1.d0, ac_mat, ao_num*ao_num) - deallocate(tmp, b_mat) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, k, l) & - !$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do l = 1, ao_num + open(unit=11, form="unformatted", file='tc_grad_square_ao', action="read") do i = 1, ao_num - do k = 1, ao_num - tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + read(11) tc_grad_square_ao(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + + else + + allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid)) + + b_mat = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint) & + !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + b_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k) enddo enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL + + tmp = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, l, ipoint) & + !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do l = 1, ao_num + tmp(l,j,ipoint) = u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) + 0.5d0 * grad12_j12(l,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + tc_grad_square_ao = 0.d0 + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & + , 1.d0, tc_grad_square_ao, ao_num*ao_num) + deallocate(tmp, b_mat) + + call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num) + + !!$OMP PARALLEL & + !!$OMP DEFAULT (NONE) & + !!$OMP PRIVATE (i, j, k, l) & + !!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num) + !!$OMP DO SCHEDULE (static) + ! do j = 1, ao_num + ! do l = 1, ao_num + ! do i = 1, ao_num + ! do k = 1, ao_num + ! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + ! enddo + ! enddo + ! enddo + ! enddo + !!$OMP END DO + !!$OMP END PARALLEL + endif - deallocate(ac_mat) + if(write_tc_integ) then + open(unit=11, form="unformatted", file='tc_grad_square_ao', action="write") + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + write(11) tc_grad_square_ao(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + endif call wall_time(time1) print*, ' Wall time for tc_grad_square_ao = ', time1 - time0 diff --git a/src/non_h_ints_mu/grad_squared_manu.irp.f b/src/non_h_ints_mu/grad_squared_manu.irp.f index 180c9588..cb9e15c4 100644 --- a/src/non_h_ints_mu/grad_squared_manu.irp.f +++ b/src/non_h_ints_mu/grad_squared_manu.irp.f @@ -11,11 +11,177 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu integer :: ipoint, i, j, k, l double precision :: weight1, ao_ik_r, ao_i_r,contrib,contrib2 double precision :: time0, time1 - double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:), tmp(:,:,:) + double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:) print*, ' providing tc_grad_square_ao_test ...' call wall_time(time0) + if(read_tc_integ) then + + open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="read") + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + read(11) tc_grad_square_ao_test(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + + else + + provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test + + allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid)) + + b_mat = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint) & + !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + b_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + tmp = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (j, l, ipoint) & + !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test) + !$OMP DO SCHEDULE (static) + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do l = 1, ao_num + tmp(l,j,ipoint) = u12sq_j1bsq_test(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b_test(l,j,ipoint) + 0.5d0 * grad12_j12_test(l,j,ipoint) + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + tc_grad_square_ao_test = 0.d0 + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid & + , 1.d0, tc_grad_square_ao_test, ao_num*ao_num) + deallocate(tmp, b_mat) + + call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num) + !do i = 1, ao_num + ! do j = 1, ao_num + ! do k = i, ao_num + + ! do l = max(j,k), ao_num + ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j)) + ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l) + ! end do + + ! !if (j.eq.k) then + ! ! do l = j+1, ao_num + ! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j)) + ! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l) + ! ! end do + ! !else + ! ! do l = j, ao_num + ! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j)) + ! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l) + ! ! enddo + ! !endif + + ! enddo + ! enddo + !enddo + !tc_grad_square_ao_test = 2.d0 * tc_grad_square_ao_test + ! !$OMP PARALLEL & + ! !$OMP DEFAULT (NONE) & + ! !$OMP PRIVATE (i, j, k, l) & + ! !$OMP SHARED (tc_grad_square_ao_test, ao_num) + ! !$OMP DO SCHEDULE (static) + ! integer :: ii + ! ii = 0 + ! do j = 1, ao_num + ! do l = 1, ao_num + ! do i = 1, ao_num + ! do k = 1, ao_num + ! if((i.lt.j) .and. (k.lt.l)) cycle + ! ii = ii + 1 + ! tc_grad_square_ao_test(k,i,l,j) = tc_grad_square_ao_test(k,i,l,j) + tc_grad_square_ao_test(l,j,k,i) + ! enddo + ! enddo + ! enddo + ! enddo + ! print *, ' ii =', ii + ! !$OMP END DO + ! !$OMP END PARALLEL + + ! !$OMP PARALLEL & + ! !$OMP DEFAULT (NONE) & + ! !$OMP PRIVATE (i, j, k, l) & + ! !$OMP SHARED (tc_grad_square_ao_test, ao_num) + ! !$OMP DO SCHEDULE (static) + ! do j = 1, ao_num + ! do l = 1, ao_num + ! do i = 1, j-1 + ! do k = 1, l-1 + ! ii = ii + 1 + ! tc_grad_square_ao_test(k,i,l,j) = tc_grad_square_ao_test(l,j,k,i) + ! enddo + ! enddo + ! enddo + ! enddo + ! print *, ' ii =', ii + ! print *, ao_num * ao_num * ao_num * ao_num + ! !$OMP END DO + ! !$OMP END PARALLEL + + endif + + if(write_tc_integ) then + open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="write") + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + write(11) tc_grad_square_ao_test(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + endif + + call wall_time(time1) + print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, ao_num, ao_num)] + + BEGIN_DOC + ! + ! tc_grad_square_ao_test_ref(k,i,l,j) = -1/2 + ! + END_DOC + + implicit none + integer :: ipoint, i, j, k, l + double precision :: weight1, ao_ik_r, ao_i_r,contrib,contrib2 + double precision :: time0, time1 + double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:), tmp(:,:,:) + + print*, ' providing tc_grad_square_ao_test_ref ...' + call wall_time(time0) + provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test allocate(ac_mat(ao_num,ao_num,ao_num,ao_num), b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid)) @@ -61,13 +227,13 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i, j, k, l) & - !$OMP SHARED (ac_mat, tc_grad_square_ao_test, ao_num) + !$OMP SHARED (ac_mat, tc_grad_square_ao_test_ref, ao_num) !$OMP DO SCHEDULE (static) do j = 1, ao_num do l = 1, ao_num do i = 1, ao_num do k = 1, ao_num - tc_grad_square_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + tc_grad_square_ao_test_ref(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) enddo enddo enddo @@ -78,7 +244,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu deallocate(ac_mat) call wall_time(time1) - print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0 + print*, ' Wall time for tc_grad_square_ao_test_ref = ', time1 - time0 END_PROVIDER diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f index 854789bd..a15f690a 100644 --- a/src/non_h_ints_mu/new_grad_tc.irp.f +++ b/src/non_h_ints_mu/new_grad_tc.irp.f @@ -25,7 +25,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_ END_DOC implicit none - integer :: ipoint, i, j + integer :: ipoint, i, j, m double precision :: time0, time1 double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 @@ -33,54 +33,76 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_ call wall_time(time0) PROVIDE j1b_type - - if(j1b_type .eq. 3) then - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) + if(read_tc_integ) then - tmp0 = 0.5d0 * v_1b(ipoint) - tmp_x = v_1b_grad(1,ipoint) - tmp_y = v_1b_grad(2,ipoint) - tmp_z = v_1b_grad(3,ipoint) - - do j = 1, ao_num - do i = 1, ao_num - - tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) - tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) - - int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x - int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y - int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z + open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="read") + do m = 1, 3 + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + read(11) int2_grad1_u12_ao(i,j,ipoint,m) + enddo + enddo enddo enddo - enddo + close(11) else - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) - - do j = 1, ao_num - do i = 1, ao_num - tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) - - int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,1) - int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,2) - int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,3) + if(j1b_type .eq. 3) then + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + tmp0 = 0.5d0 * v_1b(ipoint) + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint) + int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x + int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y + int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z + enddo enddo enddo - enddo + else + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) - int2_grad1_u12_ao *= 0.5d0 + int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,1) + int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,2) + int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,3) + enddo + enddo + enddo + int2_grad1_u12_ao *= 0.5d0 + endif endif + if(write_tc_integ) then + open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="write") + do m = 1, 3 + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + write(11) int2_grad1_u12_ao(i,j,ipoint,m) + enddo + enddo + enddo + enddo + close(11) + endif + call wall_time(time1) print*, ' Wall time for int2_grad1_u12_ao = ', time1 - time0 @@ -290,65 +312,95 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num, integer :: ipoint, i, j, k, l, m double precision :: weight1, ao_k_r, ao_i_r double precision :: time0, time1 - double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:,:) + double precision, allocatable :: b_mat(:,:,:,:) print*, ' providing tc_grad_and_lapl_ao ...' call wall_time(time0) - allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num)) + if(read_tc_integ) then - b_mat = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & - !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, & - !$OMP ao_num, n_points_final_grid, final_weight_at_r_vector) - !$OMP DO SCHEDULE (static) - do i = 1, ao_num - do k = 1, ao_num - do ipoint = 1, n_points_final_grid - - weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) - ao_i_r = aos_in_r_array_transp(ipoint,i) - ao_k_r = aos_in_r_array_transp(ipoint,k) - - b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) - b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) - b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - ac_mat = 0.d0 - do m = 1, 3 - call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & - , int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & - , 1.d0, ac_mat, ao_num*ao_num) - - enddo - deallocate(b_mat) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, k, l) & - !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao, ao_num) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do l = 1, ao_num + open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao', action="read") do i = 1, ao_num - do k = 1, ao_num - tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) - !tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + read(11) tc_grad_and_lapl_ao(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + + else + + allocate(b_mat(n_points_final_grid,ao_num,ao_num,3)) + + b_mat = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & + !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, & + !$OMP ao_num, n_points_final_grid, final_weight_at_r_vector) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + + weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) + ao_i_r = aos_in_r_array_transp(ipoint,i) + ao_k_r = aos_in_r_array_transp(ipoint,k) + + b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) + b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) + b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) enddo enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL + + tc_grad_and_lapl_ao = 0.d0 + do m = 1, 3 + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & + , 1.d0, tc_grad_and_lapl_ao, ao_num*ao_num) + + enddo + deallocate(b_mat) + + call sum_A_At(tc_grad_and_lapl_ao(1,1,1,1), ao_num*ao_num) + ! !$OMP PARALLEL & + ! !$OMP DEFAULT (NONE) & + ! !$OMP PRIVATE (i, j, k, l) & + ! !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao, ao_num) + ! !$OMP DO SCHEDULE (static) + ! do j = 1, ao_num + ! do l = 1, ao_num + ! do i = 1, ao_num + ! do k = 1, ao_num + ! tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + ! enddo + ! enddo + ! enddo + ! enddo + ! !$OMP END DO + ! !$OMP END PARALLEL - deallocate(ac_mat) + endif + + if(write_tc_integ) then + open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao', action="write") + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + write(11) tc_grad_and_lapl_ao(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + endif call wall_time(time1) print*, ' Wall time for tc_grad_and_lapl_ao = ', time1 - time0 diff --git a/src/non_h_ints_mu/new_grad_tc_manu.irp.f b/src/non_h_ints_mu/new_grad_tc_manu.irp.f index 4d85e061..47b05e52 100644 --- a/src/non_h_ints_mu/new_grad_tc_manu.irp.f +++ b/src/non_h_ints_mu/new_grad_tc_manu.irp.f @@ -24,7 +24,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po END_DOC implicit none - integer :: ipoint, i, j + integer :: ipoint, i, j, m double precision :: time0, time1 double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2 @@ -32,52 +32,73 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po call wall_time(time0) PROVIDE j1b_type - - if(j1b_type .eq. 3) then - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) + if(read_tc_integ) then - tmp0 = 0.5d0 * v_1b(ipoint) - tmp_x = v_1b_grad(1,ipoint) - tmp_y = v_1b_grad(2,ipoint) - tmp_z = v_1b_grad(3,ipoint) - - do j = 1, ao_num - do i = 1, ao_num - - tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint) - tmp2 = v_ij_u_cst_mu_j1b_test(i,j,ipoint) - - int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1) - tmp2 * tmp_x - int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2) - tmp2 * tmp_y - int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3) - tmp2 * tmp_z + open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="read") + do m = 1, 3 + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + read(11) int2_grad1_u12_ao_test(i,j,ipoint,m) + enddo + enddo enddo enddo - enddo + close(11) else - - do ipoint = 1, n_points_final_grid - x = final_grid_points(1,ipoint) - y = final_grid_points(2,ipoint) - z = final_grid_points(3,ipoint) - - do j = 1, ao_num - do i = 1, ao_num - tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) - - int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,1) - int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,2) - int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,3) + + if(j1b_type .eq. 3) then + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + tmp0 = 0.5d0 * v_1b(ipoint) + tmp_x = v_1b_grad(1,ipoint) + tmp_y = v_1b_grad(2,ipoint) + tmp_z = v_1b_grad(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint) + tmp2 = v_ij_u_cst_mu_j1b_test(i,j,ipoint) + int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1) - tmp2 * tmp_x + int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2) - tmp2 * tmp_y + int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3) - tmp2 * tmp_z + enddo enddo enddo - enddo + else + do ipoint = 1, n_points_final_grid + x = final_grid_points(1,ipoint) + y = final_grid_points(2,ipoint) + z = final_grid_points(3,ipoint) + do j = 1, ao_num + do i = 1, ao_num + tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint) + int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,1) + int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,2) + int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,3) + enddo + enddo + enddo + int2_grad1_u12_ao_test *= 0.5d0 + endif - int2_grad1_u12_ao_test *= 0.5d0 + endif + if(write_tc_integ) then + open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="write") + do m = 1, 3 + do ipoint = 1, n_points_final_grid + do j = 1, ao_num + do i = 1, ao_num + write(11) int2_grad1_u12_ao_test(i,j,ipoint,m) + enddo + enddo + enddo + enddo + close(11) endif call wall_time(time1) @@ -109,61 +130,93 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_ print*, ' providing tc_grad_and_lapl_ao_test ...' call wall_time(time0) - provide int2_grad1_u12_ao_test - - allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num)) - - b_mat = 0.d0 - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & - !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, & - !$OMP ao_num, n_points_final_grid, final_weight_at_r_vector) - !$OMP DO SCHEDULE (static) - do i = 1, ao_num - do k = 1, ao_num - do ipoint = 1, n_points_final_grid - - weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) - ao_i_r = aos_in_r_array_transp(ipoint,i) - ao_k_r = aos_in_r_array_transp(ipoint,k) - - b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) - b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) - b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL - - ac_mat = 0.d0 - do m = 1, 3 - call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & - , int2_grad1_u12_ao_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & - , 1.d0, ac_mat, ao_num*ao_num) - - enddo - deallocate(b_mat) - - !$OMP PARALLEL & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i, j, k, l) & - !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num) - !$OMP DO SCHEDULE (static) - do j = 1, ao_num - do l = 1, ao_num + if(read_tc_integ) then + + open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="read") do i = 1, ao_num - do k = 1, ao_num - tc_grad_and_lapl_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + read(11) tc_grad_and_lapl_ao_test(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + + else + + provide int2_grad1_u12_ao_test + + allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num)) + + b_mat = 0.d0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) & + !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, & + !$OMP ao_num, n_points_final_grid, final_weight_at_r_vector) + !$OMP DO SCHEDULE (static) + do i = 1, ao_num + do k = 1, ao_num + do ipoint = 1, n_points_final_grid + + weight1 = 0.5d0 * final_weight_at_r_vector(ipoint) + ao_i_r = aos_in_r_array_transp(ipoint,i) + ao_k_r = aos_in_r_array_transp(ipoint,k) + + b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) + b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) + b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3)) enddo enddo enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + !$OMP END DO + !$OMP END PARALLEL + + ac_mat = 0.d0 + do m = 1, 3 + call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 & + , int2_grad1_u12_ao_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid & + , 1.d0, ac_mat, ao_num*ao_num) + + enddo + deallocate(b_mat) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, k, l) & + !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num) + !$OMP DO SCHEDULE (static) + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + tc_grad_and_lapl_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + deallocate(ac_mat) - deallocate(ac_mat) + endif + + if(write_tc_integ) then + open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="write") + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + write(11) tc_grad_and_lapl_ao_test(l,k,j,i) + enddo + enddo + enddo + enddo + close(11) + endif call wall_time(time1) print*, ' Wall time for tc_grad_and_lapl_ao_test = ', time1 - time0 diff --git a/src/non_h_ints_mu/total_tc_int.irp.f b/src/non_h_ints_mu/total_tc_int.irp.f index 81747553..c1e010c7 100644 --- a/src/non_h_ints_mu/total_tc_int.irp.f +++ b/src/non_h_ints_mu/total_tc_int.irp.f @@ -1,6 +1,44 @@ ! --- +BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)] + + implicit none + integer :: i, j, k, l + double precision :: wall1, wall0 + + print *, ' providing ao_vartc_int_chemist ...' + call wall_time(wall0) + + if(test_cycle_tc) then + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_vartc_int_chemist(k,i,l,j) = tc_grad_square_ao_test(k,i,l,j) + ao_two_e_coul(k,i,l,j) + enddo + enddo + enddo + enddo + else + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_vartc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) + enddo + enddo + enddo + enddo + endif + + call wall_time(wall1) + print *, ' wall time for ao_vartc_int_chemist ', wall1 - wall0 + +END_PROVIDER + +! --- + BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao_num)] implicit none @@ -11,17 +49,17 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao call wall_time(wall0) if(test_cycle_tc)then - ao_tc_int_chemist = ao_tc_int_chemist_test + ao_tc_int_chemist = ao_tc_int_chemist_test else - do j = 1, ao_num - do l = 1, ao_num - do i = 1, ao_num - do k = 1, ao_num - ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) - enddo - enddo - enddo - enddo + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j) + enddo + enddo + enddo + enddo endif call wall_time(wall1) From ec082f641b36711f7af8c3f9f84d0ea0485f1629 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 02:24:16 +0100 Subject: [PATCH 3/6] tc_keywords: combined --- src/tc_keywords/EZFIO.cfg | 56 ++++++++++++++++++++++++++------------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 5d5477bc..8830a86e 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -86,7 +86,7 @@ default: False type: Threshold doc: Threshold on the convergence of the Hartree Fock energy. interface: ezfio,provider,ocaml -default: 1.e-12 +default: 1.e-10 [n_it_tcscf_max] type: Strictly_positive_int @@ -94,6 +94,12 @@ doc: Maximum number of SCF iterations interface: ezfio,provider,ocaml default: 100 +[selection_tc] +type: integer +doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative +interface: ezfio,provider,ocaml +default: 0 + [j1b_pen] type: double precision doc: exponents of the 1-body Jastrow @@ -130,12 +136,30 @@ doc: nb of Gaussians used to fit Jastrow fcts interface: ezfio,provider,ocaml default: 20 +[max_dim_diis_tcscf] +type: integer +doc: Maximum size of the DIIS extrapolation procedure +interface: ezfio,provider,ocaml +default: 15 + +[level_shift_tcscf] +type: Positive_float +doc: Energy shift on the virtual MOs to improve TCSCF convergence +interface: ezfio,provider,ocaml +default: 0. + [tcscf_algorithm] type: character*(32) doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] interface: ezfio,provider,ocaml default: Simple +[im_thresh_tcscf] +type: Threshold +doc: Thresholds on the Imag part of energy +interface: ezfio,provider,ocaml +default: 1.e-7 + [test_cycle_tc] type: logical doc: If |true|, the integrals of the three-body jastrow are computed with cycles @@ -154,29 +178,23 @@ doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph interface: ezfio,provider,ocaml default: 1.e-6 -[max_dim_diis_tcscf] -type: integer -doc: Maximum size of the DIIS extrapolation procedure +[var_tc] +type: logical +doc: If |true|, use VAR-TC interface: ezfio,provider,ocaml -default: 15 +default: False -[threshold_diis_tcscf] -type: Threshold -doc: Threshold on the convergence of the DIIS error vector during a TCSCF calculation. If 0. is chosen, the square root of thresh_tcscf will be used. +[read_tc_integ] +type: logical +doc: If |true|, read integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao interface: ezfio,provider,ocaml -default: 0. +default: False -[level_shift_tcscf] -type: Positive_float -doc: Energy shift on the virtual MOs to improve TCSCF convergence +[write_tc_integ] +type: logical +doc: If |true|, write integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao interface: ezfio,provider,ocaml -default: 0. - -[im_thresh_tcscf] -type: Threshold -doc: Thresholds on the Imag part of energy -interface: ezfio,provider,ocaml -default: 1.e-7 +default: False [debug_tc_pt2] type: integer From 6e7ca02ed161ef5ae5de667cf33971b8ba117ea1 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 02:31:04 +0100 Subject: [PATCH 4/6] bi_ort_ints: combined --- src/bi_ort_ints/total_twoe_pot.irp.f | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f index e74c6d2a..78047d1b 100644 --- a/src/bi_ort_ints/total_twoe_pot.irp.f +++ b/src/bi_ort_ints/total_twoe_pot.irp.f @@ -1,4 +1,25 @@ + +! --- + +BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, ao_num) ] + + integer :: i, j, k, l + + provide j1b_type + provide mo_r_coef mo_l_coef + + do j = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do k = 1, ao_num + ao_two_e_vartc_tot(k,i,l,j) = ao_vartc_int_chemist(k,i,l,j) + enddo + enddo + enddo + enddo + +END_PROVIDER ! --- BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ] From dac3215a652c2142f5993a88c738567453c8b576 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 02:43:33 +0100 Subject: [PATCH 5/6] tc_scf v1 of combin --- src/tc_scf/rh_tcscf.irp.f | 336 -------------------------------------- src/utils/util.irp.f | 35 ++++ 2 files changed, 35 insertions(+), 336 deletions(-) delete mode 100644 src/tc_scf/rh_tcscf.irp.f diff --git a/src/tc_scf/rh_tcscf.irp.f b/src/tc_scf/rh_tcscf.irp.f deleted file mode 100644 index 0312df5f..00000000 --- a/src/tc_scf/rh_tcscf.irp.f +++ /dev/null @@ -1,336 +0,0 @@ -! --- - -subroutine rh_tcscf() - - BEGIN_DOC - ! - ! Roothaan-Hall algorithm for TC-SCF calculation - ! - END_DOC - - implicit none - - integer :: i, j - integer :: iteration_TCSCF, dim_DIIS, index_dim_DIIS - double precision :: energy_TCSCF, energy_TCSCF_1e, energy_TCSCF_2e, energy_TCSCF_3e, gradie_TCSCF - double precision :: energy_TCSCF_previous, delta_energy_TCSCF - double precision :: gradie_TCSCF_previous, delta_gradie_TCSCF - double precision :: max_error_DIIS_TCSCF - double precision :: level_shift_save - double precision :: delta_energy_tmp, delta_gradie_tmp - double precision, allocatable :: F_DIIS(:,:,:), e_DIIS(:,:,:) - double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:) - - logical, external :: qp_stop - - - !PROVIDE ao_md5 mo_occ - PROVIDE level_shift_TCSCF - - allocate( mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num) & - , F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), e_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF) ) - - F_DIIS = 0.d0 - e_DIIS = 0.d0 - mo_l_coef_save = 0.d0 - mo_r_coef_save = 0.d0 - - call write_time(6) - - ! --- - ! Initialize energies and density matrices - - energy_TCSCF_previous = TC_HF_energy - energy_TCSCF_1e = TC_HF_one_e_energy - energy_TCSCF_2e = TC_HF_two_e_energy - energy_TCSCF_3e = 0.d0 - if(three_body_h_tc) then - energy_TCSCF_3e = diag_three_elem_hf - endif - gradie_TCSCF_previous = grad_non_hermit - delta_energy_TCSCF = 1.d0 - delta_gradie_TCSCF = 1.d0 - iteration_TCSCF = 0 - dim_DIIS = 0 - max_error_DIIS_TCSCF = 1.d0 - - ! --- - - ! Start of main SCF loop - - PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot - - do while( (max_error_DIIS_TCSCF > threshold_DIIS_nonzero_TCSCF) .or. & - !(dabs(delta_energy_TCSCF) > thresh_TCSCF) .or. & - (dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF)) ) - - iteration_TCSCF += 1 - if(iteration_TCSCF > n_it_TCSCF_max) then - print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max - stop - endif - - dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF) - - ! --- - - if((tcscf_algorithm == 'DIIS') .and. (dabs(delta_energy_TCSCF) > 1.d-6)) then - - ! store Fock and error matrices at each iteration - index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1 - do j = 1, ao_num - do i = 1, ao_num - F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j) - e_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao(i,j) - enddo - enddo - - call extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), iteration_TCSCF, dim_DIIS) - - Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot - Fock_matrix_tc_ao_beta = 0.5d0 * Fock_matrix_tc_ao_tot - !TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta - - call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) & - , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) ) - call 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) ) - TOUCH Fock_matrix_tc_mo_alpha Fock_matrix_tc_mo_beta - endif - - ! --- - - mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) - mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) - TOUCH mo_l_coef mo_r_coef - - ! --- - - ! calculate error vectors - max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo)) - - ! --- - - delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous - delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous - - ! --- - - do while((delta_gradie_tmp > 1.d-7) .and. (iteration_TCSCF > 1)) - !do while((dabs(delta_energy_tmp) > 0.5d0) .and. (iteration_TCSCF > 1)) - print *, ' very big or bad step : ', delta_energy_tmp, delta_gradie_tmp - print *, ' TC level shift = ', level_shift_TCSCF - - 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) - - if(level_shift_TCSCF <= .1d0) then - level_shift_TCSCF = 1.d0 - else - level_shift_TCSCF = level_shift_TCSCF * 3.0d0 - endif - TOUCH mo_l_coef mo_r_coef level_shift_TCSCF - - mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) - mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) - TOUCH mo_l_coef mo_r_coef - - delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous - delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous - - if(level_shift_TCSCF - level_shift_save > 40.d0) then - level_shift_TCSCF = level_shift_save * 4.d0 - SOFT_TOUCH level_shift_TCSCF - exit - endif - - dim_DIIS = 0 - enddo -! print *, ' very big step : ', delta_energy_tmp -! print *, ' TC level shift = ', level_shift_TCSCF - - ! --- - - level_shift_TCSCF = 0.d0 - !level_shift_TCSCF = level_shift_TCSCF * 0.5d0 - SOFT_TOUCH level_shift_TCSCF - - gradie_TCSCF = grad_non_hermit - energy_TCSCF = TC_HF_energy - energy_TCSCF_1e = TC_HF_one_e_energy - energy_TCSCF_2e = TC_HF_two_e_energy - energy_TCSCF_3e = 0.d0 - if(three_body_h_tc) then - energy_TCSCF_3e = diag_three_elem_hf - endif - delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous - delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous - - energy_TCSCF_previous = energy_TCSCF - gradie_TCSCF_previous = gradie_TCSCF - - - level_shift_save = level_shift_TCSCF - 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) - - - print *, ' iteration = ', iteration_TCSCF - print *, ' total TC energy = ', energy_TCSCF - print *, ' 1-e TC energy = ', energy_TCSCF_1e - print *, ' 2-e TC energy = ', energy_TCSCF_2e - print *, ' 3-e TC energy = ', energy_TCSCF_3e - print *, ' |delta TC energy| = ', dabs(delta_energy_TCSCF) - print *, ' TC gradient = ', gradie_TCSCF - print *, ' delta TC gradient = ', delta_gradie_TCSCF - print *, ' max TC DIIS error = ', max_error_DIIS_TCSCF - print *, ' TC DIIS dim = ', dim_DIIS - print *, ' TC level shift = ', level_shift_TCSCF - print *, ' ' - - call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) - call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) - - if(qp_stop()) exit - enddo - - ! --- - - print *, ' TCSCF DIIS converged !' - call print_energy_and_mos() - - call write_time(6) - - deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS) - -end - -! --- - -subroutine extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, F_ao, size_F_ao, iteration_TCSCF, dim_DIIS) - - BEGIN_DOC - ! - ! Compute the extrapolated Fock matrix using the DIIS procedure - ! - ! e = \sum_i c_i e_i and \sum_i c_i = 1 - ! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1) - ! - END_DOC - - implicit none - - integer, intent(in) :: iteration_TCSCF, size_F_ao - integer, intent(inout) :: dim_DIIS - double precision, intent(in) :: F_DIIS(ao_num,ao_num,dim_DIIS) - double precision, intent(in) :: e_DIIS(ao_num,ao_num,dim_DIIS) - double precision, intent(inout) :: F_ao(size_F_ao,ao_num) - - double precision, allocatable :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:) - - integer :: i, j, k, l, i_DIIS, j_DIIS - integer :: lwork - double precision :: rcond, ferr, berr - integer, allocatable :: iwork(:) - double precision, allocatable :: scratch(:,:) - - if(dim_DIIS < 1) then - return - endif - - allocate( B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), X_vector_DIIS(dim_DIIS+1) & - , C_vector_DIIS(dim_DIIS+1), scratch(ao_num,ao_num) ) - - ! Compute the matrices B and X - B_matrix_DIIS(:,:) = 0.d0 - do j = 1, dim_DIIS - j_DIIS = min(dim_DIIS, mod(iteration_TCSCF-j, max_dim_DIIS_TCSCF)+1) - - do i = 1, dim_DIIS - i_DIIS = min(dim_DIIS, mod(iteration_TCSCF-i, max_dim_DIIS_TCSCF)+1) - - ! Compute product of two errors vectors - do l = 1, ao_num - do k = 1, ao_num - B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + e_DIIS(k,l,i_DIIS) * e_DIIS(k,l,j_DIIS) - enddo - enddo - - enddo - enddo - - ! Pad B matrix and build the X matrix - - C_vector_DIIS(:) = 0.d0 - do i = 1, dim_DIIS - B_matrix_DIIS(i,dim_DIIS+1) = -1.d0 - B_matrix_DIIS(dim_DIIS+1,i) = -1.d0 - enddo - C_vector_DIIS(dim_DIIS+1) = -1.d0 - - deallocate(scratch) - - ! Estimate condition number of B - integer :: info - double precision :: anorm - integer, allocatable :: ipiv(:) - double precision, allocatable :: AF(:,:) - double precision, external :: dlange - - lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5) - allocate(AF(dim_DIIS+1,dim_DIIS+1)) - allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) ) - allocate(scratch(lwork,1)) - scratch(:,1) = 0.d0 - - anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1)) - - AF(:,:) = B_matrix_DIIS(:,:) - call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info) - if(info /= 0) then - dim_DIIS = 0 - return - endif - - call dgecon('1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info) - if(info /= 0) then - dim_DIIS = 0 - return - endif - - if(rcond < 1.d-14) then - dim_DIIS = 0 - return - endif - - ! solve the linear system C = B x X - - X_vector_DIIS = C_vector_DIIS - call dgesv(dim_DIIS+1, 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv , X_vector_DIIS, size(X_vector_DIIS, 1), info) - - deallocate(scratch, AF, iwork) - if(info < 0) then - stop ' bug in TC-DIIS' - endif - - ! Compute extrapolated Fock matrix - - !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) - do j = 1, ao_num - do i = 1, ao_num - F_ao(i,j) = 0.d0 - enddo - do k = 1, dim_DIIS - if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle - do i = 1,ao_num - ! FPE here - F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1) - enddo - enddo - enddo - !$OMP END PARALLEL DO - -end - -! --- - diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index e7f00ce2..41e7cad6 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -459,3 +459,38 @@ subroutine v2_over_x(v,x,res) res = 0.5d0 * (tmp - delta_E) end + +subroutine sum_A_At(A, N) + + !BEGIN_DOC + ! useful for symmetrizing a tensor without a temporary tensor + !END_DOC + + implicit none + integer, intent(in) :: N + double precision, intent(inout) :: A(N,N) + integer :: i, j + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j) & + !$OMP SHARED (A, N) + !$OMP DO + do j = 1, N + do i = j, N + A(i,j) += A(j,i) + enddo + enddo + !$OMP END DO + + !$OMP DO + do j = 2, N + do i = 1, j-1 + A(i,j) = A(j,i) + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + +end + From e436e22f2a51fa56fed0cfebcaf46b1c6e7cac5e Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Sat, 4 Mar 2023 14:28:29 +0100 Subject: [PATCH 6/6] remove Gauss_Prod when expo = 0 --- src/ao_many_one_e_ints/grad2_jmu_manu.irp.f | 143 ++++++++++++------ .../grad_lapl_jmu_manu.irp.f | 79 ++++++---- src/ao_many_one_e_ints/listj1b.irp.f | 20 +-- 3 files changed, 156 insertions(+), 86 deletions(-) diff --git a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f index f01ed5ba..8e253d75 100644 --- a/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad2_jmu_manu.irp.f @@ -18,6 +18,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n double precision :: int_gauss, dsqpi_3_2, int_j1b double precision :: factor_ij_1s, beta_ij, center_ij_1s(3), sq_pi_3_2 double precision, allocatable :: int_fit_v(:) + double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao_with1s print*, ' providing int2_grad1u2_grad2u2_j1b2_test ...' @@ -39,48 +40,61 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, & !$OMP ao_overlap_abs,sq_pi_3_2) !$OMP DO SCHEDULE(dynamic) - do ipoint = 1, n_points_final_grid - r(1) = final_grid_points(1,ipoint) - r(2) = final_grid_points(2,ipoint) - r(3) = final_grid_points(3,ipoint) - do i = 1, ao_num - do j = i, ao_num - if(ao_overlap_abs(j,i) .lt. 1.d-12) then - cycle - endif - - do i_1s = 1, List_comb_thr_b3_size(j,i) + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + do i = 1, ao_num + do j = i, ao_num + if(ao_overlap_abs(j,i) .lt. 1.d-12) then + cycle + endif - coef = List_comb_thr_b3_coef (i_1s,j,i) - beta = List_comb_thr_b3_expo (i_1s,j,i) - int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) - B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) - B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) - B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b3_j1b(1,j,i) + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_1_erf_x_2(i_fit) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) + if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle + int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j) + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- - do i_fit = 1, ng_fit_jast + do i_1s = 2, List_comb_thr_b3_size(j,i) + + coef = List_comb_thr_b3_coef (i_1s,j,i) + beta = List_comb_thr_b3_expo (i_1s,j,i) + int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) + B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) + B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) + B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) - expo_fit = expo_gauss_1_erf_x_2(i_fit) - !DIR$ FORCEINLINE - call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) - coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_1_erf_x_2(i_fit) + !DIR$ FORCEINLINE + call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) + coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef ! if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version - if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle - + if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle ! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & ! expo_fit, i, j, int_fit_v, n_points_final_grid) - int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss - - enddo + int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss + enddo enddo - enddo - enddo - enddo - !$OMP END DO - !$OMP END PARALLEL + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL do ipoint = 1, n_points_final_grid do i = 1, ao_num @@ -239,9 +253,27 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ do i = 1, ao_num do j = i, ao_num - tmp = 0.d0 - do i_1s = 1, List_comb_thr_b3_size(j,i) + + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b3_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_j_mu_x_2(i_fit) + coef_fit = coef_gauss_j_mu_x_2(i_fit) + if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle + int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) + tmp += coef_fit * int_fit + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b3_size(j,i) coef = List_comb_thr_b3_coef (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i) @@ -252,23 +284,15 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_ B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) do i_fit = 1, ng_fit_jast - expo_fit = expo_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit) !DIR$ FORCEINLINE call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) ! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle - - ! --- - - int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - - tmp += coef * coef_fit * int_fit + int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) + tmp += coef * coef_fit * int_fit enddo - - ! --- - enddo int2_u2_j1b2_test(j,i,ipoint) = tmp @@ -451,13 +475,34 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p do ipoint = 1, n_points_final_grid do i = 1, ao_num do j = i, ao_num - if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle + + if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10) cycle + r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) r(3) = final_grid_points(3,ipoint) tmp = 0.d0 - do i_1s = 1, List_comb_thr_b3_size(j,i) + + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b3_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_j_mu_1_erf(i_fit) + if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle + coef_fit = coef_gauss_j_mu_1_erf(i_fit) + int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r) + tmp += coef_fit * int_fit + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b3_size(j,i) coef = List_comb_thr_b3_coef (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i) @@ -469,9 +514,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & + (B_center(2) - r(2)) * (B_center(2) - r(2)) & + (B_center(3) - r(3)) * (B_center(3) - r(3)) - do i_fit = 1, ng_fit_jast - expo_fit = expo_gauss_j_mu_1_erf(i_fit) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.1.d-15)cycle diff --git a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f index a6a55810..5c9f81e9 100644 --- a/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f +++ b/src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f @@ -192,9 +192,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po double precision :: coef, beta, B_center(3) double precision :: tmp double precision :: wall0, wall1 + double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot + double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b + double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao_with1s - double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b print*, ' providing v_ij_u_cst_mu_j1b_test ...' @@ -208,15 +210,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po !$OMP PARALLEL DEFAULT (NONE) & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & - !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & - !$OMP SHARED (n_points_final_grid, ao_num, & - !$OMP final_grid_points, ng_fit_jast, & + !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP final_grid_points, ng_fit_jast, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & - !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & - !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & + !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP DO - !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -227,8 +228,26 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle tmp = 0.d0 - do i_1s = 1, List_comb_thr_b2_size(j,i) + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b2_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + do i_fit = 1, ng_fit_jast + expo_fit = expo_gauss_j_mu_x(i_fit) + coef_fit = coef_gauss_j_mu_x(i_fit) + if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle + int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) + tmp += coef_fit * int_fit + enddo + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b2_size(j,i) coef = List_comb_thr_b2_coef (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) @@ -236,18 +255,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) - do i_fit = 1, ng_fit_jast - expo_fit = expo_gauss_j_mu_x(i_fit) coef_fit = coef_gauss_j_mu_x(i_fit) coeftot = coef * coef_fit if(dabs(coeftot).lt.1.d-15)cycle - double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp += coef * coef_fit * int_fit enddo enddo @@ -288,9 +303,12 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, double precision :: coef, beta, B_center(3) double precision :: tmp double precision :: wall0, wall1 + double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot + double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b + double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao_with1s - double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b + dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) provide mu_erf final_grid_points j1b_pen @@ -299,17 +317,16 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, v_ij_u_cst_mu_j1b_ng_1_test = 0.d0 !$OMP PARALLEL DEFAULT (NONE) & - !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & + !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & - !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & - !$OMP SHARED (n_points_final_grid, ao_num, & - !$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, & - !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & - !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & - !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & + !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & + !$OMP SHARED (n_points_final_grid, ao_num, & + !$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, & + !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & + !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & + !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2) !$OMP DO - !do ipoint = 1, 10 do ipoint = 1, n_points_final_grid r(1) = final_grid_points(1,ipoint) r(2) = final_grid_points(2,ipoint) @@ -320,8 +337,22 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle tmp = 0.d0 - do i_1s = 1, List_comb_thr_b2_size(j,i) + ! --- --- --- + ! i_1s = 1 + ! --- --- --- + + int_j1b = ao_abs_comb_b2_j1b(1,j,i) + if(dabs(int_j1b).lt.1.d-10) cycle + expo_fit = expo_good_j_mu_1gauss + int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) + tmp += int_fit + + ! --- --- --- + ! i_1s > 1 + ! --- --- --- + + do i_1s = 2, List_comb_thr_b2_size(j,i) coef = List_comb_thr_b2_coef (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i) int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) @@ -329,18 +360,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) - ! do i_fit = 1, ng_fit_jast - expo_fit = expo_good_j_mu_1gauss coef_fit = 1.d0 coeftot = coef * coef_fit if(dabs(coeftot).lt.1.d-15)cycle - double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u) if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) - tmp += coef * coef_fit * int_fit ! enddo enddo diff --git a/src/ao_many_one_e_ints/listj1b.irp.f b/src/ao_many_one_e_ints/listj1b.irp.f index e27bf723..4698cb27 100644 --- a/src/ao_many_one_e_ints/listj1b.irp.f +++ b/src/ao_many_one_e_ints/listj1b.irp.f @@ -102,11 +102,11 @@ END_PROVIDER List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i)) enddo - print *, ' coeff, expo & cent of list b2' - do i = 1, List_all_comb_b2_size - print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i) - print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i) - enddo + !print *, ' coeff, expo & cent of list b2' + !do i = 1, List_all_comb_b2_size + ! print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i) + ! print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i) + !enddo END_PROVIDER @@ -225,11 +225,11 @@ END_PROVIDER List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i)) enddo - print *, ' coeff, expo & cent of list b3' - do i = 1, List_all_comb_b3_size - print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i) - print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i) - enddo + !print *, ' coeff, expo & cent of list b3' + !do i = 1, List_all_comb_b3_size + ! print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i) + ! print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i) + !enddo END_PROVIDER