PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

fixed confl

Browse Source
This commit is contained in:
Abdallah Ammar 2023-11-17 10:55:59 +01:00
commit 6d485d364d
10 changed files with 200 additions and 250 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
README.md
View File

@ -6,6 +6,7 @@
- [Pierre-Francois Loos](https://pfloos.github.io/WEB_LOOS) - [Pierre-Francois Loos](https://pfloos.github.io/WEB_LOOS)
- [Enzo Monino](https://enzomonino.github.io) - [Enzo Monino](https://enzomonino.github.io)
- [Antoine Marie](https://antoine-marie.github.io) - [Antoine Marie](https://antoine-marie.github.io)
- [Abdallah Ammar](https://scholar.google.com/citations?user=y437T5sAAAAJ&hl=en)
- [Anthony Scemama](https://scemama.github.io) - [Anthony Scemama](https://scemama.github.io)
# What is it? # What is it?

88
src/HF/GHF.f90
View File

@ -182,52 +182,7 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
if(nSCF > 1) Conv = maxval(abs(err)) if(nSCF > 1) Conv = maxval(abs(err))
! DIIS extrapolation ! Kinetic energy
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,F)
end if
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) then
call level_shifting(level_shift,nBas,nO,S,C,F)
end if
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X),matmul(F,X))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
Cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas2,Cp,eHF)
! Back-transform eigenvectors in non-orthogonal basis
C(:,:) = matmul(X,Cp)
! Mix guess for UHF solution in singlet states
! if(nSCF == 1) call mix_guess(nBas,nO,mix,c)
! Form super density matrix
P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO)))
! Compute individual density matrices
Paa(:,:) = P( 1:nBas , 1:nBas )
Pab(:,:) = P( 1:nBas ,nBas+1:nBas2)
Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------
! Kinetic energy
ETaa = trace_matrix(nBas,matmul(Paa,T)) ETaa = trace_matrix(nBas,matmul(Paa,T))
ETbb = trace_matrix(nBas,matmul(Pbb,T)) ETbb = trace_matrix(nBas,matmul(Pbb,T))
@ -263,9 +218,48 @@ subroutine GHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
EGHF = ET + EV + EJ + EK EGHF = ET + EV + EJ + EK
! DIIS extrapolation
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis,F_diis,err,F)
end if
! Level-shifting
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,C,F)
! Transform Fock matrix in orthogonal basis
Fp(:,:) = matmul(transpose(X),matmul(F,X))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
Cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas2,Cp,eHF)
! Back-transform eigenvectors in non-orthogonal basis
C(:,:) = matmul(X,Cp)
! Mix guess for UHF solution in singlet states
! if(nSCF == 1) call mix_guess(nBas,nO,mix,c)
! Form super density matrix
P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO)))
! Compute individual density matrices
Paa(:,:) = P( 1:nBas , 1:nBas )
Pab(:,:) = P( 1:nBas ,nBas+1:nBas2)
Pba(:,:) = P(nBas+1:nBas2, 1:nBas )
Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
! Dump results ! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',EGHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|' '|',nSCF,'|',EGHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
end do end do

2
src/HF/GHF_search.f90
View File

@ -211,6 +211,8 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
call matrix_exponential(nBas2,R,ExpR) call matrix_exponential(nBas2,R,ExpR)
c = matmul(c,ExpR) c = matmul(c,ExpR)
! call matout(nBas2,nBas2,matmul(transpose(ExpR),ExpR))
else else
write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!' write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!'

72
src/HF/RHF.f90
View File

@ -42,7 +42,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
double precision :: dipole(ncart) double precision :: dipole(ncart)
double precision :: Conv double precision :: Conv
double precision :: Gap
double precision :: rcond double precision :: rcond
double precision,external :: trace_matrix double precision,external :: trace_matrix
double precision,allocatable :: error(:,:) double precision,allocatable :: error(:,:)
@ -82,26 +81,25 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
call mo_guess(nBas,guess_type,S,Hc,X,c) call mo_guess(nBas,guess_type,S,Hc,X,c)
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Initialization ! Initialization
F_diis(:,:) = 0d0 F_diis(:,:) = 0d0
error_diis(:,:) = 0d0 error_diis(:,:) = 0d0
Conv = 1d0
n_diis = 0
nSCF = 0
rcond = 0d0 rcond = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Main SCF loop ! Main SCF loop
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
write(*,*)'| RHF calculation |' write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
write(*,*)'----------------------------------------------------' '|','#','|','E(RHF)','|','EJ(RHF)','|','EK(RHF)','|','Conv','|'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & write(*,*)'-----------------------------------------------------------------------------'
'|','#','|','HF energy','|','Conv','|','HL Gap','|'
write(*,*)'----------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF) do while(Conv > thresh .and. nSCF < maxSCF)
@ -121,6 +119,26 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error)) Conv = maxval(abs(error))
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
! Hartree energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
! Exchange energy
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
! Total energy
ERHF = ET + EV + EJ + EK
! DIIS extrapolation ! DIIS extrapolation
if(max_diis > 1) then if(max_diis > 1) then
@ -144,32 +162,14 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
! Density matrix ! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! Compute HF energy
ERHF = trace_matrix(nBas,matmul(P,Hc)) &
+ 0.5d0*trace_matrix(nBas,matmul(P,J)) &
+ 0.25d0*trace_matrix(nBas,matmul(P,K))
! Compute HOMO-LUMO gap
if(nBas > nO) then
Gap = eHF(nO+1) - eHF(nO)
else
Gap = 0d0
end if
! Dump results ! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',ERHF+ENuc,'|',Conv,'|',Gap,'|' '|',nSCF,'|',ERHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
end do end do
write(*,*)'----------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! End of SCF loop ! End of SCF loop
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -188,14 +188,6 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN
end if end if
! Compute HF energy
ET = trace_matrix(nBas,matmul(P,T))
EV = trace_matrix(nBas,matmul(P,V))
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
EK = 0.25d0*trace_matrix(nBas,matmul(P,K))
ERHF = ET + EV + EJ + EK
! Compute dipole moments ! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)

99
src/HF/ROHF.f90
View File

@ -1,5 +1,5 @@
subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EHF,e,c,Ptot) nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot)
! Perform restricted open-shell Hartree-Fock calculation ! Perform restricted open-shell Hartree-Fock calculation
@ -42,7 +42,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
double precision :: ET(nspin) double precision :: ET(nspin)
double precision :: EV(nspin) double precision :: EV(nspin)
double precision :: EJ(nsp) double precision :: EJ(nsp)
double precision :: Ex(nspin) double precision :: EK(nspin)
double precision :: dipole(ncart) double precision :: dipole(ncart)
double precision,allocatable :: cp(:,:) double precision,allocatable :: cp(:,:)
@ -61,9 +61,9 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Output variables ! Output variables
double precision,intent(out) :: EHF double precision,intent(out) :: EROHF
double precision,intent(out) :: e(nBas) double precision,intent(out) :: eHF(nBas)
double precision,intent(out) :: c(nBas,nBas) double precision,intent(inout):: c(nBas,nBas)
double precision,intent(out) :: Ptot(nBas,nBas) double precision,intent(out) :: Ptot(nBas,nBas)
! Hello world ! Hello world
@ -80,8 +80,8 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Memory allocation ! Memory allocation
allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), & allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), &
P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), & P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), &
err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis))
! Guess coefficients and demsity matrices ! Guess coefficients and demsity matrices
@ -90,6 +90,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
do ispin=1,nspin do ispin=1,nspin
P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin)))) P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin))))
end do end do
Ptot(:,:) = P(:,:,1) + P(:,:,2)
! Initialization ! Initialization
@ -106,10 +107,10 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(ROHF)','|','Ex(ROHF)','|','Conv','|' '|','#','|','E(ROHF)','|','EJ(ROHF)','|','EK(ROHF)','|','Conv','|'
write(*,*)'----------------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF) do while(Conv > thresh .and. nSCF < maxSCF)
@ -139,15 +140,43 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Check convergence ! Check convergence
err(:,:) = matmul(Ftot(:,:),matmul(Ptot(:,:),S(:,:))) - matmul(matmul(S(:,:),Ptot(:,:)),Ftot(:,:)) err(:,:) = matmul(Ftot,matmul(Ptot,S)) - matmul(matmul(S,Ptot),Ftot)
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(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
end do
! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
end do
! Total energy
EROHF = sum(ET) + sum(EV) + sum(EJ) + sum(EK)
! DIIS extrapolation ! DIIS extrapolation
if(max_diis > 1) then if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis) n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,Ftot) call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,Ftot)
end if end if
@ -168,7 +197,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Diagonalize Fock matrix to get eigenvectors and eigenvalues ! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:) cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e) call diagonalize_matrix(nBas,cp,eHF)
! Back-transform eigenvectors in non-orthogonal basis ! Back-transform eigenvectors in non-orthogonal basis
@ -181,45 +210,13 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
end do end do
Ptot(:,:) = P(:,:,1) + P(:,:,2) Ptot(:,:) = P(:,:,1) + P(:,:,2)
!------------------------------------------------------------------------
! Compute ROHF energy
!------------------------------------------------------------------------
! Kinetic energy
do ispin=1,nspin
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(:,:)))
end do
! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
end do
! Total energy
EHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
! Dump results ! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',EHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|' '|',nSCF,'|',EROHF + ENuc,'|',sum(EJ),'|',sum(EK),'|',Conv,'|'
end do end do
write(*,*)'----------------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! End of SCF loop ! End of SCF loop
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -241,13 +238,13 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN
! Compute final UHF energy ! Compute final UHF energy
call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_ROHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole) call print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole)
! Print test values ! Print test values
if(dotest) then if(dotest) then
call dump_test_value('R','ROHF energy',EHF) call dump_test_value('R','ROHF energy',EROHF)
end if end if

88
src/HF/UHF.f90
View File

@ -42,7 +42,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
double precision :: ET(nspin) double precision :: ET(nspin)
double precision :: EV(nspin) double precision :: EV(nspin)
double precision :: EJ(nsp) double precision :: EJ(nsp)
double precision :: Ex(nspin) double precision :: EK(nspin)
double precision :: dipole(ncart) double precision :: dipole(ncart)
double precision,allocatable :: cp(:,:,:) double precision,allocatable :: cp(:,:,:)
@ -91,22 +91,22 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Initialization ! Initialization
nSCF = 0
Conv = 1d0
n_diis = 0 n_diis = 0
F_diis(:,:,:) = 0d0 F_diis(:,:,:) = 0d0
err_diis(:,:,:) = 0d0 err_diis(:,:,:) = 0d0
nSCF = 0
Conv = 1d0
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! Main SCF loop ! Main SCF loop
!------------------------------------------------------------------------ !------------------------------------------------------------------------
write(*,*) write(*,*)
write(*,*)'----------------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') & write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','E(UHF)','|','Ex(UHF)','|','Conv','|' '|','#','|','E(UHF)','|','EJ(UHF)','|','EK(UHF)','|','Conv','|'
write(*,*)'----------------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF) do while(Conv > thresh .and. nSCF < maxSCF)
@ -139,7 +139,35 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
end do 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(:,:)))
end do
! Potential energy
do ispin=1,nspin
EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:)))
end do
! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
end do
! Total energy
EUHF = sum(ET) + sum(EV) + sum(EJ) + sum(EK)
! DIIS extrapolation ! DIIS extrapolation
if(max_diis > 1) then if(max_diis > 1) then
@ -174,7 +202,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
do ispin=1,nspin do ispin=1,nspin
call diagonalize_matrix(nBas,cp(:,:,ispin),eHF(:,ispin)) call diagonalize_matrix(nBas,cp(:,:,ispin),eHF(:,ispin))
end do end do
! Back-transform eigenvectors in non-orthogonal basis ! Back-transform eigenvectors in non-orthogonal basis
do ispin=1,nspin do ispin=1,nspin
@ -191,45 +219,13 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin)))
end do end do
!------------------------------------------------------------------------
! Compute UHF energy
!------------------------------------------------------------------------
! Kinetic energy
do ispin=1,nspin
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(:,:)))
end do
! Hartree energy
EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1)))
EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2)))
EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2)))
! Exchange energy
do ispin=1,nspin
Ex(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin)))
end do
! Total energy
EUHF = sum(ET(:)) + sum(EV(:)) + sum(EJ(:)) + sum(Ex(:))
! Dump results ! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X)') & write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') &
'|',nSCF,'|',EUHF + ENuc,'|',sum(Ex(:)),'|',Conv,'|' '|',nSCF,'|',EUHF + ENuc,'|',sum(EJ),'|',sum(EK),'|',Conv,'|'
end do end do
write(*,*)'----------------------------------------------------------' write(*,*)'-----------------------------------------------------------------------------'
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! End of SCF loop ! End of SCF loop
!------------------------------------------------------------------------ !------------------------------------------------------------------------
@ -251,7 +247,7 @@ subroutine UHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu
! Compute final UHF energy ! Compute final UHF energy
call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole) 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,Ex,EUHF,dipole) call print_UHF(nBas,nO,S,eHF,c,ENuc,ET,EV,EJ,EK,EUHF,dipole)
! Print test values ! Print test values

66
src/HF/print_GHF.f90
View File

@ -33,11 +33,10 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
integer :: HOMO integer :: HOMO
integer :: LUMO integer :: LUMO
double precision :: Gap double precision :: Gap
double precision :: Sx ,Sy ,Sz double precision :: Sx,Sy,Sz
double precision :: Sx2,Sy2,Sz2 double precision :: SmSp,SpSm,Sz2,S2
double precision :: SmSp,SpSm,S2 ! double precision :: na, nb
double precision :: na, nb ! double precision :: nonco_z, contam_uhf, xy_perp, contam_ghf
double precision :: nonco_z, contam_uhf, xy_perp, contam_ghf
double precision,allocatable :: Ca(:,:) double precision,allocatable :: Ca(:,:)
double precision,allocatable :: Cb(:,:) double precision,allocatable :: Cb(:,:)
@ -80,7 +79,6 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
Pba = matmul(transpose(Cb),matmul(S,Ca)) Pba = matmul(transpose(Cb),matmul(S,Ca))
Pbb = matmul(transpose(Cb),matmul(S,Cb)) Pbb = matmul(transpose(Cb),matmul(S,Cb))
! Compute components of S = (Sx,Sy,Sz) ! Compute components of S = (Sx,Sy,Sz)
Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba)) Sx = 0.5d0*(trace_matrix(nO,Pab) + trace_matrix(nO,Pba))
@ -89,17 +87,13 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
! Compute <S^2> = <Sx^2> + <Sy^2> + <Sz^2> ! Compute <S^2> = <Sx^2> + <Sy^2> + <Sz^2>
SpSm = 0d0 SpSm = 0d0
do i=1,nO do i=1,nO
do j=1,nO do j=1,nO
SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i) SpSm = SpSm + Pab(i,i)*Pba(j,j) - Pab(i,j)*Pba(j,i)
end do end do
end do end do
SpSm = 0.5d0*(trace_matrix(nO,Paa) + SpSm) SpSm = trace_matrix(nO,Paa) + SpSm
! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
SmSp = 0d0 SmSp = 0d0
do i=1,nO do i=1,nO
@ -107,10 +101,7 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i) SmSp = SmSp + Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i)
end do end do
end do end do
SmSp = 0.5d0*(trace_matrix(nO,Pbb) + SmSp) SmSp = trace_matrix(nO,Pbb) + SmSp
! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab))
Sz2 = 0d0 Sz2 = 0d0
do i=1,nO do i=1,nO
@ -119,47 +110,10 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
end do end do
end do end do
Sz2 = 0.25d0*(dble(nO) + Sz2) Sz2 = 0.25d0*(dble(nO) + Sz2)
print*,'<Sz^2> = ',Sz2
Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) &
- 0.25d0*trace_matrix(nO,Paa-Pbb)**2 &
- 0.25d0*trace_matrix(nO,matmul(Paa-Pbb,Paa-Pbb))
S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb)
do i=1,nO
do j=1,nO
S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 &
+ (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i))
end do
end do
print*,'<S^2> = ',S2
Sx2 = 0.5d0*(S2 - Sz2 + 0.5d0*(SmSp + SpSm))
print*,'<Sx^2> = ',Sx2
Sy2 = 0.5d0*(S2 - Sz2 - 0.5d0*(SmSp + SpSm))
print*,'<Sy^2> = ',Sy2
! Sx2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Pab+Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) + matmul(Pab,Pab))
! Sy2 = 0.25d0*trace_matrix(nO,Paa+Pbb) - 0.25d0*trace_matrix(nO,Pab-Pba)**2 &
! - 0.5d0*trace_matrix(nO,matmul(Paa,Pbb) - matmul(Pab,Pab))
! Sz2 = 0.25d0*trace_matrix(nO,Paa+Pbb) + 0.25d0*trace_matrix(nO,Paa-Pbb)**2 &
! - 0.25d0*trace_matrix(nO,matmul(Paa,Paa) + matmul(Pbb,Pbb)) &
! + 0.25d0*trace_matrix(nO,matmul(Pab,Pba) + matmul(Pba,Pab))
! S2 = Sz*(Sz+1d0) + trace_matrix(nO,Pbb) + 0.25d0*trace_matrix(nO,Paa+Pbb)
! do i=1,nO
! do j=1,nO
! S2 = S2 - 0.25d0*(Paa(i,j) - Pbb(i,j))**2 &
! + (Pba(i,i)*Pab(j,j) - Pba(i,j)*Pab(j,i))
! end do
! end do
! print*,'<S^2> = ',S2
! Compute <S^2> from Sz^2, S^+S^- and S^-S^+
S2 = Sz2 + 0.5d0*(SpSm + SmSp)
@ -229,10 +183,8 @@ subroutine print_GHF(nBas,nBas2,nO,eHF,C,P,S,ENuc,ET,EV,EJ,EK,EGHF,dipole)
write(*,'(A33,1X,F16.6)') ' <Sy> = ',Sy write(*,'(A33,1X,F16.6)') ' <Sy> = ',Sy
write(*,'(A33,1X,F16.6)') ' <Sz> = ',Sz write(*,'(A33,1X,F16.6)') ' <Sz> = ',Sz
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.6)') ' <Sx**2> = ',Sx2 write(*,'(A33,1X,F16.6)') ' <Sx^2+Sy^2> = ',S2 - Sz2
write(*,'(A33,1X,F16.6)') ' <Sy**2> = ',Sy2
write(*,'(A33,1X,F16.6)') ' <Sz**2> = ',Sz2 write(*,'(A33,1X,F16.6)') ' <Sz**2> = ',Sz2
write(*,'(A33,1X,F16.6)') ' <S**2> = ',Sx2+Sy2+Sz2
write(*,'(A33,1X,F16.6)') ' <S**2> = ',S2 write(*,'(A33,1X,F16.6)') ' <S**2> = ',S2
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A36)') ' Dipole moment (Debye) ' write(*,'(A36)') ' Dipole moment (Debye) '

14
src/HF/print_ROHF.f90
View File

@ -24,6 +24,8 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
double precision :: Gap(nspin) double precision :: Gap(nspin)
double precision :: S,S2 double precision :: S,S2
logical :: dump_orb = .false.
! HOMO and LUMO ! HOMO and LUMO
do ispin=1,nspin do ispin=1,nspin
@ -96,11 +98,13 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole)
! Print results ! Print results
write(*,'(A50)') '-----------------------------------------' if(dump_orb) then
write(*,'(A50)') 'ROHF orbital coefficients ' write(*,'(A50)') '-----------------------------------------'
write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') 'ROHF orbital coefficients '
call matout(nBas,nBas,c) write(*,'(A50)') '-----------------------------------------'
write(*,*) call matout(nBas,nBas,c)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' ROHF orbital energies (au) ' write(*,'(A50)') ' ROHF orbital energies (au) '
write(*,'(A50)') '---------------------------------------' write(*,'(A50)') '---------------------------------------'

9
src/utils/DIIS_extrapolation.f90
View File

@ -9,11 +9,15 @@ subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout)
! Input variables ! Input variables
integer,intent(in) :: n_err,n_e integer,intent(in) :: n_err,n_e
double precision,intent(in) :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis) double precision,intent(in) :: error_in(n_err)
double precision,intent(in) :: error(n_err,n_diis)
double precision,intent(in) :: e(n_e,n_diis)
! Local variables ! Local variables
double precision,allocatable :: A(:,:),b(:),w(:) double precision,allocatable :: A(:,:)
double precision,allocatable :: b(:)
double precision,allocatable :: w(:)
! Output variables ! Output variables
@ -45,7 +49,6 @@ subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout)
! Solve linear system ! Solve linear system
! call easy_linear_solve(n_diis+1,A,b,w)
call linear_solve(n_diis+1,A,b,w,rcond) call linear_solve(n_diis+1,A,b,w,rcond)
! Extrapolate ! Extrapolate

11
src/utils/wrap_lapack.f90
View File

@ -175,7 +175,16 @@ subroutine linear_solve(N,A,b,x,rcond)
integer,allocatable :: ipiv(:),iwork(:) integer,allocatable :: ipiv(:),iwork(:)
double precision,allocatable :: AF(:,:),work(:) double precision,allocatable :: AF(:,:),work(:)
lwork = 3*N ! Find optimal size for temporary arrays
allocate(work(1))
lwork = -1
call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)
lwork = int(work(1))
deallocate(work)
allocate(AF(N,N),ipiv(N),work(lwork),iwork(N)) allocate(AF(N,N),ipiv(N),work(lwork),iwork(N))
call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info) call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)