Added module DensityFit

scripts/ezfio_interface/ei_handler.py

@@ -368,7 +368,10 @@ def create_ezfio_stuff(dict_ezfio_cfg, config_or_default="config"):
                 except ValueError:
                     a_size_raw.append(dim)
                 else:
-                    a_size_raw.append("{0}+{1}+1".format(begin, end))
+                    if begin[0] == '-':
+                        a_size_raw.append("{0}+{1}+1".format(end, begin[1:]))
+                    else:
+                       a_size_raw.append("{0}-{1}+1".format(end, begin))
 
             size_raw = ",".join(a_size_raw)
 

src/DensityFit/Makefile

@@ -0,0 +1,6 @@
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/DensityFit/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+AOs

src/DensityFit/README.rst

@@ -0,0 +1,24 @@
+=================
+DensityFit Module
+=================
+
+In this module, the basis of all the products of atomic orbitals is built.
+
+Documentation
+=============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+
+
+Needed Modules
+==============
+
+.. Do not edit this section. It was auto-generated from the
+.. NEEDED_MODULES file.
+
+.. image:: tree_dependancy.png
+
+* `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+

src/DensityFit/aux_basis.irp.f

@@ -0,0 +1,54 @@
+ BEGIN_PROVIDER [ integer, aux_basis_num ]
+&BEGIN_PROVIDER [ integer, aux_basis_num_8 ]
+ implicit none
+ BEGIN_DOC
+ ! Number of auxiliary basis functions
+ END_DOC
+ integer :: align_double
+ aux_basis_num = ao_num * (ao_num+1)/2
+ aux_basis_num_8 = align_double(aux_basis_num)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, aux_basis_idx, (2,aux_basis_num) ]
+ implicit none
+ BEGIN_DOC
+! aux_basis_idx(k) -> i,j
+ END_DOC
+ integer :: i,j,k
+ k=0
+ do j=1,ao_num
+   do i=1,j
+     k = k+1
+     aux_basis_idx(1,k) = i
+     aux_basis_idx(2,k) = j
+   enddo
+ enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, aux_basis_overlap_matrix, (aux_basis_num_8,aux_basis_num) ]
+ implicit none
+ BEGIN_DOC  
+! Auxiliary basis set
+ END_DOC
+ integer :: m,n,i,j,k,l
+ double precision :: ao_four_overlap
+
+ aux_basis_overlap_matrix(1,1) = ao_four_overlap(1,1,1,1)
+ !$OMP PARALLEL DO PRIVATE(i,j,k,l,m,n) SCHEDULE(GUIDED)
+ do m=1,aux_basis_num
+   i = aux_basis_idx(1,m)
+   j = aux_basis_idx(2,m)
+   do n=1,m
+     k = aux_basis_idx(1,n)
+     l = aux_basis_idx(2,n)
+     aux_basis_overlap_matrix(m,n) = ao_four_overlap(i,j,k,l)
+     aux_basis_overlap_matrix(n,m) = aux_basis_overlap_matrix(m,n) 
+   enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+
+
+

src/DensityFit/overlap.irp.f

@@ -0,0 +1,66 @@
+double precision function ao_four_overlap(i,j,k,l)
+  implicit none
+  BEGIN_DOC
+! \int \chi_i(r) \chi_j(r) \chi_k(r) \chi_l(r) dr
+  END_DOC
+  integer,intent(in)             :: i,j,k,l
+  integer                        :: p,q,r,s
+  double precision               :: I_center(3),J_center(3),K_center(3),L_center(3)
+  integer                        :: num_i,num_j,num_k,num_l,dim1,I_power(3),J_power(3),K_power(3),L_power(3)
+  double precision               :: overlap_x,overlap_y,overlap_z, overlap
+  include 'include/constants.F'
+  double precision               :: P_new(0:max_dim,3),P_center(3),fact_p,pp
+  double precision               :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
+  integer                        :: iorder_p(3), iorder_q(3)
+  
+  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_four_overlap = 0.d0
+  
+  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
+  
+  do p = 1, ao_prim_num(i)
+    double precision               :: coef1
+    coef1 = ao_coef_normalized_ordered_transp(p,i)
+    do q = 1, ao_prim_num(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)
+      double precision               :: coef2
+      coef2 = coef1*ao_coef_normalized_ordered_transp(q,j)*fact_p
+      do r = 1, ao_prim_num(k)
+        double precision               :: coef3
+        coef3 = coef2*ao_coef_normalized_ordered_transp(r,k)
+        do s = 1, ao_prim_num(l)
+          double precision               :: general_primitive_integral
+          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)
+          double precision               :: coef4
+          coef4 = coef3*ao_coef_normalized_ordered_transp(s,l)*fact_q
+          call overlap_gaussian_xyz(P_center,Q_center,pp,qq,iorder_p,iorder_q,overlap_x,overlap_y,overlap_z,overlap,dim1)
+          ao_four_overlap +=  coef4 * overlap 
+        enddo ! s
+      enddo  ! r
+    enddo   ! q
+  enddo    ! p
+    
+end
+
+! TODO : Schwartz acceleration
+
+
+

src/FCIdump/README.rst

@@ -10,9 +10,6 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`fcidump <http://github.com/LCPQ/quantum_package/tree/master/src/FCIdump/fcidump.irp.f#L1>`_
-  Undocumented
-
 
 
 Needed Modules

src/Integrals_Monoelec/README.rst

@@ -103,6 +103,12 @@ Documentation
 `ao_pseudo_integral_non_local <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L119>`_
   Local pseudo-potential
 
+`pseudo_grid <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f#L223>`_
+  Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
+  .br
+  <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
+  title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
+
 `mo_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/Integrals_Monoelec/pot_mo_ints.irp.f#L1>`_
   interaction nuclear electron on the MO basis
 

src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f

@@ -4,7 +4,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
 ! Pseudo-potential
   END_DOC
   if (do_pseudo) then
-    ao_pseudo_integral = ao_pseudo_integral_local + ao_pseudo_integral_non_local 
+    ao_pseudo_integral = ao_pseudo_integral_local !+ ao_pseudo_integral_non_local 
   else
     ao_pseudo_integral = 0.d0
   endif
@@ -220,4 +220,67 @@ END_PROVIDER
 
  END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, pseudo_grid, (pseudo_grid_size,ao_num,-pseudo_lmax:pseudo_lmax,0:pseudo_lmax,nucl_num)  ]
+ implicit none
+ BEGIN_DOC
+! Grid points for f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C
+!
+! <img src="http://latex.codecogs.com/gif.latex?f(|r-r_A|)&space;=&space;\int&space;Y_{lm}^{C}&space;(|r-r_C|,&space;\Omega_C)&space;\chi_i^{A}&space;(r-r_A)&space;d\Omega_C"
+! title="f(|r-r_A|) = \int Y_{lm}^{C} (|r-r_C|, \Omega_C) \chi_i^{A} (r-r_A) d\Omega_C" />
+ END_DOC
+ ! l,m    : Y(l,m) parameters
+ ! c(3)   : pseudopotential center
+ ! a(3)   : Atomic Orbital center
+ ! n_a(3) : Powers of x,y,z in the Atomic Orbital
+ ! g_a    : Atomic Orbital exponent      
+ ! r      : Distance between the Atomic Orbital center and the considered point
+ double precision, external :: ylm_orb
+ integer :: n_a(3)
+ double precision :: a(3), c(3), g_a
+ integer :: i,j,k,l,m,n,p
+ double precision :: r(pseudo_grid_size), dr, Ulc
+ double precision, parameter :: rmax= 10.d0
+
+ dr = rmax/dble(pseudo_grid_size)
+ r(1) = 0.d0
+ do j=2,pseudo_grid_size
+   r(j) = r(j-1) + dr
+ enddo
+
+ pseudo_grid = 0.d0
+ do k=1,nucl_num
+   c(1:3) = nucl_coord(k,1:3)
+   do l=0,pseudo_lmax
+     do i=1,ao_num
+       a(1:3) = nucl_coord(ao_nucl(i),1:3)
+       n_a(1:3) = ao_power(i,1:3)
+       do j=1,pseudo_grid_size
+         do p=1,ao_prim_num(i)
+           g_a = ao_expo_ordered_transp(p,i)
+           do m=-l,l
+             double precision               :: y
+             !            y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
+             !            if (y /= 0.d0) then
+             !              print *, 'y = ', y
+             !              print *, 'l = ', l
+             !              print *, 'm = ', m
+             !              print *, 'c = ', c
+             !              print *, 'a = ', a
+             !              print *, 'n = ', n_a
+             !              print *, 'g = ', g_a
+             !              print *, 'r = ', r(j)
+             !              print *,  ''
+             !            endif
+             y = ylm_orb(l,m,c,a,n_a,g_a,r(j))
+             pseudo_grid(j,i,m,l,k) = pseudo_grid(j,i,m,l,k) +           &
+                 ao_coef_normalized_ordered_transp(p,i)*y
+           enddo
+         enddo
+       enddo
+     enddo
+   enddo
+ enddo
+
+END_PROVIDER
+
 

src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f

@@ -28,5 +28,7 @@ BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_n
     enddo
   enddo
   !$OMP END PARALLEL DO
+  call ezfio_set_pseudo_pseudo_matrix(mo_pseudo_integral(1:mo_tot_num,1:mo_tot_num))
 END_PROVIDER
 
+

src/Integrals_Monoelec/pseudopot.f90

@@ -201,7 +201,7 @@ double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax)
 ! |_ (_) (_ (_| | (/_ 
 !                     
 
-double precision :: fourpi,f,prod,prodp,binom,accu,bigR,bigI,ylm
+double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm
 double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big
 double precision :: areal,freal,breal,t1,t2,int_prod_bessel
 double precision :: arg
@@ -327,7 +327,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
  do k1=0,n_a(1)
  do k2=0,n_a(2)
  do k3=0,n_a(3)
-  array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3)  &
+  array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
                           *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)
  enddo
  enddo
@@ -336,7 +336,7 @@ else if(ac.ne.0.d0.and.bc.ne.0.d0)then
  do k1p=0,n_b(1)
  do k2p=0,n_b(2)
  do k3p=0,n_b(3)
-  array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p)  &
+  array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)  &
                              *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
  enddo
  enddo
@@ -448,7 +448,7 @@ else if(ac.eq.0.d0.and.bc.ne.0.d0)then
  do k2p=0,n_b(2)
  do k3p=0,n_b(3)
 
-  array_coefs_B(k1p,k2p,k3p)=binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p)  &
+  array_coefs_B(k1p,k2p,k3p)=binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)  &
                              *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
  enddo
  enddo
@@ -534,7 +534,7 @@ else if(ac.ne.0.d0.and.bc.eq.0.d0)then
  do k2=0,n_a(2)
  do k3=0,n_a(3)
 
-  array_coefs_A(k1,k2,k3)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3)  &
+  array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
                           *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)
 
  enddo
@@ -705,7 +705,7 @@ integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p
 double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0
 double precision,allocatable :: array_R_loc(:,:,:)
 double precision,allocatable :: array_coefs(:,:,:,:,:,:)
-double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg
+double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg
 
  fourpi=4.d0*dacos(-1.d0)
  f=fourpi**1.5d0
@@ -762,9 +762,9 @@ allocate (array_coefs(0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:ntot_max,0:n
  do k1p=0,n_b(1)
  do k2p=0,n_b(2)
  do k3p=0,n_b(3)
-  array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom(n_a(1),k1)*binom(n_a(2),k2)*binom(n_a(3),k3)  &
+  array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
   *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) &
-  *binom(n_b(1),k1p)*binom(n_b(2),k2p)*binom(n_b(3),k3p)  &
+  *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)  &
   *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p)
  enddo
  enddo
@@ -823,7 +823,7 @@ double precision function bigI(lambda,mu,l,m,k1,k2,k3)
 implicit none
 integer lambda,mu,l,m,k1,k2,k3
 integer k,i,kp,ip
-double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom,fact,coef_pm
+double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm
 pi=dacos(-1.d0)
 
 if(mu.gt.0.and.m.gt.0)then
@@ -834,8 +834,8 @@ do k=0,mu/2
  do i=0,lambda-mu
   do kp=0,m/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3)
    enddo
   enddo
@@ -868,7 +868,7 @@ do i=0,lambda
  do kp=0,m/2
   do ip=0,l-m
    cylm=factor1*coef_pm(lambda,i)
-   cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+   cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
    sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3)
   enddo
  enddo
@@ -884,7 +884,7 @@ factor2=dsqrt((2*l+1)/(4.d0*pi))
 do k=0,mu/2
  do i=0,lambda-mu
   do ip=0,l
-   cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+   cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
    cylmp=factor2*coef_pm(l,ip)
    sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3)
   enddo
@@ -904,8 +904,8 @@ do k=0,(mu-1)/2
  do i=0,lambda-mu
   do kp=0,(m-1)/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3)
    enddo
   enddo
@@ -926,7 +926,7 @@ do i=0,lambda
  do kp=0,(m-1)/2
   do ip=0,l-m
    cylm=factor1*coef_pm(lambda,i)
-   cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+   cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
    sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3)
   enddo
  enddo
@@ -944,7 +944,7 @@ factor2=dsqrt((2*l+1)/(4.d0*pi))
 do k=0,(mu-1)/2
  do i=0,lambda-mu
    do ip=0,l
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
     cylmp=factor2*coef_pm(l,ip)
     sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3)
    enddo
@@ -964,8 +964,8 @@ do k=0,mu/2
  do i=0,lambda-mu
   do kp=0,(m-1)/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3)
    enddo
   enddo
@@ -985,8 +985,8 @@ do k=0,(mu-1)/2
  do i=0,lambda-mu
   do kp=0,m/2
    do ip=0,l-m
-    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
+    cylm=(-1.d0)**k*factor1*dsqrt(2.d0)*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu)
-    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
+    cylmp=(-1.d0)**kp*factor2*dsqrt(2.d0)*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m)
     sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3)
    enddo
   enddo
@@ -1485,6 +1485,7 @@ end
         a=bessel_mod_exp(n,x)
         return
        endif
+       print *,  n,x
        if(n.eq.0)a=dsinh(x)/x
        if(n.eq.1)a=(x*dcosh(x)-dsinh(x))/x**2
        if(n.ge.2)a=bessel_mod_recur(n-2,x)-(2*n-1)/x*bessel_mod_recur(n-1,x)
@@ -1699,14 +1700,14 @@ end
 double precision function coef_pm(n,k)
 implicit none
 integer n,k
-double precision arg,binom,binom_gen
+double precision arg,binom_func,binom_gen
 if(n.eq.0.and.k.ne.0)stop 'coef_pm not defined'
 if(n.eq.0.and.k.eq.0)then
 coef_pm=1.d0
 return
 endif
 arg=0.5d0*dfloat(n+k-1)
-coef_pm=2.d0**n*binom(n,k)*binom_gen(arg,n)
+coef_pm=2.d0**n*binom_func(n,k)*binom_gen(arg,n)
 end
 
 !! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k
@@ -1735,7 +1736,7 @@ end
 double precision function ylm_bis(l,m,theta,phi)
 implicit none
 integer l,m,k,i
-double precision x,y,z,theta,phi,sum,factor,pi,binom,fact,coef_pm,cylm
+double precision x,y,z,theta,phi,sum,factor,pi,binom_func,fact,coef_pm,cylm
 pi=dacos(-1.d0)
 x=dsin(theta)*dcos(phi)
 y=dsin(theta)*dsin(phi)
@@ -1745,7 +1746,7 @@ if(m.gt.0)then
 sum=0.d0
 do k=0,m/2
  do i=0,l-m
-  cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m)
+  cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m)
   sum=sum+cylm*x**(m-2*k)*y**(2*k)*z**i
  enddo
 enddo
@@ -1765,7 +1766,7 @@ m=-m
 sum=0.d0
 do k=0,(m-1)/2
  do i=0,l-m
-  cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m)
+  cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m)
   sum=sum+cylm*x**(m-(2*k+1))*y**(2*k+1)*z**i
  enddo
 enddo
@@ -1800,13 +1801,6 @@ end
 !!
 !!   with  a_k= 2^n  binom(n,k) binom( (n+k-1)/2, n )
 
-double precision function binom(i,j)
- implicit none
- integer :: i,j
- double precision :: fact
- binom = fact(i)/(fact(j)*fact(i-j))
-end
-
 double precision function binom_gen(alpha,n)
  implicit none
  integer :: n,k
@@ -2087,4 +2081,90 @@ end
       end
 
 
+! l,m    : Y(l,m) parameters
+! c(3)   : pseudopotential center
+! a(3)   : Atomic Orbital center
+! n_a(3) : Powers of x,y,z in the Atomic Orbital
+! g_a    : Atomic Orbital exponent
+! r      : Distance between the Atomic Orbital center and the considered point
+double precision function ylm_orb(l,m,c,a,n_a,g_a,r)
+implicit none
+integer lmax_max,ntot_max
+parameter (lmax_max=2)
+parameter (ntot_max=14)
+integer l,m
+double precision a(3),g_a,c(3)
+double precision prod,binom_func,accu,bigI,ylm,bessel_mod
+double precision theta_AC0,phi_AC0,ac,factor,fourpi,arg,r,areal
+integer ntotA,mu,k1,k2,k3,lambda
+integer n_a(3)
+double precision &
+array_I_A(0:lmax_max+ntot_max,-(lmax_max+ntot_max):lmax_max+ntot_max,0:ntot_max,0:ntot_max,0:ntot_max)
+double precision array_coefs_A(0:ntot_max,0:ntot_max,0:ntot_max), y
 
+ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2)
+arg=g_a*(ac**2+r**2)
+fourpi=4.d0*dacos(-1.d0)
+factor=fourpi*dexp(-arg)
+areal=2.d0*g_a*ac
+ntotA=n_a(1)+n_a(2)+n_a(3)
+
+if(ntotA.gt.ntot_max)stop 'increase ntot_max'
+
+if(ac.eq.0.d0)then
+ ylm_orb=dsqrt(fourpi)*r**ntotA*dexp(-g_a*r**2)*bigI(0,0,l,m,n_a(1),n_a(2),n_a(3))
+ return
+else
+
+ theta_AC0=dacos( (a(3)-c(3))/ac )
+ phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac)
+
+ do k1=0,n_a(1)
+  do k2=0,n_a(2)
+   do k3=0,n_a(3)
+  array_coefs_A(k1,k2,k3)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)  &
+  *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3) &
+  *r**(k1+k2+k3)
+   enddo
+  enddo
+ enddo
+
+ do lambda=0,l+ntotA
+  do mu=-lambda,lambda
+   do k1=0,n_a(1)
+   do k2=0,n_a(2)
+   do k3=0,n_a(3)
+    array_I_A(lambda,mu,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3)
+   enddo
+   enddo
+   enddo
+  enddo
+ enddo
+
+ accu=0.d0
+ do lambda=0,l+ntotA
+  do mu=-lambda,lambda
+   y = ylm(lambda,mu,theta_AC0,phi_AC0)
+   if (y == 0.d0) then
+     cycle
+   endif
+   do k1=0,n_a(1)
+   do k2=0,n_a(2)
+   do k3=0,n_a(3)
+    prod=y*array_coefs_A(k1,k2,k3)*array_I_A(lambda,mu,k1,k2,k3)
+    if (prod == 0.d0) then
+      cycle
+    endif
+    if (areal*r < 100.d0) then ! overflow!
+      accu=accu+prod*bessel_mod(areal*r,lambda)
+    endif
+   enddo
+   enddo
+   enddo
+  enddo
+ enddo
+ ylm_orb=factor*accu
+ return
+endif
+
+end

src/NEEDED_MODULES

@@ -1 +1 @@
-AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils
+AOs Integrals_Bielec Bitmask CAS_SD CID CID_SC2_selected CID_selected CIS CISD CISD_SC2_selected CISD_selected DDCI_selected Determinants Electrons Ezfio_files FCIdump Full_CI Generators_CAS Generators_full Hartree_Fock MOGuess Molden Integrals_Monoelec MOs MP2 MRCC Nuclei Pseudo Selectors_full Utils DensityFit

src/Pseudo/EZFIO.cfg

@@ -55,3 +55,22 @@ doc:  Using pseudo potential integral of not
 interface: input
 default: False
 
+[pseudo_grid_size]
+type: integer
+doc: Size of the QMC grid 
+interface: input
+default: 100
+
+[pseudo_grid]
+type: double precision
+doc: QMC grid 
+interface: output
+size: (pseudo.pseudo_grid_size,ao_basis.ao_num,-pseudo.pseudo_lmax:pseudo.pseudo_lmax,0:pseudo.pseudo_lmax,nuclei.nucl_num)
+
+[pseudo_matrix]
+type: double precision
+doc: QMC grid 
+interface: output
+size: (mo_basis.mo_tot_num,mo_basis.mo_tot_num)
+
+

src/Utils/one_e_integration.irp.f

@@ -131,13 +131,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
   integer                        :: iorder_p(3)
   
   call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
-  if(fact_p.lt.1d-10)then
+!  if(fact_p.lt.1d-20)then
-    overlap_x = 0.d0
+!    overlap_x = 0.d0
-    overlap_y = 0.d0
+!    overlap_y = 0.d0
-    overlap_z = 0.d0
+!    overlap_z = 0.d0
-    overlap = 0.d0
+!    overlap = 0.d0
-    return
+!    return
-  endif
+!  endif
   integer                        :: nmax
   double precision               :: F_integral
   nmax = maxval(iorder_p)