Choose between with and without Cartesian prod.

SHCI/SHCI.cpp

@@ -360,26 +360,64 @@ int main(int argc, char* argv[]) {
       }
     }
 
+    // Condition checking only if nonsymmetric determinants
+    if (nAlphaDets != nBetaDets && schd.alpha_beta_cartesian_product == false){
+      cout << "Number of alpha dets not equal to beta dets ("<<nAlphaDets<<") and ("<< nBetaDets <<")"<<endl;
+        exit(0);
+    }
+
     // Make HF determinant
-    Dets.resize(nAlphaDets*nBetaDets);
-    for (int i=0; i<nAlphaDets; i++) {
+    if (schd.alpha_beta_cartesian_product == true){
+      cout<< endl;
+      cout << "Running diagonalization with Cartesian product symmetrization" << endl; 
+      cout << "Subspace dimension: " << nAlphaDets * nBetaDets << endl; 
+      cout<< endl;
+      Dets.resize(nAlphaDets * nBetaDets);
+      for (int i=0; i<nAlphaDets; i++) {
       for (int j=0; j<nBetaDets; j++) {
-        Determinant& d = Dets[i*nBetaDets+j];
-        for (int a=0; a<nalpha; a++)
-          d.setocc(occAlpha[i][a]*2, true);
-        for (int a=0; a<nbeta; a++) 
-          d.setocc(occBeta[j][a]*2+1, true);
-
-
-        double E = d.Energy(I1, I2, coreE);
-        //pout << d << " Given Ref. Energy:    "
-        //     << format("%18.10f") % (E) << endl;
-        if (E < lowestEnergy) {
-          lowestEnergy = E;
-          lowestEnergyDet = i*nBetaDets+j;
+          Determinant& d = Dets[i*nBetaDets+j];
+          for (int a=0; a<nalpha; a++)
+            d.setocc(occAlpha[i][a]*2, true);
+          for (int a=0; a<nbeta; a++) 
+            d.setocc(occBeta[j][a]*2+1, true);
+
+
+          double E = d.Energy(I1, I2, coreE);
+          //pout << d << " Given Ref. Energy:    "
+          //     << format("%18.10f") % (E) << endl;
+          if (E < lowestEnergy) {
+            lowestEnergy = E;
+            lowestEnergyDet = i*nBetaDets+j;
+          }
         }
       }
     }
+
+    if (schd.alpha_beta_cartesian_product == false){
+      cout<< endl;
+      cout << "Running diagonalization WITHOUT Cartesian product symmetrization" << endl; 
+      cout << "Subspace dimension: " << nAlphaDets << endl; 
+      cout<< endl;
+      // Make HF determinant
+      Dets.resize(nAlphaDets);
+      for (int i=0; i<nAlphaDets; i++) {
+          Determinant& d = Dets[i];
+          for (int a=0; a<nalpha; a++)
+            d.setocc(occAlpha[i][a]*2, true);
+          for (int a=0; a<nbeta; a++) 
+            d.setocc(occBeta[i][a]*2+1, true);
+
+          double E = d.Energy(I1, I2, coreE);
+          //pout << d << " Given Ref. Energy:    "
+          //     << format("%18.10f") % (E) << endl;
+          if (E < lowestEnergy) {
+            lowestEnergy = E;
+            lowestEnergyDet = i;
+          }
+
+      }
+    }
+
   }
 
   HFoccupied.resize(1);

SHCI/input.cpp

@@ -121,6 +121,7 @@ void readInput(string input, std::vector<std::vector<int> >& occupied,
   schd.sampleNewDets = false;
   schd.precondition = false;
   schd.cdfciTol = 1e-12;
+  schd.alpha_beta_cartesian_product = false;
 
   schd.Bvalue = 0;
   schd.Bdirection.resize(0);
@@ -176,6 +177,8 @@ void readInput(string input, std::vector<std::vector<int> >& occupied,
       schd.ncore = atoi(tok[1].c_str());
     else if (boost::iequals(ArgName, "noio"))
       schd.io = false;
+    else if (boost::iequals(ArgName, "alpha_beta_cartesian_product"))
+      schd.alpha_beta_cartesian_product = true;
     else if (boost::iequals(ArgName, "dolcc"))
       schd.doLCC = true;
     else if (boost::iequals(ArgName, "io"))

SHCI/input.h

@@ -108,11 +108,13 @@ struct schedule {
     & cdfci_on                                \
     & cdfciIter                               \
     & z_threshold                             \
+    & alpha_beta_cartesian_product            \
     & max_determinants;
   }
   // clang-format on
 
  public:
+  bool alpha_beta_cartesian_product;
   double davidsonTol;
   double davidsonTolLoose;
   rdmType RdmType;