LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-03-09 18:22:25 +01:00
Code Releases Activity

Starting HFB scf

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-01-27 15:24:43 +01:00
parent 12044b800d
commit cf1ad337e5
4 changed files with 396 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

229
src/HF/HFB.f90 Normal file
View File

@ -0,0 +1,229 @@
subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,F)
! Perform Hartree-Fock Bogoliubov calculation
implicit none
include 'parameters.h'
! Input variables
logical,intent(in) :: dotest
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
integer,intent(in) :: nO
integer,intent(in) :: nNuc
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int(nBas,nBas,ncart)
! Local variables
integer :: nSCF
integer :: nBas_Sq
integer :: n_diis
double precision :: ET
double precision :: EV
double precision :: EJ
double precision :: EK
double precision :: dipole(ncart)
double precision :: Conv
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: err(:,:)
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: Fp(:,:)
! Output variables
double precision,intent(out) :: EHFB
double precision,intent(out) :: eHF(nOrb)
double precision,intent(inout):: c(nBas,nOrb)
double precision,intent(out) :: P(nBas,nBas)
double precision,intent(out) :: F(nBas,nBas)
! Hello world
write(*,*)
write(*,*)'*****************************'
write(*,*)'* HF Bogoliubov Calculation *'
write(*,*)'*****************************'
write(*,*)
! Useful quantities
nBas_Sq = nBas*nBas
! Memory allocation
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
allocate(err(nBas,nBas))
allocate(cp(nOrb,nOrb))
allocate(Fp(nOrb,nOrb))
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Guess coefficients and density matrix
P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO)))
! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
! c(1,1), nBas, c(1,1), nBas, &
! 0.d0, P(1,1), nBas)
! Initialization
n_diis = 0
F_diis(:,:) = 0d0
err_diis(:,:) = 0d0
rcond = 0d0
Conv = 1d0
nSCF = 0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'-----------------------------------------------------------------------------'
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(HFB)','|','EJ(HFB)','|','EK(HFB)','|','Conv','|'
write(*,*)'-----------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Build Fock matrix
call Hartree_matrix_AO_basis(nBas,P,ERI,J)
call exchange_matrix_AO_basis(nBas,P,ERI,K)
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
! Check convergence
err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
if(nSCF > 1) Conv = maxval(abs(err))
! 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
EHFB = ET + EV + EJ + EK
! DIIS extrapolation
if(max_diis > 1) then
n_diis = min(n_diis+1,max_diis)
call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F)
end if
! Level shift
if(level_shift > 0d0 .and. Conv > thresh) then
call level_shifting(level_shift,nBas,nOrb,nO,S,c,F)
endif
! Diagonalize Fock matrix
Fp = matmul(transpose(X),matmul(F,X))
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nOrb,cp,eHF)
c = matmul(X,cp)
! Density matrix
P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
! call dgemm('N', 'T', nBas, nBas, nO, 2.d0, &
! c(1,1), nBas, c(1,1), nBas, &
! 0.d0, P(1,1), nBas)
! 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,E10.2,1X,A1,1X)') &
'|',nSCF,'|',EHFB + ENuc,'|',EJ,'|',EK,'|',Conv,'|'
end do
write(*,*)'-----------------------------------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
stop
end if
! Compute dipole moments
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_HFB(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EHFB,dipole)
! Testing zone
if(dotest) then
call dump_test_value('R','HFB energy',EHFB)
call dump_test_value('R','HFB HOMO energy',eHF(nO))
call dump_test_value('R','HFB LUMO energy',eHF(nO+1))
call dump_test_value('R','HFB dipole moment',norm2(dipole))
end if
! Memory deallocation
deallocate(J,K,err,cp,Fp,err_diis,F_diis)
end subroutine

90
src/HF/print_HFB.f90 Normal file
View File

@ -0,0 +1,90 @@
! ---
subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole)
! Print one-electron energies and other stuff for G0W0
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: eHF(nOrb)
double precision,intent(in) :: cHF(nBas,nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
double precision,intent(in) :: EV
double precision,intent(in) :: EJ
double precision,intent(in) :: EK
double precision,intent(in) :: ERHF
double precision,intent(in) :: dipole(ncart)
! Local variables
integer :: ixyz
integer :: HOMO
integer :: LUMO
double precision :: Gap
double precision :: S,S2
logical :: dump_orb = .false.
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = eHF(LUMO)-eHF(HOMO)
S2 = 0d0
S = 0d0
! Dump results
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A33)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' One-electron energy = ',ET + EV,' au'
write(*,'(A33,1X,F16.10,A3)') ' Kinetic energy = ',ET,' au'
write(*,'(A33,1X,F16.10,A3)') ' Potential energy = ',EV,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Two-electron energy = ',EJ + EK,' au'
write(*,'(A33,1X,F16.10,A3)') ' Hartree energy = ',EJ,' au'
write(*,'(A33,1X,F16.10,A3)') ' Exchange energy = ',EK,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.10,A3)') ' Electronic energy = ',ERHF,' au'
write(*,'(A33,1X,F16.10,A3)') ' Nuclear repulsion = ',ENuc,' au'
write(*,'(A33,1X,F16.10,A3)') ' HFB energy = ',ERHF + ENuc,' au'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.6,A3)') ' HFB HOMO energy = ',eHF(HOMO)*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HFB LUMO energy = ',eHF(LUMO)*HaToeV,' eV'
write(*,'(A33,1X,F16.6,A3)') ' HFB HOMO-LUMO gap = ',Gap*HaToeV,' eV'
write(*,'(A50)') '---------------------------------------'
write(*,'(A33,1X,F16.6)') ' <Sz> = ',S
write(*,'(A33,1X,F16.6)') ' <S^2> = ',S2
write(*,'(A50)') '---------------------------------------'
write(*,'(A36)') ' Dipole moment (Debye) '
write(*,'(10X,4A10)') 'X','Y','Z','Tot.'
write(*,'(10X,4F10.4)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD
write(*,'(A50)') '---------------------------------------'
write(*,*)
! Print results
if(dump_orb) then
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' HFB orbital coefficients '
write(*,'(A50)') '---------------------------------------'
call matout(nBas, nOrb, cHF)
write(*,*)
end if
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' HFB orbital energies (au) '
write(*,'(A50)') '---------------------------------------'
call vecout(nOrb, eHF)
write(*,*)
end subroutine

80
src/QuAcK/BQuAcK.f90
View File

@ -1,4 +1,6 @@
subroutine BQuAcK(working_dir,dotest,doHFB)
subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, &
guess_type,mix)
! Restricted branch of QuAcK
@ -11,10 +13,40 @@ subroutine BQuAcK(working_dir,dotest,doHFB)
logical,intent(in) :: doHFB
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
integer,intent(in) :: nO
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nOrb)
double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart)
integer,intent(in) :: maxSCF_HF,max_diis_HF
double precision,intent(in) :: thresh_HF,level_shift,mix
integer,intent(in) :: guess_type
! Local variables
double precision :: start_HF ,end_HF ,t_HF
double precision :: start_int, end_int, t_int
double precision,allocatable :: eHF(:)
double precision,allocatable :: cHF(:,:)
double precision,allocatable :: PHF(:,:)
double precision,allocatable :: FHF(:,:)
double precision :: ERHF,EHFB
! double precision,allocatable :: dipole_int_MO(:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
! double precision,allocatable :: ERI_MO(:,:,:,:)
write(*,*)
write(*,*) '******************************'
write(*,*) '* Bogoliubov Branch of QuAcK *'
@ -25,14 +57,42 @@ subroutine BQuAcK(working_dir,dotest,doHFB)
! Memory allocation !
!-------------------!
!---------------------!
! Hartree-Fock module !
!---------------------!
allocate(eHF(nOrb))
allocate(cHF(nBas,nOrb))
allocate(PHF(nBas,nBas))
allocate(FHF(nBas,nBas))
! allocate(dipole_int_MO(nOrb,nOrb,ncart))
! allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
allocate(ERI_AO(nBas,nBas,nBas,nBas))
call wall_time(start_int)
call read_2e_integrals(working_dir,nBas,ERI_AO)
call wall_time(end_int)
t_int = end_int - start_int
write(*,*)
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for reading 2e-integrals =',t_int,' seconds'
write(*,*)
!--------------------------------!
! Hartree-Fock Bogoliubov module !
!--------------------------------!
if(doHFB) then
! Run first a RHF calculation
call wall_time(start_HF)
! call HFB(dotest)
call RHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,FHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds'
write(*,*)
! Continue with a HFB calculation
call wall_time(start_HF)
call HFB(dotest,maxSCF_HF,thresh_HF,max_diis_HF,level_shift,nNuc,ZNuc,rNuc,ENuc, &
nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,FHF)
call wall_time(end_HF)
t_HF = end_HF - start_HF
@ -41,4 +101,14 @@ subroutine BQuAcK(working_dir,dotest,doHFB)
end if
! Memory deallocation
deallocate(eHF)
deallocate(cHF)
deallocate(PHF)
deallocate(FHF)
! deallocate(dipole_int_MO)
! deallocate(ERI_MO)
deallocate(ERI_AO)
end subroutine

3
src/QuAcK/QuAcK.f90
View File

@ -289,7 +289,8 @@ program QuAcK
! Bogoliubov QuAcK branch !
!--------------------------!
if(doBQuAcK) &
call BQuAcK(working_dir,doGtest,doHFB)
call BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, &
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix)
!-----------!
! Stop Test !