From e864eb1cf383d69c80f9edacaa93211ccfb7b8e7 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 11 May 2020 16:04:16 +0200 Subject: [PATCH] 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