A few changes 1. Heisenberg model has been added, 2. btas library has been added 3. USE_BTAS flas is added to makefile to disable btas library

Author: sanshar

Operators.C

@@ -739,7 +739,7 @@ void SpinAdapted::CreCreDesComp::build(const SpinBlock& b)
       SpinAdapted::operatorfunctions::TensorTrace(rightBlock, *op, &b, &(b.get_stateInfo()), *this, 1.0);
     }  
 
-  if (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY){
+  if (dmrginp.hamiltonian() != HUBBARD){
     if (loopBlock->has(CRE_DESCOMP))
       {
 
@@ -904,7 +904,7 @@ void SpinAdapted::Ham::build(const SpinBlock& b)
   f = boost::bind(&opxop::cxcddcomp, rightBlock, _1, &b, op_add); 
   for_all_multithread(leftBlock->get_op_array(CRE), f);  
 
-  if (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY) {    
+  if (dmrginp.hamiltonian() != HUBBARD) {    
     op_add =  (otherBlock->get_op_array(CRE_DESCOMP).is_local() && loopBlock->get_op_array(CRE_DES).is_local())? op_array : op_distributed;
     f = boost::bind(&opxop::cdxcdcomp, otherBlock, _1, &b, op_add);
     for_all_multithread(loopBlock->get_op_array(CRE_DES), f);

StateInfo.h

@@ -62,14 +62,20 @@ class StateInfo
 
   std::vector<SpinQuantum> quanta; //the quantas present
   std::vector<int> quantaStates;  //the number of each quanta
-  ObjectMatrix<char> allowedQuanta; //
-  ObjectMatrix< std::vector<int> > quantaMap;
-  std::vector<int> leftUnMapQuanta;
+  ObjectMatrix<char> allowedQuanta; //tensor product of left and right stateinfo states
+
+  //takes i and j quanta of the left and the right state and gives a vector of all resulting quantas. (for abelian symmetry the vector would be of length 1)
+  ObjectMatrix< std::vector<int> > quantaMap; 
+
+  //for the uncollectedstateinfo each quanta is made up of a left and a right quanta, these two vectors unmap an uncollected quata back to the left and the right quanta.
+  std::vector<int> leftUnMapQuanta;  
   std::vector<int> rightUnMapQuanta;
 
   int totalStates;
+
+  //the next three are concerned with blocking and unblocking the states
   std::vector<int> unBlockedIndex;
-  std::vector< std::vector<int> > oldToNewState;
+  std::vector< std::vector<int> > oldToNewState; //quanta[I] is the same as quanta[Ij] where Ij are the indices inside the Ith vector.
   std::vector<int> newQuantaMap;
   bool initialised;
  public:

btas

@@ -0,0 +1 @@
+Subproject commit 41f057fc97409503653d19abb60853c33a599a94

csf.C

@@ -336,9 +336,21 @@ std::vector<SpinAdapted::Csf > SpinAdapted::CSFUTIL::spinfockstrings(const std::
     occ_rep2[dmrginp.spatial_to_spin()[I]+2*irrepsize-2] = 1;
     Slater s1(occ_rep1, 1), s2(occ_rep2, 1); map<Slater, double > m1, m2;
     m1[s1]= 1.0; m2[s2] = 1.0;
+
+    if(dmrginp.hamiltonian() == HEISENBERG) {
+      thisSiteCsf.push_back( Csf(m2, 1, SpinSpace(1), 1, IrrepVector(Irrep.getirrep(), irrepsize-1))); //1,1,L    
+      for (int i=tensorops[0].Szops.size(); i> 0; i--)
+	ladderentry.push_back(applyTensorOp(tensorops[0], i-1));
+      singleSiteLadder.push_back(ladderentry);
+
+      numcsfs.push_back(thisSiteCsf.size());    
+      for (int i=0; i<thisSiteCsf.size(); i++)
+	singleSiteCsf.push_back( thisSiteCsf[i]);
+      continue;
+    }
+
     thisSiteCsf.push_back( Csf(m1, 0, SpinSpace(0), 0, IrrepVector(0,0))); //0,0,0
-    thisSiteCsf.push_back( Csf(m2, 1, SpinSpace(1), 1, IrrepVector(Irrep.getirrep(), irrepsize-1))); //1,1,L
-    
+    thisSiteCsf.push_back( Csf(m2, 1, SpinSpace(1), 1, IrrepVector(Irrep.getirrep(), irrepsize-1))); //1,1,L    
     ladderentry.push_back(Csf(m1, 0, SpinSpace(0), 0, IrrepVector(0,0))); singleSiteLadder.push_back(ladderentry);
     ladderentry.clear();
 
@@ -525,13 +537,15 @@ std::vector< SpinAdapted::Csf > SpinAdapted::Csf::distribute (const int n, const
   na = (n + sp) / 2;
   nb = (n - sp) / 2;
 
+  if(dmrginp.hamiltonian() == HEISENBERG && nb != 0) return s;
   // let's count how many left and right orbitals there actually are                                                                      
   int actualOrbs = dmrginp.spatial_to_spin()[right] - dmrginp.spatial_to_spin()[left];
   int actualNA, actualNB;
   actualNA = actualOrbs/2;
   actualNB = actualOrbs/2;
 
-  // cannot form this combination                                                                                                         
+  // cannot form this combination 
+
   if (na > actualNA || nb > actualNB || na < 0 || nb < 0)
     {
       return s;

density.C

@@ -256,7 +256,7 @@ void DensityMatrix::add_onedot_noise(const std::vector<Wavefunction>& wave_solut
       for_all_multithread(leftBlock->get_op_array(CRE), onedot_noise);
     }
 
-    if (dmrginp.hamiltonian() == QUANTUM_CHEMISTRY) {
+    if (dmrginp.hamiltonian() != HUBBARD) {
     if (leftBlock->has(CRE_CRE))
     {
       onedot_noise.set_opType(CRE_CRE);

dmrg.C

@@ -50,6 +50,7 @@ Sandeep Sharma and Garnet K.-L. Chan
 
 void dmrg(double sweep_tol);
 void restart(double sweep_tol, bool reset_iter);
+void dmrg_stateSpecific(double sweep_tol, int targetState);
 void ReadInput(char* conf);
 void fullrestartGenblock();
 void license() {
@@ -107,8 +108,15 @@ int calldmrg(char* input, char* output)
 
   SweepParams sweep_copy;
   bool direction_copy; int restartsize_copy;
-
-
+  
+  /*
+  {
+   SweepParams sweepParams;
+   sweepParams.ls_dw.resize(0);
+   sweepParams.ls_energy.resize(0);
+   Sweep::do_overlap(sweepParams, false, true, false, 0);
+  }
+  */  
   switch(dmrginp.calc_type()) {
 
   case (DMRG):
@@ -292,6 +300,7 @@ int calldmrg(char* input, char* output)
 
     break;
   }
+  //case (EXCITEDDMRG) :
 
   return 0;
 
@@ -339,6 +348,7 @@ void restart(double sweep_tol, bool reset_iter)
   int domoreIter = 2;
 
   sweepParams.restorestate(direction, restartsize);
+
   if(reset_iter) { //this is when you restart from the start of the sweep
     sweepParams.set_sweep_iter() = 0;
     sweepParams.set_restart_iter() = 0;
@@ -422,11 +432,12 @@ void dmrg(double sweep_tol)
   bool dodiis = false;
 
   int domoreIter = 0;
+  bool direction;
 
   //initialize array of size m_maxiter or dmrginp.max_iter() for dw and energy
-
-
+  sweepParams.current_root() = 0;
   last_fe = Sweep::do_one(sweepParams, true, true, false, 0);
+  direction = false;
   while ((fabs(last_fe - old_fe) > sweep_tol) || (fabs(last_be - old_be) > sweep_tol) || 
 	 (dmrginp.algorithm_method() == TWODOT_TO_ONEDOT && dmrginp.twodot_to_onedot_iter()+1 >= sweepParams.get_sweep_iter()) )
     {
@@ -435,6 +446,7 @@ void dmrg(double sweep_tol)
       if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
 	break;
       last_be = Sweep::do_one(sweepParams, false, false, false, 0);
+      direction = true;
       if (dmrginp.outputlevel() > 0) 
          pout << "Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
 
@@ -458,6 +470,7 @@ void dmrg(double sweep_tol)
       }
 
       last_fe = Sweep::do_one(sweepParams, false, true, false, 0);
+      direction = false;
 
       new_states=sweepParams.get_keep_states();
 
@@ -486,5 +499,73 @@ void dmrg(double sweep_tol)
   }
 
   const int nroots = dmrginp.nroots(sweepParams.get_sweep_iter());
+
+  if(nroots > 1 && dmrginp.doStateSpecific()) {
+    dmrginp.setStateSpecific() = true;
+    for (int i=0; i<nroots; i++) {
+      if (i != 0) 
+	SweepGenblock::do_one(sweepParams, false, direction, false, 0, i);
+      dmrg_stateSpecific(sweep_tol, i);
+    }
+  }
+}
+
+void dmrg_stateSpecific(double sweep_tol, int targetState)
+{
+  double last_fe = 10.e6;
+  double last_be = 10.e6;
+  double old_fe = 0.;
+  double old_be = 0.;
+  int ls_count=0;
+  SweepParams sweepParams;
+  int old_states=sweepParams.get_keep_states();
+  int new_states;
+  double old_error=0.0;
+  double old_energy=0.0;
+  // warm up sweep ...
+
+  bool direction;
+  int restartsize;
+  sweepParams.restorestate(direction, restartsize);
+
+  //initialize array of size m_maxiter or dmrginp.max_iter() for dw and energy
+  sweepParams.current_root() = targetState;
+
+  sweepParams.set_sweep_iter() = 0;
+  sweepParams.set_restart_iter() = 0;
+
+  last_fe = Sweep::do_one(sweepParams, false, direction, true, 0);
+
+  while ((fabs(last_fe - old_fe) > sweep_tol) || (fabs(last_be - old_be) > sweep_tol)  )
+    {
+      old_fe = last_fe;
+      old_be = last_be;
+      if(dmrginp.max_iter() <= 10) //CHANGE THIS TO SOME INPUT PARAMETER
+	break;
+      last_be = Sweep::do_one(sweepParams, false, !direction, false, 0);
+      if (dmrginp.outputlevel() > 0) 
+         pout << "Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
+
+      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
+         break;
+
+
+      last_fe = Sweep::do_one(sweepParams, false, direction, false, 0);
+
+      new_states=sweepParams.get_keep_states();
+
+
+      if (dmrginp.outputlevel() > 0)
+         pout << "Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
+
+    }
+  if(dmrginp.max_iter() <= sweepParams.get_sweep_iter()) {
+    
+    pout << "Maximum sweep iterations achieved " << std::endl;
+  }
+
+  //SweepGenblock::makeRotationsAndOverlaps(sweepParams, false, !direction, targetState);
+  //SweepGenblock::makeRotationsAndOverlaps(sweepParams, false, direction, targetState);
+  const int nroots = dmrginp.nroots(sweepParams.get_sweep_iter());
 }
 

dmrg_tests/heisenberg_2d/FCIDUMP

@@ -0,0 +1,68 @@
+ &FCI NORB= 16,NELEC= 16,MS2= 0,
+  ORBSYM=1,1,1,1, 1,1,1,1, 1,1,1,1, 1,1,1,1,
+  ISYM=1
+ &END
+-0.25 1 1 13 13
+-0.5 1 13 13 1
+-0.25 1 1 5 5
+-0.5 1 5 5 1
+-0.25 1 1 4 4
+-0.5 1 4 4 1
+-0.25 1 1 2 2
+-0.5 1 2 2 1
+-0.25 2 2 14 14
+-0.5 2 14 14 2
+-0.25 2 2 6 6
+-0.5 2 6 6 2
+-0.25 2 2 3 3
+-0.5 2 3 3 2
+-0.25 3 3 15 15
+-0.5 3 15 15 3
+-0.25 3 3 7 7
+-0.5 3 7 7 3
+-0.25 3 3 4 4
+-0.5 3 4 4 3
+-0.25 4 4 16 16
+-0.5 4 16 16 4
+-0.25 4 4 8 8
+-0.5 4 8 8 4
+-0.25 5 5 9 9
+-0.5 5 9 9 5
+-0.25 5 5 8 8
+-0.5 5 8 8 5
+-0.25 5 5 6 6
+-0.5 5 6 6 5
+-0.25 6 6 10 10
+-0.5 6 10 10 6
+-0.25 6 6 7 7
+-0.5 6 7 7 6
+-0.25 7 7 11 11
+-0.5 7 11 11 7
+-0.25 7 7 8 8
+-0.5 7 8 8 7
+-0.25 8 8 12 12
+-0.5 8 12 12 8
+-0.25 9 9 13 13
+-0.5 9 13 13 9
+-0.25 9 9 12 12
+-0.5 9 12 12 9
+-0.25 9 9 10 10
+-0.5 9 10 10 9
+-0.25 10 10 14 14
+-0.5 10 14 14 10
+-0.25 10 10 11 11
+-0.5 10 11 11 10
+-0.25 11 11 15 15
+-0.5 11 15 15 11
+-0.25 11 11 12 12
+-0.5 11 12 12 11
+-0.25 12 12 16 16
+-0.5 12 16 16 12
+-0.25 13 13 16 16
+-0.5 13 16 16 13
+-0.25 13 13 14 14
+-0.5 13 14 14 13
+-0.25 14 14 15 15
+-0.5 14 15 15 14
+-0.25 15 15 16 16
+-0.5 15 16 16 15

dmrg_tests/heisenberg_2d/dmrg.conf

@@ -0,0 +1,18 @@
+nelec 16
+spin 0
+
+
+hf_occ integral
+schedule 
+0 100  1.0e-6 1e-05 
+4 500 1e-6 1.0e-6
+end
+maxiter 20
+twodot_to_onedot 10
+sweep_tol 1e-9
+
+orbitals FCIDUMP
+
+heisenberg
+
+outputlevel 2

dmrg_tests/runtest

@@ -32,7 +32,11 @@ python ../test_energy.py 1  1.0e-6 -6.56819216
 python ../test_onepdm.py spatial_onepdm.0.0 spat_onepdm.0.0 1e-7
 cd ../
 
-
+echo "performing heisenberg test..."
+cd heisenberg_2d
+ ../../block.spin_adapted dmrg.conf >dmrg.out
+python ../test_energy.py 1  1.0e-6 -11.2284830996
+cd ../
 
 #this first calculated twopdm with a small M for two lowest states (so the energy is not converged)
 #checks if the twopdm traced with integrals matches the energies

fci.C

@@ -75,7 +75,9 @@ void SpinAdapted::Sweep::fullci(double sweep_tol)
   double tol = sweepParams.get_davidson_tol();
 
   pout << "\t\t\t Solving the Wavefunction "<<endl;
-  Solver::solve_wavefunction(solution, energies, big, tol, BASIC, false, true, false, sweepParams.get_additional_noise());
+  int currentState = 0;
+  std::vector<Wavefunction> lowerStates;
+  Solver::solve_wavefunction(solution, energies, big, tol, BASIC, false, true, false, sweepParams.get_additional_noise(), currentState, lowerStates);
   for (int i=0; i<nroots; i++) {
     pout << "fullci energy "<< energies[i]+dmrginp.get_coreenergy()<<endl;
   }