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

TODO print output info add DIIS

Browse Source
This commit is contained in:
Mauricio Rodriguez-Mayorga 2025-01-30 16:05:29 +01:00
parent b69451ba38
commit 2744c97f2e
2 changed files with 36 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
src/HF/fermi_dirac_occ.f90
View File

@ -12,12 +12,14 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
! Local variables
logical :: backward
integer :: iorb
integer :: isteps
double precision :: delta_chem_pot
double precision :: chem_pot_change
double precision :: grad_electrons
double precision :: trace_1rdm
double precision :: trace_old
! Output variables
@ -27,10 +29,12 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
! Initialize variables
isteps=0
delta_chem_pot = 1.0d-3
trace_1rdm = -1.0d0
chem_pot_change = 0.0d0
backward = .false.
isteps = 0
delta_chem_pot = 1.0d-1
chem_pot_change = 0d0
grad_electrons = 1d0
trace_1rdm = -1d0
write(*,*)
write(*,*)' Fermi-Dirac distribution for the occ numbers'
@ -40,13 +44,34 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
'|','Tr[1D]','|','Chem. Pot.','|'
write(*,*)'-------------------------------------'
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_1rdm=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
! First approach close the value with an error lower than 1
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_old=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_old,'|',chem_pot,'|'
do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + delta_chem_pot
Occ(:) = fermi_dirac(eHF,chem_pot,temperature)
trace_1rdm=sum(Occ(:))
write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
'|',trace_1rdm,'|',chem_pot,'|'
if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then
backward=.true.
chem_pot = chem_pot - 2d0*delta_chem_pot
delta_chem_pot=-delta_chem_pot
endif
enddo
! Do final search
write(*,*)'-------------------------------------'
isteps=0
delta_chem_pot = 1.0d-1
do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 )
isteps = isteps + 1
chem_pot = chem_pot + chem_pot_change

7
src/HF/fix_chem_pot.f90
View File

@ -17,13 +17,13 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB
double precision :: delta_chem_pot
double precision :: chem_pot_change
double precision :: grad_electrons
double precision :: trace_up
double precision :: trace_down
double precision :: trace_old
! Output variables
double precision :: trace_1rdm
double precision :: trace_up
double precision :: trace_down
double precision :: trace_old
double precision,intent(inout):: chem_pot
double precision,intent(inout):: cp(nOrb2,nOrb2)
double precision,intent(inout):: R(nOrb2,nOrb2)
@ -72,6 +72,7 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB
! Do final search
write(*,*)'------------------------------------------------------'
isteps = 0
delta_chem_pot = 1.0d-3
do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 )