diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90 index c35bda5..b6b2264 100644 --- a/src/HF/HFB.f90 +++ b/src/HF/HFB.f90 @@ -111,7 +111,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, ! Initialization - thrs_N = 1d-6 + thrs_N = 1d-8 n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 @@ -122,7 +122,7 @@ 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)>1e-4) then ! TODO + if(abs(temperature)>1e-4) then allocate(Occ(nOrb)) Occ(:) = 0d0 Occ(1:nO) = 1d0 @@ -152,8 +152,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, write(*,*) write(*,*) 'Enterning HFB SCF procedure' write(*,*) - do while(Conv > thresh .and. nSCF < 2) - !do while(Conv > thresh .and. nSCF < maxSCF) + do while(Conv > thresh .and. nSCF < maxSCF) ! Increment diff --git a/src/HF/fermi_dirac_occ.f90 b/src/HF/fermi_dirac_occ.f90 index e34ccd6..e48ce87 100644 --- a/src/HF/fermi_dirac_occ.f90 +++ b/src/HF/fermi_dirac_occ.f90 @@ -47,14 +47,14 @@ subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF) '|',trace_1rdm,'|',chem_pot,'|' - do while( abs(trace_1rdm-nO) > 1.0d-8 .and. isteps <= 100 ) + do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 ) isteps = isteps + 1 chem_pot = chem_pot + chem_pot_change Occ(:) = fermi_dirac(eHF,chem_pot,temperature) trace_1rdm = sum(Occ(:)) grad_electrons = ( sum(fermi_dirac(eHF,chem_pot+delta_chem_pot,temperature)) & - sum(fermi_dirac(eHF,chem_pot-delta_chem_pot,temperature)) )/(2.0d0*delta_chem_pot) - chem_pot_change = -(trace_1rdm-nO)/grad_electrons + chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10) write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & '|',trace_1rdm,'|',chem_pot,'|' enddo diff --git a/src/HF/fix_chem_pot.f90 b/src/HF/fix_chem_pot.f90 index 0f4f755..3847e94 100644 --- a/src/HF/fix_chem_pot.f90 +++ b/src/HF/fix_chem_pot.f90 @@ -11,35 +11,33 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB ! Local variables - logical :: is_up + logical :: backward integer :: iorb integer :: isteps - double precision :: trace_curr,trace_down,trace_up - double precision :: chem_pot_curr,chem_pot_down,chem_pot_up double precision :: delta_chem_pot - double precision :: golden_ratio - double precision :: a - double precision :: c + double precision :: chem_pot_change + double precision :: grad_electrons ! 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) double precision,intent(inout):: eHFB_(nOrb2) double precision,intent(inout):: H_hfb(nOrb2,nOrb2) - ! Initialize delta_chem_pot + ! Initialize - is_up=.false. - isteps=0 - golden_ratio = 1.618033988 - delta_chem_pot = 5d-1 - trace_down = 0d0 - trace_up = 0d0 - chem_pot_down = 0d0 - chem_pot_up = 0d0 + backward = .false. + isteps = 0 + delta_chem_pot = 1.0d-1 + chem_pot_change = 0d0 + grad_electrons = 1d0 + trace_1rdm = -1d0 ! Set H_HFB to its non-chemical potential dependent contribution @@ -49,93 +47,47 @@ subroutine fix_chem_pot(nO,nOrb,nOrb2,thrs_N,trace_1rdm,chem_pot,H_hfb,cp,R,eHFB enddo write(*,*) - write(*,*)'-------------------------------------' - write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') & - '|','Tr[1D]','|','Chem. Pot.','|' - write(*,*)'-------------------------------------' + write(*,*)'------------------------------------------------------' + write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1A15,2X,A1)') & + '|','Tr[1D]','|','Chem. Pot.','|','Grad N','|' + write(*,*)'------------------------------------------------------' - ! Set interval to search (the good chem_pot is in [chem_pot_down,chem_pot_up]) - call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_curr,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot,'|' - chem_pot_curr=chem_pot - if(trace_curr=nO+1 .or. isteps>100) exit ! max 50 au steps for mu (is a lot) - enddo - chem_pot_up=chem_pot - else - ! Decrease chem_pot to occupy less orbs. - do - isteps = isteps+1 - chem_pot = chem_pot - delta_chem_pot - call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_down,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_down,'|',chem_pot,'|' - if(trace_down<=nO-1 .or. isteps>100) exit ! max 50 au steps for mu (is a lot) - enddo - chem_pot_down=chem_pot - endif - if(is_up) then - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot_curr,'|' - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_up ,'|',chem_pot_up ,'|' - trace_down=trace_curr - chem_pot_down=chem_pot_curr - else - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot_curr,'|' - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_down,'|',chem_pot_down,'|' - trace_up=trace_curr - chem_pot_up=chem_pot_curr - endif + ! First approach close the value with an error lower than 1 - ! Use Golden-section search algorithm to find chem_pot - isteps=0 - do - isteps = isteps+1 - a=(chem_pot_up-chem_pot_down)/(1d0+golden_ratio) - c=a/golden_ratio - chem_pot_curr=chem_pot_down+a - chem_pot=chem_pot_curr+c - call diag_H_hfb(nOrb,nOrb2,chem_pot_curr,trace_curr,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_curr,'|',chem_pot_curr,'|' - call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) - write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') & - '|',trace_1rdm,'|',chem_pot,'|' - if(abs(trace_1rdm-nO)1000) exit ! 1000 steps for finding mu - if(is_up) then - if(trace_1rdm>=trace_curr .or. abs(trace_1rdm-trace_curr)<1d-8 ) then - chem_pot_down=chem_pot_curr - trace_down=trace_curr - else - chem_pot_up=chem_pot - trace_up=trace_1rdm - endif - else - if(trace_1rdm>=trace_curr .or. abs(trace_1rdm-trace_curr)<1d-8 ) then - chem_pot_up=chem_pot - trace_up=trace_1rdm - else - chem_pot_down=chem_pot_curr - trace_down=trace_curr - endif - endif + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_old,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & + '|',trace_old,'|',chem_pot,'|',grad_electrons,'|' + do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 ) + isteps = isteps + 1 + chem_pot = chem_pot + delta_chem_pot + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|' + 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 - write(*,*)'-------------------------------------' + ! Do final search + + isteps = 0 + delta_chem_pot = 1.0d-3 + do while( abs(trace_1rdm-nO) > thrs_N .and. isteps <= 100 ) + isteps = isteps + 1 + chem_pot = chem_pot + chem_pot_change + call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot+delta_chem_pot,trace_up,H_hfb,cp,R,eHFB_) + call diag_H_hfb(nOrb,nOrb2,chem_pot-delta_chem_pot,trace_down,H_hfb,cp,R,eHFB_) + grad_electrons = (trace_up-trace_down)/(2.0d0*delta_chem_pot) + chem_pot_change = -(trace_1rdm-nO)/(grad_electrons+1d-10) + write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') & + '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|' + enddo + write(*,*)'------------------------------------------------------' write(*,*) - + ! Reset H_HFB to its chemical potential version do iorb=1,nOrb