From 92ad3766ebef835e2e1f9888215699b3d4007bcf Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 27 Apr 2020 11:31:24 +0200 Subject: [PATCH 1/6] added two_body_dens_rout.irp.f --- src/cas_based_on_top/two_body_dens_rout.irp.f | 189 ++++++++++++++++++ src/two_body_rdm/act_2_rdm.irp.f | 2 + 2 files changed, 191 insertions(+) create mode 100644 src/cas_based_on_top/two_body_dens_rout.irp.f diff --git a/src/cas_based_on_top/two_body_dens_rout.irp.f b/src/cas_based_on_top/two_body_dens_rout.irp.f new file mode 100644 index 00000000..be17fcaf --- /dev/null +++ b/src/cas_based_on_top/two_body_dens_rout.irp.f @@ -0,0 +1,189 @@ + +subroutine give_n2_ii_val_ab(r1,r2,two_bod_dens) + implicit none + BEGIN_DOC +! contribution from purely inactive orbitals to n2_{\Psi^B}(r_1,r_2) for a CAS wave function + END_DOC + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: two_bod_dens + integer :: i,j,m,n,i_m,i_n + integer :: i_i,i_j + double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) + double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + ! You extract the inactive orbitals + allocate(mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) + do i_m = 1, n_inact_orb + mos_array_inact_r1(i_m) = mos_array_r1(list_inact(i_m)) + enddo + do i_m = 1, n_inact_orb + mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + allocate(mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) + do i_m = 1, n_basis_orb + mos_array_basis_r1(i_m) = mos_array_r1(list_basis(i_m)) + mos_array_basis_r2(i_m) = mos_array_r2(list_basis(i_m)) + enddo + + two_bod_dens = 0.d0 + ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space + do m = 1, n_inact_orb ! electron 1 + ! You browse all OCCUPIED BETA electrons in the \mathcal{A} space + do n = 1, n_inact_orb ! electron 2 + ! two_bod_dens(r_1,r_2) = n_alpha(r_1) * n_beta(r_2) + two_bod_dens += mos_array_inact_r1(m) * mos_array_inact_r1(m) * mos_array_inact_r2(n) * mos_array_inact_r2(n) + enddo + enddo +end + + +subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate) + BEGIN_DOC +! contribution from inactive and active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function + END_DOC + implicit none + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: two_bod_dens + integer :: i,orb_i,a,orb_a,n,m,b + double precision :: rho + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) + double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) + double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) + + two_bod_dens = 0.d0 + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + + ! You extract the inactive orbitals + allocate( mos_array_inact_r1(n_inact_orb) , mos_array_inact_r2(n_inact_orb) ) + do i = 1, n_inact_orb + mos_array_inact_r1(i) = mos_array_r1(list_inact(i)) + enddo + do i= 1, n_inact_orb + mos_array_inact_r2(i) = mos_array_r2(list_inact(i)) + enddo + + ! You extract the active orbitals + allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) + do i= 1, n_act_orb + mos_array_act_r1(i) = mos_array_r1(list_act(i)) + enddo + do i= 1, n_act_orb + mos_array_act_r2(i) = mos_array_r2(list_act(i)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) + do i= 1, n_basis_orb + mos_array_basis_r1(i) = mos_array_r1(list_basis(i)) + enddo + do i= 1, n_basis_orb + mos_array_basis_r2(i) = mos_array_r2(list_basis(i)) + enddo + + ! Contracted density : intermediate quantity + ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) + allocate(rho_tilde(n_inact_orb,n_act_orb)) + two_bod_dens = 0.d0 + do a = 1, n_act_orb + do i = 1, n_inact_orb + rho_tilde(i,a) = 0.d0 + do b = 1, n_act_orb + rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) + two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho + rho_tilde(i,a) += rho * mos_array_inact_r1(i) * mos_array_act_r2(b) + enddo + enddo + enddo +end + + +subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate) + BEGIN_DOC +! contribution from purely active orbitals to n2_{\Psi^B}(r_1,r_2) for the "istate" state of a CAS wave function + END_DOC + implicit none + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out):: two_bod_dens + integer :: i,orb_i,a,orb_a,n,m,b,c,d + double precision :: rho + double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) + double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) + double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) + + two_bod_dens = 0.d0 + ! You get all orbitals in r_1 and r_2 + allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) + call give_all_mos_at_r(r1,mos_array_r1) + call give_all_mos_at_r(r2,mos_array_r2) + + ! You extract the active orbitals + allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) + do i= 1, n_act_orb + mos_array_act_r1(i) = mos_array_r1(list_act(i)) + enddo + do i= 1, n_act_orb + mos_array_act_r2(i) = mos_array_r2(list_act(i)) + enddo + + ! You extract the orbitals belonging to the space \mathcal{B} + allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) + do i= 1, n_basis_orb + mos_array_basis_r1(i) = mos_array_r1(list_basis(i)) + enddo + do i= 1, n_basis_orb + mos_array_basis_r2(i) = mos_array_r2(list_basis(i)) + enddo + + ! Contracted density : intermediate quantity + ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) + allocate(rho_tilde(n_act_orb,n_act_orb)) + two_bod_dens = 0.d0 + rho_tilde = 0.d0 + do a = 1, n_act_orb ! 1 + do b = 1, n_act_orb ! 2 + do c = 1, n_act_orb ! 1 + do d = 1, n_act_orb ! 2 + rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) + rho_tilde(b,a) += rho + two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) + enddo + enddo + enddo + enddo + +end + +subroutine give_n2_cas(r1,r2,istate,n2_psi) + implicit none + BEGIN_DOC +! returns mu(r), f_psi, n2_psi for a general cas wave function + END_DOC + integer, intent(in) :: istate + double precision, intent(in) :: r1(3),r2(3) + double precision, intent(out) :: n2_psi + double precision :: two_bod_dens_ii + double precision :: two_bod_dens_ia + double precision :: two_bod_dens_aa + ! inactive-inactive part of n2_psi(r1,r2) + call give_n2_ii_val_ab(r,r,two_bod_dens_ii) + ! inactive-active part of n2_psi(r1,r2) + call give_n2_ia_val_ab(r,r,two_bod_dens_ia,istate) + ! active-active part of n2_psi(r1,r2) + call give_n2_aa_val_ab(r,r,two_bod_dens_aa,istate) + + n2_psi = n2_ii_val_ab + n2_ia_val_ab + n2_aa_val_ab + n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa + +end diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 2473d3c8..e3265572 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -2,6 +2,8 @@ BEGIN_PROVIDER [double precision, act_2_rdm_ab_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)] implicit none BEGIN_DOC +! 12 12 +! 1 2 1 2 == ! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons ! ! From 2047abcdb0dbf36ce7876cfae227df235511f6f3 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 29 Apr 2020 14:48:28 +0200 Subject: [PATCH 2/6] modifs --- src/cas_based_on_top/two_body_dens_rout.irp.f | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/cas_based_on_top/two_body_dens_rout.irp.f b/src/cas_based_on_top/two_body_dens_rout.irp.f index be17fcaf..19d7632f 100644 --- a/src/cas_based_on_top/two_body_dens_rout.irp.f +++ b/src/cas_based_on_top/two_body_dens_rout.irp.f @@ -92,16 +92,12 @@ subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate) enddo ! Contracted density : intermediate quantity - ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) - allocate(rho_tilde(n_inact_orb,n_act_orb)) two_bod_dens = 0.d0 do a = 1, n_act_orb do i = 1, n_inact_orb - rho_tilde(i,a) = 0.d0 do b = 1, n_act_orb rho = one_e_act_dm_beta_mo_for_dft(b,a,istate) + one_e_act_dm_alpha_mo_for_dft(b,a,istate) two_bod_dens += mos_array_inact_r1(i) * mos_array_inact_r1(i) * mos_array_act_r2(a) * mos_array_act_r2(b) * rho - rho_tilde(i,a) += rho * mos_array_inact_r1(i) * mos_array_act_r2(b) enddo enddo enddo @@ -147,16 +143,12 @@ subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate) enddo ! Contracted density : intermediate quantity - ! rho_tilde(i,a) = \sum_b rho(b,a) * phi_i(1) * phi_j(2) - allocate(rho_tilde(n_act_orb,n_act_orb)) two_bod_dens = 0.d0 - rho_tilde = 0.d0 do a = 1, n_act_orb ! 1 do b = 1, n_act_orb ! 2 do c = 1, n_act_orb ! 1 do d = 1, n_act_orb ! 2 rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate) - rho_tilde(b,a) += rho two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) enddo enddo @@ -177,13 +169,12 @@ subroutine give_n2_cas(r1,r2,istate,n2_psi) double precision :: two_bod_dens_ia double precision :: two_bod_dens_aa ! inactive-inactive part of n2_psi(r1,r2) - call give_n2_ii_val_ab(r,r,two_bod_dens_ii) + call give_n2_ii_val_ab(r1,r2,two_bod_dens_ii) ! inactive-active part of n2_psi(r1,r2) - call give_n2_ia_val_ab(r,r,two_bod_dens_ia,istate) + call give_n2_ia_val_ab(r1,r2,two_bod_dens_ia,istate) ! active-active part of n2_psi(r1,r2) - call give_n2_aa_val_ab(r,r,two_bod_dens_aa,istate) + call give_n2_aa_val_ab(r1,r2,two_bod_dens_aa,istate) - n2_psi = n2_ii_val_ab + n2_ia_val_ab + n2_aa_val_ab n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa end From 9737de21b782309ae497dab0a7ffac854ae05819 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Wed, 29 Apr 2020 15:11:48 +0200 Subject: [PATCH 3/6] removed spurious dependency --- src/cas_based_on_top/two_body_dens_rout.irp.f | 32 ++----------------- 1 file changed, 2 insertions(+), 30 deletions(-) diff --git a/src/cas_based_on_top/two_body_dens_rout.irp.f b/src/cas_based_on_top/two_body_dens_rout.irp.f index 19d7632f..4a57a868 100644 --- a/src/cas_based_on_top/two_body_dens_rout.irp.f +++ b/src/cas_based_on_top/two_body_dens_rout.irp.f @@ -9,7 +9,6 @@ subroutine give_n2_ii_val_ab(r1,r2,two_bod_dens) integer :: i,j,m,n,i_m,i_n integer :: i_i,i_j double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) - double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) ! You get all orbitals in r_1 and r_2 allocate(mos_array_r1(mo_num) , mos_array_r2(mo_num) ) @@ -24,13 +23,6 @@ subroutine give_n2_ii_val_ab(r1,r2,two_bod_dens) mos_array_inact_r2(i_m) = mos_array_r2(list_inact(i_m)) enddo - ! You extract the orbitals belonging to the space \mathcal{B} - allocate(mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) - do i_m = 1, n_basis_orb - mos_array_basis_r1(i_m) = mos_array_r1(list_basis(i_m)) - mos_array_basis_r2(i_m) = mos_array_r2(list_basis(i_m)) - enddo - two_bod_dens = 0.d0 ! You browse all OCCUPIED ALPHA electrons in the \mathcal{A} space do m = 1, n_inact_orb ! electron 1 @@ -55,7 +47,6 @@ subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate) double precision :: rho double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) double precision, allocatable :: mos_array_inact_r1(:),mos_array_inact_r2(:) - double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) two_bod_dens = 0.d0 @@ -74,7 +65,7 @@ subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate) enddo ! You extract the active orbitals - allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) + allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) ) do i= 1, n_act_orb mos_array_act_r1(i) = mos_array_r1(list_act(i)) enddo @@ -82,15 +73,6 @@ subroutine give_n2_ia_val_ab(r1,r2,two_bod_dens,istate) mos_array_act_r2(i) = mos_array_r2(list_act(i)) enddo - ! You extract the orbitals belonging to the space \mathcal{B} - allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) - do i= 1, n_basis_orb - mos_array_basis_r1(i) = mos_array_r1(list_basis(i)) - enddo - do i= 1, n_basis_orb - mos_array_basis_r2(i) = mos_array_r2(list_basis(i)) - enddo - ! Contracted density : intermediate quantity two_bod_dens = 0.d0 do a = 1, n_act_orb @@ -115,7 +97,6 @@ subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate) integer :: i,orb_i,a,orb_a,n,m,b,c,d double precision :: rho double precision, allocatable :: mos_array_r1(:) , mos_array_r2(:) - double precision, allocatable :: mos_array_basis_r1(:),mos_array_basis_r2(:) double precision, allocatable :: mos_array_act_r1(:),mos_array_act_r2(:) two_bod_dens = 0.d0 @@ -125,7 +106,7 @@ subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate) call give_all_mos_at_r(r2,mos_array_r2) ! You extract the active orbitals - allocate( mos_array_act_r1(n_basis_orb) , mos_array_act_r2(n_basis_orb) ) + allocate( mos_array_act_r1(n_act_orb) , mos_array_act_r2(n_act_orb) ) do i= 1, n_act_orb mos_array_act_r1(i) = mos_array_r1(list_act(i)) enddo @@ -133,15 +114,6 @@ subroutine give_n2_aa_val_ab(r1,r2,two_bod_dens,istate) mos_array_act_r2(i) = mos_array_r2(list_act(i)) enddo - ! You extract the orbitals belonging to the space \mathcal{B} - allocate( mos_array_basis_r1(n_basis_orb) , mos_array_basis_r2(n_basis_orb) ) - do i= 1, n_basis_orb - mos_array_basis_r1(i) = mos_array_r1(list_basis(i)) - enddo - do i= 1, n_basis_orb - mos_array_basis_r2(i) = mos_array_r2(list_basis(i)) - enddo - ! Contracted density : intermediate quantity two_bod_dens = 0.d0 do a = 1, n_act_orb ! 1 From b0823fe00341810439df9d1af3b18668d519ac5b Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 30 Apr 2020 19:14:50 +0200 Subject: [PATCH 4/6] added shank --- src/utils/shank.irp.f | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 src/utils/shank.irp.f diff --git a/src/utils/shank.irp.f b/src/utils/shank.irp.f new file mode 100644 index 00000000..4ec10bbe --- /dev/null +++ b/src/utils/shank.irp.f @@ -0,0 +1,39 @@ +double precision function shank3_f(array,n,nmax) + implicit none + integer, intent(in) :: n,nmax + double precision, intent(in) :: array(0:nmax) ! array of the partial sums + integer :: ntmp + double precision :: shank1(0:nmax),shank2(0:nmax),shank3(0:nmax) + ntmp = n + call shank(array,ntmp,nmax,shank1) + ntmp = ntmp - 2 + call shank(shank1,ntmp,nmax,shank2) + ntmp = ntmp - 2 + call shank(shank2,ntmp,nmax,shank3) + ntmp = ntmp - 2 + shank3_f = shank3(ntmp) +end + + +subroutine shank(array,n,nmax,shank1) + implicit none + integer, intent(in) :: n,nmax + double precision, intent(in) :: array(0:nmax) + double precision, intent(out) :: shank1(0:nmax) + integer :: i,j + double precision :: shank_function + do i = 1, n-1 + shank1(i-1) = shank_function(array,i,nmax) + enddo +end + +double precision function shank_function(array,i,n) + implicit none + integer, intent(in) :: i,n + double precision, intent(in) :: array(0:n) + double precision :: b_n, b_n1 + b_n = array(i) - array(i-1) + b_n1 = array(i+1) - array(i) + shank_function = array(i+1) - b_n1*b_n1/(b_n1-b_n) +end + From 61df4e01dfebc7614284351a338ac70310c92a22 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 30 Apr 2020 19:35:21 +0200 Subject: [PATCH 5/6] added shank --- src/utils/shank.irp.f | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/utils/shank.irp.f b/src/utils/shank.irp.f index 4ec10bbe..238538f3 100644 --- a/src/utils/shank.irp.f +++ b/src/utils/shank.irp.f @@ -34,6 +34,12 @@ double precision function shank_function(array,i,n) double precision :: b_n, b_n1 b_n = array(i) - array(i-1) b_n1 = array(i+1) - array(i) - shank_function = array(i+1) - b_n1*b_n1/(b_n1-b_n) + if(dabs(b_n1-b_n).lt.1.d-12)then + shank_function = array(i+1) + else + shank_function = array(i+1) - b_n1*b_n1/(b_n1-b_n) + endif + end + From e864eb1cf383d69c80f9edacaa93211ccfb7b8e7 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 11 May 2020 16:04:16 +0200 Subject: [PATCH 6/6] added the possibility to use no_vvvv integrals from EZFIO --- src/becke_numerical_grid/EZFIO.cfg | 19 +++++++++++++ src/becke_numerical_grid/grid_becke.irp.f | 7 ++++- src/becke_numerical_grid/list_angular_grid | 32 ++++++++++++++++++++++ src/casscf/casscf.irp.f | 1 + src/mo_two_e_ints/EZFIO.cfg | 6 ++++ src/mo_two_e_ints/four_idx_novvvv.irp.f | 16 ++++++----- src/utils/integration.irp.f | 1 - src/utils/shank.irp.f | 29 +++++++++++++------- 8 files changed, 92 insertions(+), 19 deletions(-) create mode 100644 src/becke_numerical_grid/list_angular_grid diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg index ca2100a1..17e8fd90 100644 --- a/src/becke_numerical_grid/EZFIO.cfg +++ b/src/becke_numerical_grid/EZFIO.cfg @@ -14,3 +14,22 @@ type: double precision doc: threshold on the weight of a given grid point interface: ezfio,provider,ocaml default: 1.e-20 + +[my_grid_becke] +type: logical +doc: if True, the number of angular and radial grid points are read from EZFIO +interface: ezfio,provider,ocaml +default: False + +[my_n_pt_r_grid] +type: integer +doc: Number of radial grid points given from input +interface: ezfio,provider,ocaml +default: 300 + +[my_n_pt_a_grid] +type: integer +doc: Number of angular grid points given from input. Warning, this number cannot be any integer. See file list_angular_grid +interface: ezfio,provider,ocaml +default: 1202 + diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index e72f6460..bedf0a84 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/src/becke_numerical_grid/grid_becke.irp.f @@ -8,7 +8,8 @@ ! ! These numbers are automatically set by setting the grid_type_sgn parameter END_DOC -select case (grid_type_sgn) +if(.not.my_grid_becke)then + select case (grid_type_sgn) case(0) n_points_radial_grid = 23 n_points_integration_angular = 170 @@ -25,6 +26,10 @@ select case (grid_type_sgn) write(*,*) '!!! Quadrature grid not available !!!' stop end select +else + n_points_radial_grid = my_n_pt_r_grid + n_points_integration_angular = my_n_pt_a_grid +endif END_PROVIDER BEGIN_PROVIDER [integer, n_points_grid_per_atom] diff --git a/src/becke_numerical_grid/list_angular_grid b/src/becke_numerical_grid/list_angular_grid new file mode 100644 index 00000000..252c7432 --- /dev/null +++ b/src/becke_numerical_grid/list_angular_grid @@ -0,0 +1,32 @@ +0006 +0014 +0026 +0038 +0050 +0074 +0086 +0110 +0146 +0170 +0194 +0230 +0266 +0302 +0350 +0434 +0590 +0770 +0974 +1202 +1454 +1730 +2030 +2354 +2702 +3074 +3470 +3890 +4334 +4802 +5294 +5810 diff --git a/src/casscf/casscf.irp.f b/src/casscf/casscf.irp.f index d83aa271..950cfd55 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -5,6 +5,7 @@ program casscf END_DOC call reorder_orbitals_for_casscf no_vvvv_integrals = .True. + touch no_vvvv_integrals pt2_max = 0.02 SOFT_TOUCH no_vvvv_integrals pt2_max call run_stochastic_cipsi diff --git a/src/mo_two_e_ints/EZFIO.cfg b/src/mo_two_e_ints/EZFIO.cfg index bec74552..ea47c51c 100644 --- a/src/mo_two_e_ints/EZFIO.cfg +++ b/src/mo_two_e_ints/EZFIO.cfg @@ -11,3 +11,9 @@ interface: ezfio,provider,ocaml default: 1.e-15 ezfio_name: threshold_mo +[no_vvvv_integrals] +type: logical +doc: If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices +interface: ezfio,provider,ocaml +default: false + diff --git a/src/mo_two_e_ints/four_idx_novvvv.irp.f b/src/mo_two_e_ints/four_idx_novvvv.irp.f index 054d0a35..a4cdfd52 100644 --- a/src/mo_two_e_ints/four_idx_novvvv.irp.f +++ b/src/mo_two_e_ints/four_idx_novvvv.irp.f @@ -1,11 +1,11 @@ -BEGIN_PROVIDER [ logical, no_vvvv_integrals ] - implicit none - BEGIN_DOC +!BEGIN_PROVIDER [ logical, no_vvvv_integrals ] +! implicit none +! BEGIN_DOC ! If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices - END_DOC - - no_vvvv_integrals = .False. -END_PROVIDER +! END_DOC +! +! no_vvvv_integrals = .False. +!END_PROVIDER BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ] implicit none @@ -56,6 +56,8 @@ subroutine four_idx_novvvv BEGIN_DOC ! Retransform MO integrals for next CAS-SCF step END_DOC + print*,'Using partial transformation' + print*,'It will not transform all integrals with at least 3 indices within the virtuals' integer :: i,j,k,l,n_integrals double precision, allocatable :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:) double precision, external :: get_ao_two_e_integral diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f index 3ff1bb42..c907e425 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -75,7 +75,6 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, P_new(0,1) = 0.d0 P_new(0,2) = 0.d0 P_new(0,3) = 0.d0 - !DIR$ FORCEINLINE call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) if (fact_k < thresh) then diff --git a/src/utils/shank.irp.f b/src/utils/shank.irp.f index 238538f3..d459399f 100644 --- a/src/utils/shank.irp.f +++ b/src/utils/shank.irp.f @@ -1,17 +1,26 @@ -double precision function shank3_f(array,n,nmax) +double precision function shank_general(array,n,nmax) implicit none integer, intent(in) :: n,nmax double precision, intent(in) :: array(0:nmax) ! array of the partial sums - integer :: ntmp - double precision :: shank1(0:nmax),shank2(0:nmax),shank3(0:nmax) + integer :: ntmp,i + double precision :: sum(0:nmax),shank1(0:nmax) + if(n.lt.3)then + print*,'You asked to Shank a sum but the order is smaller than 3 ...' + print*,'n = ',n + print*,'stopping ....' + stop + endif ntmp = n - call shank(array,ntmp,nmax,shank1) - ntmp = ntmp - 2 - call shank(shank1,ntmp,nmax,shank2) - ntmp = ntmp - 2 - call shank(shank2,ntmp,nmax,shank3) - ntmp = ntmp - 2 - shank3_f = shank3(ntmp) + sum = array + i = 0 + do while(ntmp.ge.2) + i += 1 +! print*,'i = ',i + call shank(sum,ntmp,nmax,shank1) + ntmp = ntmp - 2 + sum = shank1 + shank_general = shank1(ntmp) + enddo end