From 889ad27c20c39214d04fe2d1ea51548dae60d972 Mon Sep 17 00:00:00 2001
From: Antoine Marie <amarie@irsamc.ups-tlse>
Date: Thu, 14 Dec 2023 10:14:00 +0100
Subject: [PATCH] fixing nuclear energy

---
 PyDuck.py                   |  7 +++++++
 src/utils/read_geometry.f90 | 15 +++++++++------
 2 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/PyDuck.py b/PyDuck.py
index 809a7a8..1f77a4b 100644
--- a/PyDuck.py
+++ b/PyDuck.py
@@ -76,6 +76,13 @@ for i in range(len(list_pos_atom)):
     f.write(list_pos_atom[i][0]+' '+str(list_pos_atom[i][1][0])+' '+str(list_pos_atom[i][1][1])+' '+str(list_pos_atom[i][1][2])+'\n')
 f.close()
 
+#Compute nuclear energy and put it in a file
+subprocess.call(['rm', working_dir + '/int/ENuc.dat'])
+f = open(working_dir+'/int/ENuc.dat','w')
+f.write(mol.energy_nuc())
+f.write(' ')
+f.close()
+
 #Compute 1e integrals
 ovlp = mol.intor('int1e_ovlp')#Overlap matrix elements
 v1e  = mol.intor('int1e_nuc') #Nuclear repulsion matrix elements
diff --git a/src/utils/read_geometry.f90 b/src/utils/read_geometry.f90
index 5c2da13..5ee3916 100644
--- a/src/utils/read_geometry.f90
+++ b/src/utils/read_geometry.f90
@@ -43,13 +43,16 @@ subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
 ! Compute nuclear repulsion energy
 
   ENuc = 0
+  open(unit=3,file='int/ENuc.dat')
+  read(3,*) ENuc
+  close(unit=3)
 
-  do i=1,nNuc-1
-    do j=i+1,nNuc
-      RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2
-      ENuc = ENuc + ZNuc(i)*ZNuc(j)/(AntoBo*sqrt(RAB))
-    end do
-  end do
+  ! do i=1,nNuc-1
+  !   do j=i+1,nNuc
+  !     RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2
+  !     ENuc = ENuc + ZNuc(i)*ZNuc(j)/(AntoBo*sqrt(RAB))
+  !   end do
+  ! end do
 
 ! Close file with geometry specification
   close(unit=10)