diff --git a/src/MonoInts/int.f90 b/src/MonoInts/int.f90 index b67e6102..640688d0 100644 --- a/src/MonoInts/int.f90 +++ b/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)) diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index f430ace9..6ad38d6a 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -8,19 +8,43 @@ 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 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 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) + !$OMP n_pt_max_integrals,ao_nucl_elec_integral,nucl_num,nucl_charge) n_pt_in = n_pt_max_integrals !$OMP DO SCHEDULE (guided) do j = 1, ao_num @@ -50,7 +74,15 @@ C_center(1) = nucl_coord(k,1) 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) + 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 diff --git a/src/MonoInts/test_michel.irp.f b/src/MonoInts/test_michel.irp.f index e290ea4a..c7e3c685 100644 --- a/src/MonoInts/test_michel.irp.f +++ b/src/MonoInts/test_michel.irp.f @@ -137,19 +137,31 @@ program compute_integrals_pseudo ! \_ (_| | (_ |_| | ! - write(*,*)'a?' - read*,a(1),a(2),a(3) +! write(*,*)'a?' +! read*,a(1),a(2),a(3) !write(*,*)'b?' !read*,b(1),b(2),b(3) - b(1)=-0.1d0 - b(2)=-0.2d0 - b(3)=0.3d0 - !write(*,*)'a?' - !read*,c(1),c(2),c(3) - c(1)=0.1d0 - c(2)=0.2d0 - c(3)=0.3d0 +! b(1)=-0.1d0 +! b(2)=-0.2d0 +! b(3)=0.3d0 +! !write(*,*)'a?' +! !read*,c(1),c(2),c(3) +! c(1)=0.1d0 +! c(2)=0.2d0 +! 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