Works for HFB + Delta

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 
 

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

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              :: chem_pot_change
-  double precision              :: a
+  double precision              :: grad_electrons
-  double precision              :: c
 
 ! 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.
+  backward = .false.
-  isteps=0
+  isteps = 0
-  golden_ratio = 1.618033988
+  delta_chem_pot  = 1.0d-1
-  delta_chem_pot = 5d-1
+  chem_pot_change = 0d0
-  trace_down = 0d0
+  grad_electrons  = 1d0
-  trace_up   = 0d0
+  trace_1rdm      = -1d0
-  chem_pot_down = 0d0 
-  chem_pot_up   = 0d0 
 
   ! 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(*,*)'------------------------------------------------------'
-  write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') &
+  write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1A15,2X,A1)') &
-          '|','Tr[1D]','|','Chem. Pot.','|'
+          '|','Tr[1D]','|','Chem. Pot.','|','Grad N','|'
-  write(*,*)'-------------------------------------'
+  write(*,*)'------------------------------------------------------'
 
-  ! Set interval to search (the good chem_pot is in [chem_pot_down,chem_pot_up])
+  ! First approach close the value with an error lower than 1
-  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) then
-   ! Increase chem_pot to occupy more orbs.
-   is_up=.true.
-   do
-    isteps = isteps+1
-    chem_pot = chem_pot + delta_chem_pot
-    call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_up,H_hfb,cp,R,eHFB_)
-    write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
-    '|',trace_up,'|',chem_pot,'|'
-    if(trace_up>=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
 
-  ! Use Golden-section search algorithm to find chem_pot
+  call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_old,H_hfb,cp,R,eHFB_)
-  isteps=0
+  write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
-  do
+  '|',trace_old,'|',chem_pot,'|',grad_electrons,'|'
-    isteps = isteps+1
+  do while( abs(trace_1rdm-nO) > 1.0d0 .and. isteps <= 100 )
-    a=(chem_pot_up-chem_pot_down)/(1d0+golden_ratio)
+   isteps = isteps + 1
-    c=a/golden_ratio
+   chem_pot = chem_pot + delta_chem_pot
-    chem_pot_curr=chem_pot_down+a
+   call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_)
-    chem_pot=chem_pot_curr+c
+   write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1F16.10,1X,A1)') &
-    call diag_H_hfb(nOrb,nOrb2,chem_pot_curr,trace_curr,H_hfb,cp,R,eHFB_)
+   '|',trace_1rdm,'|',chem_pot,'|',grad_electrons,'|'
-    write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
+   if( abs(trace_1rdm-nO) > abs(trace_old-nO) .and. .not.backward ) then
-     '|',trace_curr,'|',chem_pot_curr,'|'
+    backward=.true.
-    call diag_H_hfb(nOrb,nOrb2,chem_pot,trace_1rdm,H_hfb,cp,R,eHFB_)
+    chem_pot = chem_pot - 2d0*delta_chem_pot
-    write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
+    delta_chem_pot=-delta_chem_pot
-     '|',trace_1rdm,'|',chem_pot,'|'
+   endif
-    if(abs(trace_1rdm-nO)<thrs_N .or. isteps>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
   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