Working on S > 0.

src/csf/cfgCI_utils.c

@@ -342,7 +342,11 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
     ************************************/
 
 
+    //printf(" --- In convet ----\n");
     convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
+    //printf(" --- done bf det basis ---- row=%d col=%d\n",rowsbftodetI,colsbftodetI);
+
+    //printRealMatrix(bftodetmatrixI,rowsbftodetI,colsbftodetI);
 
     int rowsI = 0;
     int colsI = 0;
@@ -350,6 +354,8 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
     //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
     getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
 
+    //printf("Overlap matrix\n");
+    //printRealMatrix(overlapMatrixI,rowsI,colsI);
 
     /***********************************
          Get Orthonormalization Matrix
@@ -359,6 +365,9 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
 
     gramSchmidt(overlapMatrixI, rowsI, colsI, orthoMatrixI);
 
+    //printf("Ortho matrix\n");
+    //printRealMatrix(orthoMatrixI,rowsI,colsI);
+
     /***********************************
          Get Final CSF to Det Matrix
     ************************************/
@@ -1367,7 +1376,9 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
     int ndets = 0;
     int NBF = 0;
     double dNSOMO = NSOMO*1.0;
+    // MS = alpha_num - beta_num
     double nalpha = (NSOMO + MS)/2.0;
+    //printf(" in convertbftodet : MS=%d nalpha=%3.2f\n",MS,nalpha);
     ndets = (int)binom(dNSOMO, nalpha);
 
     Tree dettree = (Tree){  .rootNode = NULL, .NBF = -1 };

src/csf/conversion.irp.f

@@ -24,7 +24,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
   double precision, intent(out)  :: psi_coef_cfg_out(n_CSF,N_st)
   integer*8                      :: Isomo, Idomo, mask
   integer(bit_kind)              :: Ialpha(N_int) ,Ibeta(N_int)
-  integer                        :: rows, cols, i, j, k
+  integer                        :: rows, cols, i, j, k, salpha
   integer                        :: startdet, enddet
   integer                        :: ndetI
   integer                        :: getNSOMO
@@ -65,13 +65,15 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
     enddo
 
     if(iand(s,1) .EQ. 0) then
-      bfIcfg = max(1,nint((binom(s,s/2)-binom(s,(s/2)+1))))
+      salpha = (s + MS)/2
+      bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
     else
-      bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
+      bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
     endif
 
     ! perhaps blocking with CFGs of same seniority
     ! can be more efficient
+    print *,"bfIcfg=",bfIcfg, " ndetI=",ndetI
     allocate(tempBuffer(bfIcfg,ndetI))
     tempBuffer = DetToCSFTransformationMatrix(s,:bfIcfg,:ndetI)
 
@@ -99,7 +101,7 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
   double precision,intent(in)    :: psi_coef_cfg_in(n_CSF,N_st)
   double precision,intent(out)   :: psi_coef_det(N_det,N_st)
   double precision               :: tmp_psi_coef_det(maxDetDimPerBF,N_st)
-  integer                        :: s, bfIcfg
+  integer                        :: s, bfIcfg, salpha
   integer                        :: countcsf
   integer(bit_kind)              :: Ialpha(N_int), Ibeta(N_int)
   integer                        :: rows, cols, i, j, k
@@ -110,6 +112,8 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
   double precision,allocatable   :: tempCoeff (:,:)
   double precision               :: phasedet
   integer                        :: idx
+  integer MS
+  MS = elec_alpha_num-elec_beta_num
 
   countcsf = 0
 
@@ -123,7 +127,8 @@ subroutine convertWFfromCSFtoDET(N_st,psi_coef_cfg_in, psi_coef_det)
       if (psi_configuration(k,1,i) == 0_bit_kind) cycle
       s = s + popcnt(psi_configuration(k,1,i))
     enddo
-    bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1))))
+    salpha = (s + MS)/2
+    bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
 
     allocate(tempCoeff(bfIcfg,N_st))
 

src/csf/sigma_vector.irp.f

@@ -20,12 +20,15 @@
   ! The maximum number of SOMOs for the current calculation.
   ! required for the calculation of prototype arrays.
   END_DOC
+  integer MS, ialpha
+  MS = elec_alpha_num-elec_beta_num
   NSOMOMax = min(elec_num, cfg_nsomo_max + 2)
   NSOMOMin = max(0,cfg_nsomo_min-2)
   ! Note that here we need NSOMOMax + 2 sizes
-  NCSFMax = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2)-binom(NSOMOMax,((NSOMOMax+1)/2)+1)))) ! TODO: NCSFs for MS=0
+  ialpha = (NSOMOMax + MS)/2
+  NCSFMax = max(1,nint((binom(NSOMOMax,ialpha)-binom(NSOMOMax,ialpha+1)))) ! TODO: NCSFs for MS=0 (CHECK)
   NBFMax = NCSFMax
-  maxDetDimPerBF = max(1,nint((binom(NSOMOMax,(NSOMOMax+1)/2))))
+  maxDetDimPerBF = max(1,nint((binom(NSOMOMax,ialpha))))
   NMO = n_act_orb
   integer i,j,k,l
   integer startdet,enddet
@@ -34,11 +37,9 @@
   integer dimcsfpercfg
   integer detDimperBF
   real*8 :: coeff, binom1, binom2
-  integer MS
   integer ncfgpersomo
   real*8, external :: logabsgamma
   detDimperBF = 0
-  MS = elec_alpha_num-elec_beta_num
   ! number of cfgs = number of dets for 0 somos
   n_CSF = cfg_seniority_index(NSOMOMin)-1
   ncfgprev = cfg_seniority_index(NSOMOMin)
@@ -105,7 +106,11 @@
   !  endif
   !enddo
   n_CSF = 0
-  ncfgprev = cfg_seniority_index(0) ! should be 1
+  print *," -9(((((((((((((( NSOMOMin=",NSOMOMin
+  ncfgprev = cfg_seniority_index(NSOMOMin) ! can be -1 
+  if(ncfgprev.eq.-1)then
+    ncfgprev=1
+  endif
   do i=NSOMOMin,NSOMOMax+2,2
     !k=0
     !do while((cfg_seniority_index(i+2+k) .eq. -1) .and. (k.le.NSOMOMax))
@@ -126,9 +131,11 @@
       dimcsfpercfg = 2
     else
       if(iand(MS,1) .EQ. 0) then
-        dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1))))
+        ialpha = (i + MS)/2
+        dimcsfpercfg = max(1,nint((binom(i,ialpha)-binom(i,ialpha+1))))
       else
-        dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2))))
+        ialpha = (i + MS)/2
+        dimcsfpercfg = max(1,nint((binom(i,ialpha)-binom(i,ialpha+1))))
       endif
     endif
     n_CSF += ncfg*dimcsfpercfg
@@ -212,7 +219,7 @@ end subroutine get_phase_qp_to_cfg
   integer(bit_kind) :: Ialpha(N_int),Ibeta(N_int)
   integer   :: rows, cols, i, j, k
   integer   :: startdet, enddet, idx
-  integer*8 MS
+  integer*8 MS, salpha
   integer ndetI
   integer :: getNSOMO
   real*8,dimension(:,:),allocatable    :: tempBuffer
@@ -231,8 +238,11 @@ end subroutine get_phase_qp_to_cfg
     Isomo = IBSET(0_8, i) - 1_8
     ! rows = Ncsfs
     ! cols = Ndets
-    bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1))))
+    salpha = (i+MS)/2
-    ndetI = max(1,nint((binom(i,(i+1)/2))))
+    bfIcfg = max(1,nint((binom(i,salpha)-binom(i,salpha+1))))
+    ndetI =  max(1,nint((binom(i,salpha))))
+    !bfIcfg = max(1,nint((binom(i,(i+1)/2)-binom(i,((i+1)/2)+1))))
+    !ndetI =  max(1,nint((binom(i,(i+1)/2))))
 
     allocate(tempBuffer(bfIcfg,ndetI))
     call getCSFtoDETTransformationMatrix(Isomo, MS, NBFMax, maxDetDimPerBF, tempBuffer)
@@ -1364,6 +1374,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   integer(omp_lock_kind), allocatable :: lock(:)
   call omp_set_max_active_levels(1)
 
+  print *," sze = ",sze
   allocate(lock(sze))
   do i=1,sze
     call omp_init_lock(lock(i))

src/csf/tree_utils.c

@@ -311,3 +311,13 @@ void callBlasMatxMat(double *A, int rowA, int colA, double *B, int rowB, int col
             break;
     }
 }
+
+void printRealMatrix(double *orthoMatrix, int rows, int cols){
+  int i,j;
+  for(i=0;i<rows;++i){
+    for(j=0;j<cols;++j){
+      printf(" %3.5f ",orthoMatrix[i*cols + j]);
+    }
+    printf("\n");
+  }
+}

src/davidson/diagonalize_ci.irp.f

@@ -71,8 +71,10 @@ END_PROVIDER
 
    if (diag_algorithm == "Davidson") then
 
-     if (do_csf) then
+     !if (do_csf) then
-       if (sigma_vector_algorithm == 'det') then
+     if (.true.) then
+       !if (sigma_vector_algorithm == 'det') then
+       if (.false.) then
          call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
            size(CI_eigenvectors,1),CI_electronic_energy,               &
            N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)