From 3b75503bc189527aaaef0f0d03934ec1ef544978 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Fri, 7 Aug 2020 11:47:51 +0200
Subject: [PATCH] minor changes

---
 src/determinants/psi_energy_mono_elec.irp.f | 30 +++++++++++++++++++++
 src/utils/integration.irp.f                 | 10 ++++++-
 2 files changed, 39 insertions(+), 1 deletion(-)

diff --git a/src/determinants/psi_energy_mono_elec.irp.f b/src/determinants/psi_energy_mono_elec.irp.f
index 74e69160..fc0ffba9 100644
--- a/src/determinants/psi_energy_mono_elec.irp.f
+++ b/src/determinants/psi_energy_mono_elec.irp.f
@@ -27,3 +27,33 @@
   psi_energy_h_core(i) = psi_energy_h_core(i) * accu
  enddo
 END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, v_ne_psi_energy, (N_states) ]
+  implicit none
+  integer :: i
+  integer :: j,k
+  double precision :: tmp(mo_num,mo_num),mono_ints(mo_num,mo_num)
+  BEGIN_DOC
+! v_ne_psi_energy = $\langle \Psi | v_ne |\Psi \rangle$
+!
+! computed using the :c:data:`one_e_dm_mo_alpha` +
+! :c:data:`one_e_dm_mo_beta` and :c:data:`mo_one_e_integrals`
+  END_DOC
+  v_ne_psi_energy = 0.d0
+  do i = 1, N_states
+   do j = 1, mo_num
+    do k = 1, mo_num
+     v_ne_psi_energy(i) += mo_integrals_n_e(k,j) * (one_e_dm_mo_alpha(k,j,i) + one_e_dm_mo_beta(k,j,i))
+    enddo
+   enddo
+  enddo
+ double precision :: accu
+ do i = 1, N_states
+  accu = 0.d0
+  do j = 1, mo_num
+   accu += one_e_dm_mo_alpha(j,j,i) + one_e_dm_mo_beta(j,j,i)
+  enddo
+  accu = (elec_alpha_num + elec_beta_num ) / accu
+  v_ne_psi_energy(i) = v_ne_psi_energy(i) * accu
+ enddo
+END_PROVIDER
diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f
index c907e425..da8120a1 100644
--- a/src/utils/integration.irp.f
+++ b/src/utils/integration.irp.f
@@ -30,7 +30,11 @@ subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alph
   ab = alpha * beta
   d_AB = (A_center - B_center) * (A_center - B_center)
   P_center = (alpha * A_center + beta * B_center) * p_inv
-  fact_k = exp(-ab*p_inv * d_AB)
+  if(dabs(ab*p_inv * d_AB).lt.50.d0)then
+   fact_k = exp(-ab*p_inv * d_AB)
+  else
+   fact_k = 0.d0
+  endif
 
   ! Recenter the polynomials P_a and P_b on x
   !DIR$ FORCEINLINE
@@ -78,6 +82,10 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   !DIR$ FORCEINLINE
   call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
   if (fact_k < thresh) then
+    P_center = 0.d0
+    p = 1.d-10
+    P_new = 0.d0
+    iorder = -1
     fact_k = 0.d0
     return
   endif