diff --git a/src/HF/HFB.f90 b/src/HF/HFB.f90
index 0279f67..518c88d 100644
--- a/src/HF/HFB.f90
+++ b/src/HF/HFB.f90
@@ -144,7 +144,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
    do iorb=1,nOrb
     P(:,:)     = P(:,:)     + Occ(iorb)                        * &
                 matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) 
-    Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb)))  * &
+    Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb))))  * &
                 matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) 
    enddo
   endif
@@ -158,7 +158,7 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
    do iorb=1,nOrb
     P(:,:)     = P(:,:)     + Occ(iorb)                        * &
                 matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) 
-    Panom(:,:) = Panom(:,:) + sqrt(Occ(iorb)*(1d0-Occ(iorb)))  * &
+    Panom(:,:) = Panom(:,:) + sqrt(abs(Occ(iorb)*(1d0-Occ(iorb))))  * &
                 matmul(c(:,iorb:iorb),transpose(c(:,iorb:iorb))) 
    enddo
   endif
@@ -343,11 +343,13 @@ subroutine HFB(dotest,maxSCF,thresh,max_diis,level_shift,nNuc,ZNuc,rNuc,ENuc,
 ! Compute dipole moments, occupation numbers and || Anomalous density||
 ! also print the restart file
 
-  eigVEC(:,:) = 0d0
+  deallocate(eigVEC,eigVAL)
+  allocate(eigVEC(nOrb,nOrb),eigVAL(nOrb))
+  eigVEC(1:nOrb,1:nOrb) = 0d0
   eigVEC(1:nOrb,1:nOrb) = R(1:nOrb,1:nOrb)
-  call diagonalize_matrix(nOrb2,eigVEC,eigVAL)
-  Occ(1:nOrb)   = eigVAL(nOrb+1:nOrb2)
-  c(1:nBas,1:nOrb) = matmul(X(1:nBas,1:nOrb),eigVEC(1:nOrb,nOrb+1:nOrb2))
+  call diagonalize_matrix(nOrb,eigVEC,eigVAL)
+  Occ(1:nOrb)   = eigVAL(1:nOrb)
+  c = matmul(X,eigVEC)
   norm_anom = trace_matrix(nOrb,matmul(transpose(R(1:nOrb,nOrb+1:nOrb2)),R(1:nOrb,nOrb+1:nOrb2)))
   call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
   call write_restart_HFB(nBas,nOrb,Occ,c,chem_pot)
diff --git a/src/HF/read_restart_HFB.f90 b/src/HF/read_restart_HFB.f90
index 8d523c2..5fb09ff 100644
--- a/src/HF/read_restart_HFB.f90
+++ b/src/HF/read_restart_HFB.f90
@@ -22,6 +22,7 @@ subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
   integer                           :: nOrb_
   double precision                  :: chem_pot_
   double precision                  :: max_diff
+  double precision                  :: val_read
   double precision,allocatable      :: eigVAL(:)
   double precision,allocatable      :: c_tmp(:,:)
   double precision,allocatable      :: S_mol(:,:)
@@ -37,11 +38,18 @@ subroutine read_restart_HFB(nBas, nOrb, Occ, c, S, chem_pot)
 
   allocate(eigVAL(nOrb),c_tmp(nBas,nOrb),S_mol(nOrb,nOrb),X_mol(nOrb,nOrb))
 
+  c_tmp=0d0
+  S_mol=0d0
+  X_mol=0d0
+
   open(unit=iunit,form='unformatted',file='hfb_bin',status='old')
   read(iunit) nBas_,nOrb_ 
   read(iunit) chem_pot_
   do iorb=1,nOrb 
-   read(iunit) c_tmp(:,iorb) 
+   do ibas=1,nBas
+    read(iunit) val_read 
+    c_tmp(ibas,iorb) = val_read
+   enddo
   enddo
   do iorb=1,nOrb 
    read(iunit) eigVAL(iorb) 
diff --git a/src/HF/write_restart_HFB.f90 b/src/HF/write_restart_HFB.f90
index bcd7ef7..d5b4f97 100644
--- a/src/HF/write_restart_HFB.f90
+++ b/src/HF/write_restart_HFB.f90
@@ -19,6 +19,7 @@ subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
 ! Local variables
 
   integer                            :: ibas,iorb,iunit=666
+  double precision                   :: val_print
 
 ! Dump results
 
@@ -26,10 +27,16 @@ subroutine write_restart_HFB(nBas, nOrb, Occ, c, chem_pot)
   write(iunit) nBas,nOrb 
   write(iunit) chem_pot
   do iorb=1,nOrb 
-   write(iunit) c(:,nOrb+1-iorb) 
+   do ibas=1,nBas
+    val_print = c(ibas,nOrb+1-iorb)
+    if(abs(val_print)<1d-8) val_print=0d0
+    write(iunit) val_print
+   enddo
   enddo
   do iorb=1,nOrb 
-   write(iunit) Occ(nOrb+1-iorb) 
+   val_print = Occ(nOrb+1-iorb)
+   if(abs(val_print)<1d-8) val_print=0d0
+   write(iunit) val_print
   enddo
   write(iunit) iunit 
   close(iunit)