LCPQ
/
quack
mirror of https://github.com/pfloos/quack
4
1
Fork 0
Code Releases Activity

fix UHF

This commit is contained in:
Pierre-Francois Loos 2023-10-24 14:39:02 +02:00
parent 83e5c9b07c
commit 95a3715a3f
8 changed files with 45 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -13,7 +13,7 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F
# G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW
T F F F F F
F F F F F F
# G0T0pp* evGTpp* qsGTpp* G0T0eh evGTeh qsGTeh
F F F F F F
# * unrestricted version available

27
src/HF/HF.f90
View File

@ -39,34 +39,15 @@ subroutine HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF,thresh,max_di
double precision :: start_HF ,end_HF ,t_HF
integer :: nSCF
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: dipole(ncart)
double precision :: Conv
double precision :: Gap
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: Fp(:,:)
! Output variables
logical,intent(out) :: unrestricted
double precision,intent(out) :: EHF
double precision,intent(out) :: epsHF(nBas)
double precision,intent(out) :: cHF(nBas,nBas)
double precision,intent(out) :: PHF(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
double precision,intent(out) :: epsHF(nBas,nspin)
double precision,intent(out) :: cHF(nBas,nBas,nspin)
double precision,intent(out) :: PHF(nBas,nBas,nspin)
double precision,intent(out) :: F(nBas,nBas,nspin)
!------------------------------------------------------------------------
! Compute RHF energy

2
src/HF/ROHF.f90
View File

@ -239,6 +239,6 @@ subroutine ROHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc
! Compute final UHF energy
call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_ROHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
call print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
end subroutine

58
src/HF/UHF.f90
View File

@ -112,35 +112,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
nSCF = nSCF + 1
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),e(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
! Mix guess for UHF solution in singlet states
if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c)
! Compute density matrix
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do
! Build Coulomb repulsion
do ispin=1,nspin
@ -189,6 +160,35 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
end if
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:,:) = Fp(:,:,:)
do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),e(:,ispin))
end do
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
end do
! Mix guess for UHF solution in singlet states
if(mix .and. nSCF == 1) call mix_guess(nBas,nO,c)
! Compute density matrix
do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------

4
src/HF/mo_guess.f90
View File

@ -12,10 +12,6 @@ subroutine mo_guess(nBas,guess_type,S,Hc,X,c)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
! Local variables
integer :: nSCF
! Output variables
double precision,intent(out) :: c(nBas,nBas)

4
src/HF/print_ROHF.f90
View File

@ -1,4 +1,4 @@
subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
subroutine print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
! Print one- and two-electron energies and other stuff for RoHF calculation
@ -7,7 +7,6 @@ subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: Ov(nBas,nBas)
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: ENuc
@ -23,7 +22,6 @@ subroutine print_ROHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
double precision :: S_exact,S2_exact
double precision :: S,S2
! HOMO and LUMO

14
src/QuAcK/QuAcK.f90
View File

@ -42,8 +42,8 @@ program QuAcK
double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: dipole_int_aa(:,:,:)
double precision,allocatable :: dipole_int_bb(:,:,:)
double precision,allocatable :: F_AO(:,:)
double precision,allocatable :: F_MO(:,:)
double precision,allocatable :: F_AO(:,:,:)
double precision,allocatable :: F_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: ixyz
@ -185,7 +185,7 @@ program QuAcK
allocate(cHF(nBas,nBas,nspin),epsHF(nBas,nspin),PHF(nBas,nBas,nspin),S(nBas,nBas),T(nBas,nBas), &
V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas),dipole_int_AO(nBas,nBas,ncart), &
dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas))
dipole_int_MO(nBas,nBas,ncart),F_AO(nBas,nBas,nspin))
! Read integrals
@ -215,8 +215,8 @@ program QuAcK
call wall_time(start_HF)
call HF(doRHF,doUHF,doROHF,doRMOM,doUMOM,unrestricted,maxSCF_HF,thresh_HF,max_diis_HF, &
guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO, &
ERI_AO,dipole_int_AO,X,EHF,epsHF,cHF,PHF)
guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F_AO,ERI_AO, &
dipole_int_AO,X,EHF,epsHF,cHF,PHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -291,7 +291,7 @@ program QuAcK
! Memory allocation
allocate(ERI_MO(nBas,nBas,nBas,nBas),F_MO(nBas,nBas))
allocate(ERI_MO(nBas,nBas,nBas,nBas),F_MO(nBas,nBas,nspin))
! Read and transform dipole-related integrals
@ -304,7 +304,7 @@ program QuAcK
call AOtoMO_integral_transform(1,1,1,1,nBas,cHF,ERI_AO,ERI_MO)
F_MO(:,:) = F_AO(:,:)
F_MO(:,:,1) = F_AO(:,:,1)
call AOtoMO_transform(nBas,cHF,F_MO)
end if

4
src/make_ninja.py
View File

@ -57,7 +57,7 @@ else:
compile_gfortran_mac = """
FC = gfortran
AR = libtool -static -o
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
CC = gcc
CXX = g++
LAPACK=-lblas -llapack
@ -68,7 +68,7 @@ FIX_ORDER_OF_LIBS=
compile_gfortran_linux = """
FC = gfortran
AR = ar crs
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
FFLAGS = -I$IDIR -J$IDIR -fbacktrace -Wall -g -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant
CC = gcc
CXX = g++
LAPACK=-lblas -llapack