New int.f90 working with A=B=C

src/MonoInts/int.f90

@@ -560,12 +560,22 @@ double precision int_prod_bessel_loc,binom,accu,prod,ylm,bigI,arg
  Vloc=0.d0
  return
  endif
- freal=dexp(-g_a*ac**2-g_b*bc**2)
 
  ntotA=n_a(1)+n_a(2)+n_a(3)
  ntotB=n_b(1)+n_b(2)+n_b(3)
  ntot=ntotA+ntotB
 
+ if(ac.eq.0.d0.and.bc.eq.0.d0)then
+  accu=0.d0
+  do k=1,klocmax 
+   accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k))
+  enddo
+  Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3))
+  return
+ endif
+
+ freal=dexp(-g_a*ac**2-g_b*bc**2)
+ 
  d2=0.d0
  do i=1,3
   d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i))

src/MonoInts/pot_ao_ints.irp.f

@@ -8,17 +8,41 @@
  double precision :: A_center(3),B_center(3),C_center(3)
  integer :: power_A(3),power_B(3)
  integer :: i,j,k,l,n_pt_in,m
- double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult
+ double precision ::overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult, Vloc
  integer :: nucl_numC
  ! Important for OpenMP
 
  ao_nucl_elec_integral = 0.d0
 
+  !                 
+  ! |   _   _  _. | 
+  ! |_ (_) (_ (_| | 
+  !                 
+
+  integer klocmax
+  integer, allocatable ::  n_k(:)
+  double precision, allocatable ::  v_k(:), dz_k(:)
+
+  call ezfio_get_pseudo_klocmax(klocmax)
+
+  allocate(n_k(klocmax),v_k(klocmax), dz_k(klocmax))
+
+  call ezfio_get_pseudo_v_k(v_k)
+  call ezfio_get_pseudo_n_k(n_k)
+  call ezfio_get_pseudo_dz_k(dz_k)
+
+  print*, "klocmax", klocmax
+
+  print*, "n_k_ezfio", n_k
+  print*, "v_k_ezfio",v_k
+  print*, "dz_k_ezfio", dz_k
+
 
  !$OMP PARALLEL &
  !$OMP DEFAULT (NONE) &
  !$OMP PRIVATE (i, j, k, l, m, alpha, beta, A_center, B_center, C_center, power_A, power_B, &
- !$OMP  num_A,num_B,Z,c,n_pt_in) &
+ !$OMP          num_A, num_B, Z, c, n_pt_in,&
+ !$OMP          v_k, n_k, dz_k, klocmax) &
  !$OMP SHARED (ao_num,ao_prim_num,ao_expo_transp,ao_power,ao_nucl,nucl_coord,ao_coef_transp, &
  !$OMP         n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge)
  n_pt_in = n_pt_max_integrals
@@ -51,6 +75,14 @@
       C_center(2) = nucl_coord(k,2)
       C_center(3) = nucl_coord(k,3)
       c = c + Z*NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in)
+
+      print*, A_center
+      print*, B_center
+      print*, C_center
+
+      print*, Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center)
+!      c = c + Z*Vloc(klocmax,v_k,n_k,dz_k,A_center,power_A,alpha,B_center,power_B,beta,C_center)
+
      enddo
      ao_nucl_elec_integral(i,j) = ao_nucl_elec_integral(i,j) - &
        ao_coef_transp(l,j)*ao_coef_transp(m,i)*c

src/MonoInts/test_michel.irp.f

@@ -137,18 +137,30 @@ program compute_integrals_pseudo
   ! \_ (_| | (_ |_| | 
   !                   
 
-  write(*,*)'a?'
+!  write(*,*)'a?'
-  read*,a(1),a(2),a(3)
+!  read*,a(1),a(2),a(3)
   !write(*,*)'b?'
   !read*,b(1),b(2),b(3)
-  b(1)=-0.1d0
+!  b(1)=-0.1d0
-  b(2)=-0.2d0
+!  b(2)=-0.2d0
-  b(3)=0.3d0
+!  b(3)=0.3d0
-  !write(*,*)'a?'
+!  !write(*,*)'a?'
-  !read*,c(1),c(2),c(3)
+!  !read*,c(1),c(2),c(3)
-  c(1)=0.1d0
+!  c(1)=0.1d0
-  c(2)=0.2d0
+!  c(2)=0.2d0
-  c(3)=0.3d0
+!  c(3)=0.3d0
+  
+  a(1)= 0.d0
+  a(2)= 0.d0
+  a(3)= 0.d0
+
+  b(1)= 0.d0
+  b(2)= 0.d0
+  b(3)= 0.d0
+
+  c(1)= 0.d0
+  c(2)= 0.d0
+  c(3)= 0.d0
 
   print*,'ntps? rmax for brute force integration'
   read*,npts,rmax