Merge branch 'dev'

OverlapHelement.C

@@ -39,7 +39,7 @@ int main(int argc, char* argv []) {
   //readMPSFromDiskAndInitializeStaticVariables();
   //initializeGlobalMPS(mpsstate);
 
-  if (rank =0) 
+  if (rank ==0) 
     printf("Reading file %s\n", argv[2]);
   test(argv[2]);
   //exit(0);

Stackspinblock.C

@@ -878,6 +878,26 @@ void StackSpinBlock::BuildSumBlock(int condition, StackSpinBlock& lBlock, StackS
   }
   dmrginp.buildsumblock -> start();
   BuildSumBlockSkeleton(condition, lBlock, rBlock, collectQuanta, compState);
+  //To reduced the number of operators in onepdm calculation, two index operators are removed. However, they are still needed for single orbital sites.
+  //For the sum block built from dummy sites (singlet embedding), there are no two index operators.
+  if(dmrginp.do_npdm_ops() && (dmrginp.calc_type() == RESTART_ONEPDM || dmrginp.calc_type() == ONEPDM))
+    if(sites.size()==1)
+    {
+      if (!is_direct()) {
+        ops[CRE_DES] = make_new_stackop(CRE_DES, true);
+        ops[CRE_CRE] = make_new_stackop(CRE_CRE, true);
+        ops[DES_CRE] = make_new_stackop(DES_CRE, true);
+        ops[DES_DES] = make_new_stackop(DES_DES, true);
+      }
+      else
+      {
+        ops[CRE_DES] = make_new_stackop(CRE_DES, false);
+        ops[CRE_CRE] = make_new_stackop(CRE_CRE, false);
+        ops[DES_CRE] = make_new_stackop(DES_CRE, false);
+        ops[DES_DES] = make_new_stackop(DES_DES, false);
+
+      }
+    }
 
   totalMemory = build_iterators();
   if (totalMemory != 0)

dmrg.C

@@ -490,6 +490,18 @@ int calldmrg(char* input, char* output)
   else if (dmrginp.calc_type() == TINYCALC) {
     Sweep::tiny(sweep_tol);
   }
+  else if (dmrginp.calc_type() == RESTART_ONEPDM) {
+    Npdm::npdm(NPDM_ONEPDM);
+  }
+  else if (dmrginp.calc_type() == RESTART_T_ONEPDM) {
+    Npdm::npdm(NPDM_ONEPDM,true);
+  }
+  else if (dmrginp.calc_type() == RESTART_TWOPDM) {
+    Npdm::npdm(NPDM_TWOPDM);
+  }
+  else if (dmrginp.calc_type() == RESTART_T_TWOPDM) {
+    Npdm::npdm(NPDM_TWOPDM,true);
+  }
   else if (dmrginp.calc_type() == RESTART_THREEPDM) {
     Npdm::npdm(NPDM_THREEPDM);
   }

fciqmchelper.C

@@ -4,6 +4,7 @@
 #include "Stackspinblock.h"
 #include "initblocks.h"
 #ifndef SERIAL
+#include "mpi.h"
 #include <boost/mpi.hpp>
 #endif
 
@@ -396,13 +397,36 @@ namespace SpinAdapted{
     p2out << system<<endl;
 
     std::vector<Matrix> Rotationa, Rotationb;
+    /*
+    int sysquanta = system.get_stateInfo().quanta.size();
+    if (mpigetrank() == 0) {
+      Rotationa = statea.getSiteTensors(0);
+      Rotationb = stateb.getSiteTensors(0);
+      Rotationa.resize(sysquanta);
+      Rotationb.resize(sysquanta);
+    }
+
+#ifndef SERIAL
+    mpi::communicator world;
+    mpi::broadcast(world, Rotationa, 0);
+    mpi::broadcast(world, Rotationb, 0);
+#endif
+
+    system.transform_operators(const_cast<std::vector<Matrix>&>(Rotationa), 
+			       const_cast<std::vector<Matrix>&>(Rotationb));
+    */
     int sys_add = true; bool direct = true; 
 
+
     std::vector<int> rotSites(2,0);
     int sweepIters = dmrginp.spinAdapted() ? dmrginp.last_site() -2 : dmrginp.last_site()/2-2;
+    int normToComp = sweepIters/2;
+
     for (int i=0; i<sweepIters-1; i++) {
-      pout << i<<" out of "<<sweepIters-1<<endl;
+      pout << i<<" out of "<<sweepIters-1<<"  ";
       StackSpinBlock newSystem;
+      if (i>=normToComp && !system.has(CRE_DESCOMP))
+	system.addAllCompOps();
       system.addAdditionalOps();
 
       StackSpinBlock dotsite(i+1, i+1, integralIndex, false);
@@ -419,8 +443,14 @@ namespace SpinAdapted{
       mpi::broadcast(calc, Rotationa, 0);
       mpi::broadcast(calc, Rotationb, 0);
 #endif
-      InitBlocks::InitNewSystemBlock(system, dotsite, newSystem, 0, statebindex, sys_add, direct, integralIndex, DISTRIBUTED_STORAGE, false, true);
-
+      if (i <normToComp )  {
+	pout << "norm ops "<<endl;
+	InitBlocks::InitNewSystemBlock(system, dotsite, newSystem, 0, statebindex, sys_add, direct, integralIndex, DISTRIBUTED_STORAGE, true, false);
+      }
+      else {
+	pout << "comp ops "<<endl;
+	InitBlocks::InitNewSystemBlock(system, dotsite, newSystem, 0, statebindex, sys_add, direct, integralIndex, DISTRIBUTED_STORAGE, false, true);
+      }
       newSystem.transform_operators(const_cast<std::vector<Matrix>&>(Rotationa), 
 				    const_cast<std::vector<Matrix>&>(Rotationb));
       {
@@ -436,6 +466,9 @@ namespace SpinAdapted{
     StackSpinBlock dotsite1(sweepIters, sweepIters, integralIndex, false);
     StackSpinBlock dotsite2(sweepIters+1, sweepIters+1, integralIndex, false);
 
+    //For molecule has at most 4 orbitals, there is at most one iteration.
+    //System does not have CompOps.
+	  system.addAllCompOps();
 
     system.addAdditionalOps();
     InitBlocks::InitNewSystemBlock(system, dotsite1, newSystem, 0, statebindex, sys_add, direct, integralIndex, DISTRIBUTED_STORAGE, false, true);
@@ -450,6 +483,22 @@ namespace SpinAdapted{
     stateaw.LoadWavefunctionInfo(s, rotSites, statea, true);
     statebw.LoadWavefunctionInfo(s, rotSites, stateb, true);
 
+#ifndef SERIAL
+    mpi::broadcast(calc, stateaw, 0);
+    mpi::broadcast(calc, statebw, 0);
+#endif
+    if (mpigetrank() != 0) {
+      double* dataa = Stackmem[omprank].allocate(stateaw.memoryUsed());
+      stateaw.set_data(dataa);
+      stateaw.allocateOperatorMatrix();
+      double* datab = Stackmem[omprank].allocate(statebw.memoryUsed());
+      statebw.set_data(datab);
+      statebw.allocateOperatorMatrix();
+    }
+    calc.barrier();
+    MPI_Bcast(stateaw.get_data(), stateaw.memoryUsed(), MPI_DOUBLE, 0, Calc);
+    MPI_Bcast(statebw.get_data(), statebw.memoryUsed(), MPI_DOUBLE, 0, Calc);
+
     StackWavefunction temp; temp.initialise(stateaw);
     temp.Clear();
 

input.C

@@ -184,6 +184,7 @@ void SpinAdapted::Input::initialize_defaults()
   m_maxM = 0;
   m_lastM = 500;
   m_startM = 250;
+  m_bra_M = 0;
 
   m_calc_ri_4pdm=false;
   m_store_ripdm_readable=false;
@@ -967,6 +968,20 @@ SpinAdapted::Input::Input(const string& config_name) {
       {
         m_npdm_multinode = false;
       }
+      else if (boost::iequals(keyword, "specificpdm"))
+      {
+        if(tok.size() ==2)
+          m_specificpdm.push_back(atoi(tok[1].c_str()));
+        else if(tok.size() ==3) 
+        {
+          m_specificpdm.push_back(atoi(tok[1].c_str()));
+          m_specificpdm.push_back(atoi(tok[2].c_str()));
+        }
+        else {
+          pout << "keyword "<<keyword<<" should be followed by one or two numbers and then an endline";
+          abort();
+        }
+      }
 
 
 
@@ -1175,6 +1190,19 @@ SpinAdapted::Input::Input(const string& config_name) {
       else if (boost::iequals(keyword,  "reset_iterations") || boost::iequals(keyword,  "reset_iter") || boost::iequals(keyword,  "reset_iters")) {
 	m_reset_iterations = true;
       }
+      else if (boost::iequals(keyword,  "bra_M"))
+      {
+        //The number of renormalized bases in bra wave function. 
+        //Currently only used in transition denstiy matrix calcualtions.
+	      if (tok.size() != 2) {
+	        pout << "keyword bra_M should be followed by a single integer and then an endline"<<endl;
+	        pout << "error found in the following line "<<endl;
+	        pout << msg<<endl;
+	        abort();
+	      }
+        m_bra_M = atoi(tok[1].c_str());
+
+      }
 
       else
       {
@@ -1246,6 +1274,8 @@ SpinAdapted::Input::Input(const string& config_name) {
   mpi::broadcast(world, m_calc_type, 0);
   mpi::broadcast(world, m_calc_procs, 0);
   mpi::broadcast(world, m_baseState, 0);
+  mpi::broadcast(world, m_useSharedMemory, 0);
+
 #endif
 
   //make the scratch files

input.h

@@ -87,6 +87,7 @@ class Input {
   int m_lastM;
   int m_startM;
   int m_maxM;
+  int m_bra_M;
   int m_integral_disk_storage_thresh;
   int m_num_Integrals;
 
@@ -114,6 +115,7 @@ class Input {
   bool m_pdm_unsorted;
   bool m_store_nonredundant_pdm;
   bool m_npdm_intermediate;
+  std::vector<int> m_specificpdm;
   bool m_npdm_multinode;
   bool m_set_Sz;
   int m_maxiter;
@@ -189,12 +191,12 @@ class Input {
     ar & m_deflation_min_size & m_deflation_max_size & m_outputlevel & m_reorderfile;
     ar & m_algorithm_type & m_twodot_to_onedot_iter & m_orbformat & m_calc_procs;
     ar & m_nquanta & m_sys_add & m_env_add & m_do_fci & m_no_transform ;
-    ar & m_do_pdm & m_do_npdm_ops & m_do_npdm_in_core & m_npdm_generate & m_new_npdm_code  & m_transition_diff_spatial_irrep & m_occupied_orbitals;
+    ar & m_do_pdm & m_do_npdm_ops & m_do_npdm_in_core & m_npdm_generate & m_new_npdm_code & m_specificpdm & m_transition_diff_spatial_irrep & m_occupied_orbitals;
     ar & m_store_spinpdm &m_spatpdm_disk_dump & m_pdm_unsorted & m_npdm_intermediate & m_npdm_multinode;
     ar & m_maxj & m_ninej & m_maxiter & m_do_deriv & m_oneindex_screen_tol & m_twoindex_screen_tol & m_quantaToKeep & m_noise_type;
     ar & m_sweep_tol & m_restart & m_backward & m_fullrestart & m_restart_warm & m_reset_iterations & m_calc_type & m_ham_type & m_warmup;
     ar & m_do_diis & m_diis_error & m_start_diis_iter & m_diis_keep_states & m_diis_error_tol & m_num_spatial_orbs;
-    ar & m_spatial_to_spin & m_spin_to_spatial & m_maxM & m_schedule_type_backward & m_schedule_type_default &m_integral_disk_storage_thresh;
+    ar & m_spatial_to_spin & m_spin_to_spatial & m_maxM & m_bra_M & m_schedule_type_backward & m_schedule_type_default &m_integral_disk_storage_thresh;
     ar & n_twodot_noise & m_twodot_noise & m_twodot_gamma & m_guessState & m_useSharedMemory ;
     ar & m_calc_ri_4pdm & m_store_ripdm_readable & m_nevpt2 & m_conventional_nevpt2 & m_kept_nevpt2_states & NevPrint;
   }
@@ -524,6 +526,8 @@ class Input {
   bool &spatpdm_disk_dump() {return m_spatpdm_disk_dump;}
   const bool &pdm_unsorted() const {return m_pdm_unsorted;}
   bool &pdm_unsorted(){return m_pdm_unsorted;}
+  const std::vector<int> &specificpdm() const {return m_specificpdm;}
+  std::vector<int> &specificpdm() {return m_specificpdm;}
   const bool &store_nonredundant_pdm() const { return m_store_nonredundant_pdm;}
   bool &store_nonredundant_pdm() { return m_store_nonredundant_pdm;}
   int slater_size() const {return m_norbs;}
@@ -533,6 +537,8 @@ class Input {
   const bool &npdm_intermediate() const { return m_npdm_intermediate; }
   bool &npdm_multinode() { return m_npdm_multinode; }
   const bool &npdm_multinode() const { return m_npdm_multinode; }
+  int bra_M() const {return m_bra_M;}
+
 };
 }
 #endif

linear.C

@@ -132,7 +132,7 @@ void SpinAdapted::Linear::block_davidson(vector<StackWavefunction>& b, DiagonalM
   int sigmasize=0, bsize= currentRoot == -1 ? dmrginp.nroots() : 1;
   int converged_roots = 0;
   int maxiter = h_diag.Ncols() - lowerStates.size();
-  while(true)
+  while(true && iter < 100)
   {
     //p3out << "\t\t\t Davidson Iteration :: " << iter << endl;
 

modules/ResponseTheory/sweepResponse.C

@@ -270,7 +270,7 @@ void SpinAdapted::SweepResponse::BlockAndDecimate (SweepParams &sweepParams, Sta
 			     sweepParams.set_lowest_error(), rotatematrix, 
 			     sweepParams.get_keep_states(), sweepParams.get_keep_qstates(), 
 			     sweepParams.get_davidson_tol(), big, sweepParams.get_guesstype(), 
-			     0.0, 0.0, //noise 
+			     sweepParams.get_noise(), sweepParams.get_additional_noise(), //noise 
 			     sweepParams.get_onedot(), system, systemDot, environment, 
 			     dot_with_sys, useSlater, sweepParams.get_sweep_iter(), targetState, 
 			     lowerStates, &bratracedMatrix);
@@ -279,12 +279,12 @@ void SpinAdapted::SweepResponse::BlockAndDecimate (SweepParams &sweepParams, Sta
 
   p1out <<"\t\t\t Performing Renormalization "<<endl;
 
-  rotatematrix.resize(0);
+  //rotatematrix.resize(0);
   pout << "**** STACK MEMORY REMAINING before renormalization***** "<<1.0*(Stackmem[0].size-Stackmem[0].memused)*sizeof(double)/1.e9<<" GB"<<endl;
 
   if(mpigetrank() == 0) {
-    ScaleAdd(sweepParams.get_noise()*(max(1.e-5, trace(branoiseMatrix))), branoiseMatrix, bratracedMatrix);
-    sweepParams.set_lowest_error() = makeRotateMatrix(bratracedMatrix, rotatematrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
+    //ScaleAdd(sweepParams.get_noise()*(max(1.e-5, trace(branoiseMatrix))), branoiseMatrix, bratracedMatrix);
+    //sweepParams.set_lowest_error() = makeRotateMatrix(bratracedMatrix, rotatematrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
     bratracedMatrix.deallocate();
   }
   branoiseMatrix.deallocate();

modules/npdm/npdm.C

@@ -143,8 +143,11 @@ void npdm_block_and_decimate( Npdm_driver_base& npdm_driver, SweepParams &sweepP
     tracedMatrixB.makedensitymatrix(solution[1], big, 1.0);
     rotateMatrixB.clear();
     if (!mpigetrank()){
-      double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
-      error = makeRotateMatrix(tracedMatrixB, rotateMatrixB, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
+      double error = makeRotateMatrix(tracedMatrixB, rotateMatrixB, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
+      if (dmrginp.bra_M() == 0)
+        error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
+      else
+        error = makeRotateMatrix(tracedMatrix, rotateMatrix, dmrginp.bra_M(), sweepParams.get_keep_qstates());
     }
     tracedMatrixB.deallocate();
     tracedMatrix.deallocate();
@@ -410,6 +413,38 @@ void npdm(NpdmOrder npdm_order, bool transitionpdm)
     }
   }
 
+  if (dmrginp.specificpdm().size()!=0)
+  {
+    Timer timer;
+    dmrginp.set_fullrestart() = true;
+    sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy;
+	  dmrginp.npdm_generate() = true;
+    if (dmrginp.specificpdm().size()==1)
+      SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[0]); //this will generate the cd operators                               
+    else if (dmrginp.specificpdm().size()==2)
+      SweepGenblock::do_one(sweepParams, false, !direction, false, 0, dmrginp.specificpdm()[0], dmrginp.specificpdm()[1]); //this will generate the cd operators                               
+    else 
+      abort();
+		dmrginp.npdm_generate() = false;
+    double ecpu = timer.elapsedcputime();
+    double ewall=timer.elapsedwalltime();
+    p3out << "\t\t\t NPDM SweepGenblock time " << ewall << " " << ecpu << endl;
+    dmrginp.set_fullrestart() = false;
+
+    sweepParams = sweep_copy; direction = direction_copy; restartsize = restartsize_copy;
+    Timer timerX;
+    npdm_driver->clear();
+    if (dmrginp.specificpdm().size()==1)
+      npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[0]);
+    else if (dmrginp.specificpdm().size()==2)
+      npdm_do_one_sweep(*npdm_driver, sweepParams, false, direction, false, 0,dmrginp.specificpdm()[0],dmrginp.specificpdm()[1]);
+    else
+      abort();
+    ecpu = timerX.elapsedcputime();ewall=timerX.elapsedwalltime();
+    p3out << "\t\t\t NPDM sweep time " << ewall << " " << ecpu << endl;
+    return;
+  }
+
 
 
   if(dmrginp.transition_diff_irrep()){

set_stackspinblock_components.C

@@ -369,6 +369,9 @@ void StackSpinBlock::default_op_components(bool direct, bool haveNormops, bool h
 	ops[CRE_DES_DESCOMP] = make_new_stackop(CRE_DES_DESCOMP, true);
     }
 
+    if(dmrginp.do_npdm_ops() && (dmrginp.calc_type() == RESTART_ONEPDM || dmrginp.calc_type() == ONEPDM))
+      return;
+
     //for hubbard model if we want to calculate twopdm we still need cd operators
     if (dmrginp.hamiltonian() != HUBBARD || dmrginp.do_npdm_ops()) {
       if (haveNormops || dmrginp.do_npdm_ops()) {
@@ -437,6 +440,10 @@ void StackSpinBlock::default_op_components(bool direct, bool haveNormops, bool h
       if (!dmrginp.do_npdm_ops())
 	ops[CRE_DES_DESCOMP] = make_new_stackop(CRE_DES_DESCOMP, false);
     }
+
+    //Only need one index operators for one pdm calculation except for single dot.
+    if(dmrginp.do_npdm_ops() && (dmrginp.calc_type() == RESTART_ONEPDM || dmrginp.calc_type() == ONEPDM))
+      return;
 
     //for hubbard model if we want to calculate twopdm we still need cd operators
     if (dmrginp.hamiltonian() != HUBBARD || dmrginp.do_npdm_ops()) {

wrapper.C

@@ -245,12 +245,12 @@ void test(char* infile)
   std::vector< std::vector<double> > Overlap(nstates, std::vector<double>(nstates, 0.0));
 
 
-  for (int i=1; i<nstates; i++) {
+  for (int i=0; i<nstates; i++) {
     if(mpigetrank() == 0)
       printf("starting row : %i\n", i);
-    for (int j=0; j<1; j++) {
+    for (int j=0; j<=i; j++) {
       double h=0,o=0;
-      calcHamiltonianAndOverlap(states[i], states[j], h, o, i==j, 0);
+      calcHamiltonianAndOverlap(states[i], states[j], h, o);
       ham[i][j] = h; ham[j][i] = h;
       Overlap[i][j] = o; Overlap[j][i] = o;
       if (mpigetrank() == 0)