fixing nuclear energy

This commit is contained in:
Antoine Marie 2023-12-14 10:14:00 +01:00
parent 21b1bbda5c
commit 889ad27c20
2 changed files with 16 additions and 6 deletions

View File

@ -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

View File

@ -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)