cleaning GHF code

src/HF/GHF.f90

@@ -22,7 +22,7 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   double precision,intent(in)   :: rNuc(nNuc,ncart)
   double precision,intent(in)   :: ENuc
 
-  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nO
   double precision,intent(in)   :: Ov(nBas,nBas)
   double precision,intent(in)   :: T(nBas,nBas)
   double precision,intent(in)   :: V(nBas,nBas)
@@ -36,9 +36,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   integer                       :: nSCF
   integer                       :: nBasSq
   integer                       :: nBas2Sq
-  integer                       :: nOa
-  integer                       :: nOb
-  integer                       :: nOcc
   integer                       :: n_diis
   double precision              :: Conv
   double precision              :: rcond
@@ -48,7 +45,6 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   double precision              :: Ex,Exaaaa,Exabba,Exbaab,Exbbbb
   double precision              :: dipole(ncart)
 
-  double precision,allocatable  :: Ca(:,:),Cb(:,:)
   double precision,allocatable  :: Jaa(:,:),Jbb(:,:)
   double precision,allocatable  :: Kaa(:,:),Kab(:,:),Kba(:,:),Kbb(:,:)
   double precision,allocatable  :: Faa(:,:),Fab(:,:),Fba(:,:),Fbb(:,:)
@@ -86,18 +82,14 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
   nBasSq  = nBas *nBas
   nBas2Sq = nBas2*nBas2
 
-  nOa  = nO(1)
-  nOb  = nO(2)
-  nOcc = nOa + nOb
-
 ! Memory allocation
 
-  allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),Jaa(nBas,nBas),Jbb(nBas,nBas), &
+  allocate(Jaa(nBas,nBas),Jbb(nBas,nBas),                                 &
            Kaa(nBas,nBas),Kab(nBas,nBas),Kba(nBas,nBas),Kbb(nBas,nBas),   &
            Faa(nBas,nBas),Fab(nBas,nBas),Fba(nBas,nBas),Fbb(nBas,nBas),   &
            Paa(nBas,nBas),Pab(nBas,nBas),Pba(nBas,nBas),Pbb(nBas,nBas),   &
-           F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2),              &
+           F(nBas2,nBas2),Fp(nBas2,nBas2),Cp(nBas2,nBas2),H(nBas2,nBas2), & 
-           H(nBas2,nBas2),S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2), &
+           S(nBas2,nBas2),X(nBas2,nBas2),err(nBas2,nBas2),                &
            err_diis(nBas2Sq,max_diis),F_diis(nBas2Sq,max_diis))
 
 ! Initialization
@@ -131,12 +123,14 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
   call mo_guess(nBas2,guess_type,S,H,X,C)
 
-  P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc)))
+! Construct super density matrix
+
+  P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO)))
 
   Paa(:,:) = P(     1:nBas ,     1:nBas )
   Pab(:,:) = P(     1:nBas ,nBas+1:nBas2)
   Pba(:,:) = P(nBas+1:nBas2,     1:nBas )
-  Pbb(:,:) = P(nBas+1:nBAs2,nBAs+1:nBas2)
+  Pbb(:,:) = P(nBas+1:nBas2,nBas+1:nBas2)
 
 !------------------------------------------------------------------------
 ! Main SCF loop
@@ -189,18 +183,14 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 !   DIIS extrapolation
 
     if(max_diis > 1) then
-
       n_diis = min(n_diis+1,max_diis)
       call DIIS_extrapolation(rcond,nBas2Sq,nBas2Sq,n_diis,err_diis(:,1:n_diis),F_diis(:,1:n_diis),err,F)
-
     end if
 
 !   Level-shifting
 
     if(level_shift > 0d0 .and. Conv > thresh) then
-
+      call level_shifting(level_shift,nBas,nO,S,C,F)
-      call level_shifting(level_shift,nBas,nOa+nOb,S,C,F)
-
     end if
 
 !  Transform Fock matrix in orthogonal basis
@@ -216,17 +206,15 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
     C(:,:) = matmul(X,Cp)
 
-!   Form individual coefficient matrices
-
-    Ca(1:nBas,1:nBas2) = C(     1:nBas ,1:nBas2) 
-    Cb(1:nBas,1:nBas2) = C(nBas+1:nBas2,1:nBas2) 
-
 !   Mix guess for UHF solution in singlet states
 
 !   if(nSCF == 1) call mix_guess(nBas,nO,mix,c)
 
+!   Form super density matrix
+
+    P(:,:) = matmul(C(:,1:nO),transpose(C(:,1:nO)))
+
 !   Compute individual density matrices
-    P(:,:) = matmul(C(:,1:nOcc),transpose(C(:,1:nOcc)))
 
     Paa(:,:) = P(     1:nBas ,     1:nBas )
     Pab(:,:) = P(     1:nBas ,nBas+1:nBas2)
@@ -298,9 +286,12 @@ subroutine GHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
 
   end if
 
-! Compute final UHF energy
+! Compute dipole moments
+
+  call dipole_moment(nBas2,P,nNuc,ZNuc,rNuc,dipole_int,dipole)
+
+! Compute final GHF energy
 
-! call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int,dipole)
 ! call print_GHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EHF,dipole)
 
 end subroutine 

src/QuAcK/GQuAcK.f90

@@ -1,5 +1,5 @@
 subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, &
-                  nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO,                 &
+                  nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO,                    &
                   maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP,                       &
                   spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF,       &
                   maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW,                              &
@@ -16,11 +16,10 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   logical,intent(in)            :: doG0W0,doevGW,doqsGW
 
   integer,intent(in)            :: nNuc,nBas
-  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nC
-  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nO
-  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nV
-  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nR
-  integer,intent(in)            :: nS(nspin)
   double precision,intent(in)   :: ENuc
 
   double precision,intent(in)   :: ZNuc(nNuc),rNuc(nNuc,ncart)
@@ -77,6 +76,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
 
   integer                       :: ixyz
   integer                       :: nBas2
+  integer                       :: nS
 
   write(*,*)
   write(*,*) '*******************************'
@@ -131,17 +131,24 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   ! 4-index transform 
   
   allocate(Ca(nBas,nBas2),Cb(nBas,nBas2),ERI_tmp(nBas2,nBas2,nBas2,nBas2))
+
   Ca(:,:) = cHF(1:nBas,1:nBas2)
   Cb(:,:) = cHF(nBas+1:nBas2,1:nBas2)
+
   call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Ca,Ca,Ca,ERI_AO,ERI_tmp)
   ERI_MO(:,:,:,:) = ERI_tmp(:,:,:,:)
+
   call AOtoMO_integral_transform_GHF(nBas,nBas2,Ca,Cb,Ca,Cb,ERI_AO,ERI_tmp)
   ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
+
   call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Ca,Cb,Ca,ERI_AO,ERI_tmp)
   ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
+
   call AOtoMO_integral_transform_GHF(nBas,nBas2,Cb,Cb,Cb,Cb,ERI_AO,ERI_tmp)
   ERI_MO(:,:,:,:) = ERI_MO(:,:,:,:) + ERI_tmp(:,:,:,:)
 
+  deallocate(Ca,Cb,ERI_tmp)
+
   call wall_time(end_AOtoMO)
 
   t_AOtoMO = end_AOtoMO - start_AOtoMO
@@ -152,10 +159,12 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
 ! Stability analysis of HF solution !
 !-----------------------------------!
 
+  nS = (nO - nC)*(nV - nR)
+
   if(dostab) then
 
     call wall_time(start_stab)
-    call GHF_stability(nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),epsHF,ERI_MO)
+    call GHF_stability(nBas2,nC,nO,nV,nR,nS,epsHF,ERI_MO)
     call wall_time(end_stab)
 
     t_stab = end_stab - start_stab
@@ -173,7 +182,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
   if(doMP) then
 
     call wall_time(start_MP)
-    call GMP(doMP2,reg_MP,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),ERI_MO,ENuc,EHF,epsHF)
+    call GMP(doMP2,reg_MP,nBas2,nC,nO,nV,nR,ERI_MO,ENuc,EHF,epsHF)
     call wall_time(end_MP)
 
     t_MP = end_MP - start_MP
@@ -192,7 +201,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
 
     call wall_time(start_RPA)
     call GRPA(dophRPA,dophRPAx,doppRPA,TDA,doACFDT,exchange_kernel, & 
-              nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV),ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
+              nBas2,nC,nO,nV,nR,nS,ENuc,EHF,ERI_MO,dipole_int_MO,epsHF,cHF,S)
     call wall_time(end_RPA)
 
     t_RPA = end_RPA - start_RPA
@@ -233,7 +242,7 @@ subroutine GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqs
     call wall_time(start_GW)
     call GGW(doG0W0,doevGW,doqsGW,maxSCF_GW,thresh_GW,max_diis_GW,doACFDT,       &
              exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA, & 
-             lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,sum(nC),sum(nO),sum(nV),sum(nR),sum(nO)*sum(nV), &
+             lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas2,nC,nO,nV,nR,nS,      &
              EHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,epsHF)
     call wall_time(end_GW)
   

src/QuAcK/QuAcK.f90

@@ -20,7 +20,6 @@ program QuAcK
   integer                       :: nO(nspin)
   integer                       :: nV(nspin)
   integer                       :: nR(nspin)
-  integer                       :: nS(nspin)
   double precision              :: ENuc
 
   double precision,allocatable  :: ZNuc(:),rNuc(:,:)
@@ -126,9 +125,6 @@ program QuAcK
 ! nV   = number of virtual orbitals (see below)  !
 ! nR   = number of Rydberg orbitals              !
 ! nBas = number of basis functions (see below)   !
-!      = nO + nV                                 !
-! nS   = number of single excitation             !
-!      = nO*nV                                   !
 !------------------------------------------------!
 
   call read_molecule(nNuc,nO,nC,nR)
@@ -143,7 +139,6 @@ program QuAcK
 !---------------------------------------!
 
   call read_basis_pyscf(nBas,nO,nV)
-  nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:))
 
 !--------------------------------------!
 ! Read one- and two-electron integrals !
@@ -193,7 +188,7 @@ program QuAcK
     call RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                               &
                 dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
                 doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,            &
-                doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,   &
+                doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,      &
                 S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
                 guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA,              &
                 maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW,  &
@@ -208,7 +203,7 @@ program QuAcK
     call UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                                      &
                 dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
                 doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,            &
-                doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,   &
+                doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,      &
                 S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
                 guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA,              &
                 maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW,  &
@@ -221,7 +216,7 @@ program QuAcK
 
   if(doGQuAcK) & 
     call GQuAcK(doGHF,dostab,doMP2,dophRPA,dophRPAx,doppRPA,doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2,     &
-                nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO,                 &
+                nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO,    &
                 maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP,                           & 
                 spin_conserved,spin_flip,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF,           & 
                 maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW,                                  &

src/QuAcK/RQuAcK.f90

@@ -1,7 +1,7 @@
 subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                               &   
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & 
                   doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,            &
-                  doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,   & 
+                  doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,      & 
                   S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
                   guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA,                       &
                   maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW,  & 
@@ -25,11 +25,10 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC
   logical,intent(in)            :: doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh
 
   integer,intent(in)            :: nNuc,nBas
-  integer,intent(in)            :: nC(nspin)
+  integer,intent(in)            :: nC
-  integer,intent(in)            :: nO(nspin)
+  integer,intent(in)            :: nO
-  integer,intent(in)            :: nV(nspin)
+  integer,intent(in)            :: nV
-  integer,intent(in)            :: nR(nspin)
+  integer,intent(in)            :: nR
-  integer,intent(in)            :: nS(nspin)
   double precision,intent(in)   :: ENuc
 
   double precision,intent(in)   :: ZNuc(nNuc),rNuc(nNuc,ncart)
@@ -93,6 +92,7 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC
   double precision,allocatable  :: dipole_int_MO(:,:,:)
   double precision,allocatable  :: ERI_MO(:,:,:,:)
   integer                       :: ixyz
+  integer                       :: nS
 
   write(*,*)
   write(*,*) '******************************'
@@ -168,6 +168,8 @@ subroutine RQuAcK(doRHF,doROHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCC
 ! Stability analysis of HF solution !
 !-----------------------------------!
 
+  nS = (nO - nC)*(nV - nR)
+
   if(dostab) then
 
     call wall_time(start_stab)

src/QuAcK/UQuAcK.f90

@@ -1,7 +1,7 @@
 subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,                                      &
                   dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, &
                   doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW,            &
-                  doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,nS,ENuc,ZNuc,rNuc,   &
+                  doG0T0pp,doevGTpp,doqsGTpp,doG0T0eh,doevGTeh,doqsGTeh,nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc,      &
                   S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift,                     &
                   guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA,              &
                   maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW,  &
@@ -27,7 +27,6 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,
   integer,intent(in)            :: nO(nspin)
   integer,intent(in)            :: nV(nspin)
   integer,intent(in)            :: nR(nspin)
-  integer,intent(in)            :: nS(nspin)
   double precision,intent(in)   :: ENuc
 
   double precision,intent(in)   :: ZNuc(nNuc),rNuc(nNuc,ncart)
@@ -91,6 +90,7 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,
   double precision,allocatable  :: dipole_int_aa(:,:,:),dipole_int_bb(:,:,:)
   double precision,allocatable  :: ERI_aaaa(:,:,:,:),ERI_aabb(:,:,:,:),ERI_bbbb(:,:,:,:)
   integer                       :: ixyz
+  integer                       :: nS(nspin)
 
   write(*,*)
   write(*,*) '********************************'
@@ -186,6 +186,8 @@ subroutine UQuAcK(doUHF,dostab,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,
 ! Stability analysis of HF solution !
 !-----------------------------------!
 
+  nS(:) = (nO(:) - nC(:))*(nV(:) - nR(:))
+
   if(dostab) then
 
     call wall_time(start_stab)