Merge pull request #187 from v1j4y/csf_verified

CSF-CIPSI with native Sigma-Vector in CFG basis

src/csf/cfgCI_interface.f90

@@ -46,6 +46,24 @@ module cfunctions
            real   (kind=C_DOUBLE ),intent(out) :: csftodetmatrix(rowsmax,colsmax)
          end subroutine getCSFtoDETTransformationMatrix
       end interface
+       interface
+         subroutine gramSchmidt(A, m, n, B) bind(C, name='gramSchmidt')
+           import C_INT32_T, C_INT64_T, C_DOUBLE
+           integer(kind=C_INT32_T),value,intent(in) :: m
+           integer(kind=C_INT32_T),value,intent(in) :: n
+           real (kind=C_DOUBLE ),intent(in) :: A(m,n)
+           real (kind=C_DOUBLE ),intent(out) :: B(m,n)
+         end subroutine gramSchmidt
+      end interface
+       interface
+         subroutine gramSchmidt_qp(A, m, n, B) bind(C, name='gramSchmidt_qp')
+           import C_INT32_T, C_INT64_T, C_DOUBLE
+           integer(kind=C_INT32_T),value,intent(in) :: m
+           integer(kind=C_INT32_T),value,intent(in) :: n
+           real (kind=C_DOUBLE ),intent(in) :: A(m,n)
+           real (kind=C_DOUBLE ),intent(out) :: B(m,n)
+         end subroutine gramSchmidt_qp
+      end interface
 end module cfunctions
 
 subroutine f_dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) &

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,19 +29,19 @@ 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){
-    int nparcoupl = (NSOMO + MS)/2;
+    int nparcoupl = (NSOMO + MS)/2; // n_alpha
-    int nparcouplp1 = ((NSOMO + MS)/2)+1;
+    int nparcouplp1 = ((NSOMO + MS)/2)+1; // n_alpha + 1
     double tmpndets=0.0;
     if(NSOMO == 0){
         (*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>
@@ -252,6 +253,26 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
     buildTreeDriver(bftree, *NSOMO, MS, NBF);
 }
 
+void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols);
+
+void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
+  int i,j;
+  //for(j=0;j<cols;++j){
+  //  for(i=0;i<rows;++i){
+  //    printf(" %3.2f ",overlapMatrix[j*rows + i]);
+  //  }
+  //  printf("\n");
+  //}
+  // Call the function ortho_qr from qp
+  ortho_qr_csf(overlapMatrix, rows, orthoMatrix, rows, cols);
+  //for(j=0;j<cols;++j){
+  //  for(i=0;i<rows;++i){
+  //    printf(" %3.2f ",orthoMatrix[j*rows + i]);
+  //  }
+  //  printf("\n");
+  //}
+}
+
 void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
 
     // vector
@@ -341,8 +362,12 @@ void convertCSFtoDetBasis(int64_t Isomo, int MS, int rowsmax, int colsmax, doubl
          Get BFtoDeterminant Matrix
     ************************************/
 
-    printf("In convertcsftodet\n");
+
+    //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 +375,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 +386,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
     ************************************/
@@ -1340,11 +1370,11 @@ void getbftodetfunction(Tree *dettree, int NSOMO, int MS, int *BF1, double *rowv
     for(int i = 0; i < npairs; i++){
         for(int j = 0; j < NSOMO; j++) {
             inpdet[j] = detslist[i*NSOMO + j];
-            printf(" %d ",inpdet[j]);
+            //printf(" %d ",inpdet[j]);
         }
-        printf("\n");
+        //printf("\n");
         findAddofDetDriver(dettree, NSOMO, inpdet, &addr);
-        printf("(%d) - addr  = %d\n",i,addr);
+        //printf("(%d) - addr  = %d\n",i,addr);
         // Calculate the phase for cfg to QP2 conversion
         //get_phase_cfg_to_qp_inpList(inpdet, NSOMO, &phase_cfg_to_qp);
         //rowvec[addr] = 1.0 * phaselist[i]*phase_cfg_to_qp/sqrt(fac);
@@ -1363,12 +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;
-    double nalpha = (NSOMO + MS)/2.0;
+    // MS = alpha_num - beta_num
-    ndets = (int)binom(dNSOMO, nalpha);
+    int nalpha = (NSOMO + MS)/2;
+    //printf(" in convertbftodet : MS=%d nalpha=%3.2f\n",MS,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));
@@ -1389,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);
@@ -1411,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;
@@ -1465,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));
@@ -1551,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){
@@ -1669,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
@@ -1676,8 +1713,14 @@ void getApqIJMatrixDriverArrayInp(int64_t Isomo, int64_t Jsomo, int32_t orbp, in
     int colsI = 0;
 
     //getOverlapMatrix(Isomo, MS, &overlapMatrixI, &rowsI, &colsI, &NSOMO);
-    //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));
 
@@ -1689,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;
@@ -1696,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));
 
@@ -1713,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
@@ -1735,13 +1790,14 @@ 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++)
         for(int j = 0; j < rowsJ; j++)
             CSFICSFJApqIJ[j*rowsI + i] = tmpCSFICSFJApqIJ[i*rowsJ + j];
 
-
     // Garbage collection
     free(overlapMatrixI);
     free(overlapMatrixJ);

src/csf/configuration_CI_sigma_helpers.irp.f

@@ -1,3 +1,592 @@
+use bitmasks
+
+ BEGIN_PROVIDER [ integer(bit_kind), alphasIcfg_list , (N_int,2,N_configuration,mo_num*(mo_num))]
+&BEGIN_PROVIDER [ integer, NalphaIcfg_list, (N_configuration) ]
+  implicit none
+  !use bitmasks
+  BEGIN_DOC
+  ! Documentation for alphasI
+  ! Returns the associated alpha's for
+  ! the input configuration Icfg.
+  END_DOC
+
+  integer                        :: idxI ! The id of the Ith CFG
+  integer(bit_kind)              :: Icfg(N_int,2)
+  integer                        :: NalphaIcfg
+  logical,dimension(:,:),allocatable :: tableUniqueAlphas
+  integer                        :: listholes(mo_num)
+  integer                        :: holetype(mo_num) ! 1-> SOMO 2->DOMO
+  integer                        :: nholes
+  integer                        :: nvmos
+  integer                        :: listvmos(mo_num)
+  integer                        :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
+  integer*8                      :: Idomo
+  integer*8                      :: Isomo
+  integer*8                      :: Jdomo
+  integer*8                      :: Jsomo
+  integer*8                      :: diffSOMO
+  integer*8                      :: diffDOMO
+  integer*8                      :: xordiffSOMODOMO
+  integer                        :: ndiffSOMO
+  integer                        :: ndiffDOMO
+  integer                        :: nxordiffSOMODOMO
+  integer                        :: ndiffAll
+  integer                        :: i,ii
+  integer                        :: j,jj
+  integer                        :: k,kk
+  integer                        :: kstart
+  integer                        :: kend
+  integer                        :: Nsomo_I
+  integer                        :: hole
+  integer                        :: p
+  integer                        :: q
+  integer                        :: countalphas
+  logical                        :: pqAlreadyGenQ
+  logical                        :: pqExistsQ
+  logical                        :: ppExistsQ
+  integer*8                      :: MS
+
+  double precision               :: t0, t1
+  call wall_time(t0)
+
+  MS = elec_alpha_num-elec_beta_num
+
+  allocate(tableUniqueAlphas(mo_num,mo_num))
+  NalphaIcfg_list = 0
+
+  do idxI = 1, N_configuration
+
+    Icfg  = psi_configuration(:,:,idxI)
+
+    Isomo = iand(act_bitmask(1,1),Icfg(1,1))
+    Idomo = iand(act_bitmask(1,1),Icfg(1,2))
+
+    ! find out all pq holes possible
+    nholes = 0
+    ! holes in SOMO
+    do ii = 1,n_act_orb
+      i = list_act(ii)
+      if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
+        nholes += 1
+        listholes(nholes) = i
+        holetype(nholes) = 1
+      endif
+    end do
+    ! holes in DOMO
+    do ii = 1,n_act_orb
+      i = list_act(ii)
+      if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
+        nholes += 1
+        listholes(nholes) = i
+        holetype(nholes) = 2
+      endif
+    end do
+
+    ! find vmos
+    listvmos = -1
+    vmotype = -1
+    nvmos = 0
+    do ii = 1,n_act_orb
+      i = list_act(ii)
+      if(IAND(Idomo,(IBSET(0_8,i-1))) .EQ. 0) then
+        if(IAND(Isomo,(IBSET(0_8,i-1))) .EQ. 0) then
+          nvmos += 1
+          listvmos(nvmos) = i
+          vmotype(nvmos) = 1
+        else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1) then
+          nvmos += 1
+          listvmos(nvmos) = i
+          vmotype(nvmos) = 2
+        end if
+      end if
+    end do
+
+    tableUniqueAlphas = .FALSE.
+
+    ! Now find the allowed (p,q) excitations
+    Isomo = iand(act_bitmask(1,1),Icfg(1,1))
+    Idomo = iand(act_bitmask(1,1),Icfg(1,2))
+    Nsomo_I = POPCNT(Isomo)
+    if(Nsomo_I .EQ. 0) then
+      kstart = 1
+    else
+      kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))
+    endif
+    kend = idxI-1
+
+    do i = 1,nholes
+      p = listholes(i)
+      do j = 1,nvmos
+        q = listvmos(j)
+        if(p .EQ. q) cycle
+        if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
+          ! SOMO -> VMO
+          Jsomo = IBCLR(Isomo,p-1)
+          Jsomo = IBSET(Jsomo,q-1)
+          Jdomo = Idomo
+          kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
+          kend = idxI-1
+        else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
+          ! SOMO -> SOMO
+          Jsomo = IBCLR(Isomo,p-1)
+          Jsomo = IBCLR(Jsomo,q-1)
+          Jdomo = IBSET(Idomo,q-1)
+          ! Check for Minimal alpha electrons (MS)
+          if(POPCNT(Jsomo).ge.MS)then
+            kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
+            kend = idxI-1
+          else
+            cycle
+          endif
+        else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
+          ! DOMO -> VMO
+          Jsomo = IBSET(Isomo,p-1)
+          Jsomo = IBSET(Jsomo,q-1)
+          Jdomo = IBCLR(Idomo,p-1)
+          kstart = cfg_seniority_index(Nsomo_I)
+          kend = idxI-1
+        else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
+          ! DOMO -> SOMO
+          Jsomo = IBSET(Isomo,p-1)
+          Jsomo = IBCLR(Jsomo,q-1)
+          Jdomo = IBCLR(Idomo,p-1)
+          Jdomo = IBSET(Jdomo,q-1)
+          kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
+          kend = idxI-1
+        else
+          print*,"Something went wrong in obtain_associated_alphaI"
+        endif
+        ! Check for Minimal alpha electrons (MS)
+        if(POPCNT(Jsomo).lt.MS)then
+          cycle
+        endif
+
+        ! Again, we don't have to search from 1
+        ! we just use seniority to find the
+        ! first index with NSOMO - 2 to NSOMO + 2
+        ! this is what is done in kstart, kend
+
+        pqAlreadyGenQ = .FALSE.
+        ! First check if it can be generated before
+        do k = kstart, kend
+          diffSOMO = IEOR(Jsomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)))
+          ndiffSOMO = POPCNT(diffSOMO)
+          if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
+          diffDOMO = IEOR(Jdomo,iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)))
+          xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+          ndiffDOMO = POPCNT(diffDOMO)
+          nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+          nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
+          !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
+          if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
+            pqAlreadyGenQ = .TRUE.
+            ppExistsQ = .TRUE.
+            EXIT
+          endif
+          if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+            pqAlreadyGenQ = .TRUE.
+            EXIT
+          endif
+        end do
+
+        if(pqAlreadyGenQ) cycle
+
+        pqExistsQ = .FALSE.
+
+        if(.NOT. pqExistsQ) then
+          tableUniqueAlphas(p,q) = .TRUE.
+        endif
+      end do
+    end do
+
+    !print *,tableUniqueAlphas(:,:)
+
+    ! prune list of alphas
+    Isomo = Icfg(1,1)
+    Idomo = Icfg(1,2)
+    Jsomo = Icfg(1,1)
+    Jdomo = Icfg(1,2)
+    NalphaIcfg = 0
+    do i = 1, nholes
+      p = listholes(i)
+      do j = 1, nvmos
+        q = listvmos(j)
+        if(p .EQ. q) cycle
+        if(tableUniqueAlphas(p,q)) then
+          if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
+            ! SOMO -> VMO
+            Jsomo = IBCLR(Isomo,p-1)
+            Jsomo = IBSET(Jsomo,q-1)
+            Jdomo = Idomo
+          else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
+            ! SOMO -> SOMO
+            Jsomo = IBCLR(Isomo,p-1)
+            Jsomo = IBCLR(Jsomo,q-1)
+            Jdomo = IBSET(Idomo,q-1)
+            if(POPCNT(Jsomo).ge.MS)then
+              kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
+              kend = idxI-1
+            else
+              cycle
+            endif
+          else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
+            ! DOMO -> VMO
+            Jsomo = IBSET(Isomo,p-1)
+            Jsomo = IBSET(Jsomo,q-1)
+            Jdomo = IBCLR(Idomo,p-1)
+          else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
+            ! DOMO -> SOMO
+            Jsomo = IBSET(Isomo,p-1)
+            Jsomo = IBCLR(Jsomo,q-1)
+            Jdomo = IBCLR(Idomo,p-1)
+            Jdomo = IBSET(Jdomo,q-1)
+          else
+            print*,"Something went wrong in obtain_associated_alphaI"
+          endif
+
+          ! SOMO
+          !print *,i,j,"|",NalphaIcfg, Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+          if(POPCNT(Jsomo) .ge. NSOMOMin) then
+            NalphaIcfg += 1
+            alphasIcfg_list(1,1,idxI,NalphaIcfg) = Jsomo
+            alphasIcfg_list(1,2,idxI,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+            NalphaIcfg_list(idxI) = NalphaIcfg
+          endif
+        endif
+      end do
+    end do
+
+    ! Check if this Icfg has been previously generated as a mono
+    ppExistsQ = .False.
+    Isomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,1))
+    Idomo = iand(reunion_of_act_virt_bitmask(1,1),Icfg(1,2))
+    kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
+    do k = kstart, idxI-1
+      diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
+      ndiffSOMO = POPCNT(diffSOMO)
+      if (ndiffSOMO /= 2) cycle
+      diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
+      xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+      ndiffDOMO = POPCNT(diffDOMO)
+      nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+      if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4)) then
+        ppExistsQ = .TRUE.
+        EXIT
+      endif
+    end do
+    ! Diagonal part (pp,qq)
+    if(nholes > 0 .AND. (.NOT. ppExistsQ))then
+      ! SOMO
+      if(POPCNT(Jsomo) .ge. NSOMOMin) then
+        NalphaIcfg += 1
+        alphasIcfg_list(1,1,idxI,NalphaIcfg) = Icfg(1,1)
+        alphasIcfg_list(1,2,idxI,NalphaIcfg) = Icfg(1,2)
+        NalphaIcfg_list(idxI) = NalphaIcfg
+      endif
+    endif
+
+    NalphaIcfg = 0
+  enddo ! end loop idxI
+  call wall_time(t1)
+  print *, 'Preparation : ', t1 - t0
+
+END_PROVIDER
+
+  subroutine obtain_associated_alphaI(idxI, Icfg, alphasIcfg, NalphaIcfg)
+  implicit none
+  use bitmasks
+  BEGIN_DOC
+  ! Documentation for alphasI
+  ! Returns the associated alpha's for
+  ! the input configuration Icfg.
+  END_DOC
+
+  integer,intent(in)                 :: idxI ! The id of the Ith CFG
+  integer(bit_kind),intent(in)       :: Icfg(N_int,2)
+  integer,intent(out)                :: NalphaIcfg
+  integer(bit_kind),intent(out)      :: alphasIcfg(N_int,2,*)
+  logical,dimension(:,:),allocatable :: tableUniqueAlphas
+  integer                            :: listholes(mo_num)
+  integer                            :: holetype(mo_num) ! 1-> SOMO 2->DOMO
+  integer                            :: nholes
+  integer                            :: nvmos
+  integer                            :: listvmos(mo_num)
+  integer                            :: vmotype(mo_num) ! 1 -> VMO 2 -> SOMO
+  integer*8                          :: Idomo
+  integer*8                          :: Isomo
+  integer*8                          :: Jdomo
+  integer*8                          :: Jsomo
+  integer*8                          :: diffSOMO
+  integer*8                          :: diffDOMO
+  integer*8                          :: xordiffSOMODOMO
+  integer                            :: ndiffSOMO
+  integer                            :: ndiffDOMO
+  integer                            :: nxordiffSOMODOMO
+  integer                            :: ndiffAll
+  integer                            :: i, ii
+  integer                            :: j, jj
+  integer                            :: k, kk
+  integer                            :: kstart
+  integer                            :: kend
+  integer                            :: Nsomo_I
+  integer                            :: hole
+  integer                            :: p
+  integer                            :: q
+  integer                            :: countalphas
+  logical                            :: pqAlreadyGenQ
+  logical                            :: pqExistsQ
+  logical                            :: ppExistsQ
+  Isomo = iand(act_bitmask(1,1),Icfg(1,1))
+  Idomo = iand(act_bitmask(1,1),Icfg(1,2))
+  !print*,"Input cfg"
+  !call debug_spindet(Isomo,1)
+  !call debug_spindet(Idomo,1)
+
+  ! find out all pq holes possible
+  nholes = 0
+  ! holes in SOMO
+  do ii = 1,n_act_orb
+    i = list_act(ii)
+     if(POPCNT(IAND(Isomo,IBSET(0_8,i-1))) .EQ. 1) then
+        nholes += 1
+        listholes(nholes) = i
+        holetype(nholes) = 1
+     endif
+  end do
+  ! holes in DOMO
+  do ii = 1,n_act_orb
+    i = list_act(ii)
+     if(POPCNT(IAND(Idomo,IBSET(0_8,i-1))) .EQ. 1) then
+        nholes += 1
+        listholes(nholes) = i
+        holetype(nholes) = 2
+     endif
+  end do
+
+  ! find vmos
+  listvmos = -1
+  vmotype = -1
+  nvmos = 0
+  do ii = 1,n_act_orb
+    i = list_act(ii)
+     !print *,i,IBSET(0,i-1),POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))), POPCNT(IAND(Idomo,(IBSET(0_8,i-1))))
+     if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 0 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0) then
+        nvmos += 1
+        listvmos(nvmos) = i
+        vmotype(nvmos) = 1
+     else if(POPCNT(IAND(Isomo,(IBSET(0_8,i-1)))) .EQ. 1 .AND. POPCNT(IAND(Idomo,(IBSET(0_8,i-1)))) .EQ. 0 ) then
+        nvmos += 1
+        listvmos(nvmos) = i
+        vmotype(nvmos) = 2
+     end if
+  end do
+
+  !print *,"Nvmo=",nvmos
+  !print *,listvmos
+  !print *,vmotype
+
+  allocate(tableUniqueAlphas(mo_num,mo_num))
+  tableUniqueAlphas = .FALSE.
+
+  ! Now find the allowed (p,q) excitations
+  Isomo = iand(act_bitmask(1,1),Icfg(1,1))
+  Idomo = iand(act_bitmask(1,1),Icfg(1,2))
+  Nsomo_I = POPCNT(Isomo)
+  if(Nsomo_I .EQ. 0) then
+    kstart = 1
+  else
+    kstart = cfg_seniority_index(max(NSOMOMin,Nsomo_I-2))
+  endif
+  kend = idxI-1
+  !print *,"Isomo"
+  !call debug_spindet(Isomo,1)
+  !call debug_spindet(Idomo,1)
+
+  !print *,"Nholes=",nholes," Nvmos=",nvmos, " idxi=",idxI
+  !do i = 1,nholes
+  !   print *,i,"->",listholes(i)
+  !enddo
+  !do i = 1,nvmos
+  !   print *,i,"->",listvmos(i)
+  !enddo
+
+  do i = 1,nholes
+     p = listholes(i)
+     do j = 1,nvmos
+        q = listvmos(j)
+        if(p .EQ. q) cycle
+        if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
+           ! SOMO -> VMO
+           Jsomo = IBCLR(Isomo,p-1)
+           Jsomo = IBSET(Jsomo,q-1)
+           Jdomo = Idomo
+           kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
+           kend = idxI-1
+        else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
+           ! SOMO -> SOMO
+           Jsomo = IBCLR(Isomo,p-1)
+           Jsomo = IBCLR(Jsomo,q-1)
+           Jdomo = IBSET(Idomo,q-1)
+           kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-4)))
+           kend = idxI-1
+        else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
+           ! DOMO -> VMO
+           Jsomo = IBSET(Isomo,p-1)
+           Jsomo = IBSET(Jsomo,q-1)
+           Jdomo = IBCLR(Idomo,p-1)
+           kstart = cfg_seniority_index(Nsomo_I)
+           kend = idxI-1
+        else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
+           ! DOMO -> SOMO
+           Jsomo = IBSET(Isomo,p-1)
+           Jsomo = IBCLR(Jsomo,q-1)
+           Jdomo = IBCLR(Idomo,p-1)
+           Jdomo = IBSET(Jdomo,q-1)
+           kstart = max(1,cfg_seniority_index(max(NSOMOMin,Nsomo_I-2)))
+           kend = idxI-1
+        else
+           print*,"Something went wrong in obtain_associated_alphaI"
+        endif
+
+        ! Again, we don't have to search from 1
+        ! we just use seniortiy to find the
+        ! first index with NSOMO - 2 to NSOMO + 2
+        ! this is what is done in kstart, kend
+
+        pqAlreadyGenQ = .FALSE.
+        ! First check if it can be generated before
+        do k = kstart, kend
+           diffSOMO = IEOR(Jsomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
+           ndiffSOMO = POPCNT(diffSOMO)
+           if((ndiffSOMO .NE. 0) .AND. (ndiffSOMO .NE. 2)) cycle
+           diffDOMO = IEOR(Jdomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
+           xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+           ndiffDOMO = POPCNT(diffDOMO)
+           nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+           nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO
+           !if(POPCNT(IEOR(diffSOMO,diffDOMO)) .LE. 1 .AND. ndiffDOMO .LT. 3) then
+           if((ndiffSOMO+ndiffDOMO) .EQ. 0) then
+              pqAlreadyGenQ = .TRUE.
+              ppExistsQ = .TRUE.
+              EXIT
+           endif
+           if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+              pqAlreadyGenQ = .TRUE.
+              !EXIT
+              !ppExistsQ = .TRUE.
+              !print *,i,k,ndiffSOMO,ndiffDOMO
+              !call debug_spindet(Jsomo,1)
+              !call debug_spindet(Jdomo,1)
+              !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,1,k)),1)
+              !call debug_spindet(iand(reunion_of_act_virt_bitmask(1,1),psi_configuration(1,2,k)),1)
+              EXIT
+           endif
+        end do
+
+        !print *,"(,",p,",",q,")",pqAlreadyGenQ
+
+        if(pqAlreadyGenQ) cycle
+
+        pqExistsQ = .FALSE.
+        ! now check if this exists in the selected list
+        !do k = idxI+1, N_configuration
+        !   diffSOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jsomo),psi_configuration(1,1,k))
+        !   diffDOMO = IEOR(OR(reunion_of_act_virt_bitmask(1,1),Jdomo),psi_configuration(1,2,k))
+        !   ndiffSOMO = POPCNT(diffSOMO)
+        !   ndiffDOMO = POPCNT(diffDOMO)
+        !   if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
+        !      pqExistsQ = .TRUE.
+        !      EXIT
+        !   endif
+        !end do
+
+        if(.NOT. pqExistsQ) then
+           tableUniqueAlphas(p,q) = .TRUE.
+           !print *,p,q
+           !call debug_spindet(Jsomo,1)
+           !call debug_spindet(Jdomo,1)
+        endif
+     end do
+  end do
+
+  !print *,tableUniqueAlphas(:,:)
+
+  ! prune list of alphas
+  Isomo = Icfg(1,1)
+  Idomo = Icfg(1,2)
+  Jsomo = Icfg(1,1)
+  Jdomo = Icfg(1,2)
+  NalphaIcfg = 0
+  do i = 1, nholes
+     p = listholes(i)
+     do j = 1, nvmos
+        q = listvmos(j)
+        if(p .EQ. q) cycle
+        if(tableUniqueAlphas(p,q)) then
+           if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 1) then
+              ! SOMO -> VMO
+              Jsomo = IBCLR(Isomo,p-1)
+              Jsomo = IBSET(Jsomo,q-1)
+              Jdomo = Idomo
+           else if(holetype(i) .EQ. 1 .AND. vmotype(j) .EQ. 2) then
+              ! SOMO -> SOMO
+              Jsomo = IBCLR(Isomo,p-1)
+              Jsomo = IBCLR(Jsomo,q-1)
+              Jdomo = IBSET(Idomo,q-1)
+           else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 1) then
+              ! DOMO -> VMO
+              Jsomo = IBSET(Isomo,p-1)
+              Jsomo = IBSET(Jsomo,q-1)
+              Jdomo = IBCLR(Idomo,p-1)
+           else if(holetype(i) .EQ. 2 .AND. vmotype(j) .EQ. 2) then
+              ! DOMO -> SOMO
+              Jsomo = IBSET(Isomo,p-1)
+              Jsomo = IBCLR(Jsomo,q-1)
+              Jdomo = IBCLR(Idomo,p-1)
+              Jdomo = IBSET(Jdomo,q-1)
+           else
+              print*,"Something went wrong in obtain_associated_alphaI"
+           endif
+
+           ! SOMO
+           NalphaIcfg += 1
+           !print *,i,j,"|",NalphaIcfg
+           alphasIcfg(1,1,NalphaIcfg) = Jsomo
+           alphasIcfg(1,2,NalphaIcfg) = IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+           !print *,"I = ",idxI, " Na=",NalphaIcfg," - ",Jsomo, IOR(Jdomo,ISHFT(1_8,n_core_orb)-1)
+        endif
+     end do
+  end do
+
+  ! Check if this Icfg has been previously generated as a mono
+  ppExistsQ = .False.
+  Isomo = iand(act_bitmask(1,1),Icfg(1,1))
+  Idomo = iand(act_bitmask(1,1),Icfg(1,2))
+  do k = 1, idxI-1
+     diffSOMO = IEOR(Isomo,iand(act_bitmask(1,1),psi_configuration(1,1,k)))
+     diffDOMO = IEOR(Idomo,iand(act_bitmask(1,1),psi_configuration(1,2,k)))
+     xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+     ndiffSOMO = POPCNT(diffSOMO)
+     ndiffDOMO = POPCNT(diffDOMO)
+     nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+     if((ndiffSOMO+ndiffDOMO+nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+        ppExistsQ = .TRUE.
+        EXIT
+     endif
+  end do
+  ! Diagonal part (pp,qq)
+  if(nholes > 0 .AND. (.NOT. ppExistsQ))then
+     ! SOMO
+     NalphaIcfg += 1
+     !print *,p,q,"|",holetype(i),vmotype(j),NalphaIcfg
+     !call debug_spindet(Idomo,1)
+     !call debug_spindet(Jdomo,1)
+     alphasIcfg(1,1,NalphaIcfg) = Icfg(1,1)
+     alphasIcfg(1,2,NalphaIcfg) = Icfg(1,2)
+  endif
+
+  end subroutine
+
   function getNSOMO(Icfg) result(NSOMO)
     implicit none
     integer(bit_kind),intent(in)   :: Icfg(N_int,2)
@@ -8,98 +597,3 @@
        NSOMO += POPCNT(Icfg(i,1))
     enddo
   end function getNSOMO
-
-subroutine convertOrbIdsToModelSpaceIds(Ialpha, Jcfg, p, q, extype, pmodel, qmodel)
-  implicit none
-  BEGIN_DOC
-  ! This function converts the orbital ids
-  ! in real space to those used in model space
-  ! in order to identify the matrices required
-  ! for the calculation of MEs.
-  !
-  ! The type of excitations are ordered as follows:
-  ! Type 1 - SOMO -> SOMO
-  ! Type 2 - DOMO -> VMO
-  ! Type 3 - SOMO -> VMO
-  ! Type 4 - DOMO -> SOMO
-  END_DOC
-  integer(bit_kind),intent(in)         :: Ialpha(N_int,2)
-  integer(bit_kind),intent(in)         :: Jcfg(N_int,2)
-  integer,intent(in)                   :: p,q
-  integer,intent(in)                   :: extype
-  integer,intent(out)                  :: pmodel,qmodel
-  integer*8                            :: Isomo
-  integer*8                            :: Idomo
-  integer*8                            :: Jsomo
-  integer*8                            :: Jdomo
-  integer*8                            :: mask
-  integer*8                            :: Isomotmp
-  integer*8                            :: Jsomotmp
-  integer                              :: pos0,pos0prev
-
-  ! TODO Flag (print) when model space indices is > 64
-  Isomo = Ialpha(1,1)
-  Idomo = Ialpha(1,2)
-  Jsomo = Jcfg(1,1)
-  Jdomo = Jcfg(1,2)
-  pos0prev = 0
-  pmodel = p
-  qmodel = q
-
-  if(p .EQ. q) then
-     pmodel = 1
-     qmodel = 1
-  else
-     !print *,"input pq=",p,q,"extype=",extype
-     !call debug_spindet(Isomo,1)
-     !call debug_spindet(Idomo,1)
-     !call debug_spindet(Jsomo,1)
-     !call debug_spindet(Jdomo,1)
-     select case(extype)
-       case (1)
-          ! SOMO -> SOMO
-          ! remove all domos
-          !print *,"type -> SOMO -> SOMO"
-          mask = ISHFT(1_8,p) - 1
-          Isomotmp = IAND(Isomo,mask)
-          pmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
-          mask = ISHFT(1_8,q) - 1
-          Isomotmp = IAND(Isomo,mask)
-          qmodel = POPCNT(mask) - POPCNT(XOR(Isomotmp,mask))
-       case (2)
-          ! DOMO -> VMO
-          ! remove all domos except one at p
-          !print *,"type -> DOMO -> VMO"
-          mask = ISHFT(1_8,p) - 1
-          Jsomotmp = IAND(Jsomo,mask)
-          pmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
-          mask = ISHFT(1_8,q) - 1
-          Jsomotmp = IAND(Jsomo,mask)
-          qmodel = POPCNT(mask) - POPCNT(XOR(Jsomotmp,mask))
-       case (3)
-          ! SOMO -> VMO
-          !print *,"type -> SOMO -> VMO"
-          !Isomo = IEOR(Isomo,Jsomo)
-          mask = ISHFT(1_8,p) - 1
-          Isomo = IAND(Isomo,mask)
-          pmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
-          mask = ISHFT(1_8,q) - 1
-          Jsomo = IAND(Jsomo,mask)
-          qmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
-       case (4)
-          ! DOMO -> SOMO
-          ! remove all domos except one at p
-          !print *,"type -> DOMO -> SOMO"
-          !Isomo = IEOR(Isomo,Jsomo)
-          mask = ISHFT(1_8,p) - 1
-          Jsomo = IAND(Jsomo,mask)
-          pmodel = POPCNT(mask) - POPCNT(XOR(Jsomo,mask))
-          mask = ISHFT(1_8,q) - 1
-          Isomo = IAND(Isomo,mask)
-          qmodel = POPCNT(mask) - POPCNT(XOR(Isomo,mask))
-       case default
-          print *,"something is wrong in convertOrbIdsToModelSpaceIds"
-     end select
-  endif
-     !print *,p,q,"model ids=",pmodel,qmodel
-end subroutine convertOrbIdsToModelSpaceIds

src/csf/configurations.irp.f

@@ -458,8 +458,9 @@ end
 
 END_PROVIDER
 
- BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num) ]
+ BEGIN_PROVIDER [ integer, cfg_seniority_index, (0:elec_num+2) ]
 &BEGIN_PROVIDER [ integer, cfg_nsomo_max ]
+&BEGIN_PROVIDER [ integer, cfg_nsomo_min ]
   implicit none
   BEGIN_DOC
  ! Returns the index in psi_configuration of the first cfg with
@@ -467,9 +468,10 @@ END_PROVIDER
  !
  ! cfg_nsomo_max : Max number of SOMO in the current wave function
  END_DOC
- integer :: i, k, s, sold
+ integer :: i, k, s, sold, soldmin
  cfg_seniority_index(:) = -1
  sold = -1
+ soldmin = 2000
  cfg_nsomo_max = 0
  do i=1,N_configuration
    s = 0
@@ -482,6 +484,10 @@ END_PROVIDER
      cfg_seniority_index(s) = i
      cfg_nsomo_max = s
    endif
+   if (soldmin .GT. s ) then
+     soldmin = s
+     cfg_nsomo_min = s
+   endif
  enddo
 END_PROVIDER
 
@@ -743,7 +749,7 @@ BEGIN_PROVIDER [ integer(bit_kind), dominant_dets_of_cfgs, (N_int,2,N_dominant_d
  enddo
 END_PROVIDER
 
-subroutine binary_search_cfg(cfgInp,addcfg)
+subroutine binary_search_cfg(cfgInp,addcfg,bit_tmp)
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -755,29 +761,100 @@ subroutine binary_search_cfg(cfgInp,addcfg)
   END_DOC
   integer(bit_kind), intent(in)  :: cfgInp(N_int,2)
   integer          , intent(out) :: addcfg
-  integer :: i,j,k,r,l
+  integer*8,          intent(in) :: bit_tmp(0:N_configuration+1)
-  integer*8 :: key, key2
-  logical :: found
-  !integer*8, allocatable :: bit_tmp(:)
-  !integer*8, external            :: configuration_search_key
 
-  !allocate(bit_tmp(0:N_configuration))
+  logical                        :: found
-  !bit_tmp(0) = 0
+  integer                        :: l, r, j, k
-  do i=1,N_configuration
+  integer*8                      :: key
-    !bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int)
+
+  integer*8, external            :: configuration_search_key
+
+  key = configuration_search_key(cfgInp,N_int)
+
+  ! Binary search
+  l = 0
+  r = N_configuration+1
+IRP_IF WITHOUT_SHIFTRL
+  j = ishft(r-l,-1)
+IRP_ELSE
+  j = shiftr(r-l,1)
+IRP_ENDIF
+  do while (j>=1)
+    j = j+l
+    if (bit_tmp(j) == key) then
+      ! Find 1st element which matches the key
+      if (j > 1) then
+        do while (j>1 .and. bit_tmp(j-1) == key)
+          j = j-1
+        enddo
+      endif
+      ! Find correct element matching the key
+      do while (bit_tmp(j) == key)
         found = .True.
         do k=1,N_int
-      found = found .and. (psi_configuration(k,1,i) == cfgInp(k,1)) &
+          found = found .and. (psi_configuration(k,1,j) == cfgInp(k,1))&
-                    .and. (psi_configuration(k,2,i) == cfgInp(k,2))
+                        .and. (psi_configuration(k,2,j) == cfgInp(k,2))
         enddo
         if (found) then
-      addcfg = i
+          addcfg = j
-      exit
+          return
         endif
+        j = j+1
       enddo
+      addcfg = -1
+      return
+    else if (bit_tmp(j) > key) then
+      r = j
+    else
+      l = j
+    endif
+IRP_IF WITHOUT_SHIFTRL
+    j = ishft(r-l,-1)
+IRP_ELSE
+    j = shiftr(r-l,1)
+IRP_ENDIF
+  enddo
+
+  addcfg = -1
+  return
 
 end subroutine
 
+!subroutine binary_search_cfg(cfgInp,addcfg)
+!  use bitmasks
+!  implicit none
+!  BEGIN_DOC
+!  ! Documentation for binary_search
+!  ! 
+!  ! Does a binary search to find 
+!  ! the address of a configuration in a list of
+!  ! configurations.
+!  END_DOC
+!  integer(bit_kind), intent(in)  :: cfgInp(N_int,2)
+!  integer          , intent(out) :: addcfg
+!  integer :: i,j,k,r,l
+!  integer*8 :: key, key2
+!  logical :: found
+!  !integer*8, allocatable :: bit_tmp(:)
+!  !integer*8, external            :: configuration_search_key
+!
+!  !allocate(bit_tmp(0:N_configuration))
+!  !bit_tmp(0) = 0
+!  do i=1,N_configuration
+!    !bit_tmp(i) = configuration_search_key(psi_configuration(1,1,i),N_int)
+!    found = .True.
+!    do k=1,N_int
+!      found = found .and. (psi_configuration(k,1,i) == cfgInp(k,1)) &
+!                    .and. (psi_configuration(k,2,i) == cfgInp(k,2))
+!    enddo
+!    if (found) then
+!      addcfg = i
+!      exit
+!    endif
+!  enddo
+!
+!end subroutine
+!
  BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
 &BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ]
 &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]

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,9 +65,11 @@ 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))))
+      salpha = (s + MS)/2
+      bfIcfg = max(1,nint((binom(s,salpha)-binom(s,salpha+1))))
     endif
 
     ! perhaps blocking with CFGs of same seniority
@@ -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/create_excitations.irp.f

@@ -226,7 +226,7 @@ subroutine generate_all_singles_cfg(cfg,singles,n_singles,Nint)
   enddo
 end
 
-subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_singles,ex_type_singles,n_singles,Nint)
+subroutine generate_all_singles_cfg_with_type(bit_tmp,cfgInp,singles,idxs_singles,pq_singles,ex_type_singles,n_singles,Nint)
   implicit none
   use bitmasks
   BEGIN_DOC
@@ -238,6 +238,7 @@ subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_sin
   ! ex_type_singles : on output contains type of excitations  :
   !
   END_DOC
+  integer*8, intent(in)          :: bit_tmp(0:N_configuration+1)
   integer, intent(in)            :: Nint
   integer, intent(inout)         :: n_singles
   integer, intent(out)           :: idxs_singles(*)
@@ -248,20 +249,26 @@ subroutine generate_all_singles_cfg_with_type(cfgInp,singles,idxs_singles,pq_sin
   integer(bit_kind)              :: Jdet(Nint,2)
 
   integer           :: i,k, n_singles_ma, i_hole, i_particle, ex_type, addcfg
+  integer           :: ii,kk
   integer(bit_kind) :: single(Nint,2)
   logical           :: i_ok
 
+
   n_singles = 0
   !TODO
   !Make list of Somo  and Domo for holes
   !Make list of Unocc and Somo for particles
-  do i_hole = 1+n_core_orb, n_core_orb + n_act_orb
+  !do i_hole = 1+n_core_orb, n_core_orb + n_act_orb
-    do i_particle = 1+n_core_orb, n_core_orb + n_act_orb
+  do ii = 1, n_act_orb
+    i_hole = list_act(ii)
+    !do i_particle = 1+n_core_orb, n_core_orb + n_act_orb
+    do kk = 1, n_act_orb
+      i_particle = list_act(kk)
       if(i_hole .EQ. i_particle) cycle
       addcfg = -1
       call do_single_excitation_cfg_with_type(cfgInp,single,i_hole,i_particle,ex_type,i_ok)
       if (i_ok) then
-        call binary_search_cfg(single,addcfg)
+        call binary_search_cfg(single,addcfg,bit_tmp)
         if(addcfg .EQ. -1) cycle
         n_singles = n_singles + 1
         do k=1,Nint

src/csf/obtain_I_foralpha.irp.f

@@ -0,0 +1,397 @@
+subroutine obtain_connected_J_givenI(idxI, givenI, connectedI, idxs_connectedI, nconnectedI,ntotalconnectedI)
+  implicit none
+  use bitmasks
+  BEGIN_DOC
+  ! Documentation for obtain_connected_I_foralpha
+  ! This function returns all those selected configurations
+  ! which are connected to the input configuration
+  ! givenI by a single excitation.
+  !
+  ! The type of excitations are ordered as follows:
+  ! Type 1 - SOMO -> SOMO
+  ! Type 2 - DOMO -> VMO
+  ! Type 3 - SOMO -> VMO
+  ! Type 4 - DOMO -> SOMO
+  !
+  ! Order of operators
+  ! \alpha> = a^\dag_p a_q |I> = E_pq |I>
+  END_DOC
+  integer          ,intent(in)             :: idxI
+  integer(bit_kind),intent(in)             :: givenI(N_int,2)
+  integer(bit_kind),intent(out)            :: connectedI(N_int,2,*)
+  integer          ,intent(out)            :: idxs_connectedI(*)
+  integer,intent(out)                      :: nconnectedI
+  integer,intent(out)                      :: ntotalconnectedI
+  integer*8                                :: Idomo
+  integer*8                                :: Isomo
+  integer*8                                :: Jdomo
+  integer*8                                :: Jsomo
+  integer*8                                :: IJsomo
+  integer*8                                :: diffSOMO
+  integer*8                                :: diffDOMO
+  integer*8                                :: xordiffSOMODOMO
+  integer                                  :: ndiffSOMO
+  integer                                  :: ndiffDOMO
+  integer                                  :: nxordiffSOMODOMO
+  integer :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes
+  integer :: listholes(mo_num)
+  integer :: holetype(mo_num)
+  integer :: end_index
+  integer :: Nsomo_I
+
+  ! 
+  ! 2 2 1 1 0 0 : 1 1 0 0 0 0
+  !               0 0 1 1 0 0
+  !
+  ! 2 1 1 1 1 0 : 1 0 0 0 0 0
+  !               0 1 1 1 1 0
+  !xorS           0 1 0 0 1 0 : 2
+  !xorD           0 1 0 0 0 0 : 1 
+  !xorSD          0 0 0 0 1 0 : 1
+  !                           -----
+  !                             4
+  ! 1 1 1 1 1 1 : 0 0 0 0 0 0
+  !               1 1 1 1 1 1
+  !               1 1 0 0 1 1 : 4
+  !               1 1 0 0 0 0 : 2 
+  !               0 0 0 0 1 1 : 2
+  !                           -----
+  !                             8
+  !
+
+  nconnectedI = 0
+  ntotalconnectedI = 0
+  end_index = N_configuration
+
+  ! Since CFGs are sorted wrt to seniority
+  ! we don't have to search the full CFG list
+  Isomo = givenI(1,1)
+  Idomo = givenI(1,2)
+  Nsomo_I = POPCNT(Isomo)
+  end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_I+6,elec_num))-1)
+  if(end_index .LT. 0) end_index= N_configuration
+  !end_index = N_configuration
+  !print *,"Start and End = ",idxI, end_index
+
+
+  p = 0
+  q = 0
+  do i=idxI,end_index
+     !if(.True.) then
+     !  nconnectedI += 1
+     !  connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+     !  idxs_connectedI(nconnectedI)=i
+     !  cycle
+     !endif
+     Isomo = givenI(1,1)
+     Idomo = givenI(1,2)
+     Jsomo = psi_configuration(1,1,i)
+     Jdomo = psi_configuration(1,2,i)
+     diffSOMO = IEOR(Isomo,Jsomo)
+     ndiffSOMO = POPCNT(diffSOMO)
+     diffDOMO = IEOR(Idomo,Jdomo)
+     xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+     ndiffDOMO = POPCNT(diffDOMO)
+     nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+     nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+     if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+       !-------
+       ! MONO |
+       !-------
+        nconnectedI += 1
+        connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+        idxs_connectedI(nconnectedI)=i
+        ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1))
+     else if((nxordiffSOMODOMO .EQ. 8) .AND. ndiffSOMO .EQ. 4) then
+       !----------------------------
+       ! DOMO -> VMO + DOMO -> VMO |
+       !----------------------------
+       nconnectedI += 1
+       connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+       idxs_connectedI(nconnectedI)=i
+        ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1))
+     else if((nxordiffSOMODOMO .EQ. 6) .AND. ndiffSOMO .EQ. 2) then
+       !----------------------------
+       ! DOUBLE
+       !----------------------------
+       nconnectedI += 1
+       connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+       idxs_connectedI(nconnectedI)=i
+        ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1))
+     else if((nxordiffSOMODOMO .EQ. 2) .AND. ndiffSOMO .EQ. 3) then
+       !-----------------
+       ! DOUBLE
+       !-----------------
+       nconnectedI += 1
+       connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+       idxs_connectedI(nconnectedI)=i
+        ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1))
+     else if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 0) then
+       !-----------------
+       ! DOUBLE
+       !-----------------
+       nconnectedI += 1
+       connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+       idxs_connectedI(nconnectedI)=i
+        ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1))
+     else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
+       !--------
+       ! I = I |
+       !--------
+       nconnectedI += 1
+       connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+       idxs_connectedI(nconnectedI)= i
+        ! find out all pq holes possible
+        nholes = 0
+        ! holes in SOMO
+        Isomo = psi_configuration(1,1,i)
+        Idomo = psi_configuration(1,2,i)
+        do iii = 1,n_act_orb
+          ii = list_act(iii)
+           if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
+              nholes += 1
+              listholes(nholes) = ii
+              holetype(nholes) = 1
+           endif
+        end do
+        ! holes in DOMO
+        do iii = 1,n_act_orb
+          ii = list_act(iii)
+           if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
+              nholes += 1
+              listholes(nholes) = ii
+              holetype(nholes) = 2
+           endif
+        end do
+        ntotalconnectedI += max(1,(psi_config_data(i,2)-psi_config_data(i,1)+1)*nholes)
+     endif
+  end do
+
+end subroutine obtain_connected_J_givenI
+
+subroutine obtain_connected_I_foralpha(idxI, Ialpha, connectedI, idxs_connectedI, nconnectedI, excitationIds, excitationTypes, diagfactors)
+  implicit none
+  use bitmasks
+  BEGIN_DOC
+  ! Documentation for obtain_connected_I_foralpha
+  ! This function returns all those selected configurations
+  ! which are connected to the input configuration
+  ! Ialpha by a single excitation.
+  !
+  ! The type of excitations are ordered as follows:
+  ! Type 1 - SOMO -> SOMO
+  ! Type 2 - DOMO -> VMO
+  ! Type 3 - SOMO -> VMO
+  ! Type 4 - DOMO -> SOMO
+  !
+  ! Order of operators
+  ! \alpha> = a^\dag_p a_q |I> = E_pq |I>
+  END_DOC
+  integer          ,intent(in)             :: idxI
+  integer(bit_kind),intent(in)             :: Ialpha(N_int,2)
+  integer(bit_kind),intent(out)            :: connectedI(N_int,2,*)
+  integer          ,intent(out)            :: idxs_connectedI(*)
+  integer,intent(out)                      :: nconnectedI
+  integer,intent(out)                      :: excitationIds(2,*)
+  integer,intent(out)                      :: excitationTypes(*)
+  real*8 ,intent(out)                      :: diagfactors(*)
+  integer*8                                :: Idomo
+  integer*8                                :: Isomo
+  integer*8                                :: Jdomo
+  integer*8                                :: Jsomo
+  integer*8                                :: IJsomo
+  integer*8                                :: diffSOMO
+  integer*8                                :: diffDOMO
+  integer*8                                :: xordiffSOMODOMO
+  integer                                  :: ndiffSOMO
+  integer                                  :: ndiffDOMO
+  integer                        :: nxordiffSOMODOMO
+  integer                        :: iii,ii,i,j,k,l,p,q,nsomoJ,nsomoalpha,starti,endi,extyp,nholes
+  integer                        :: listholes(mo_num)
+  integer                        :: holetype(mo_num)
+  integer                        :: end_index
+  integer                        :: Nsomo_alpha
+  integer*8                      :: MS
+  MS = elec_alpha_num-elec_beta_num
+
+  nconnectedI = 0
+  end_index = N_configuration
+
+  ! Since CFGs are sorted wrt to seniority
+  ! we don't have to search the full CFG list
+  Isomo = Ialpha(1,1)
+  Idomo = Ialpha(1,2)
+  Nsomo_alpha = POPCNT(Isomo)
+  end_index = min(N_configuration,cfg_seniority_index(min(Nsomo_alpha+4,elec_num))-1)
+  if(end_index .LT. 0) end_index= N_configuration
+  end_index = N_configuration
+
+
+  p = 0
+  q = 0
+  if (N_int > 1) stop 'obtain_connected_i_foralpha : N_int > 1'
+  do i=idxI,end_index
+     Isomo = Ialpha(1,1)
+     Idomo = Ialpha(1,2)
+     Jsomo = psi_configuration(1,1,i)
+     Jdomo = psi_configuration(1,2,i)
+     ! Check for Minimal alpha electrons (MS)
+     if(POPCNT(Isomo).lt.MS)then
+       cycle
+     endif
+     diffSOMO = IEOR(Isomo,Jsomo)
+     ndiffSOMO = POPCNT(diffSOMO)
+     !if(idxI.eq.1)then
+     !  print *," \t idxI=",i," diffS=",ndiffSOMO," popJs=", POPCNT(Jsomo)," popIs=",POPCNT(Isomo)
+     !endif
+     diffDOMO = IEOR(Idomo,Jdomo)
+     xordiffSOMODOMO = IEOR(diffSOMO,diffDOMO)
+     ndiffDOMO = POPCNT(diffDOMO)
+     nxordiffSOMODOMO = POPCNT(xordiffSOMODOMO)
+     nxordiffSOMODOMO += ndiffSOMO + ndiffDOMO 
+     if((nxordiffSOMODOMO .EQ. 4) .AND. ndiffSOMO .EQ. 2) then
+        select case(ndiffDOMO)
+        case (0)
+           ! SOMO -> VMO
+           !print *,"obt SOMO -> VMO"
+           extyp = 3
+           IJsomo = IEOR(Isomo, Jsomo)
+!IRP_IF WITHOUT_TRAILZ
+!           p = (popcnt(ieor( IAND(Isomo,IJsomo) , IAND(Isomo,IJsomo) -1))-1) + 1
+!IRP_ELSE
+           p = TRAILZ(IAND(Isomo,IJsomo)) + 1
+!IRP_ENDIF
+           IJsomo = IBCLR(IJsomo,p-1)
+!IRP_IF WITHOUT_TRAILZ
+!           q = (popcnt(ieor(IJsomo,IJsomo-1))-1) + 1
+!IRP_ELSE
+           q = TRAILZ(IJsomo) + 1
+!IRP_ENDIF
+        case (1)
+           ! DOMO -> VMO
+           ! or
+           ! SOMO -> SOMO
+           nsomoJ = POPCNT(Jsomo)
+           nsomoalpha = POPCNT(Isomo)
+           if(nsomoJ .GT. nsomoalpha) then
+              ! DOMO -> VMO
+              !print *,"obt DOMO -> VMO"
+              extyp = 2
+!IRP_IF WITHOUT_TRAILZ
+!              p = (popcnt(ieor( IEOR(Idomo,Jdomo),IEOR(Idomo,Jdomo) -1))-1) + 1
+!IRP_ELSE
+              p = TRAILZ(IEOR(Idomo,Jdomo)) + 1
+!IRP_ENDIF
+              Isomo = IEOR(Isomo, Jsomo)
+              Isomo = IBCLR(Isomo,p-1)
+!IRP_IF WITHOUT_TRAILZ
+!              q = (popcnt(ieor(Isomo,Isomo-1))-1) + 1
+!IRP_ELSE
+              q = TRAILZ(Isomo) + 1
+!IRP_ENDIF
+           else
+              ! SOMO -> SOMO
+              !print *,"obt SOMO -> SOMO"
+              extyp = 1
+!IRP_IF WITHOUT_TRAILZ
+!              q = (popcnt(ieor( IEOR(Idomo,Jdomo), IEOR(Idomo,Jdomo)-1))-1) + 1
+!IRP_ELSE
+              q = TRAILZ(IEOR(Idomo,Jdomo)) + 1
+!IRP_ENDIF
+              Isomo = IEOR(Isomo, Jsomo)
+              Isomo = IBCLR(Isomo,q-1)
+!IRP_IF WITHOUT_TRAILZ
+!              p = (popcnt(ieor(Isomo,Isomo-1))-1) + 1
+!IRP_ELSE
+              p = TRAILZ(Isomo) + 1
+!IRP_ENDIF
+              ! Check for Minimal alpha electrons (MS)
+              !if(POPCNT(Isomo).lt.MS)then
+              !  cycle
+              !endif
+           end if
+        case (2)
+           ! DOMO -> SOMO
+           !print *,"obt DOMO -> SOMO"
+           extyp = 4
+           IJsomo = IEOR(Isomo, Jsomo)
+!IRP_IF WITHOUT_TRAILZ
+!           p = (popcnt(ieor( IAND(Jsomo,IJsomo), IAND(Jsomo,IJsomo)-1))-1) + 1
+!IRP_ELSE
+           p = TRAILZ(IAND(Jsomo,IJsomo)) + 1
+!IRP_ENDIF
+           IJsomo = IBCLR(IJsomo,p-1)
+!IRP_IF WITHOUT_TRAILZ
+!           q = (popcnt(ieor( IJsomo , IJsomo -1))-1) + 1
+!IRP_ELSE
+           q = TRAILZ(IJsomo) + 1
+!IRP_ENDIF
+        case default
+           print *,"something went wront in get connectedI"
+        end select
+        starti = psi_config_data(i,1)
+        endi   = psi_config_data(i,2)
+        nconnectedI += 1
+        do k=1,N_int
+          connectedI(k,1,nconnectedI) = psi_configuration(k,1,i)
+          connectedI(k,2,nconnectedI) = psi_configuration(k,2,i)
+        enddo
+        idxs_connectedI(nconnectedI)=starti
+        excitationIds(1,nconnectedI)=p
+        excitationIds(2,nconnectedI)=q
+        excitationTypes(nconnectedI) = extyp
+        diagfactors(nconnectedI) = 1.0d0
+     else if((ndiffSOMO + ndiffDOMO) .EQ. 0) then
+        ! find out all pq holes possible
+        nholes = 0
+        ! holes in SOMO
+        Isomo = psi_configuration(1,1,i)
+        Idomo = psi_configuration(1,2,i)
+        do iii = 1,n_act_orb
+          ii = list_act(iii)
+           if(POPCNT(IAND(Isomo,IBSET(0_8,ii-1))) .EQ. 1) then
+              nholes += 1
+              listholes(nholes) = ii
+              holetype(nholes) = 1
+           endif
+        end do
+        ! holes in DOMO
+        do iii = 1,n_act_orb
+          ii = list_act(iii)
+           if(POPCNT(IAND(Idomo,IBSET(0_8,ii-1))) .EQ. 1) then
+              nholes += 1
+              listholes(nholes) = ii
+              holetype(nholes) = 2
+           endif
+        end do
+
+        do k=1,nholes
+           p = listholes(k)
+           q = p
+           extyp = 1
+           if(holetype(k) .EQ. 1) then
+              starti = psi_config_data(i,1)
+              endi   = psi_config_data(i,2)
+              nconnectedI += 1
+              connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+              idxs_connectedI(nconnectedI)=starti
+              excitationIds(1,nconnectedI)=p
+              excitationIds(2,nconnectedI)=q
+              excitationTypes(nconnectedI) = extyp
+              diagfactors(nconnectedI) = 1.0d0
+           else
+              starti = psi_config_data(i,1)
+              endi   = psi_config_data(i,2)
+              nconnectedI += 1
+              connectedI(:,:,nconnectedI) = psi_configuration(:,:,i)
+              idxs_connectedI(nconnectedI)=starti
+              excitationIds(1,nconnectedI)=p
+              excitationIds(2,nconnectedI)=q
+              excitationTypes(nconnectedI) = extyp
+              diagfactors(nconnectedI) = 2.0d0
+           endif
+        enddo
+     endif
+  end do
+
+end subroutine obtain_connected_I_foralpha

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;
 
 }
 
@@ -311,3 +315,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/csf/tree_utils.h

@@ -47,6 +47,7 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
 void getSetBits(int64_t n, int *nsetbits);
 void getOverlapMatrix(int64_t Isomo, int64_t MS, double **overlapMatrixptr, int *rows, int *cols, int *NSOMOout);
 void getOverlapMatrix_withDet(double *bftodetmatrixI, int rowsbftodetI, int colsbftodetI, int64_t Isomo, int64_t MS, double **overlapMatrixI, int *rowsI, int *colsI, int *NSOMO);
+void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix);
 void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix);
 
 

src/davidson/diagonalization_hcfg.irp.f

@@ -0,0 +1,624 @@
+subroutine davidson_diag_h_cfg(dets_in,u_in,dim_in,energies,sze,sze_csf,N_st,N_st_diag,Nint,dressing_state,converged)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Davidson diagonalization.
+  !
+  ! dets_in : bitmasks corresponding to determinants
+  !
+  ! u_in : guess coefficients on the various states. Overwritten
+  !   on exit
+  !
+  ! dim_in : leftmost dimension of u_in
+  !
+  ! sze : Number of determinants
+  !
+  ! N_st : Number of eigenstates
+  !
+  END_DOC
+  integer, intent(in)            :: dim_in, sze, sze_csf, N_st, N_st_diag, Nint
+  integer(bit_kind), intent(in)  :: dets_in(Nint,2,sze)
+  double precision, intent(inout) :: u_in(dim_in,N_st_diag)
+  double precision, intent(out)  :: energies(N_st_diag)
+  integer, intent(in)            :: dressing_state
+  logical, intent(out)           :: converged
+  double precision, allocatable  :: H_jj(:)
+
+  double precision, external     :: diag_H_mat_elem, diag_S_mat_elem
+  integer                        :: i,k
+  ASSERT (N_st > 0)
+  ASSERT (sze > 0)
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  PROVIDE mo_two_e_integrals_in_map
+  allocate(H_jj(sze))
+
+  H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
+  !$OMP PARALLEL DEFAULT(NONE)                                       &
+      !$OMP  SHARED(sze,H_jj, dets_in,Nint)                    &
+      !$OMP  PRIVATE(i)
+  !$OMP DO SCHEDULE(static)
+  do i=2,sze
+    H_jj(i)  = diag_H_mat_elem(dets_in(1,1,i),Nint)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  if (dressing_state > 0) then
+    do k=1,N_st
+      do i=1,sze
+        H_jj(i)  += u_in(i,k) * dressing_column_h(i,k)
+      enddo
+    enddo
+  endif
+
+  call davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N_st,N_st_diag,Nint,dressing_state,converged)
+  deallocate(H_jj)
+end
+
+
+subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N_st,N_st_diag_in,Nint,dressing_state,converged)
+  use bitmasks
+  use mmap_module
+  implicit none
+  BEGIN_DOC
+  ! Davidson diagonalization with specific diagonal elements of the H matrix
+  !
+  ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson
+  !
+  ! dets_in : bitmasks corresponding to determinants
+  !
+  ! u_in : guess coefficients on the various states. Overwritten
+  !   on exit
+  !
+  ! dim_in : leftmost dimension of u_in
+  !
+  ! sze : Number of determinants
+  !
+  ! N_st : Number of eigenstates
+  !
+  ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze
+  !
+  END_DOC
+  integer, intent(in)            :: dim_in, sze, sze_csf, N_st, N_st_diag_in, Nint
+  integer(bit_kind), intent(in)  :: dets_in(Nint,2,sze)
+  double precision,  intent(in)  :: H_jj(sze)
+  integer, intent(in)            :: dressing_state
+  double precision, intent(inout) :: u_in(dim_in,N_st_diag_in)
+  double precision, intent(out)  :: energies(N_st_diag_in)
+
+  integer                        :: iter, N_st_diag
+  integer                        :: i,j,k,l,m,kk,ii,ll
+  logical, intent(inout)         :: converged
+
+  double precision, external     :: u_dot_v, u_dot_u
+
+  integer                        :: k_pairs, kl
+
+  integer                        :: iter2, itertot
+  double precision, allocatable  :: y(:,:), h(:,:), lambda(:)
+  double precision, allocatable  :: s_tmp(:,:)
+  double precision               :: diag_h_mat_elem
+  double precision, allocatable  :: residual_norm(:)
+  character*(16384)              :: write_buffer
+  double precision               :: to_print(2,N_st)
+  double precision               :: cpu, wall
+  integer                        :: shift, shift2, itermax, istate
+  double precision               :: r1, r2, alpha
+  logical                        :: state_ok(N_st_diag_in*davidson_sze_max)
+  integer                        :: nproc_target
+  integer                        :: order(N_st_diag_in)
+  double precision               :: cmax
+  double precision, allocatable  :: U(:,:), U_csf(:,:), overlap(:,:)
+  double precision, allocatable  :: tmpU(:,:), tmpW(:,:)
+  double precision, pointer      :: W(:,:), W_csf(:,:)
+  logical                        :: disk_based
+  double precision               :: energy_shift(N_st_diag_in*davidson_sze_max)
+
+  include 'constants.include.F'
+
+  N_st_diag = N_st_diag_in
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda
+  if (N_st_diag*3 > sze) then
+    print *,  'error in Davidson :'
+    print *,  'Increase n_det_max_full to ', N_st_diag*3
+    stop -1
+  endif
+
+  itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1
+  itertot = 0
+
+  if (state_following) then
+    allocate(overlap(N_st_diag*itermax, N_st_diag*itermax))
+  else
+    allocate(overlap(1,1))  ! avoid 'if' for deallocate
+  endif
+  overlap = 0.d0
+
+  PROVIDE nuclear_repulsion expected_s2 psi_bilinear_matrix_order psi_bilinear_matrix_order_reverse threshold_davidson_pt2 threshold_davidson_from_pt2
+
+  call write_time(6)
+  write(6,'(A)') ''
+  write(6,'(A)') 'Davidson Diagonalization'
+  write(6,'(A)') '------------------------'
+  write(6,'(A)') ''
+
+  ! Find max number of cores to fit in memory
+  ! -----------------------------------------
+
+  nproc_target = nproc
+  double precision :: rss
+  integer :: maxab
+  maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
+
+  m=1
+  disk_based = .False.
+  call resident_memory(rss)
+  do
+    r1 = 8.d0 *                                   &! bytes
+         ( dble(sze)*(N_st_diag)                  &! U
+         + dble(sze_csf)*(N_st_diag*itermax)      &! U_csf
+         + dble(sze)*(N_st_diag)                  &! W
+         + dble(sze_csf)*(N_st_diag*itermax)      &! W_csf
+         + 3.0d0*(N_st_diag*itermax)**2           &! h,y,s_tmp
+         + 1.d0*(N_st_diag*itermax)               &! lambda
+         + 1.d0*(N_st_diag)                       &! residual_norm
+                                                   ! In H_u_0_nstates_zmq
+         + 2.d0*(N_st_diag*N_det)                 &! u_t, v_t, on collector
+         + 2.d0*(N_st_diag*N_det)                 &! u_t, v_t, on slave
+         + 0.5d0*maxab                            &! idx0 in H_u_0_nstates_openmp_work_*
+         + nproc_target *                         &! In OMP section
+           ( 1.d0*(N_int*maxab)                   &! buffer
+           + 3.5d0*(maxab) )                      &! singles_a, singles_b, doubles, idx
+         ) / 1024.d0**3
+
+    if (nproc_target == 0) then
+      call check_mem(r1,irp_here)
+      nproc_target = 1
+      exit
+    endif
+
+    if (r1+rss < qp_max_mem) then
+      exit
+    endif
+
+    if (itermax > 4) then
+      itermax = itermax - 1
+    else if (m==1.and.disk_based_davidson) then
+      m=0
+      disk_based = .True.
+      itermax = 6
+    else
+      nproc_target = nproc_target - 1
+    endif
+
+  enddo
+  nthreads_davidson = nproc_target
+  TOUCH nthreads_davidson
+  call write_int(6,N_st,'Number of states')
+  call write_int(6,N_st_diag,'Number of states in diagonalization')
+  call write_int(6,sze,'Number of determinants')
+  call write_int(6,sze_csf,'Number of CSFs')
+  call write_int(6,nproc_target,'Number of threads for diagonalization')
+  call write_double(6, r1, 'Memory(Gb)')
+  if (disk_based) then
+    print *, 'Using swap space to reduce RAM'
+  endif
+
+  !---------------
+
+  write(6,'(A)') ''
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = 'Iter'
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//'       Energy        Residual '
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+  write_buffer = '====='
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') write_buffer(1:6+41*N_st)
+
+
+  if (disk_based) then
+    ! Create memory-mapped files for W and S
+    type(c_ptr) :: ptr_w, ptr_s
+    integer :: fd_s, fd_w
+    call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
+        8, fd_w, .False., ptr_w)
+    call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/))
+  else
+    allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax))
+  endif
+
+  allocate(                                                          &
+      ! Large
+      U(sze,N_st_diag),                                              &
+      U_csf(sze_csf,N_st_diag*itermax),                              &
+
+      ! Small
+      h(N_st_diag*itermax,N_st_diag*itermax),                        &
+      y(N_st_diag*itermax,N_st_diag*itermax),                        &
+      s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
+      residual_norm(N_st_diag),                                      &
+      lambda(N_st_diag*itermax))
+
+
+  h = 0.d0
+  U = 0.d0
+  y = 0.d0
+  s_tmp = 0.d0
+
+
+  ASSERT (N_st > 0)
+  ASSERT (N_st_diag >= N_st)
+  ASSERT (sze > 0)
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+
+  ! Davidson iterations
+  ! ===================
+
+  converged = .False.
+  call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),U_csf(1,1))
+  do k=N_st+1,N_st_diag
+    do i=1,sze_csf
+        call random_number(r1)
+        call random_number(r2)
+        r1 = dsqrt(-2.d0*dlog(r1))
+        r2 = dtwo_pi*r2
+        U_csf(i,k) = r1*dcos(r2) * u_csf(i,k-N_st)
+    enddo
+    U_csf(k,k) = u_csf(k,k) + 10.d0
+  enddo
+  do k=1,N_st_diag
+    call normalize(U_csf(1,k),sze_csf)
+  enddo
+  call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1))
+
+  do while (.not.converged)
+    itertot = itertot+1
+    if (itertot == 8) then
+      exit
+    endif
+
+    do iter=1,itermax-1
+
+      shift  = N_st_diag*(iter-1)
+      shift2 = N_st_diag*iter
+
+!      if ((iter > 1).or.(itertot == 1)) then
+        ! Compute |W_k> = \sum_i |i><i|H|u_k>
+        ! -----------------------------------
+
+        !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift+1),U)
+        PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
+        if ((sze > 100000).and.distributed_davidson) then
+            !call H_u_0_nstates_zmq   (W,U,N_st_diag,sze)
+            allocate(tmpW(N_st_diag,sze_csf))
+            allocate(tmpU(N_st_diag,sze_csf))
+            do kk=1,N_st_diag
+              do ii=1,sze_csf
+                tmpU(kk,ii) = U_csf(ii,shift+kk)
+              enddo
+            enddo
+            call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1)
+            do kk=1,N_st_diag
+              do ii=1,sze_csf
+                W_csf(ii,shift+kk)=tmpW(kk,ii)
+              enddo
+            enddo
+            deallocate(tmpW)
+            deallocate(tmpU)
+        else
+            !call H_u_0_nstates_openmp(W,U,N_st_diag,sze)
+            allocate(tmpW(N_st_diag,sze_csf))
+            allocate(tmpU(N_st_diag,sze_csf))
+            do kk=1,N_st_diag
+              do ii=1,sze_csf
+                tmpU(kk,ii) = U_csf(ii,shift+kk)
+              enddo
+            enddo
+            !tmpU     =0.0d0
+            !tmpU(1,2)=1.0d0
+            double precision               :: irp_rdtsc
+            double precision               :: ticks_0, ticks_1
+            integer*8                      :: irp_imax
+            irp_imax = 1
+            !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()
+            !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----"
+            do kk=1,N_st_diag
+              do ii=1,sze_csf
+                W_csf(ii,shift+kk)=tmpW(kk,ii)
+              enddo
+            enddo
+
+            !U_csf = 0.0d0
+            !U_csf(1,1) = 1.0d0
+            !u_in = 0.0d0
+            !call convertWFfromCSFtoDET(N_st_diag,tmpU,U2)
+            !call H_u_0_nstates_openmp(u_in,U2,N_st_diag,sze)
+            !call convertWFfromDETtoCSF(N_st_diag,u_in(1,1),W_csf2(1,1))
+            !do i=1,sze_csf
+            !  print *,"I=",i," qp=",W_csf2(i,1)," my=",W_csf(i,1)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1))
+            !  if(dabs(dabs(W_csf2(i,1))-dabs(W_csf(i,1))) .gt. 1.0e-10)then
+            !    print *,"somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)," diff=",dabs(W_csf2(i,1))-dabs(W_csf(i,1))
+            !  endif
+            !end do
+            !stop
+            deallocate(tmpW)
+            deallocate(tmpU)
+        endif
+!      else
+!         ! Already computed in update below
+!         continue
+!      endif
+
+      if (dressing_state > 0) then
+
+        if (N_st == 1) then
+
+          l = dressed_column_idx(1)
+          double precision :: f
+          f = 1.0d0/psi_coef(l,1)
+          do istate=1,N_st_diag
+            do i=1,sze
+              W(i,istate) += dressing_column_h(i,1) *f * U(l,istate)
+              W(l,istate) += dressing_column_h(i,1) *f * U(i,istate)
+            enddo
+
+          enddo
+
+        else
+
+          call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+            psi_coef, size(psi_coef,1), &
+            U(1,1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+
+          call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, &
+            dressing_column_h, size(dressing_column_h,1), s_tmp, size(s_tmp,1), &
+            1.d0, W(1,1), size(W,1))
+
+
+          call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, &
+            dressing_column_h, size(dressing_column_h,1), &
+            U(1,1), size(U,1), 0.d0, s_tmp, size(s_tmp,1))
+
+          call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, &
+            psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), &
+            1.d0, W(1,1), size(W,1))
+
+        endif
+      endif
+
+      !call convertWFfromDETtoCSF(N_st_diag,W,W_csf(1,shift+1))
+
+      ! Compute h_kl = <u_k | W_l> = <u_k| H |u_l>
+      ! -------------------------------------------
+
+      call dgemm('T','N', shift2, shift2, sze_csf,                       &
+          1.d0, U_csf, size(U_csf,1), W_csf, size(W_csf,1),                          &
+          0.d0, h, size(h,1))
+      call dgemm('T','N', shift2, shift2, sze_csf,                       &
+          1.d0, U_csf, size(U_csf,1), U_csf, size(U_csf,1),                          &
+          0.d0, s_tmp, size(s_tmp,1))
+
+      ! Diagonalize h
+      ! ---------------
+
+       integer :: lwork, info
+       double precision, allocatable :: work(:)
+
+       y = h
+       lwork = -1
+       allocate(work(1))
+       call dsygv(1,'V','U',shift2,y,size(y,1), &
+          s_tmp,size(s_tmp,1), lambda, work,lwork,info)
+       lwork = int(work(1))
+       deallocate(work)
+       allocate(work(lwork))
+       call dsygv(1,'V','U',shift2,y,size(y,1), &
+          s_tmp,size(s_tmp,1), lambda, work,lwork,info)
+       deallocate(work)
+       if (info /= 0) then
+         stop 'DSYGV Diagonalization failed'
+       endif
+
+      ! Compute Energy for each eigenvector
+      ! -----------------------------------
+
+      call dgemm('N','N',shift2,shift2,shift2,                       &
+          1.d0, h, size(h,1), y, size(y,1),                          &
+          0.d0, s_tmp, size(s_tmp,1))
+
+      call dgemm('T','N',shift2,shift2,shift2,                       &
+          1.d0, y, size(y,1), s_tmp, size(s_tmp,1),                  &
+          0.d0, h, size(h,1))
+
+      do k=1,shift2
+        lambda(k) = h(k,k)
+      enddo
+
+      if (state_following) then
+
+        overlap = -1.d0
+        do i=1,shift2
+          do k=1,shift2
+            overlap(k,i) = dabs(y(k,i))
+          enddo
+        enddo
+        do k=1,N_st
+          cmax = -1.d0
+          do i=1,N_st
+            if (overlap(i,k) > cmax) then
+              cmax = overlap(i,k)
+              order(k) = i
+            endif
+          enddo
+          do i=1,N_st_diag
+            overlap(order(k),i) = -1.d0
+          enddo
+        enddo
+        overlap = y
+        do k=1,N_st
+          l = order(k)
+          if (k /= l) then
+            y(1:shift2,k) = overlap(1:shift2,l)
+          endif
+        enddo
+        do k=1,N_st
+          overlap(k,1) = lambda(k)
+        enddo
+
+      endif
+
+
+      ! Express eigenvectors of h in the csf basis
+      ! ------------------------------------------
+
+      call dgemm('N','N', sze_csf, N_st_diag, shift2,                    &
+          1.d0, U_csf, size(U_csf,1), y, size(y,1), 0.d0, U_csf(1,shift2+1), size(U_csf,1))
+      call convertWFfromCSFtoDET(N_st_diag,U_csf(1,shift2+1),U)
+
+      call dgemm('N','N', sze_csf, N_st_diag, shift2,                    &
+          1.d0, W_csf, size(W_csf,1), y, size(y,1), 0.d0, W_csf(1,shift2+1), size(W_csf,1))
+      call convertWFfromCSFtoDET(N_st_diag,W_csf(1,shift2+1),W)
+
+      ! Compute residual vector and davidson step
+      ! -----------------------------------------
+
+      !if (without_diagonal) then
+      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
+      do k=1,N_st_diag
+        do i=1,sze
+           U(i,k) =  (lambda(k) * U(i,k) - W(i,k) )      &
+              /max(H_jj(i) - lambda (k),1.d-2)
+        enddo
+      enddo
+      !$OMP END PARALLEL DO
+      !else
+      !  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
+      !  do k=1,N_st_diag
+      !    do i=1,sze
+      !       U(i,k) =  (lambda(k) * U(i,k) - W(i,k) )  
+      !    enddo
+      !  enddo
+      !  !$OMP END PARALLEL DO
+      !endif
+
+      do k=1,N_st
+        residual_norm(k) = u_dot_u(U(1,k),sze)
+        to_print(1,k) = lambda(k) + nuclear_repulsion
+        to_print(2,k) = residual_norm(k)
+      enddo
+      call convertWFfromDETtoCSF(N_st_diag,U,U_csf(1,shift2+1))
+
+      if ((itertot>1).and.(iter == 1)) then
+        !don't print
+        continue
+      else
+        write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st)
+      endif
+
+      ! Check convergence
+      if (iter > 1) then
+        if (threshold_davidson_from_pt2) then
+          converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson_pt2
+        else
+          converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson
+        endif
+      endif
+
+      do k=1,N_st
+        if (residual_norm(k) > 1.d8) then
+          print *, 'Davidson failed'
+          stop -1
+        endif
+      enddo
+      if (converged) then
+        exit
+      endif
+
+      logical, external :: qp_stop
+      if (qp_stop()) then
+        converged = .True.
+        exit
+      endif
+
+
+    enddo
+
+    ! Re-contract U
+    ! -------------
+
+    call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0,      &
+        W_csf, size(W_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
+    do k=1,N_st_diag
+      do i=1,sze_csf
+        W_csf(i,k) = u_in(i,k)
+      enddo
+    enddo
+    call convertWFfromCSFtoDET(N_st_diag,W_csf,W)
+
+    call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0,      &
+        U_csf, size(U_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1))
+    do k=1,N_st_diag
+      do i=1,sze_csf
+        U_csf(i,k) = u_in(i,k)
+      enddo
+    enddo
+    call convertWFfromCSFtoDET(N_st_diag,U_csf,U)
+
+  enddo
+
+
+  call nullify_small_elements(sze,N_st_diag,U,size(U,1),threshold_davidson_pt2)
+  do k=1,N_st_diag
+    do i=1,sze
+      u_in(i,k) = U(i,k)
+    enddo
+  enddo
+
+  do k=1,N_st_diag
+    energies(k) = lambda(k)
+  enddo
+  write_buffer = '======'
+  do i=1,N_st
+    write_buffer = trim(write_buffer)//' ================ ==========='
+  enddo
+  write(6,'(A)') trim(write_buffer)
+  write(6,'(A)') ''
+  call write_time(6)
+
+  if (disk_based)then
+    ! Remove temp files
+    integer, external :: getUnitAndOpen
+    call munmap( (/int(sze,8),int(N_st_diag*itermax,8)/), 8, fd_w, ptr_w )
+    fd_w = getUnitAndOpen(trim(ezfio_work_dir)//'davidson_w','r')
+    close(fd_w,status='delete')
+  else
+    deallocate(W, W_csf)
+  endif
+
+  deallocate (                                                       &
+      residual_norm,                                                 &
+      U, U_csf, overlap,                                             &
+      h, y, s_tmp,                                                   &
+      lambda                                                         &
+      )
+  FREE nthreads_davidson
+end
+
+
+
+
+
+
+

src/davidson/diagonalization_hcsf_dressed.irp.f

@@ -89,7 +89,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
   double precision, intent(out)  :: energies(N_st_diag_in)
 
   integer                        :: iter, N_st_diag
-  integer                        :: i,j,k,l,m
+  integer                        :: i,j,k,l,m,kk
   logical, intent(inout)         :: converged
 
   double precision, external     :: u_dot_v, u_dot_u

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -149,7 +149,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   character*(16384)              :: write_buffer
   double precision               :: to_print(3,N_st)
   double precision               :: cpu, wall
-  integer                        :: shift, shift2, itermax, istate
+  integer                        :: shift, shift2, itermax, istate, ii
   double precision               :: r1, r2, alpha
   logical                        :: state_ok(N_st_diag_in*davidson_sze_max)
   integer                        :: nproc_target
@@ -356,7 +356,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         if ((sze > 100000).and.distributed_davidson) then
             call H_S2_u_0_nstates_zmq   (W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze)
         else
+            double precision               :: irp_rdtsc
+            double precision               :: ticks_0, ticks_1
+            integer*8                      :: irp_imax
+            irp_imax = 1
+            !ticks_0 = irp_rdtsc()
             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()
+            !print *,' ----Cycles:',(ticks_1-ticks_0)/dble(irp_imax)," ----"
         endif
         S(1:sze,shift+1:shift+N_st_diag) = real(S_d(1:sze,1:N_st_diag))
 !      else

src/davidson/diagonalize_ci.irp.f

@@ -5,7 +5,8 @@ BEGIN_PROVIDER [ character*(3), sigma_vector_algorithm ]
  !
  ! If 'cfg', use <Psi_csf|H|Psi_csf> in Davidson
  END_DOC
- sigma_vector_algorithm = 'det'
+ !sigma_vector_algorithm = 'det'
+ sigma_vector_algorithm = 'cfg'
 END_PROVIDER
 
 BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ]
@@ -71,7 +72,7 @@ END_PROVIDER
    if (diag_algorithm == 'Davidson') then
 
      if (do_csf) then
-!      if (sigma_vector_algorithm == 'det') then
+       if (sigma_vector_algorithm == 'det') then
          call davidson_diag_H_csf (psi_det,                  &
                                    CI_eigenvectors,          &
                                    size(CI_eigenvectors,1),  &
@@ -83,14 +84,14 @@ END_PROVIDER
                                    N_int,                    &
                                    0,                        &
                                    converged)
-!       else if (sigma_vector_algorithm == 'cfg') then
+       else if (sigma_vector_algorithm == 'cfg') then
-!       call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
+          call davidson_diag_H_cfg(psi_det,CI_eigenvectors, &
-!         size(CI_eigenvectors,1),CI_electronic_energy,               &
+          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)
+          N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)
-!       else
+       else
-!         print *, irp_here
+          print *, irp_here
-!         stop 'bug'
+          stop 'bug'
-!      endif
+       endif
      else
        call davidson_diag_HS2(psi_det,                  &
                               CI_eigenvectors,          &

src/mo_two_e_ints/core_quantities.irp.f

@@ -59,3 +59,45 @@ BEGIN_PROVIDER [ double precision, h_core_ri, (mo_num, mo_num) ]
  enddo
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ double precision, h_act_ri, (mo_num, mo_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Active Hamiltonian with 3-index exchange integrals:
+ !
+ ! $\tilde{h}{pq} = h_{pq} - \frac{1}{2}\sum_{k} g(pk,kq)$
+ END_DOC
+
+ integer :: i,j, k
+ integer :: p,q, r
+ ! core-core contribution
+ h_act_ri = core_fock_operator
+ !print *,' Bef----hact(1,14)=',h_act_ri(4,14)
+ ! act-act contribution
+ do p=1,n_act_orb
+   j=list_act(p)
+   do q=1,n_act_orb
+     i=list_act(q)
+     h_act_ri(i,j) = mo_one_e_integrals(i,j)
+   enddo
+   do r=1,n_act_orb
+     k=list_act(r)
+     do q=1,n_act_orb
+       i=list_act(q)
+       h_act_ri(i,j) = h_act_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j)
+     enddo
+   enddo
+ enddo
+ ! core-act contribution
+ !do p=1,n_act_orb
+ !  j=list_core(p)
+ !  do k=1,n_core_orb
+ !    do q=1,n_act_orb
+ !      i=list_act(q)
+ !      h_act_ri(i,j) = h_act_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j)
+ !    enddo
+ !  enddo
+ !enddo
+ !print *,' Aft----hact(1,14)=',h_act_ri(4,14), mo_one_e_integrals(4,14)
+END_PROVIDER
+

src/utils/linear_algebra.irp.f

@@ -1133,6 +1133,103 @@ subroutine ortho_svd(A,LDA,m,n)
 
 end
 
+subroutine ortho_qr_withB(A,LDA,B,m,n)
+  implicit none
+  BEGIN_DOC
+  ! Orthogonalization using Q.R factorization
+  !
+  ! A : Overlap Matrix
+  !
+  ! LDA : leftmost dimension of A
+  !
+  ! m : Number of rows of A
+  !
+  ! n : Number of columns of A
+  !
+  ! B : Output orthogonal basis
+  !
+  END_DOC
+  integer, intent(in)            :: m,n, LDA
+  double precision, intent(inout) :: A(LDA,n)
+  double precision, intent(inout) :: B(LDA,n)
+
+  integer                        :: LWORK, INFO
+  integer, allocatable           :: jpvt(:)
+  double precision, allocatable  :: TAU(:), WORK(:)
+  double precision, allocatable  :: C(:,:)
+  double precision :: norm
+  integer :: i,j
+
+  allocate (TAU(min(m,n)), WORK(1))
+  allocate (jpvt(n))
+  !print *," In function ortho"
+  B = A
+
+  jpvt(1:n)=1
+
+  LWORK=-1
+  call dgeqp3( m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
+
+  ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+  LWORK=max(n,int(WORK(1)))
+
+  deallocate(WORK)
+  allocate(WORK(LWORK))
+  call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
+  print *,A
+  print *,jpvt
+  deallocate(WORK,TAU)
+  !stop
+
+  !LWORK=-1
+  !call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
+  !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+  !LWORK=max(n,int(WORK(1)))
+
+  !deallocate(WORK)
+  !allocate(WORK(LWORK))
+  !call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
+
+  !LWORK=-1
+  !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
+  !! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
+  !LWORK=max(n,int(WORK(1)))
+
+  !deallocate(WORK)
+  !allocate(WORK(LWORK))
+  !call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
+  !
+  !allocate(C(LDA,n))
+  !call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA)
+  !norm = 0.0d0
+  !B = 0.0d0
+  !!print *,C
+  !do i=1,m
+  !  norm = 0.0d0
+  !  do j=1,n
+  !    norm = norm + C(j,i)*C(j,i)
+  !  end do
+  !  norm = 1.0d0/dsqrt(norm)
+  !  do j=1,n
+  !    B(j,i) = C(j,i)
+  !  end do
+  !end do
+  !print *,B
+
+
+  !deallocate(WORK,TAU)
+end
+
+subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
+  use iso_c_binding
+  integer(c_int32_t), value      :: LDA
+  integer(c_int32_t), value      :: m
+  integer(c_int32_t), value      :: n
+  integer(c_int16_t)             :: A(LDA,n)
+  integer(c_int16_t)             :: B(LDA,n)
+  call ortho_qr_withB(A,LDA,B,m,n)
+end subroutine ortho_qr_csf
+
 subroutine ortho_qr(A,LDA,m,n)
   implicit none
   BEGIN_DOC