From b1673f66a844cfe6dc5f5b558368d596f037fc56 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 21 Apr 2020 01:02:48 +0200 Subject: [PATCH 01/11] istep in davidson_parallel --- src/davidson/davidson_parallel.irp.f | 5 +++++ src/zmq/utils.irp.f | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/davidson/davidson_parallel.irp.f b/src/davidson/davidson_parallel.irp.f index c0d94b35..aed81063 100644 --- a/src/davidson/davidson_parallel.irp.f +++ b/src/davidson/davidson_parallel.irp.f @@ -438,6 +438,11 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze) ipos=1 do imin=1,N_det,tasksize imax = min(N_det,imin-1+tasksize) + if (imin==1) then + istep = 2 + else + istep = 1 + endif do ishift=0,istep-1 write(task(ipos:ipos+50),'(4(I11,1X),1X,1A)') imin, imax, ishift, istep, '|' ipos = ipos+50 diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f index fc6f2ba6..67b386e5 100644 --- a/src/zmq/utils.irp.f +++ b/src/zmq/utils.irp.f @@ -585,7 +585,7 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in) stop 'Wrong end of job' endif - do i=3600,1,-1 + do i=360,1,-1 rc = f77_zmq_send(zmq_to_qp_run_socket, 'end_job '//trim(zmq_state),8+len(trim(zmq_state)),0) rc = f77_zmq_recv(zmq_to_qp_run_socket, message, 512, 0) if (trim(message(1:13)) == 'error waiting') then From 92ad3766ebef835e2e1f9888215699b3d4007bcf Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 27 Apr 2020 11:31:24 +0200 Subject: [PATCH 02/11] 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 03/11] 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 04/11] 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 05/11] 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 06/11] 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 07/11] 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 From a62456b238464145fd834fa8192050a7aaafb330 Mon Sep 17 00:00:00 2001 From: Michel Caffarel Date: Tue, 12 May 2020 21:57:05 +0200 Subject: [PATCH 08/11] Avoid to provide ao basis when reading integrals --- src/ao_basis/aos.irp.f | 13 +++++++++++++ src/ao_one_e_ints/ao_overlap.irp.f | 1 + src/ao_two_e_ints/two_e_integrals.irp.f | 7 +++++-- 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index e95ac711..f747b1de 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -3,6 +3,7 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ] BEGIN_DOC ! Max number of primitives. END_DOC +print *, 'XXXX', irp_here ao_prim_num_max = maxval(ao_prim_num) END_PROVIDER @@ -19,6 +20,7 @@ END_PROVIDER C_A(1) = 0.d0 C_A(2) = 0.d0 C_A(3) = 0.d0 +print *, 'XXXX', irp_here ao_coef_normalized = 0.d0 do i=1,ao_num @@ -65,6 +67,7 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalization_libint_factor, (ao_num) integer :: l, powA(3), nz integer :: i,j,k nz=100 +print *, 'XXXX', irp_here C_A(1) = 0.d0 C_A(2) = 0.d0 C_A(3) = 0.d0 @@ -99,6 +102,7 @@ END_PROVIDER integer :: iorder(ao_prim_num_max) double precision :: d(ao_prim_num_max,2) integer :: i,j +print *, 'XXXX', irp_here do i=1,ao_num do j=1,ao_prim_num(i) iorder(j) = j @@ -121,6 +125,7 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered_transp, (ao_prim_n ! Transposed :c:data:`ao_coef_normalized_ordered` END_DOC integer :: i,j +print *, 'XXXX', irp_here do j=1, ao_num do i=1, ao_prim_num_max ao_coef_normalized_ordered_transp(i,j) = ao_coef_normalized_ordered(j,i) @@ -135,6 +140,7 @@ BEGIN_PROVIDER [ double precision, ao_expo_ordered_transp, (ao_prim_num_max,ao_n ! Transposed :c:data:`ao_expo_ordered` END_DOC integer :: i,j +print *, 'XXXX', irp_here do j=1, ao_num do i=1, ao_prim_num_max ao_expo_ordered_transp(i,j) = ao_expo_ordered(j,i) @@ -151,6 +157,7 @@ END_PROVIDER ! :math:`l` value of the |AO|: :math`a+b+c` in :math:`x^a y^b z^c` END_DOC integer :: i +print *, 'XXXX', irp_here do i=1,ao_num ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) ao_l_char(i) = l_to_character(ao_l(i)) @@ -167,6 +174,7 @@ integer function ao_power_index(nx,ny,nz) ! :math:`\frac{1}{2} (l-n_x) (l-n_x+1) + n_z + 1` END_DOC integer :: l +print *, 'XXXX', irp_here l = nx + ny + nz ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1 end @@ -177,6 +185,7 @@ BEGIN_PROVIDER [ character*(128), l_to_character, (0:7)] ! Character corresponding to the "l" value of an |AO| END_DOC implicit none +print *, 'XXXX', irp_here l_to_character(0)='s' l_to_character(1)='p' l_to_character(2)='d' @@ -195,6 +204,7 @@ END_PROVIDER ! Number of |AOs| per atom END_DOC integer :: i +print *, 'XXXX', irp_here Nucl_N_Aos = 0 do i = 1, ao_num Nucl_N_Aos(ao_nucl(i)) +=1 @@ -209,6 +219,7 @@ END_PROVIDER END_DOC integer :: i integer, allocatable :: nucl_tmp(:) +print *, 'XXXX', irp_here allocate(nucl_tmp(nucl_num)) nucl_tmp = 0 Nucl_Aos = 0 @@ -229,6 +240,7 @@ END_PROVIDER ! By convention, for p,d,f and g |AOs|, we take the index ! of the |AO| with the the corresponding power in the x axis END_DOC +print *, 'XXXX', irp_here do i = 1, nucl_num Nucl_num_shell_Aos(i) = 0 @@ -276,6 +288,7 @@ BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] END_DOC integer :: i character*(4) :: give_ao_character_space +print *, 'XXXX', irp_here do i=1,ao_num if(ao_l(i)==0)then diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index d7300936..d7a4d640 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -109,6 +109,7 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ] double precision :: A_center(3), B_center(3) integer :: power_A(3), power_B(3) double precision :: lower_exp_val, dx +print *, "XXX---", irp_here if (is_periodic) then do j=1,ao_num do i= 1,ao_num diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index acc9ab7a..973a9011 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -343,8 +343,6 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] integer :: kk, m, j1, i1, lmax character*(64) :: fmt - integral = ao_two_e_integral(1,1,1,1) - double precision :: map_mb PROVIDE read_ao_two_e_integrals io_ao_two_e_integrals if (read_ao_two_e_integrals) then @@ -360,6 +358,11 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] call wall_time(wall_1) call cpu_time(cpu_1) + if (.True.) then + ! Avoid openMP + integral = ao_two_e_integral(1,1,1,1) + endif + integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') From 329bcf805bd03768243d41fe2c16f9822796f22d Mon Sep 17 00:00:00 2001 From: Michel Caffarel Date: Tue, 12 May 2020 22:46:39 +0200 Subject: [PATCH 09/11] Fixed schwartz screening when integrals are read --- src/ao_basis/aos.irp.f | 13 -- src/ao_one_e_ints/ao_overlap.irp.f | 1 - src/ao_one_e_ints/pot_ao_ints.irp.f | 10 +- src/ao_one_e_ints/screening.irp.f | 5 +- src/ao_two_e_ints/screening.irp.f | 6 +- src/ao_two_e_ints/two_e_integrals.irp.f | 253 ++++++++++++------------ src/mo_two_e_ints/mo_bi_integrals.irp.f | 2 - 7 files changed, 137 insertions(+), 153 deletions(-) diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index f747b1de..e95ac711 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -3,7 +3,6 @@ BEGIN_PROVIDER [ integer, ao_prim_num_max ] BEGIN_DOC ! Max number of primitives. END_DOC -print *, 'XXXX', irp_here ao_prim_num_max = maxval(ao_prim_num) END_PROVIDER @@ -20,7 +19,6 @@ END_PROVIDER C_A(1) = 0.d0 C_A(2) = 0.d0 C_A(3) = 0.d0 -print *, 'XXXX', irp_here ao_coef_normalized = 0.d0 do i=1,ao_num @@ -67,7 +65,6 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalization_libint_factor, (ao_num) integer :: l, powA(3), nz integer :: i,j,k nz=100 -print *, 'XXXX', irp_here C_A(1) = 0.d0 C_A(2) = 0.d0 C_A(3) = 0.d0 @@ -102,7 +99,6 @@ END_PROVIDER integer :: iorder(ao_prim_num_max) double precision :: d(ao_prim_num_max,2) integer :: i,j -print *, 'XXXX', irp_here do i=1,ao_num do j=1,ao_prim_num(i) iorder(j) = j @@ -125,7 +121,6 @@ BEGIN_PROVIDER [ double precision, ao_coef_normalized_ordered_transp, (ao_prim_n ! Transposed :c:data:`ao_coef_normalized_ordered` END_DOC integer :: i,j -print *, 'XXXX', irp_here do j=1, ao_num do i=1, ao_prim_num_max ao_coef_normalized_ordered_transp(i,j) = ao_coef_normalized_ordered(j,i) @@ -140,7 +135,6 @@ BEGIN_PROVIDER [ double precision, ao_expo_ordered_transp, (ao_prim_num_max,ao_n ! Transposed :c:data:`ao_expo_ordered` END_DOC integer :: i,j -print *, 'XXXX', irp_here do j=1, ao_num do i=1, ao_prim_num_max ao_expo_ordered_transp(i,j) = ao_expo_ordered(j,i) @@ -157,7 +151,6 @@ END_PROVIDER ! :math:`l` value of the |AO|: :math`a+b+c` in :math:`x^a y^b z^c` END_DOC integer :: i -print *, 'XXXX', irp_here do i=1,ao_num ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3) ao_l_char(i) = l_to_character(ao_l(i)) @@ -174,7 +167,6 @@ integer function ao_power_index(nx,ny,nz) ! :math:`\frac{1}{2} (l-n_x) (l-n_x+1) + n_z + 1` END_DOC integer :: l -print *, 'XXXX', irp_here l = nx + ny + nz ao_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1 end @@ -185,7 +177,6 @@ BEGIN_PROVIDER [ character*(128), l_to_character, (0:7)] ! Character corresponding to the "l" value of an |AO| END_DOC implicit none -print *, 'XXXX', irp_here l_to_character(0)='s' l_to_character(1)='p' l_to_character(2)='d' @@ -204,7 +195,6 @@ END_PROVIDER ! Number of |AOs| per atom END_DOC integer :: i -print *, 'XXXX', irp_here Nucl_N_Aos = 0 do i = 1, ao_num Nucl_N_Aos(ao_nucl(i)) +=1 @@ -219,7 +209,6 @@ END_PROVIDER END_DOC integer :: i integer, allocatable :: nucl_tmp(:) -print *, 'XXXX', irp_here allocate(nucl_tmp(nucl_num)) nucl_tmp = 0 Nucl_Aos = 0 @@ -240,7 +229,6 @@ END_PROVIDER ! By convention, for p,d,f and g |AOs|, we take the index ! of the |AO| with the the corresponding power in the x axis END_DOC -print *, 'XXXX', irp_here do i = 1, nucl_num Nucl_num_shell_Aos(i) = 0 @@ -288,7 +276,6 @@ BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ] END_DOC integer :: i character*(4) :: give_ao_character_space -print *, 'XXXX', irp_here do i=1,ao_num if(ao_l(i)==0)then diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f index d7a4d640..d7300936 100644 --- a/src/ao_one_e_ints/ao_overlap.irp.f +++ b/src/ao_one_e_ints/ao_overlap.irp.f @@ -109,7 +109,6 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ] double precision :: A_center(3), B_center(3) integer :: power_A(3), power_B(3) double precision :: lower_exp_val, dx -print *, "XXX---", irp_here if (is_periodic) then do j=1,ao_num do i= 1,ao_num diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 486ff534..4108ce71 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -3,6 +3,8 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] ! Nucleus-electron interaction, in the |AO| basis set. ! ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + ! + ! These integrals also contain the pseudopotential integrals. END_DOC implicit none double precision :: alpha, beta, gama, delta @@ -75,11 +77,11 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] !$OMP END DO !$OMP END PARALLEL - endif + IF (DO_PSEUDO) THEN + ao_integrals_n_e += ao_pseudo_integrals + ENDIF - IF (DO_PSEUDO) THEN - ao_integrals_n_e += ao_pseudo_integrals - ENDIF + endif if (write_ao_integrals_n_e) then diff --git a/src/ao_one_e_ints/screening.irp.f b/src/ao_one_e_ints/screening.irp.f index aa23a9b8..1bbe3c73 100644 --- a/src/ao_one_e_ints/screening.irp.f +++ b/src/ao_one_e_ints/screening.irp.f @@ -3,14 +3,11 @@ logical function ao_one_e_integral_zero(i,k) integer, intent(in) :: i,k ao_one_e_integral_zero = .False. - if (.not.(read_ao_one_e_integrals.or.is_periodic)) then + if (.not.((io_ao_integrals_overlap/='None').or.is_periodic)) then if (ao_overlap_abs(i,k) < ao_integrals_threshold) then ao_one_e_integral_zero = .True. return endif endif - if (ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then - ao_one_e_integral_zero = .True. - endif end diff --git a/src/ao_two_e_ints/screening.irp.f b/src/ao_two_e_ints/screening.irp.f index 5095deec..d3230370 100644 --- a/src/ao_two_e_ints/screening.irp.f +++ b/src/ao_two_e_ints/screening.irp.f @@ -8,8 +8,8 @@ logical function ao_two_e_integral_zero(i,j,k,l) ao_two_e_integral_zero = .True. return endif - endif - if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then - ao_two_e_integral_zero = .True. + if (ao_two_e_integral_schwartz(j,l)*ao_two_e_integral_schwartz(i,k) < ao_integrals_threshold) then + ao_two_e_integral_zero = .True. + endif endif end diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index 973a9011..b6e959d7 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -18,89 +18,89 @@ double precision function ao_two_e_integral(i,j,k,l) if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l) - return - endif + else - dim1 = n_pt_max_integrals + dim1 = n_pt_max_integrals - num_i = ao_nucl(i) - num_j = ao_nucl(j) - num_k = ao_nucl(k) - num_l = ao_nucl(l) - ao_two_e_integral = 0.d0 + num_i = ao_nucl(i) + num_j = ao_nucl(j) + num_k = ao_nucl(k) + num_l = ao_nucl(l) + ao_two_e_integral = 0.d0 - if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - I_center(p) = nucl_coord(num_i,p) - J_center(p) = nucl_coord(num_j,p) - K_center(p) = nucl_coord(num_k,p) - L_center(p) = nucl_coord(num_l,p) - enddo + if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + I_center(p) = nucl_coord(num_i,p) + J_center(p) = nucl_coord(num_j,p) + K_center(p) = nucl_coord(num_k,p) + L_center(p) = nucl_coord(num_l,p) + enddo - double precision :: coef1, coef2, coef3, coef4 - double precision :: p_inv,q_inv - double precision :: general_primitive_integral + double precision :: coef1, coef2, coef3, coef4 + double precision :: p_inv,q_inv + double precision :: general_primitive_integral - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & - I_power,J_power,I_center,J_center,dim1) - p_inv = 1.d0/pp - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& - ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & - K_power,L_power,K_center,L_center,dim1) - q_inv = 1.d0/qq - integral = general_primitive_integral(dim1, & - P_new,P_center,fact_p,pp,p_inv,iorder_p, & - Q_new,Q_center,fact_q,qq,q_inv,iorder_q) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p + do p = 1, ao_prim_num(i) + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) + call give_explicit_poly_and_gaussian(P_new,P_center,pp,fact_p,iorder_p,& + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j), & + I_power,J_power,I_center,J_center,dim1) + p_inv = 1.d0/pp + do r = 1, ao_prim_num(k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + call give_explicit_poly_and_gaussian(Q_new,Q_center,qq,fact_q,iorder_q,& + ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l), & + K_power,L_power,K_center,L_center,dim1) + q_inv = 1.d0/qq + integral = general_primitive_integral(dim1, & + P_new,P_center,fact_p,pp,p_inv,iorder_p, & + Q_new,Q_center,fact_q,qq,q_inv,iorder_q) + ao_two_e_integral = ao_two_e_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p - else + else - do p = 1, 3 - I_power(p) = ao_power(i,p) - J_power(p) = ao_power(j,p) - K_power(p) = ao_power(k,p) - L_power(p) = ao_power(l,p) - enddo - double precision :: ERI + do p = 1, 3 + I_power(p) = ao_power(i,p) + J_power(p) = ao_power(j,p) + K_power(p) = ao_power(k,p) + L_power(p) = ao_power(l,p) + enddo + double precision :: ERI - do p = 1, ao_prim_num(i) - coef1 = ao_coef_normalized_ordered_transp(p,i) - do q = 1, ao_prim_num(j) - coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) - do r = 1, ao_prim_num(k) - coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) - do s = 1, ao_prim_num(l) - coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) - integral = ERI( & - ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& - I_power(1),J_power(1),K_power(1),L_power(1), & - I_power(2),J_power(2),K_power(2),L_power(2), & - I_power(3),J_power(3),K_power(3),L_power(3)) - ao_two_e_integral = ao_two_e_integral + coef4 * integral - enddo ! s - enddo ! r - enddo ! q - enddo ! p + do p = 1, ao_prim_num(i) + coef1 = ao_coef_normalized_ordered_transp(p,i) + do q = 1, ao_prim_num(j) + coef2 = coef1*ao_coef_normalized_ordered_transp(q,j) + do r = 1, ao_prim_num(k) + coef3 = coef2*ao_coef_normalized_ordered_transp(r,k) + do s = 1, ao_prim_num(l) + coef4 = coef3*ao_coef_normalized_ordered_transp(s,l) + integral = ERI( & + ao_expo_ordered_transp(p,i),ao_expo_ordered_transp(q,j),ao_expo_ordered_transp(r,k),ao_expo_ordered_transp(s,l),& + I_power(1),J_power(1),K_power(1),L_power(1), & + I_power(2),J_power(2),K_power(2),L_power(2), & + I_power(3),J_power(3),K_power(3),L_power(3)) + ao_two_e_integral = ao_two_e_integral + coef4 * integral + enddo ! s + enddo ! r + enddo ! q + enddo ! p + + endif endif - end double precision function ao_two_e_integral_schwartz_accel(i,j,k,l) @@ -350,71 +350,72 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_in_map ] call map_load_from_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) print*, 'AO integrals provided' ao_two_e_integrals_in_map = .True. - return - endif + else - print*, 'Providing the AO integrals' - call wall_time(wall_0) - call wall_time(wall_1) - call cpu_time(cpu_1) + print*, 'Providing the AO integrals' + call wall_time(wall_0) + call wall_time(wall_1) + call cpu_time(cpu_1) - if (.True.) then - ! Avoid openMP - integral = ao_two_e_integral(1,1,1,1) - endif - - integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull - call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') - - character(len=:), allocatable :: task - allocate(character(len=ao_num*12) :: task) - write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' - do l=1,ao_num - write(task,fmt) (i,l, i=1,l) - integer, external :: add_task_to_taskserver - if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then - stop 'Unable to add task to server' + if (.True.) then + ! Avoid openMP + integral = ao_two_e_integral(1,1,1,1) endif - enddo - deallocate(task) - integer, external :: zmq_set_running - if (zmq_set_running(zmq_to_qp_run_socket) == -1) then - print *, irp_here, ': Failed in zmq_set_running' - endif + integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull + call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'ao_integrals') - PROVIDE nproc - !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1) - i = omp_get_thread_num() - if (i==0) then - call ao_two_e_integrals_in_map_collector(zmq_socket_pull) - else - call ao_two_e_integrals_in_map_slave_inproc(i) + character(len=:), allocatable :: task + allocate(character(len=ao_num*12) :: task) + write(fmt,*) '(', ao_num, '(I5,X,I5,''|''))' + do l=1,ao_num + write(task,fmt) (i,l, i=1,l) + integer, external :: add_task_to_taskserver + if (add_task_to_taskserver(zmq_to_qp_run_socket,trim(task)) == -1) then + stop 'Unable to add task to server' endif - !$OMP END PARALLEL + enddo + deallocate(task) - call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals') + integer, external :: zmq_set_running + if (zmq_set_running(zmq_to_qp_run_socket) == -1) then + print *, irp_here, ': Failed in zmq_set_running' + endif + + PROVIDE nproc + !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+1) + i = omp_get_thread_num() + if (i==0) then + call ao_two_e_integrals_in_map_collector(zmq_socket_pull) + else + call ao_two_e_integrals_in_map_slave_inproc(i) + endif + !$OMP END PARALLEL + + call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'ao_integrals') - print*, 'Sorting the map' - call map_sort(ao_integrals_map) - call cpu_time(cpu_2) - call wall_time(wall_2) - integer(map_size_kind) :: get_ao_map_size, ao_map_size - ao_map_size = get_ao_map_size() + print*, 'Sorting the map' + call map_sort(ao_integrals_map) + call cpu_time(cpu_2) + call wall_time(wall_2) + integer(map_size_kind) :: get_ao_map_size, ao_map_size + ao_map_size = get_ao_map_size() - print*, 'AO integrals provided:' - print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' - print*, ' Number of AO integrals :', ao_map_size - print*, ' cpu time :',cpu_2 - cpu_1, 's' - print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' + print*, 'AO integrals provided:' + print*, ' Size of AO map : ', map_mb(ao_integrals_map) ,'MB' + print*, ' Number of AO integrals :', ao_map_size + print*, ' cpu time :',cpu_2 - cpu_1, 's' + print*, ' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1+tiny(1.d0)), ' )' - ao_two_e_integrals_in_map = .True. + ao_two_e_integrals_in_map = .True. + + if (write_ao_two_e_integrals.and.mpi_master) then + call ezfio_set_work_empty(.False.) + call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) + call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') + endif - if (write_ao_two_e_integrals.and.mpi_master) then - call ezfio_set_work_empty(.False.) - call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints',ao_integrals_map) - call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals('Read') endif END_PROVIDER diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index a9983e51..b926d688 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -189,7 +189,6 @@ subroutine add_integrals_to_map(mask_ijkl) two_e_tmp_2 = 0.d0 do j1 = 1,ao_num call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) - ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) enddo do j1 = 1,ao_num kmax = 0 @@ -747,7 +746,6 @@ subroutine add_integrals_to_map_no_exit_34(mask_ijkl) two_e_tmp_2 = 0.d0 do j1 = 1,ao_num call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) - ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) enddo do j1 = 1,ao_num kmax = 0 From ceeef553b1e31c09a5cf8b7f1ffeca0601c4183b Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 15 May 2020 14:53:55 +0200 Subject: [PATCH 10/11] Increased size of strings --- src/ezfio_files/ezfio.irp.f | 6 +++--- src/ezfio_files/get_unit_and_open.irp.f | 2 +- src/ezfio_files/qp_stop.irp.f | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/ezfio_files/ezfio.irp.f b/src/ezfio_files/ezfio.irp.f index e7a7a159..4f53b173 100644 --- a/src/ezfio_files/ezfio.irp.f +++ b/src/ezfio_files/ezfio.irp.f @@ -1,4 +1,4 @@ -BEGIN_PROVIDER [ character*(128), ezfio_filename ] +BEGIN_PROVIDER [ character*(1024), ezfio_filename ] implicit none BEGIN_DOC ! Name of EZFIO file. It is obtained from the QPACKAGE_INPUT environment @@ -34,7 +34,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ] ! Adjust out-of-memory killer flag such that the current process will be ! killed first by the OOM killer, allowing compute nodes to survive integer :: getpid - character*(64) :: command, pidc + character*(1024) :: command, pidc write(pidc,*) getpid() write(command,*) 'echo 15 > /proc//'//trim(adjustl(pidc))//'/oom_adj' call system(command) @@ -43,7 +43,7 @@ BEGIN_PROVIDER [ character*(128), ezfio_filename ] END_PROVIDER -BEGIN_PROVIDER [ character*(128), ezfio_work_dir ] +BEGIN_PROVIDER [ character*(1024), ezfio_work_dir ] implicit none BEGIN_DOC ! EZFIO/work/ diff --git a/src/ezfio_files/get_unit_and_open.irp.f b/src/ezfio_files/get_unit_and_open.irp.f index 71bbc76c..6440579f 100644 --- a/src/ezfio_files/get_unit_and_open.irp.f +++ b/src/ezfio_files/get_unit_and_open.irp.f @@ -17,7 +17,7 @@ integer function getUnitAndOpen(f,mode) END_DOC character*(*) :: f - character*(128) :: new_f + character*(256) :: new_f integer :: iunit logical :: is_open, exists character :: mode diff --git a/src/ezfio_files/qp_stop.irp.f b/src/ezfio_files/qp_stop.irp.f index b2618246..c1039dc3 100644 --- a/src/ezfio_files/qp_stop.irp.f +++ b/src/ezfio_files/qp_stop.irp.f @@ -1,5 +1,5 @@ - BEGIN_PROVIDER [ character*(128), qp_stop_filename ] -&BEGIN_PROVIDER [ character*(128), qp_kill_filename ] + BEGIN_PROVIDER [ character*(256), qp_stop_filename ] +&BEGIN_PROVIDER [ character*(256), qp_kill_filename ] &BEGIN_PROVIDER [ integer, qp_stop_variable ] implicit none BEGIN_DOC From c614cb39226e8c03badac1e782e0bb56ec21add0 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 15 May 2020 15:03:06 +0200 Subject: [PATCH 11/11] Fixed conversion --- bin/qp_convert_output_to_ezfio | 2 +- src/nuclei/nuclei.irp.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 5f0c4771..cbc81032 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -335,7 +335,7 @@ def write_ezfio(res, filename): def get_full_path(file_path): file_path = os.path.expanduser(file_path) file_path = os.path.expandvars(file_path) - file_path = os.path.abspath(file_path) +# file_path = os.path.abspath(file_path) return file_path diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index 88ee24e8..39e6b399 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -211,7 +211,7 @@ END_PROVIDER END_DOC integer :: iunit, i integer, external :: getUnitAndOpen - character*(128) :: filename + character*(1024) :: filename if (mpi_master) then call getenv('QP_ROOT',filename) filename = trim(filename)//'/data/list_element.txt'