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

Allow restart in HFB

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-02-06 14:52:12 +01:00
parent 2b8c9c8bd2
commit de42425683
7 changed files with 195 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
src/HF/HFB.f90
View File

@ -1,6 +1,6 @@
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,Panom,F,Delta, &
temperature,sigma,chem_pot_hf)
temperature,sigma,chem_pot_hf,restart_hfb)
! Perform Hartree-Fock Bogoliubov calculation
@ -35,6 +35,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
! Local variables
logical :: chem_pot_hf
logical :: restart_hfb
integer :: nBas2
integer :: iorb
integer :: nSCF
@ -97,6 +98,8 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
! Memory allocation
allocate(Occ(nOrb))
allocate(J(nBas,nBas))
allocate(K(nBas,nBas))
@ -132,7 +135,6 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot
if(abs(temperature)>1d-4) then
allocate(Occ(nOrb))
Occ(:) = 0d0
Occ(1:nO) = 1d0
call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
@ -145,7 +147,20 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo
deallocate(Occ)
endif
! Read restart file
if(restart_hfb) then
call read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
P(:,:) = 0d0
Panom(:,:) = 0d0
do iorb=1,nOrb
P(:,:) = P(:,:) + Occ(iorb) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb))) * &
matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb)))
enddo
endif
P(:,:) = 2d0*P(:,:)
@ -318,7 +333,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis)
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ)
deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao)
stop
@ -326,14 +341,17 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
end if
! Compute dipole moments, occupation numbers and || Anomalous density||
! also print the restart file
eigVEC(:,:) = 0d0
eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb)
call diagonalize_matrix(nOrb2,eigVEC,eigVAL)
eigVAL(:) = 2d0*eigVAL(:)
Occ(1:nOrb) = eigVAL(nOrb+1:nOrb2)
c(1:nBas,1:nOrb) = matmul(X(1:nBas,1:nOrb),eigVEC(1:nOrb,nOrb+1:nOrb2))
norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2)))
call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
call print_HFB(nBas,nOrb,nO,norm_anom,eigVAL,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole)
call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot)
call print_HFB(nBas,nOrb,nO,norm_anom,Occ,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole)
! Testing zone
@ -348,7 +366,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
! Memory deallocation
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis)
deallocate(J,K,eigVEC,H_HFB,R,eigVAL,err_diis,H_HFB_diis,Occ)
deallocate(err_ao,S_ao,X_ao,R_ao_old,H_HFB_ao)
end subroutine

12
src/HF/print_HFB.f90
View File

@ -3,7 +3,7 @@
subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole)
! Print one-electron energies and other stuff for G0W0
! Print one-electron energies and other stuff
implicit none
include 'parameters.h'
@ -12,7 +12,7 @@ subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF
integer,intent(in) :: nBas, nOrb
integer,intent(in) :: nO
double precision,intent(in) :: Occ(2*nOrb)
double precision,intent(in) :: Occ(nOrb)
double precision,intent(in) :: ENuc
double precision,intent(in) :: ET
double precision,intent(in) :: EV
@ -66,11 +66,11 @@ subroutine print_HFB(nBas, nOrb, nO, N_anom, Occ, ENuc, ET, EV, EJ, EK, EL, ERHF
write(*,'(A50)') ' HFB occupation numbers '
write(*,'(A50)') '---------------------------------------'
trace_occ=0d0
do iorb=1,2*nOrb
if(abs(Occ(2*nOrb-iorb))>1d-8) then
write(*,'(I7,10F15.8)') iorb,Occ(2*nOrb-iorb)
do iorb=1,nOrb
if(abs(Occ(nOrb+1-iorb))>1d-8) then
write(*,'(I7,10F15.8)') iorb,2d0*Occ(nOrb+1-iorb)
endif
trace_occ=trace_occ+Occ(iorb)
trace_occ=trace_occ+2d0*Occ(iorb)
enddo
write(*,*)
write(*,'(A33,1X,F16.10,A3)') ' Trace [ 1D ] = ',trace_occ,' '

116
src/HF/read_restart_HFB.f90 Normal file
View File

@ -0,0 +1,116 @@
! ---
subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
! Write the binary file used to restart calcs
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: S(nBas,nBas)
! Local variables
integer :: ibas,iorb,iorb1,iunit=667
integer :: nBas_
integer :: nOrb_
double precision :: chem_pot_
double precision :: max_diff
double precision,allocatable :: eigVAL(:)
double precision,allocatable :: c_tmp(:,:)
double precision,allocatable :: S_mol(:,:)
double precision,allocatable :: X_mol(:,:)
! Output variables
double precision,intent(out) :: chem_pot
double precision,intent(out) :: Occ(nOrb)
double precision,intent(out) :: c(nBas,nOrb)
! Dump results
allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb))
open(unit=iunit,form='unformatted',file='hfb_bin',status='old')
read(iunit) nBas_,nOrb_
read(iunit) chem_pot_
do iorb=1,nOrb
read(iunit) c_tmp(:,iorb)
enddo
do iorb=1,nOrb
read(iunit) eigVAL(iorb)
enddo
close(iunit)
if(nBas==nBas_ .and. nOrb==nOrb_) then
write(*,*)
write(*,*)' Reading restart file'
write(*,*)
chem_pot = chem_pot_
Occ(:) = eigVAL(:)
c(:,:) = c_tmp(:,:)
! Check the orthonormality
max_diff=-1d0
S_mol = matmul(transpose(c),matmul(S,c))
do iorb=1,nOrb
do iorb1=1,iorb
if(iorb==iorb1) then
if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0)
else
if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1))
endif
enddo
enddo
if(max_diff>1d-8) then
write(*,*) ' '
write(*,'(a)') ' Orthonormality violations. Applying Lowdin orthonormalization.'
write(*,*) ' '
X_mol = S_mol
S_mol = 0d0
call diagonalize_matrix(nOrb,X_mol,eigVAL)
do iorb=1,nOrb
S_mol(iorb,iorb) = 1d0/(sqrt(eigVAL(iorb)) + 1d-10)
enddo
X_mol = matmul(X_mol,matmul(S_mol,transpose(X_mol)))
c = matmul(c,X_mol)
max_diff=-1d0
S_mol = matmul(transpose(c),matmul(S,c))
do iorb=1,nOrb
do iorb1=1,iorb
if(iorb==iorb1) then
if(abs(S_mol(iorb,iorb)-1d0)>1d-8) max_diff=abs(S_mol(iorb,iorb)-1d0)
else
if(abs(S_mol(iorb,iorb1))>1d-8) max_diff=abs(S_mol(iorb,iorb1))
endif
enddo
enddo
if(max_diff<=1d-8) then
write(*,*) ' '
write(*,'(a)') ' No orthonormality violations.'
write(*,*) ' '
else
write(*,*) ' '
write(*,'(a)') ' Error in Lowdin orthonormalization.'
write(*,*) ' '
stop
endif
else
write(*,*) ' '
write(*,'(a)') ' No orthonormality violations.'
write(*,*) ' '
endif
endif
deallocate(eigVAL,c_tmp,S_mol,X_mol)
end subroutine

37
src/HF/write_restart_HFB.f90 Normal file
View File

@ -0,0 +1,37 @@
! ---
subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
! Write the binary file used to restart calcs
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nOrb
double precision,intent(in) :: chem_pot
double precision,intent(in) :: Occ(nOrb)
double precision,intent(in) :: c(nBas,nOrb)
! Local variables
integer :: ibas,iorb,iunit=666
! Dump results
open(unit=iunit,form='unformatted',file='hfb_bin')
write(iunit) nBas,nOrb
write(iunit) chem_pot
do iorb=1,nOrb
write(iunit) c(:,nOrb+1-iorb)
enddo
do iorb=1,nOrb
write(iunit) Occ(nOrb+1-iorb)
enddo
write(iunit) iunit
close(iunit)
end subroutine

5
src/QuAcK/BQuAcK.f90
View File

@ -1,6 +1,6 @@
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,temperature,sigma,chem_pot_hf)
guess_type,mix,temperature,sigma,chem_pot_hf,restart_hfb)
! Restricted branch of QuAcK
@ -13,6 +13,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
logical,intent(in) :: doHFB
logical,intent(in) :: restart_hfb
logical,intent(in) :: chem_pot_hf
integer,intent(in) :: nNuc,nBas,nOrb
integer,intent(in) :: nC
@ -95,7 +96,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
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,PanomHF, &
FHF,Delta,temperature,sigma,chem_pot_hf)
FHF,Delta,temperature,sigma,chem_pot_hf,restart_hfb)
call wall_time(end_HF)
t_HF = end_HF - start_HF

5
src/QuAcK/QuAcK.f90
View File

@ -74,6 +74,7 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest
logical :: chem_pot_hf
logical :: restart_hfb
double precision :: temperature,sigma
character(len=256) :: working_dir
@ -138,7 +139,7 @@ program QuAcK
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA, &
temperature,sigma,chem_pot_hf)
temperature,sigma,chem_pot_hf,restart_hfb)
!------------------!
! Hardware !
@ -295,7 +296,7 @@ program QuAcK
if(doBQuAcK) &
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, &
temperature,sigma,chem_pot_hf)
temperature,sigma,chem_pot_hf,restart_hfb)
!-----------!
! Stop Test !

7
src/QuAcK/read_options.f90
View File

@ -8,7 +8,7 @@ subroutine read_options(working_dir,
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA, &
temperature,sigma,chem_pot_hf)
temperature,sigma,chem_pot_hf,restart_hfb)
! Read desired methods
@ -74,6 +74,7 @@ subroutine read_options(working_dir,
logical,intent(out) :: dTDA
logical,intent(out) :: chem_pot_hf
logical,intent(out) :: restart_hfb
double precision,intent(out) :: temperature
double precision,intent(out) :: sigma
@ -227,11 +228,13 @@ subroutine read_options(working_dir,
temperature = 0d0
sigma = 1d0
chem_pot_hf = .false.
restart_hfb = .false.
read(1,*)
read(1,*) temperature,sigma,ans1
read(1,*) temperature,sigma,ans1,ans2
if(ans1 == 'T') chem_pot_hf = .true.
if(ans2 == 'T') restart_hfb = .true.
endif