LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

fix bug in PyDuck

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-12-14 11:16:44 +01:00
parent 04887bc569
commit 1c7a98ba0e
5 changed files with 29 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PyDuck.py
View File

@ -79,7 +79,7 @@ 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(str(mol.energy_nuc()))
f.write(' ')
f.close()

8
src/GT/ufG0T0pp.f90
View File

@ -469,8 +469,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
write(*,*)'-------------------------------------------'
do s=1,nH
! if(eGT(s) < eF .and. eGT(s) > eF - window) then
if(Z(s) > cutoff1) then
if(eGT(s) < eF .and. eGT(s) > eF - window) then
! if(Z(s) > cutoff1) then
write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') &
'|',s,'|',eGT(s)*HaToeV,'|',Z(s),'|'
end if
@ -545,8 +545,8 @@ subroutine ufG0T0pp(dotest,TDA_T,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
do s=1,nH
! if(eGT(s) < eF .and. eGT(s) > eF - window) then
if(Z(s) > cutoff2) then
if(eGT(s) < eF .and. eGT(s) > eF - window) then
! if(Z(s) > cutoff2) then
write(*,*)'-------------------------------------------------------------'
write(*,'(1X,A7,1X,I3,A6,I3,A1,1X,A7,F12.6,A13,F6.4,1X)') &

2
src/HF/UHF.f90
View File

@ -248,7 +248,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Compute final UHF energy
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_UHF(nBas,nO,S,eHF,c,ENuc,ET,EV,EJ,EK,EUHF,dipole)
call print_UHF(nBas,nO,S,eHF,c,P,ENuc,ET,EV,EJ,EK,EUHF,dipole)
! Print test values

2
src/HF/print_RHF.f90
View File

@ -27,7 +27,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
double precision :: Gap
double precision :: S,S2
logical :: dump_orb = .false.
logical :: dump_orb = .true.
! HOMO and LUMO

27
src/HF/print_UHF.f90
View File

@ -1,4 +1,4 @@
subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
subroutine print_UHF(nBas,nO,S,eHF,c,P,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
! Print one- and two-electron energies and other stuff for UHF calculation
@ -9,9 +9,10 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: Ov(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: c(nBas,nBas,nspin)
double precision,intent(in) :: P(nBas,nBas,nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET(nspin)
double precision,intent(in) :: EV(nspin)
@ -29,8 +30,10 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
double precision :: Gap(nspin)
double precision :: Sz
double precision :: Sx2,Sy2,Sz2
integer :: mu,nu
double precision,allocatable :: qa(:),qb(:)
logical :: dump_orb = .false.
logical :: dump_orb = .true.
! HOMO and LUMO
@ -51,8 +54,8 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
end do
Sz = 0.5d0*dble(nO(1) - nO(2))
Sx2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sy2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(Ov,c(:,1:nO(2),2)))**2)
Sx2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(S,c(:,1:nO(2),2)))**2)
Sy2 = 0.25d0*dble(nO(1) - nO(2)) + 0.5d0*nO(2) - 0.5d0*sum(matmul(transpose(c(:,1:nO(1),1)),matmul(S,c(:,1:nO(2),2)))**2)
Sz2 = 0.25d0*dble(nO(1) - nO(2))**2
! Dump results
@ -129,4 +132,18 @@ subroutine print_UHF(nBas,nO,Ov,eHF,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole)
call vecout(nBas,eHF(:,2))
write(*,*)
allocate(qa(nBas),qb(nBas))
qa(:) = 0d0
qb(:) = 0d0
do mu=1,nBas
do nu=1,nBas
qa(mu) = qa(mu) + P(mu,nu,1)*S(nu,mu)
qb(mu) = qb(mu) + P(mu,nu,2)*S(nu,mu)
end do
end do
call vecout(nBas,qa)
call vecout(nBas,qb)
end subroutine