Incl FD distrib occ

src/HF/HFB.f90

@@ -126,7 +126,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
    allocate(Occ(nOrb))
    Occ(:)     = 0d0
    Occ(1:nO)  = 1d0
-!   call fermi_dirac_occ(nOrb,chem_pot,Occ,eHF)
+   call fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
    P(:,:)      = 0d0
    Panom(:,:)  = 0d0
    do iorb=1,nOrb
@@ -148,8 +148,12 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
 !------------------------------------------------------------------------
 ! Main SCF loop
 !------------------------------------------------------------------------
-  
-  do while(Conv > thresh .and. nSCF < maxSCF)
+
+  write(*,*)
+  write(*,*) 'Enterning HFB SCF procedure'  
+  write(*,*)
+  do while(Conv > thresh .and. nSCF < 2)
+  !do while(Conv > thresh .and. nSCF < maxSCF)
 
     ! Increment 
 
@@ -217,7 +221,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
 
     ! Anomalous energy
 
-    EL = -0.25d0*trace_matrix(nBas,matmul(Panom,Delta))
+    EL = 0.25d0*trace_matrix(nBas,matmul(Panom,Delta))
 
     ! Total energy
 

src/HF/fermi_dirac_occ.f90

@@ -0,0 +1,85 @@
+subroutine fermi_dirac_occ(nO,nOrb,thrs_N,temperature,chem_pot,Occ,eHF)
+
+! Use Fermi Dirac distribution to set up fractional Occs numbers and adjust the chemical potential to integrate to the N for 2N electron systems
+
+  implicit none
+
+! Input variables
+
+  integer,intent(in)            :: nO,nOrb
+  double precision,intent(in)   :: thrs_N
+  double precision,intent(in)   :: temperature
+
+! Local variables
+
+  integer                       :: iorb
+  integer                       :: isteps
+  double precision              :: delta_chem_pot
+  double precision              :: chem_pot_change
+  double precision              :: grad_electrons
+  double precision              :: trace_1rdm
+
+! Output variables
+
+  double precision,intent(inout):: chem_pot
+  double precision,intent(inout):: eHF(nOrb) 
+  double precision,intent(inout):: Occ(nOrb) 
+
+  !  Initialize variables
+
+  isteps=0
+  delta_chem_pot = 1.0d-3
+  trace_1rdm     = -1.0d0
+  chem_pot_change = 0.0d0
+
+  write(*,*)
+  write(*,*)' Fermi-Dirac distribution for the occ numbers'
+  write(*,*)
+  write(*,*)'-------------------------------------'
+  write(*,'(1X,A1,1X,A15,1X,A1,1X,A15,1X,A1)') &
+          '|','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,'|'
+
+
+  do while( abs(trace_1rdm-nO) > 1.0d-8 .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
+    write(*,'(1X,A1,F16.10,1X,A1,F16.10,1X,A1)') &
+    '|',trace_1rdm,'|',chem_pot,'|'
+  enddo
+  write(*,*)'-------------------------------------'
+  write(*,*)
+
+  write(*,*)
+  write(*,*) ' Initial occ. numbers '
+  write(*,*)
+  do iorb=1,nOrb
+   write(*,'(3X,F16.10)') Occ(iorb)
+  enddo
+
+contains 
+
+function fermi_dirac(eHF,chem_pot,temperature)
+  implicit none
+  double precision,intent(in) :: eHF(nOrb)
+  double precision,intent(in) :: chem_pot
+  double precision,intent(in) :: temperature
+  double precision            :: fermi_dirac(nOrb)
+
+  fermi_dirac(:) = 1d0 / ( 1d0 + exp((eHF(:) - chem_pot ) / temperature ) )
+
+end function fermi_dirac
+
+end subroutine 
+