Starting chem pot

src/HF/HFB.f90

@@ -1,5 +1,5 @@
 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)
+               nBas,nOrb,nO,S,T,V,Hc,ERI,dipole_int,X,EHFB,eHF,c,P,Panom,F,Delta)
 
 ! Perform Hartree-Fock Bogoliubov calculation
 
@@ -40,6 +40,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
   double precision              :: EJ
   double precision              :: EK
   double precision              :: EL
+  double precision              :: chem_pot
   double precision              :: dipole(ncart)
 
   double precision              :: Conv
@@ -62,6 +63,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
   double precision,intent(out)  :: P(nBas,nBas)
   double precision,intent(out)  :: Panom(nBas,nBas)
   double precision,intent(out)  :: F(nBas,nBas)
+  double precision,intent(out)  :: Delta(nBas,nBas)
 
 ! Hello world
 
@@ -89,10 +91,11 @@ 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 and density matrix
+! Guess coefficients, density matrix, and chem. pot
 
+  chem_pot = 0d0
   P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO)))
-  Panom(:,:) = -P(:,:)
+  Panom(:,:) = -P(:,:) ! Do sth TODO
 
 ! Initialization
 
@@ -127,6 +130,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
     call anomalous_matrix_AO_basis(nBas,Panom,ERI,L)
     
     F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:)
+    Delta(:,:) = L(:,:)
 
     ! Check convergence 
 
@@ -151,11 +155,11 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
 
     ! Anomalous energy
 
-    EL = 0.25d0*trace_matrix(nBas,matmul(-Panom,L))
+    EL = -0.25d0*trace_matrix(nBas,matmul(Panom,L))
 
     ! Total energy
 
-    EHFB = ET + EV + EJ + EK
+    EHFB = ET + EV + EJ + EK + EL
 
     ! DIIS extrapolation
 
@@ -182,6 +186,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
     ! Density matrix
 
     P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO)))
+    Panom(:,:) = -P(:,:) ! Do sth TODO
 
     ! Dump results
 
@@ -213,7 +218,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc, &
 ! Compute dipole moments
 
   call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
-  call print_HFB(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EL,EHFB,dipole)
+  call print_HFB(nBas,nOrb,nO,eHF,c,ENuc,ET,EV,EJ,EK,EL,EHFB,chem_pot,dipole)
 
 ! Testing zone
 

src/HF/print_HFB.f90

@@ -1,7 +1,7 @@
 
 ! ---
 
-subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, dipole)
+subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, chem_pot, dipole)
 
 ! Print one-electron energies and other stuff for G0W0
 
@@ -21,6 +21,7 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, d
   double precision,intent(in)        :: EK
   double precision,intent(in)        :: EL
   double precision,intent(in)        :: ERHF
+  double precision,intent(in)        :: chem_pot
   double precision,intent(in)        :: dipole(ncart)
 
 ! Local variables
@@ -61,6 +62,8 @@ subroutine print_HFB(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, EL, ERHF, d
   write(*,'(A33,1X,F16.6,A3)')  ' HFB LUMO     energy = ',eHF(LUMO)*HaToeV,' eV'
   write(*,'(A33,1X,F16.6,A3)')  ' HFB HOMO-LUMO gap   = ',Gap*HaToeV,' eV'
   write(*,'(A50)')           '---------------------------------------'
+  write(*,'(A33,1X,F16.10,A3)') ' Chemical potential  = ',chem_pot,' au'
+  write(*,'(A50)')           '---------------------------------------'
   write(*,'(A36)')           ' Dipole moment (Debye)    '
   write(*,'(10X,4A10)')      'X','Y','Z','Tot.'
   write(*,'(10X,4F10.4)')    (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD

src/QuAcK/BQuAcK.f90

@@ -43,10 +43,9 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
   double precision,allocatable  :: PHF(:,:)
   double precision,allocatable  :: PanomHF(:,:)
   double precision,allocatable  :: FHF(:,:)
+  double precision,allocatable  :: Delta(:,:)
   double precision              :: ERHF,EHFB
-!  double precision,allocatable  :: dipole_int_MO(:,:,:)
   double precision,allocatable  :: ERI_AO(:,:,:,:)
-!  double precision,allocatable  :: ERI_MO(:,:,:,:)
 
   write(*,*)
   write(*,*) '******************************'
@@ -63,8 +62,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
   allocate(PHF(nBas,nBas))
   allocate(PanomHF(nBas,nBas))
   allocate(FHF(nBas,nBas))
-!  allocate(dipole_int_MO(nOrb,nOrb,ncart))
+  allocate(Delta(nBas,nBas))
-!  allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb))
 
   allocate(ERI_AO(nBas,nBas,nBas,nBas))
   call wall_time(start_int)
@@ -94,7 +92,8 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
     ! Continue with a HFB calculation
     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)
+             nBas,nOrb,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EHFB,eHF,cHF,PHF,PanomHF,  &
+             FHF,Delta)
     call wall_time(end_HF)
 
     t_HF = end_HF - start_HF
@@ -110,8 +109,7 @@ subroutine BQuAcK(working_dir,dotest,doHFB,nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,
   deallocate(PHF)
   deallocate(PanomHF)
   deallocate(FHF)
-!  deallocate(dipole_int_MO)
+  deallocate(Delta)
-!  deallocate(ERI_MO)
   deallocate(ERI_AO)
 
 end subroutine