PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-23 04:43:42 +01:00
Code Releases Activity

small modifs in HF

Browse Source
This commit is contained in:
Pierre-Francois Loos 2023-11-20 14:57:45 +01:00
parent 9388a0d290
commit 3c1553fd57
3 changed files with 18 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
src/HF/RHF.f90
View File

@ -44,8 +44,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision :: Conv
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: error(:,:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: err(:,:)
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
@ -74,8 +74,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Memory allocation
allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
Fp(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
allocate(J(nBas,nBas),K(nBas,nBas),err(nBas,nBas),cp(nBas,nBas),F(nBas,nBas), &
Fp(nBas,nBas),err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and density matrix
@ -84,9 +84,10 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Initialization
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
rcond = 0d0
n_diis = 0
F_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
Conv = 1d0
nSCF = 0
@ -116,8 +117,8 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Check convergence
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error))
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
if(nSCF > 1) Conv = maxval(abs(err))
! Kinetic energy
@ -144,7 +145,7 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F)
end if

11
src/HF/UHF.f90
View File

@ -94,6 +94,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
n_diis = 0
F_diis(:,:,:) = 0d0
err_diis(:,:,:) = 0d0
rcond(:) = 0d0
nSCF = 0
Conv = 1d0
@ -138,18 +139,18 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
err(:,:,ispin) = matmul(F(:,:,ispin),matmul(P(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),P(:,:,ispin)),F(:,:,ispin))
end do
if(nSCF > 1) Conv = maxval(abs(err(:,:,:)))
if(nSCF > 1) Conv = maxval(abs(err))
! Kinetic energy
do ispin=1,nspin
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:)))
ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V))
end do
! Hartree energy
@ -193,7 +194,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Transform Fock matrix in orthogonal basis
do ispin=1,nspin
Fp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(F(:,:,ispin),X(:,:)))
Fp(:,:,ispin) = matmul(transpose(X),matmul(F(:,:,ispin),X))
end do
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
@ -206,7 +207,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin
c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin))
c(:,:,ispin) = matmul(X,cp(:,:,ispin))
end do
! Mix guess for UHF solution in singlet states

2
src/HF/print_RHF.f90
View File

@ -60,7 +60,7 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole)
write(*,'(A33,1X,F16.6,A3)') ' HF LUMO energy = ',eHF(LUMO)*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HF HOMO-LUMO gap = ',Gap*HaToeV,' eV'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.6)') ' S = ',2d0*S+ 1d0
write(*,'(A33,1X,F16.6)') ' S = ',S
write(*,'(A33,1X,F16.6)') ' <S**2> = ',S2
write(*,'(A50)') '---------------------------------------'
write(*,'(A36)') ' Dipole moment (Debye) '