PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2025-04-30 04:04:50 +02:00
Code Releases Activity

TODO FD occupancies and chem pot

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-01-29 23:03:54 +01:00
parent e105665f8e
commit d0bdf0661e
4 changed files with 44 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
src/HF/HFB.f90
View File

@ -1,5 +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)
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)
! Perform Hartree-Fock Bogoliubov calculation
@ -22,6 +23,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
double precision,intent(in) :: ZNuc(nNuc)
double precision,intent(in) :: rNuc(nNuc,ncart)
double precision,intent(in) :: ENuc
double precision,intent(in) :: temperature
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
@ -51,6 +53,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
double precision :: thrs_N
double precision,external :: trace_matrix
double precision,allocatable :: eHFB_(:)
double precision,allocatable :: Occ(:)
double precision,allocatable :: err(:,:)
double precision,allocatable :: err_diis(:,:)
double precision,allocatable :: F_diis(:,:)
@ -102,7 +105,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
allocate(err_diis(nBas_Sq,max_diis))
allocate(F_diis(nBas_Sq,max_diis))
! Guess coefficients chem. pot
! Guess chem. pot.
chem_pot = 0.5d0*(eHF(nO)+eHF(nO+1))
@ -115,7 +118,26 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
rcond = 0d0
P(:,:) = matmul(c(:,1:nO), transpose(c(:,1:nO)))
Panom(:,:) = 0d0 ! Do sth TODO
Panom(:,:) = 0d0
! Use Fermi-Dirac occupancies to compute P, Panom, and chem_pot
if(abs(temperature)>1e-4) then ! TODO
allocate(Occ(nOrb))
Occ(:) = 0d0
Occ(1:nO) = 1d0
! call fermi_dirac_occ(nOrb,chem_pot,Occ,eHF)
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
deallocate(Occ)
endif
P_old(:,:) = P(:,:)
Panom_old(:,:) = Panom(:,:)
P(:,:) = 2d0*P(:,:)
@ -166,9 +188,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
! Adjust the chemical potential
if( abs(trace_1rdm-nO) > thrs_N ) then
call fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB_)
endif
! Extract P and Panom from R
@ -289,15 +309,3 @@ end subroutine
! if(level_shift > 0d0 .and. Conv > thresh) then
! call level_shifting(level_shift,nBas,nOrb,nO,S,c,F)
! endif
!block
!integer::iorb1
!write(*,*) 'Tr[1D] and chem_pot',trace_1rdm,chem_pot
!write(*,*) 'iter',nSCF
!do iorb1=1,nOrb2
! write(*,*) eHFB_(iorb1)
!enddo
!do iorb1=1,nOrb
! write(*,'(7f10.5)') P(iorb1,:)
!enddo
!end block

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)
guess_type,mix,temperature)
! Restricted branch of QuAcK
@ -19,6 +19,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
integer,intent(in) :: nV
integer,intent(in) :: nR
double precision,intent(in) :: ENuc
double precision,intent(in) :: temperature
double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,ncart)
@ -93,7 +94,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)
FHF,Delta,temperature)
call wall_time(end_HF)
t_HF = end_HF - start_HF

7
src/QuAcK/QuAcK.f90
View File

@ -73,6 +73,8 @@ program QuAcK
logical :: dotest,doRtest,doUtest,doGtest
double precision :: temperature
character(len=256) :: working_dir
! Check if the right number of arguments is provided
@ -134,7 +136,7 @@ program QuAcK
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,temperature)
!------------------!
! Hardware !
@ -290,7 +292,8 @@ 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)
S,T,V,Hc,X,dipole_int_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix, &
temperature)
!-----------!
! Stop Test !

11
src/QuAcK/read_options.f90
View File

@ -7,7 +7,7 @@ subroutine read_options(working_dir,
maxSCF_GW,thresh_GW,max_diis_GW,lin_GW,eta_GW,reg_GW,TDA_W, &
maxSCF_GT,thresh_GT,max_diis_GT,lin_GT,eta_GT,reg_GT,TDA_T, &
doACFDT,exchange_kernel,doXBS, &
dophBSE,dophBSE2,doppBSE,dBSE,dTDA)
dophBSE,dophBSE2,doppBSE,dBSE,dTDA,temperature)
! Read desired methods
@ -72,6 +72,8 @@ subroutine read_options(working_dir,
logical,intent(out) :: dBSE
logical,intent(out) :: dTDA
double precision,intent(out) :: temperature
! Local variables
character(len=1) :: ans1,ans2,ans3,ans4,ans5
@ -217,6 +219,13 @@ subroutine read_options(working_dir,
if(ans4 == 'T') dBSE = .true.
if(ans5 == 'F') dTDA = .false.
! Options for dynamical Fermi-Dirac occupancies
temperature = 0d0
read(1,*)
read(1,*) temperature
endif
! Close file with options