added the possibility to use no_vvvv integrals from EZFIO

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
+

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]

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

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

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
+

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

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

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)
+ sum = array
+ i = 0
+ do while(ntmp.ge.2)
+  i += 1
+!  print*,'i = ',i
+  call shank(sum,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 = shank1
+  shank_general = shank1(ntmp)
+ enddo
 end