From 1c7a98ba0e37a80a69c78bccb368459690a22aa1 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 14 Dec 2023 11:16:44 +0100 Subject: [PATCH] fix bug in PyDuck --- PyDuck.py | 2 +- src/GT/ufG0T0pp.f90 | 8 ++++---- src/HF/UHF.f90 | 2 +- src/HF/print_RHF.f90 | 2 +- src/HF/print_UHF.f90 | 27 ++++++++++++++++++++++----- 5 files changed, 29 insertions(+), 12 deletions(-) diff --git a/PyDuck.py b/PyDuck.py index 1f77a4b..c74dc22 100644 --- a/PyDuck.py +++ b/PyDuck.py @@ -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() diff --git a/src/GT/ufG0T0pp.f90 b/src/GT/ufG0T0pp.f90 index 68a1db6..df1c63f 100644 --- a/src/GT/ufG0T0pp.f90 +++ b/src/GT/ufG0T0pp.f90 @@ -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)') & diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 61f2991..919cb0f 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -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 diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 3a06d31..89b93ec 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -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 diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 466667f..d985b76 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -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