Triplet seems to be working. Doublet needs work.

src/csf/cfgCI_utils.c

@@ -1,5 +1,6 @@
 #include <stdint.h>
 #include <stdio.h>
+#include <assert.h>
 #include "tree_utils.h"
 
 void int_to_bin_digit(int64_t in, int count, int* out)
@@ -28,7 +29,7 @@ void getncsfs1(int *inpnsomo, int *inpms, int *outncsfs){
     int nsomo = *inpnsomo;
     int ms = *inpms;
     int nparcoupl = (nsomo + ms)/2;
-    *outncsfs = binom(nsomo, nparcoupl);
+    *outncsfs = binom((double)nsomo, (double)nparcoupl);
 }
 
 void getncsfs(int NSOMO, int MS, int *outncsfs){
@@ -39,8 +40,8 @@ void getncsfs(int NSOMO, int MS, int *outncsfs){
         (*outncsfs) = 1;
         return;
     }
-    tmpndets = binom(NSOMO, nparcoupl);
+    tmpndets = binom((double)NSOMO, (double)nparcoupl);
-    (*outncsfs) = round(tmpndets - binom(NSOMO, nparcouplp1));
+    (*outncsfs) = round(tmpndets - binom((double)NSOMO, (double)nparcouplp1));
 }
 
 #include <stdint.h>
@@ -1392,14 +1393,23 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
 void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *rows, int *cols){
 
     int NSOMO=0;
+    //printf("before getSetBits Isomo=%ld, NSOMO=%ld\n",Isomo,NSOMO);
     getSetBits(Isomo, &NSOMO);
+    //printf("Isomo=%ld, NSOMO=%ld\n",Isomo,NSOMO);
     int ndets = 0;
     int NBF = 0;
-    double dNSOMO = NSOMO*1.0;
+    //double dNSOMO = NSOMO*1.0;
     // MS = alpha_num - beta_num
-    double nalpha = (NSOMO + MS)/2.0;
+    int nalpha = (NSOMO + MS)/2;
     //printf(" in convertbftodet : MS=%d nalpha=%3.2f\n",MS,nalpha);
-    ndets = (int)binom(dNSOMO, nalpha);
+    //ndets = (int)binom(dNSOMO, nalpha);
+    if(NSOMO > 0){
+      ndets = (int)binom((double)NSOMO, (double)nalpha);
+    }
+    else if(NSOMO == 0){
+      ndets = 1;
+    }
+    else printf("Something is wrong in calcMEdetpair\n");
 
     Tree dettree = (Tree){  .rootNode = NULL, .NBF = -1 };
     dettree.rootNode = malloc(sizeof(Node));
@@ -1420,16 +1430,6 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
     }
     else{
 
-    //int addr = -1;
-    //int inpdet[NSOMO];
-    //inpdet[0] = 1;
-    //inpdet[1] = 1;
-    //inpdet[2] = 1;
-    //inpdet[3] = 0;
-    //inpdet[4] = 0;
-    //inpdet[5] = 0;
-
-    //findAddofDetDriver(&dettree, NSOMO, inpdet, &addr);
 
     int detlist[ndets];
     getDetlistDriver(&dettree, NSOMO, detlist);
@@ -1442,6 +1442,9 @@ void convertBFtoDetBasis(int64_t Isomo, int MS, double **bftodetmatrixptr, int *
     generateAllBFs(Isomo, MS, &bftree, &NBF, &NSOMO);
 
     // Initialize transformation matrix
+    //printf("MS=%d NBF=%d ndets=%d NSOMO=%d\n",MS,NBF,ndets,NSOMO);
+    assert( NBF > 0);
+    assert( ndets > 0);
     (*bftodetmatrixptr) = malloc(NBF*ndets*sizeof(double));
     (*rows) = NBF;
     (*cols) = ndets;
@@ -1496,9 +1499,10 @@ void convertBFtoDetBasisWithArrayDims(int64_t Isomo, int MS, int rowsmax, int co
     getSetBits(Isomo, &NSOMO);
     int ndets = 0;
     int NBF = 0;
-    double dNSOMO = NSOMO*1.0;
+    //double dNSOMO = NSOMO*1.0;
-    double nalpha = (NSOMO + MS)/2.0;
+    //double nalpha = (NSOMO + MS)/2.0;
-    ndets = (int)binom(dNSOMO, nalpha);
+    int nalpha = (NSOMO + MS)/2;
+    ndets = (int)binom((double)NSOMO, (double)nalpha);
 
     Tree dettree = (Tree){  .rootNode = NULL, .NBF = -1 };
     dettree.rootNode = malloc(sizeof(Node));
@@ -1582,6 +1586,7 @@ void getApqIJMatrixDims(int64_t Isomo, int64_t Jsomo, int64_t MS, int32_t *rowso
     getncsfs(NSOMOJ, MS, &NBFJ);
     (*rowsout) = NBFI;
     (*colsout) = NBFJ;
+  //exit(0);
 }
 
 void getApqIJMatrixDriver(int64_t Isomo, int64_t Jsomo, int orbp, int orbq, int64_t MS, int64_t NMO, double **CSFICSFJApqIJptr, int *rowsout, int *colsout){
@@ -1700,6 +1705,7 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
 
     int rowsbftodetI, colsbftodetI;
 
+    //printf(" 1Calling convertBFtoDetBasis Isomo=%ld MS=%ld\n",Isomo,MS);
     convertBFtoDetBasis(Isomo, MS, &bftodetmatrixI, &rowsbftodetI, &colsbftodetI);
 
     // Fill matrix
@@ -1707,7 +1713,14 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
     int colsI = 0;
 
     //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+    //printf("Isomo=%ld\n",Isomo);
     getOverlapMatrix_withDet(bftodetmatrixI, rowsbftodetI, colsbftodetI, Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
+    if(Isomo == 0){
+      rowsI = 1;
+      colsI = 1;
+    }
+
+    //printf("Isomo=%ld\n",Isomo);
 
     orthoMatrixI = malloc(rowsI*colsI*sizeof(double));
 
@@ -1719,6 +1732,7 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
 
     int rowsbftodetJ, colsbftodetJ;
 
+    //printf(" 2Calling convertBFtoDetBasis Jsomo=%ld MS=%ld\n",Jsomo,MS);
     convertBFtoDetBasis(Jsomo, MS, &bftodetmatrixJ, &rowsbftodetJ, &colsbftodetJ);
 
     int rowsJ = 0;
@@ -1726,6 +1740,10 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
     // Fill matrix
     //getOverlapMatrix(Jsomo, MS, &overlapMatrixJ, &rowsJ, &colsJ, &NSOMO);
     getOverlapMatrix_withDet(bftodetmatrixJ, rowsbftodetJ, colsbftodetJ, Jsomo, MS, &overlapMatrixJ, &rowsJ, &colsJ, &NSOMO);
+    if(Jsomo == 0){
+      rowsJ = 1;
+      colsJ = 1;
+    }
 
     orthoMatrixJ = malloc(rowsJ*colsJ*sizeof(double));
 
@@ -1743,18 +1761,25 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
 
     int transA=false;
     int transB=false;
+    //printf("1Calling blas\n");
+    //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsbftodetI,colsbftodetI,rowsA,colsA);
     callBlasMatxMat(bftodetmatrixI, rowsbftodetI, colsbftodetI, ApqIJ, rowsA, colsA, bfIApqIJ, transA, transB);
+    //printf("done\n");
 
     // now transform I in csf basis
     double *CSFIApqIJ = malloc(rowsI*colsA*sizeof(double));
     transA = false;
     transB = false;
+    //printf("2Calling blas\n");
+    //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,colsI,colsI,colsA);
     callBlasMatxMat(orthoMatrixI, rowsI, colsI, bfIApqIJ, colsI, colsA, CSFIApqIJ, transA, transB);
 
     // now transform J in BF basis
     double *CSFIbfJApqIJ = malloc(rowsI*rowsbftodetJ*sizeof(double));
     transA = false;
     transB = true;
+    //printf("3Calling blas\n");
+    //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,colsA,rowsbftodetJ,colsbftodetJ);
     callBlasMatxMat(CSFIApqIJ, rowsI, colsA, bftodetmatrixJ, rowsbftodetJ, colsbftodetJ, CSFIbfJApqIJ, transA, transB);
 
     // now transform J in CSF basis
@@ -1765,6 +1790,8 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
     double *tmpCSFICSFJApqIJ = malloc(rowsI*rowsJ*sizeof(double));
     transA = false;
     transB = true;
+    //printf("4Calling blas\n");
+    //printf("rowsA=%d colsA=%d\nrowB=%d colB=%d\n",rowsI,rowsbftodetJ,rowsJ,colsJ);
     callBlasMatxMat(CSFIbfJApqIJ, rowsI, rowsbftodetJ, orthoMatrixJ, rowsJ, colsJ, tmpCSFICSFJApqIJ, transA, transB);
     // Transfer to actual buffer in Fortran order
     for(int i = 0; i < rowsI; i++)

src/csf/conversion.irp.f

@@ -68,6 +68,7 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out)
       salpha = (s + MS)/2
       bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
     else
+      salpha = (s + MS)/2
       bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
     endif
 

src/csf/sigma_vector.irp.f

@@ -23,7 +23,13 @@
   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)
+  print *,'cfg_nsomo_min=',cfg_nsomo_min
+  print *,'cfg_nsomo_max=',cfg_nsomo_max
+  if(AND(cfg_nsomo_min , 1) .eq. 0)then
+    NSOMOMin = max(0,cfg_nsomo_min-2)
+  else
+    NSOMOMin = max(1,cfg_nsomo_min-2)
+  endif
   ! Note that here we need NSOMOMax + 2 sizes
   ialpha = (NSOMOMax + MS)/2
   NCSFMax = max(1,nint((binom(NSOMOMax,ialpha)-binom(NSOMOMax,ialpha+1)))) ! TODO: NCSFs for MS=0 (CHECK)
@@ -552,7 +558,7 @@ end subroutine get_phase_qp_to_cfg
   rows = -1
   cols = -1
   integer*8 MS
-  MS = 0
+  MS = elec_alpha_num-elec_beta_num
   real*8,dimension(:,:),allocatable :: meMatrix
   integer maxdim
 
@@ -931,7 +937,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
         noccp = holetype(k)
 
 
-        ! core-active
+        ! core-virtual
         do l = 1, n_core_orb
          jj = list_core(l)
          core_act_contrib += noccp * (2.d0 * mo_two_e_integrals_jj(jj,p) - mo_two_e_integrals_jj_exchange(jj,p))
@@ -1387,6 +1393,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
 
   allocate(diag_energies(n_CSF))
   call calculate_preconditioner_cfg(diag_energies)
+  !print *," diag energy =",diag_energies(1)
 
   MS = 0
   norm_coef_cfg=0.d0
@@ -1549,12 +1556,15 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
           !meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* h_core_ri(p,q)
           core_act_contrib = 0.0d0
           if(p.ne.q)then
-          do pp=1,n_core_orb
+            do pp=1,n_core_orb
-            n=list_core(pp)
+              n=list_core(pp)
-            core_act_contrib += 2.d0 * get_two_e_integral(p,n,q,n,mo_integrals_map) - get_two_e_integral(p,n,n,q,mo_integrals_map)
+              core_act_contrib += 2.d0 * get_two_e_integral(p,n,q,n,mo_integrals_map) - get_two_e_integral(p,n,n,q,mo_integrals_map)
-          end do
+            end do
           endif
           meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* (h_act_ri(p,q) + core_act_contrib)
+          if(jj.eq.1.and.ii.eq.1)then
+            print *,"CC=",AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI), " p=",p," q=",q
+          endif
           call omp_set_lock(lock(jj))
           do kk = 1,n_st
             psi_out(kk,jj) = psi_out(kk,jj) + meCC1 * psi_in(kk,ii)
@@ -1578,7 +1588,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   deallocate(excitationIds_single)
   deallocate(excitationTypes_single)
 
-  !print *," singles part psi(1,5)=",psi_out(1,5)
+  !print *," singles part psi(1,1)=",psi_out(1,1)
   
   allocate(listconnectedJ(N_INT,2,max(sze,10000)))
   allocate(alphas_Icfg(N_INT,2,max(sze,10000)))
@@ -1751,6 +1761,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   deallocate(diagfactors)
 
   !print *," psi(1,823)=",psi_out(1,823), " g(1 8, 3 15)=",mo_two_e_integral(1,8,3,15), " ncore=",n_core_orb
+  !print *," psi(1,1)=",psi_out(1,1)
 
   ! Add the diagonal contribution
   !$OMP DO

src/csf/tree_utils.c

@@ -1,3 +1,4 @@
+#include <assert.h>
 #include "tree_utils.h"
 
 void buildTree(Tree *bftree,
@@ -52,6 +53,7 @@ void buildTreeDriver(Tree *bftree, int NSOMO, int MS, int *NBF){
     int icpl  = 0; // keep track of the ith ms (cannot be -ve)
     int addr  = 0; // Counts the total BF's
 
+    assert(bftree->rootNode->addr == 0);
     buildTree(bftree, &(bftree->rootNode), isomo, izeros, icpl, NSOMO, MS);
 
     *NBF = bftree->rootNode->addr;
@@ -264,6 +266,8 @@ void genDetBasis(Tree *dettree, int Isomo, int MS, int *ndets){
     int NSOMO=0;
     getSetBits(Isomo, &NSOMO);
     genDetsDriver(dettree, NSOMO, MS, ndets);
+    // Closed shell case
+    if(NSOMO==0) (*ndets) = 1;
 
 }
 

src/davidson/diagonalization_hcfg.irp.f

@@ -329,10 +329,12 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
             double precision               :: ticks_0, ticks_1
             integer*8                      :: irp_imax
             irp_imax = 1
-            ticks_0 = irp_rdtsc()
+            !ticks_0 = irp_rdtsc()
             call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1)
-            ticks_1 = irp_rdtsc()
+            !ticks_1 = irp_rdtsc()
-            print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----"
+            !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----"
+            !print *,' tmpW(1,1)=',tmpW(1,1)
+            !stop
             do kk=1,N_st_diag
               do ii=1,sze_csf
                 W_csf(ii,shift+kk)=tmpW(kk,ii)

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -360,10 +360,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
             double precision               :: ticks_0, ticks_1
             integer*8                      :: irp_imax
             irp_imax = 1
-            ticks_0 = irp_rdtsc()
+            !ticks_0 = irp_rdtsc()
+            !U = 0d0
+            !U(1,1)=1.0d0
             call H_S2_u_0_nstates_openmp(W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze)
-            ticks_1 = irp_rdtsc()
+            !ticks_1 = irp_rdtsc()
-            print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----"
+            !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----"
+            !print *,W(1,1)
+            !stop
         endif
         S(1:sze,shift+1:shift+N_st_diag) = real(S_d(1:sze,1:N_st_diag))
 !      else