C++ SpinBlock::get_sites方法代码示例

  template<> void Op_component<CreCreComp>::build_iterators(SpinBlock& b)
    {
      if (b.get_sites().size () == 0) return; // blank construction (used in unset_initialised() Block copy construction, for use with STL)
      const double screen_tol = dmrginp.twoindex_screen_tol();
      vector< pair<int, int> > screened_dd_ix = (dmrginp.hamiltonian() == BCS) ?
        screened_dd_indices(b.get_complementary_sites(), b.get_sites(), *b.get_twoInt(), v_cc, v_cccc, v_cccd, screen_tol) :
        screened_dd_indices(b.get_complementary_sites(), b.get_sites(), *b.get_twoInt(), screen_tol);
      m_op.set_pair_indices(screened_dd_ix, dmrginp.last_site());      
      
      std::vector<int> orbs(2);
      for (int i = 0; i < m_op.local_nnz(); ++i)
	{
	  orbs = m_op.unmap_local_index(i);
	  std::vector<boost::shared_ptr<CreCreComp> >& vec = m_op.get_local_element(i);
	  SpinQuantum spin1 = getSpinQuantum(orbs[0]);
	  SpinQuantum spin2 = getSpinQuantum(orbs[1]);
	  std::vector<SpinQuantum> spinvec = spin1+spin2;
	  vec.resize(spinvec.size());
	  for (int j=0; j<spinvec.size(); j++) {
	    vec[j]=boost::shared_ptr<CreCreComp>(new CreCreComp);
	    SparseMatrix& op = *vec[j];
	    op.set_orbs() = orbs;
	    op.set_initialised() = true;
	    op.set_fermion() = false;
        
	      op.set_deltaQuantum(1, spinvec[j]);

	  }
	}
      
    }

void build_3index_ops( const opTypes& optype, SpinBlock& big, 
                       const opTypes& lhsType1, const opTypes& lhsType2,
                       const opTypes& rhsType1, const opTypes& rhsType2,
                       const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  // 3-index output file
//pout << "build_3index_op, ofs =" <<  big.get_op_array(optype).get_filename() << endl;
  std::ofstream ofs;
  if ( ! dmrginp.do_npdm_in_core() ) ofs.open( big.get_op_array(optype).get_filename().c_str(), std::ios::binary );

  SpinBlock* sysBlock = big.get_leftBlock();
  SpinBlock* dotBlock = big.get_rightBlock();

  // All 3 orbitals on sys or dot block
  do_3index_tensor_trace( optype, big, sysBlock, ofs, rotateMatrix, stateinfo );
  do_3index_tensor_trace( optype, big, dotBlock, ofs, rotateMatrix, stateinfo );

  bool forwards = ! ( sysBlock->get_sites().at(0) > dotBlock->get_sites().at(0) );

  // 2,1 partitioning
  if ( forwards ) {
    do_3index_1_2_tensor_products( forwards, optype, lhsType1, rhsType2, big, dotBlock, sysBlock, ofs, rotateMatrix, stateinfo );
    do_3index_2_1_tensor_products( forwards, optype, lhsType2, rhsType1, big, dotBlock, sysBlock, ofs, rotateMatrix, stateinfo );
  } else {
    do_3index_1_2_tensor_products( forwards, optype, lhsType1, rhsType2, big, sysBlock, dotBlock, ofs, rotateMatrix, stateinfo );
    do_3index_2_1_tensor_products( forwards, optype, lhsType2, rhsType1, big, sysBlock, dotBlock, ofs, rotateMatrix, stateinfo );
  }

  if ( ofs.is_open() ) ofs.close();

}

double SweepOnepdm::do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, const bool &restart, const int &restartSize)
{

  SpinBlock system;
  const int nroots = dmrginp.nroots();
  std::vector<double> finalEnergy(nroots,0.);
  std::vector<double> finalEnergy_spins(nroots,0.);
  double finalError = 0.;

  Matrix onepdm(2*dmrginp.last_site(), 2*dmrginp.last_site());onepdm=0.0;
  for (int i=0; i<nroots; i++)
    for (int j=0; j<=i; j++)      
      save_onepdm_binary(onepdm, i ,j);

  sweepParams.set_sweep_parameters();
  // a new renormalisation sweep routine
  pout << ((forward) ? "\t\t\t Starting renormalisation sweep in forwards direction" : "\t\t\t Starting renormalisation sweep in backwards direction") << endl;
  pout << "\t\t\t ============================================================================ " << endl;
  
  InitBlocks::InitStartingBlock (system,forward, sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp);

  sweepParams.set_block_iter() = 0;
 
  pout << "\t\t\t Starting block is :: " << endl << system << endl;

  SpinBlock::store (forward, system.get_sites(), system); // if restart, just restoring an existing block --
  sweepParams.savestate(forward, system.get_sites().size());
  bool dot_with_sys = true;

  sweepParams.set_guesstype() = TRANSPOSE;
	  
  SpinBlock newSystem;
  BlockAndDecimate (sweepParams, system, newSystem, warmUp, dot_with_sys);
  pout.precision(12);
  pout << "\t\t\t The lowest sweep energy : "<< sweepParams.get_lowest_energy()[0]+dmrginp.get_coreenergy()<<endl;
  pout << "\t\t\t ============================================================================ " << endl;

  for (int i=0; i<nroots; i++)
    for (int j=0; j<=i; j++) {
      load_onepdm_binary(onepdm, i ,j);
      accumulate_onepdm(onepdm);
      save_onepdm_spatial_text(onepdm, i ,j);
      save_onepdm_text(onepdm, i ,j);
      save_onepdm_spatial_binary(onepdm, i ,j);
    }
  return sweepParams.get_lowest_energy()[0];
}

//-------------------------------------------------------------------------------------------------------------------------------------------------------------
//  (Cre,Cre,Cre,Cre)
//-------------------------------------------------------------------------------------------------------------------------------------------------------------
void SpinAdapted::CreCreCreCre::build(const SpinBlock& b) { 
  dmrginp.makeopsT -> start();
  built = true;
  allocate(b.get_braStateInfo(), b.get_ketStateInfo());

  const int i = get_orbs()[0];
  const int j = get_orbs()[1];
  const int k = get_orbs()[2];
  const int l = get_orbs()[3];
  SpinBlock* leftBlock = b.get_leftBlock();
  SpinBlock* rightBlock = b.get_rightBlock();

  if (leftBlock->get_op_array(CRE_CRE_CRE_CRE).has(i,j,k,l))
  {      
    const boost::shared_ptr<SparseMatrix>& op = leftBlock->get_op_rep(CRE_CRE_CRE_CRE, quantum_ladder, i,j,k,l);
    if (rightBlock->get_sites().size() == 0) 
      SpinAdapted::operatorfunctions::TensorTrace(leftBlock, *op, &b, &(b.get_stateInfo()), *this);
    dmrginp.makeopsT -> stop();
    return;
  }
  assert(false && "Only build CRECRECRECRE in the starting block when spin-embeding is used");
}

void SweepOnepdm::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys)
{
  //mcheck("at the start of block and decimate");
  // figure out if we are going forward or backwards
  dmrginp.guessgenT -> start();
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  SpinBlock envDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  if (forward)
  {
    systemDotStart = *system.get_sites().rbegin () + 1;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = system.get_sites() [0] - 1;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  systemDot = SpinBlock(systemDotStart, systemDotEnd);

  SpinBlock environment, environmentDot, newEnvironment;

  int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();

  system.addAdditionalCompOps();
  InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, sweepParams.get_sys_add(), dmrginp.direct(), DISTRIBUTED_STORAGE, true, true);
  
  InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot,
				      sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(),
				      sweepParams.get_onedot(), nexact, useSlater, true, true, true);
  SpinBlock big;
  newSystem.set_loopblock(true);
  system.set_loopblock(false);
  newEnvironment.set_loopblock(false);
  InitBlocks::InitBigBlock(newSystem, newEnvironment, big); 

  const int nroots = dmrginp.nroots();
  std::vector<Wavefunction> solutions(nroots);

  
  for(int i=0;i<nroots;++i)
    {
      StateInfo newInfo;
      solutions[i].LoadWavefunctionInfo (newInfo, newSystem.get_sites(), i);
    }
  
#ifndef SERIAL
  mpi::communicator world;
  mpi::broadcast(world,solutions,0);
#endif

#ifdef SERIAL
  const int numprocs = 1;
#endif
#ifndef SERIAL
  const int numprocs = world.size();
#endif


  compute_onepdm(solutions, system, systemDot, newSystem, newEnvironment, big, numprocs);

}

void SpinAdapted::InitBlocks::InitNewOverlapEnvironmentBlock(SpinBlock &environment, SpinBlock& environmentDot, SpinBlock &newEnvironment, 
							     const SpinBlock &system, SpinBlock &systemDot, int leftState, int rightState,
							     const int &sys_add, const int &env_add, const bool &forward, int integralIndex,
							     const bool &onedot, const bool& dot_with_sys, int constraint)
{
  // now initialise environment Dot
  int systemDotStart, systemDotEnd, environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;
  int systemDotSize = sys_add - 1;
  int environmentDotSize = env_add - 1;
  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
    environmentDotStart = systemDotEnd + 1;
    environmentDotEnd = environmentDotStart + environmentDotSize;
    environmentStart = environmentDotEnd + 1;
    environmentEnd = dmrginp.spinAdapted() ? dmrginp.last_site() - 1 : dmrginp.last_site()/2 - 1;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
    environmentDotStart = systemDotEnd - 1;
    environmentDotEnd = environmentDotStart - environmentDotSize;
    environmentStart = environmentDotEnd - 1;
    environmentEnd = 0;
  }

  std::vector<int> environmentSites;
  environmentSites.resize(abs(environmentEnd - environmentStart) + 1);
  for (int i = 0; i < abs(environmentEnd - environmentStart) + 1; ++i) *(environmentSites.begin () + i) = min(environmentStart,environmentEnd) + i;


  p2out << "\t\t\t Restoring block of size " << environmentSites.size () << " from previous iteration" << endl;

  if(dot_with_sys && onedot) {
    newEnvironment.set_integralIndex() = integralIndex;
    SpinBlock::restore (!forward, environmentSites, newEnvironment, leftState, rightState);
  }
  else {
    environment.set_integralIndex() = integralIndex;
    SpinBlock::restore (!forward, environmentSites, environment, leftState, rightState);
  }
  if (dmrginp.outputlevel() > 0)
    mcheck("");

  // now initialise newEnvironment
  if (!dot_with_sys || !onedot)
  {
    newEnvironment.set_integralIndex() = integralIndex;
    newEnvironment.initialise_op_array(OVERLAP, false);
    //newEnvironment.set_op_array(OVERLAP) = boost::shared_ptr<Op_component<Overlap> >(new Op_component<Overlap>(false));
    newEnvironment.setstoragetype(DISTRIBUTED_STORAGE);
      
    newEnvironment.BuildSumBlock (constraint, environment, environmentDot);
    p2out << "\t\t\t Environment block " << endl << environment << endl;
    environment.printOperatorSummary();
    p2out << "\t\t\t NewEnvironment block " << endl << newEnvironment << endl;
    newEnvironment.printOperatorSummary();
  }
  else {
    p2out << "\t\t\t Environment block " << endl << newEnvironment << endl;
    newEnvironment.printOperatorSummary();
  }

}

void SpinAdapted::InitBlocks::InitNewEnvironmentBlock(SpinBlock &environment, SpinBlock& environmentDot, SpinBlock &newEnvironment, 
						      const SpinBlock &system, SpinBlock &systemDot, int leftState, int rightState,
						      const int &sys_add, const int &env_add, const bool &forward, const bool &direct, 
						      const bool &onedot, const bool &nexact, const bool &useSlater, int integralIndex, 
						      bool haveNormops, bool haveCompops, const bool& dot_with_sys, int constraint, const std::vector<SpinQuantum>& braquanta, const std::vector<SpinQuantum>& ketquanta) {
  // now initialise environment Dot
  int systemDotStart, systemDotEnd, environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;
  int systemDotSize = sys_add - 1;
  int environmentDotSize = env_add - 1;
  if (forward) {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
    environmentDotStart = systemDotEnd + 1;
    environmentDotEnd = environmentDotStart + environmentDotSize;
    environmentStart = environmentDotEnd + 1;
    environmentEnd = dmrginp.spinAdapted() ? dmrginp.last_site() - 1 : dmrginp.last_site()/2 - 1;
  } else {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
    environmentDotStart = systemDotEnd - 1;
    environmentDotEnd = environmentDotStart - environmentDotSize;
    environmentStart = environmentDotEnd - 1;
    environmentEnd = 0;
  }

  std::vector<int> environmentSites;
  environmentSites.resize(abs(environmentEnd - environmentStart) + 1);
  for (int i = 0; i < abs(environmentEnd - environmentStart) + 1; ++i) *(environmentSites.begin () + i) = min(environmentStart,environmentEnd) + i;


  // now initialise environment
  if (useSlater) { // for FCI
    StateInfo system_stateinfo = system.get_stateInfo();
    StateInfo sysdot_stateinfo = systemDot.get_stateInfo();
    StateInfo tmp;
    TensorProduct (system_stateinfo, sysdot_stateinfo, tmp, NO_PARTICLE_SPIN_NUMBER_CONSTRAINT);
    // tmp has the system+dot quantum numbers
    tmp.CollectQuanta ();
    
    // exact environment
    if (dmrginp.do_fci() || environmentSites.size() == nexact) {
      if ((!dot_with_sys && onedot) || !onedot) { // environment has dot
	environment.set_integralIndex() = integralIndex;
	environment.default_op_components(!forward, leftState==rightState);
	environment.setstoragetype(DISTRIBUTED_STORAGE);
	environment.BuildTensorProductBlock(environmentSites); // exact block
	SpinBlock::store (true, environmentSites, environment, leftState, rightState);	
      } 
      else { // environment has no dot, so newEnv = Env
	newEnvironment.set_integralIndex() = integralIndex;
	newEnvironment.default_op_components(!forward, leftState==rightState);
	newEnvironment.setstoragetype(DISTRIBUTED_STORAGE);
	newEnvironment.BuildTensorProductBlock(environmentSites);
	SpinBlock::store (true, environmentSites, newEnvironment, leftState, rightState);	
      }
    } else if (dmrginp.warmup() == LOCAL2 || dmrginp.warmup() == LOCAL3 || dmrginp.warmup() == LOCAL4) {
      int nactiveSites, ncoreSites;
      if (dmrginp.warmup() == LOCAL2) {
        nactiveSites = 1;
      } else if (dmrginp.warmup() == LOCAL3) {
        nactiveSites = 2;
      } else if (dmrginp.warmup() == LOCAL4) {
        nactiveSites = 3;
      }
      if (dot_with_sys && onedot) {
        nactiveSites += 1;
      }

      if (nactiveSites > environmentSites.size()) {
        nactiveSites = environmentSites.size();
      }
      ncoreSites = environmentSites.size() - nactiveSites;

      // figure out what sites are in the active and core sites
      int environmentActiveEnd = forward ? environmentStart + nactiveSites - 1 : environmentStart - nactiveSites + 1;
      int environmentCoreStart = forward ? environmentActiveEnd + 1 : environmentActiveEnd - 1;
      
      std::vector<int> activeSites(nactiveSites), coreSites(ncoreSites);
      for (int i = 0; i < nactiveSites; ++i) {
        activeSites[i] = min(environmentStart,environmentActiveEnd) + i;
      }
      for (int i = 0; i < ncoreSites; ++i) {
        coreSites[i] = min(environmentCoreStart,environmentEnd) + i;
      }

      SpinBlock environmentActive, environmentCore;
      environmentActive.nonactive_orb() = system.nonactive_orb();
      environmentCore.nonactive_orb() = system.nonactive_orb();
      if (coreSites.size() > 0) {
	environmentActive.set_integralIndex() = integralIndex;
	environmentCore.set_integralIndex() = integralIndex;
        environmentActive.default_op_components(!forward, leftState==rightState);
        environmentActive.setstoragetype(DISTRIBUTED_STORAGE);
        environmentCore.default_op_components(!forward, leftState==rightState);      
        environmentCore.setstoragetype(DISTRIBUTED_STORAGE);

        environmentActive.BuildTensorProductBlock(activeSites);
        environmentCore.BuildSingleSlaterBlock(coreSites);

        dmrginp.datatransfer -> start();
//.........这里部分代码省略.........

void SweepGenblock::do_one(SweepParams &sweepParams, const bool &forward, int stateA, int stateB)
{
  Timer sweeptimer;
  int integralIndex = 0;
  SpinBlock system;

  sweepParams.set_sweep_parameters();
  // a new renormalisation sweep routine
  pout << ((forward) ? "\t\t\t Starting renormalisation sweep in forwards direction" : "\t\t\t Starting renormalisation sweep in backwards direction") << endl;
  pout << "\t\t\t ============================================================================ " << endl;
  
  InitBlocks::InitStartingBlock (system,forward, stateA, stateB, sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), 0, false, false, integralIndex);

  sweepParams.set_block_iter() = 0;

  p2out << "\t\t\t Starting block is :: " << endl << system << endl;

  bool dot_with_sys = true;

  for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); )
    {
      pout << "\n\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl;
      pout << "\t\t\t ----------------------------" << endl;
      if (forward)
	{ p1out << "\t\t\t Current direction is :: Forwards " << endl; }
      else
	{ p1out << "\t\t\t Current direction is :: Backwards " << endl; }

  
      if (dmrginp.no_transform())
	      sweepParams.set_guesstype() = BASIC;
      else if ( sweepParams.get_block_iter() != 0) 
  	    sweepParams.set_guesstype() = TRANSFORM;
      else if ( sweepParams.get_block_iter() == 0 )
        sweepParams.set_guesstype() = TRANSPOSE;
      else
        sweepParams.set_guesstype() = BASIC;
      
      p1out << "\t\t\t Blocking and Decimating " << endl;
	  
      SpinBlock newSystem;

      BlockAndDecimate (sweepParams, system, newSystem, false, dot_with_sys, stateA, stateB);

      system = newSystem;

      SpinBlock::store(forward, system.get_sites(), system, stateA, stateB);

      //system size is going to be less than environment size
      if (forward && system.get_complementary_sites()[0] >= dmrginp.last_site()/2)
	dot_with_sys = false;
      if (!forward && system.get_sites()[0]-1 < dmrginp.last_site()/2)
	dot_with_sys = false;

      ++sweepParams.set_block_iter();
    }
  pout << "\t\t\t Finished Generate-Blocks Sweep. " << endl;
  pout << "\t\t\t ============================================================================ " << endl;

  // update the static number of iterations

  ++sweepParams.set_sweep_iter();

  ecpu = sweeptimer.elapsedcputime(); ewall = sweeptimer.elapsedwalltime();
  pout << "\t\t\t Elapsed Sweep CPU  Time (seconds): " << setprecision(3) << ecpu << endl;
  pout << "\t\t\t Elapsed Sweep Wall Time (seconds): " << setprecision(3) << ewall << endl;

}

double SweepTwopdm::do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, const bool &restart, const int &restartSize, int state)
{
  Timer sweeptimer;
  int integralIndex = 0;
  if (dmrginp.hamiltonian() == BCS) {
    pout << "Two PDM with BCS calculations is not implemented" << endl;
    exit(0);
  }
  pout.precision(12);
  SpinBlock system;
  const int nroots = dmrginp.nroots();
  std::vector<double> finalEnergy(nroots,0.);
  std::vector<double> finalEnergy_spins(nroots,0.);
  double finalError = 0.;

  sweepParams.set_sweep_parameters();
  // a new renormalisation sweep routine
  pout << ((forward) ? "\t\t\t Starting renormalisation sweep in forwards direction" : "\t\t\t Starting renormalisation sweep in backwards direction") << endl;
  pout << "\t\t\t ============================================================================ " << endl;
  
  InitBlocks::InitStartingBlock (system,forward, sweepParams.current_root(), sweepParams.current_root(), sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp, integralIndex);
  if(!restart)
    sweepParams.set_block_iter() = 0;
 
  pout << "\t\t\t Starting block is :: " << endl << system << endl;
  if (!restart) 
    SpinBlock::store (forward, system.get_sites(), system, sweepParams.current_root(), sweepParams.current_root()); // if restart, just restoring an existing block --
  sweepParams.savestate(forward, system.get_sites().size());
  bool dot_with_sys = true;

  array_4d<double> twopdm(2*dmrginp.last_site(), 2*dmrginp.last_site(), 2*dmrginp.last_site(), 2*dmrginp.last_site());
  twopdm.Clear();

  save_twopdm_binary(twopdm, state, state); 


  for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); )
    {
      pout << "\n\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl;
      pout << "\t\t\t ----------------------------" << endl;
      if (forward)
	p1out << "\t\t\t Current direction is :: Forwards " << endl;
      else
	p1out << "\t\t\t Current direction is :: Backwards " << endl;

      //if (SHOW_MORE) pout << "system block" << endl << system << endl;
  
      if (dmrginp.no_transform())
	      sweepParams.set_guesstype() = BASIC;
      else if (!warmUp && sweepParams.get_block_iter() != 0) 
  	    sweepParams.set_guesstype() = TRANSFORM;
      else if (!warmUp && sweepParams.get_block_iter() == 0 && 
                ((dmrginp.algorithm_method() == TWODOT_TO_ONEDOT && dmrginp.twodot_to_onedot_iter() != sweepParams.get_sweep_iter()) ||
                  dmrginp.algorithm_method() != TWODOT_TO_ONEDOT))
        sweepParams.set_guesstype() = TRANSPOSE;
      else
        sweepParams.set_guesstype() = BASIC;
      
      p1out << "\t\t\t Blocking and Decimating " << endl;
	  
      SpinBlock newSystem;

      BlockAndDecimate (sweepParams, system, newSystem, warmUp, dot_with_sys, state);

      for(int j=0;j<nroots;++j)
        pout << "\t\t\t Total block energy for State [ " << j << 
	  " ] with " << sweepParams.get_keep_states()<<" :: " << sweepParams.get_lowest_energy()[j] <<endl;              

      finalEnergy_spins = ((sweepParams.get_lowest_energy()[0] < finalEnergy[0]) ? sweepParams.get_lowest_energy_spins() : finalEnergy_spins);
      finalEnergy = ((sweepParams.get_lowest_energy()[0] < finalEnergy[0]) ? sweepParams.get_lowest_energy() : finalEnergy);
      finalError = max(sweepParams.get_lowest_error(),finalError);

      system = newSystem;

      pout << system<<endl;
      
      SpinBlock::store (forward, system.get_sites(), system, sweepParams.current_root(), sweepParams.current_root());	 	

      p1out << "\t\t\t saving state " << system.get_sites().size() << endl;
      ++sweepParams.set_block_iter();
      //sweepParams.savestate(forward, system.get_sites().size());
    }
  //for(int j=0;j<nroots;++j)
  {int j = state;
    pout << "\t\t\t Finished Sweep with " << sweepParams.get_keep_states() << " states and sweep energy for State [ " << j 
	 << " ] with Spin [ " << dmrginp.molecule_quantum().get_s()  << " ] :: " << finalEnergy[j] << endl;
  }
  pout << "\t\t\t Largest Error for Sweep with " << sweepParams.get_keep_states() << " states is " << finalError << endl;
  pout << "\t\t\t ============================================================================ " << endl;

  int i = state, j = state;
  //for (int j=0; j<=i; j++) {
  load_twopdm_binary(twopdm, i, j); 
  //calcenergy(twopdm, i);
  save_twopdm_text(twopdm, i, j);
  save_spatial_twopdm_text(twopdm, i, j);
  save_spatial_twopdm_binary(twopdm, i, j);
  

  // update the static number of iterations
//.........这里部分代码省略.........

void SpinAdapted::mps_nevpt::type1::Startup(const SweepParams &sweepParams, const bool &forward, perturber& pb, int baseState) {

#ifndef SERIAL
  mpi::communicator world;
#endif
  assert(forward);
  SpinBlock system;
  system.nonactive_orb() =pb.orb();
  bool restart=false, warmUp = false;
  int forward_starting_size=1, backward_starting_size=0, restartSize =0;
  InitBlocks::InitStartingBlock(system, forward, pb.wavenumber(), baseState, forward_starting_size, backward_starting_size, restartSize, restart, warmUp, 0,pb.braquanta, pb.ketquanta); 

  SpinBlock::store (forward, system.get_sites(), system, pb.wavenumber(), baseState); // if restart, just restoring an existing block --

  for (int i=0; i<mps_nevpt::sweepIters; i++) {
    SpinBlock newSystem;
    SpinBlock dotSystem(i+1,i+1,pb.orb(),false);

    system.addAdditionalCompOps();
    //newSystem.default_op_components(true, system, dotSystem, true, true, false);
    newSystem.perturb_op_components(false, system, dotSystem, pb);
    newSystem.setstoragetype(DISTRIBUTED_STORAGE);
    newSystem.BuildSumBlock(LessThanQ, system, dotSystem, pb.braquanta, pb.ketquanta);
    newSystem.printOperatorSummary();
    //SpinBlock Environment, big;
    //SpinBlock::restore (!forward, newSystem.get_complementary_sites() , Environment, baseState, baseState);
    //TODO
    //SpinBlock::restore (!forward, newSystem.get_complementary_sites() , Environment,sweepParams.current_root(),sweepParams.current_root());

    //big.BuildSumBlock(PARTICLE_SPIN_NUMBER_CONSTRAINT, newSystem, Environment, pb.braquanta, pb.ketquanta);

    //StateInfo envStateInfo;
    StateInfo ketStateInfo;
    StateInfo braStateInfo;
    StateInfo halfbraStateInfo;// It has the same left and right StateInfo as braStateInfo. However, its total quanta is pb.ketquanta.
    // It is used to project solution into to braStateInfo.

    std::vector<Wavefunction> solution; solution.resize(1);
    std::vector<Wavefunction> outputState; outputState.resize(1);
    std::vector<Wavefunction> solutionprojector; solutionprojector.resize(1);
    solution[0].LoadWavefunctionInfo(ketStateInfo, newSystem.get_sites(), baseState);
    #ifndef SERIAL
      broadcast(world, ketStateInfo, 0);
      broadcast(world, solution, 0);
    #endif
    outputState[0].AllowQuantaFor(newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta);
    outputState[0].set_onedot(solution[0].get_onedot());
    outputState[0].Clear();
    solutionprojector[0].AllowQuantaFor(newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.ketquanta);
    solutionprojector[0].set_onedot(solution[0].get_onedot());
    solutionprojector[0].Clear();
    //TensorProduct (newSystem.get_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta[0], EqualQ, braStateInfo);
    //TODO
    //TensorProduct do not support const StateInfo&
    TensorProduct (newSystem.set_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.braquanta[0], EqualQ, braStateInfo);
    TensorProduct (newSystem.set_braStateInfo(), *(ketStateInfo.rightStateInfo), pb.ketquanta[0], EqualQ, halfbraStateInfo);

    //StateInfo::restore(forward, environmentsites, envStateInfo, baseState);

    //DiagonalMatrix e;
    //if(i == 0)
    //  GuessWave::guess_wavefunctions(solution, e, big, TRANSPOSE, true, true, 0.0, baseState); 
    //else
    //  GuessWave::guess_wavefunctions(solution, e, big, TRANSFORM, true, true, 0.0, baseState); 


    //SpinAdapted::operatorfunctions::Product(&newSystem, ccd, solution[0], &ketStateInfo, stateb.getw(), temp, SpinQuantum(0, SpinSpace(0), IrrepSpace(0)), true, 1.0);

    

    boost::shared_ptr<SparseMatrix> O;
    if (pb.type() == TwoPerturbType::Va)
      O = newSystem.get_op_array(CDD_SUM).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
    if (pb.type() == TwoPerturbType::Vi)
      O = newSystem.get_op_array(CCD_SUM).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
    boost::shared_ptr<SparseMatrix> overlap = newSystem.get_op_array(OVERLAP).get_local_element(0)[0]->getworkingrepresentation(&newSystem);
    SpinAdapted::operatorfunctions::TensorMultiply(*O, &braStateInfo, &ketStateInfo , solution[0], outputState[0], pb.delta, true, 1.0);
    SpinAdapted::operatorfunctions::TensorMultiply(*overlap, &halfbraStateInfo, &ketStateInfo , solution[0], solutionprojector[0], overlap->get_deltaQuantum(0), true, 1.0);
    DensityMatrix bratracedMatrix(newSystem.get_braStateInfo());
    bratracedMatrix.allocate(newSystem.get_braStateInfo());
    double norm = DotProduct(outputState[0], outputState[0]);
    if(norm > NUMERICAL_ZERO)
      SpinAdapted::operatorfunctions::MultiplyProduct(outputState[0], Transpose(const_cast<Wavefunction&> (outputState[0])), bratracedMatrix, 0.5/norm);
    SpinAdapted::operatorfunctions::MultiplyProduct(solutionprojector[0], Transpose(const_cast<Wavefunction&> (solutionprojector[0])), bratracedMatrix, 0.5);
    std::vector<Matrix> brarotateMatrix, ketrotateMatrix;
    LoadRotationMatrix (newSystem.get_sites(), ketrotateMatrix, baseState);
    double error;
    if (!mpigetrank())
      error = makeRotateMatrix(bratracedMatrix, brarotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
    #ifndef SERIAL
      broadcast(world, ketrotateMatrix, 0);
      broadcast(world, brarotateMatrix, 0);
    #endif

    SaveRotationMatrix (newSystem.get_sites(), brarotateMatrix, pb.wavenumber());
    newSystem.transform_operators(brarotateMatrix,ketrotateMatrix);
    SpinBlock::store (forward, newSystem.get_sites(), newSystem, pb.wavenumber(), baseState); // if restart, just restoring an existing block --
    system=newSystem;
  }
  //TODO
//.........这里部分代码省略.........

double SpinAdapted::mps_nevpt::type1::do_one(SweepParams &sweepParams, const bool &warmUp, const bool &forward, const bool &restart, const int &restartSize, perturber& pb, int baseState)
{
  int integralIndex = 0;
  SpinBlock system;
  system.nonactive_orb() = pb.orb();
  const int nroots = dmrginp.nroots(sweepParams.get_sweep_iter());

  std::vector<double> finalEnergy(nroots,-1.0e10);
  std::vector<double> finalEnergy_spins(nroots,0.);
  double finalError = 0.;

  sweepParams.set_sweep_parameters();
  // a new renormalisation sweep routine
  if (forward)
    if (dmrginp.outputlevel() > 0)
      pout << "\t\t\t Starting sweep "<< sweepParams.set_sweep_iter()<<" in forwards direction"<<endl;
  else
    if (dmrginp.outputlevel() > 0)
    {
      pout << "\t\t\t Starting sweep "<< sweepParams.set_sweep_iter()<<" in backwards direction" << endl;
      pout << "\t\t\t ============================================================================ " << endl;
    }

  InitBlocks::InitStartingBlock (system,forward, baseState, pb.wavenumber(), sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp, integralIndex, pb.braquanta, pb.ketquanta);
  if(!restart)
    sweepParams.set_block_iter() = 0;

 
  if (dmrginp.outputlevel() > 0)
    pout << "\t\t\t Starting block is :: " << endl << system << endl;

  SpinBlock::store (forward, system.get_sites(), system, pb.wavenumber(), baseState); // if restart, just restoring an existing block --
  sweepParams.savestate(forward, system.get_sites().size());
  bool dot_with_sys = true;
  vector<int> syssites = system.get_sites();

  if (restart)
  {
    if (forward && system.get_complementary_sites()[0] >= dmrginp.last_site()/2)
      dot_with_sys = false;
    if (!forward && system.get_sites()[0]-1 < dmrginp.last_site()/2)
      dot_with_sys = false;
  }
  if (dmrginp.outputlevel() > 0)
    mcheck("at the very start of sweep");  // just timer

  for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); ) // get_n_iters() returns the number of blocking iterations needed in one sweep
    {
      if (dmrginp.outputlevel() > 0)
      {
        pout << "\t\t\t Block Iteration :: " << sweepParams.get_block_iter() << endl;
        pout << "\t\t\t ----------------------------" << endl;
      }
      if (dmrginp.outputlevel() > 0) {
	    if (forward)
      {
	      pout << "\t\t\t Current direction is :: Forwards " << endl;
      }
	    else
      {
	      pout << "\t\t\t Current direction is :: Backwards " << endl;
      }
      }

      if (sweepParams.get_block_iter() != 0) 
	sweepParams.set_guesstype() = TRANSFORM;
      else
        sweepParams.set_guesstype() = TRANSPOSE;


      
      if (dmrginp.outputlevel() > 0)
         pout << "\t\t\t Blocking and Decimating " << endl;
	  
      SpinBlock newSystem; // new system after blocking and decimating
      newSystem.nonactive_orb() = pb.orb();

      //Need to substitute by:
     // if (warmUp )
     //   Startup(sweepParams, system, newSystem, dot_with_sys, pb.wavenumber(), baseState);
     // else {
     //   BlockDecimateAndCompress (sweepParams, system, newSystem, false, dot_with_sys, pb.wavenumber(), baseState);
     // }
      
        BlockDecimateAndCompress (sweepParams, system, newSystem, warmUp, dot_with_sys,pb, baseState);
      //Need to substitute by?


      system = newSystem;
      if (dmrginp.outputlevel() > 0){
	    pout << system<<endl;
	    pout << system.get_braStateInfo()<<endl;
	    system.printOperatorSummary();
      }
      
      //system size is going to be less than environment size
      if (forward && system.get_complementary_sites()[0] >= dmrginp.last_site()/2)
	    dot_with_sys = false;
      if (!forward && system.get_sites()[0]-1 < dmrginp.last_site()/2)
	    dot_with_sys = false;
//.........这里部分代码省略.........

void SpinAdapted::mps_nevpt::type1::BlockDecimateAndCompress (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, perturber& pb, int baseState)
{
  int sweepiter = sweepParams.get_sweep_iter();
  if (dmrginp.outputlevel() > 0) {
    mcheck("at the start of block and decimate");
    pout << "\t\t\t dot with system "<<dot_with_sys<<endl;
    pout <<endl<< "\t\t\t Performing Blocking"<<endl;
  }
  // figure out if we are going forward or backwards
  dmrginp.guessgenT -> start();
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  SpinBlock environment, environmentDot, newEnvironment;
  SpinBlock big;
  environment.nonactive_orb() = pb.orb();
  newEnvironment.nonactive_orb() = pb.orb();
  int systemDotStart, systemDotEnd;
  int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  int environmentDotSize = sweepParams.get_env_add() -1;
  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
    environmentDotStart = systemDotEnd + 1;
    environmentDotEnd = environmentDotStart + environmentDotSize;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
    environmentDotStart = systemDotEnd - 1;
    environmentDotEnd = environmentDotStart - environmentDotSize;
  }
  systemDot = SpinBlock(systemDotStart, systemDotEnd, pb.orb());
  environmentDot = SpinBlock(environmentDotStart, environmentDotEnd, pb.orb());

  Sweep::makeSystemEnvironmentBigBlocks(system, systemDot, newSystem, environment, environmentDot, newEnvironment, big, sweepParams, dot_with_sys, useSlater, system.get_integralIndex(), pb.wavenumber(), baseState,pb.braquanta,pb.ketquanta);


  //analyse_operator_distribution(big);
  dmrginp.guessgenT -> stop();
  dmrginp.multiplierT -> start();
  std::vector<Matrix> rotatematrix;

  if (dmrginp.outputlevel() > 0)
    mcheck(""); 
  if (dmrginp.outputlevel() > 0) {
    if (!dot_with_sys && sweepParams.get_onedot())  { pout << "\t\t\t System  Block"<<system;    }
    else pout << "\t\t\t System  Block"<<newSystem;
    pout << "\t\t\t Environment Block"<<newEnvironment<<endl;
    pout << "\t\t\t Solving wavefunction "<<endl;
  }

  std::vector<Wavefunction> solution; solution.resize(1);
  std::vector<Wavefunction> outputState; outputState.resize(1);

  DiagonalMatrix e;


  //read the 0th wavefunction which we keep on the ket side because by default the ket stateinfo is used to initialize wavefunction
  //also when you use spinblock operators to multiply a state, it does so from the ket side i.e.  H|ket>
  //GuessWave::guess_wavefunctions(solution, e, big, sweepParams.set_guesstype(), sweepParams.get_onedot(), dot_with_sys, 0.0, baseState); 
  GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.set_guesstype(), sweepParams.get_onedot(), baseState, dot_with_sys, 0.0); 

#ifndef SERIAL
  mpi::communicator world;
  broadcast(world, solution, 0);
#endif
  
  outputState[0].AllowQuantaFor(big.get_leftBlock()->get_braStateInfo(), big.get_rightBlock()->get_braStateInfo(),pb.braquanta);
  outputState[0].set_onedot(sweepParams.get_onedot());
  outputState[0].Clear();
  if (pb.type() == TwoPerturbType::Va)
    big.multiplyCDD_sum(solution[0],&(outputState[0]),MAX_THRD);
  if (pb.type() == TwoPerturbType::Vi)
    big.multiplyCCD_sum(solution[0],&(outputState[0]),MAX_THRD);

  //davidson_f(solution[0], outputState[0]);
  SpinBlock newbig;

  if (sweepParams.get_onedot() && !dot_with_sys)
  {
    InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, baseState, pb.wavenumber(), systemDot.size(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, false, true,NO_PARTICLE_SPIN_NUMBER_CONSTRAINT,pb.braquanta,pb.ketquanta);
    InitBlocks::InitBigBlock(newSystem, environment, newbig,pb.braquanta,pb.ketquanta); 

    Wavefunction tempwave = outputState[0];
    GuessWave::onedot_shufflesysdot(big.get_braStateInfo(), newbig.get_braStateInfo(), outputState[0], tempwave);  
    outputState[0] = tempwave;

    tempwave = solution[0];
    GuessWave::onedot_shufflesysdot(big.get_ketStateInfo(), newbig.get_ketStateInfo(), solution[0], tempwave);  
    solution[0] = tempwave;

    big.get_rightBlock()->clear();
    big.clear();
  }
  else
    newbig = big;
  
//.........这里部分代码省略.........

void SweepOnepdm::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int state)
{
  //mcheck("at the start of block and decimate");
  // figure out if we are going forward or backwards
  dmrginp.guessgenT -> start();
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  SpinBlock envDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;

  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), true);

  SpinBlock environment, environmentDot, newEnvironment;
  int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();
  
  newSystem.set_integralIndex() = system.get_integralIndex();
  newSystem.default_op_components(dmrginp.direct(), system, systemDot, false, false, true);
  newSystem.erase(CRE_CRE_DESCOMP);
  newSystem.erase(CRE_CRE);
  newSystem.erase(HAM);
  newSystem.setstoragetype(DISTRIBUTED_STORAGE_FOR_ONEPDM);
  newSystem.BuildSumBlock (NO_PARTICLE_SPIN_NUMBER_CONSTRAINT, system, systemDot);
  if (dmrginp.outputlevel() > 0) {
    pout << "\t\t\t NewSystem block " << endl << newSystem << endl;
    newSystem.printOperatorSummary();
  }

  
  InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot, sweepParams.current_root(), sweepParams.current_root(),
				      sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(),
				      sweepParams.get_onedot(), nexact, useSlater, system.get_integralIndex(), false, false, true);
  SpinBlock big;
  newSystem.set_loopblock(true);
  system.set_loopblock(false);
  newEnvironment.set_loopblock(false);
  InitBlocks::InitBigBlock(newSystem, newEnvironment, big); 

  const int nroots = dmrginp.nroots();
  std::vector<Wavefunction> solution(1);

  DiagonalMatrix e;
  GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0); 

#ifndef SERIAL
  mpi::communicator world;
  mpi::broadcast(world, solution, 0);
#endif

  std::vector<Matrix> rotateMatrix;
  DensityMatrix tracedMatrix(newSystem.get_stateInfo());
  tracedMatrix.allocate(newSystem.get_stateInfo());
  tracedMatrix.makedensitymatrix(solution, big, std::vector<double>(1,1.0), 0.0, 0.0, false);
  rotateMatrix.clear();
  if (!mpigetrank())
    double error = makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
  

#ifndef SERIAL
  mpi::broadcast(world,rotateMatrix,0);
#endif
#ifdef SERIAL
  const int numprocs = 1;
#endif
#ifndef SERIAL
  const int numprocs = world.size();
#endif

  Matrix onepdm;
  load_onepdm_binary(onepdm, state ,state);
  Matrix pairmat;
  if (dmrginp.hamiltonian() == BCS)
    load_pairmat_binary(pairmat, state ,state);

  if (sweepParams.get_block_iter() == 0) {
    //this is inface a combination of  2_0_0, 1_1_0 and 0_2_0
    p2out << "\t\t\t compute 2_0_0"<<endl;
    compute_one_pdm_2_0_0(solution[0], solution[0], big, onepdm);
    if (dmrginp.hamiltonian() == BCS)
      compute_pair_2_0_0(solution[0], solution[0], big, pairmat);
    p2out << "\t\t\t compute 1_1_0"<<endl;
    compute_one_pdm_1_1_0(solution[0], solution[0], big, onepdm);
    if (dmrginp.hamiltonian() == BCS)    
      compute_pair_1_1_0(solution[0], solution[0], big, pairmat);
  }
//.........这里部分代码省略.........

void SweepGenblock::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int stateA, int stateB)
{
  if (dmrginp.outputlevel() > 0) 
    mcheck("at the start of block and decimate");
  p1out << "\t\t\t Performing Blocking"<<endl;
  dmrginp.guessgenT -> start();
  // figure out if we are going forward or backwards  
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  dmrginp.sysdotmake->start();
  systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), stateA==stateB);
  dmrginp.sysdotmake->stop();

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();

  dmrginp.guessgenT -> stop();
  dmrginp.datatransfer -> start();
  system.addAdditionalCompOps();
  dmrginp.datatransfer -> stop();
  dmrginp.initnewsystem->start();
  InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, stateA, stateB, sweepParams.get_sys_add(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, dot_with_sys, true);
  dmrginp.initnewsystem->stop();

  pout << "\t\t\t System  Block"<<newSystem;
  newSystem.printOperatorSummary();

  std::vector<Matrix> leftrotateMatrix, rightrotateMatrix;

  LoadRotationMatrix (newSystem.get_sites(), leftrotateMatrix, stateA);
  LoadRotationMatrix (newSystem.get_sites(), rightrotateMatrix, stateB);

#ifndef SERIAL
  mpi::communicator world;
  broadcast(world, leftrotateMatrix, 0);
  broadcast(world, rightrotateMatrix, 0);
#endif

  p1out <<"\t\t\t Performing Renormalization "<<endl<<endl;

  dmrginp.operrotT->start();
  if (stateB == stateA)
    newSystem.transform_operators(leftrotateMatrix);
  else
    newSystem.transform_operators(leftrotateMatrix, rightrotateMatrix);
  dmrginp.operrotT->stop();


  if (dmrginp.outputlevel() > 0) 
    //mcheck("after rotation and transformation of block");
  p2out <<newSystem<<endl;
  newSystem.printOperatorSummary();
  //mcheck("After renorm transform");

  p2out << *dmrginp.guessgenT<<" "<<*dmrginp.multiplierT<<" "<<*dmrginp.operrotT<< "  "<<globaltimer.totalwalltime()<<" timer "<<endl;
  p2out << *dmrginp.makeopsT<<"  "<<*dmrginp.initnewsystem<<"  "<<*dmrginp.sysdotmake<<"  "<<*dmrginp.buildcsfops<<" makeops "<<endl;
  p2out << *dmrginp.datatransfer<<" datatransfer "<<endl;
  p2out <<"oneindexopmult   twoindexopmult   Hc  couplingcoeff"<<endl;  
  p2out << *dmrginp.oneelecT<<" "<<*dmrginp.twoelecT<<" "<<*dmrginp.hmultiply<<" "<<*dmrginp.couplingcoeff<<" hmult"<<endl;
  p2out << *dmrginp.buildsumblock<<" "<<*dmrginp.buildblockops<<" build block"<<endl;
  p2out << *dmrginp.blockintegrals<<"  "<<*dmrginp.blocksites<<"  "<<*dmrginp.statetensorproduct<<"  "<<*dmrginp.statecollectquanta<<"  "<<*dmrginp.buildsumblock<<" "<<*dmrginp.builditeratorsT<<"  "<<*dmrginp.diskio<<" build sum block"<<endl;
  p2out << "addnoise  S_0_opxop  S_1_opxop   S_2_opxop"<<endl;
  p3out << *dmrginp.addnoise<<" "<<*dmrginp.s0time<<" "<<*dmrginp.s1time<<" "<<*dmrginp.s2time<<endl;

}

void SweepGenblock::BlockAndDecimate (SweepParams &sweepParams, SpinBlock& system, SpinBlock& newSystem, const bool &useSlater, const bool& dot_with_sys, int stateA, int stateB)
{
  if (dmrginp.outputlevel() > 0) 
    mcheck("at the start of block and decimate");
  pout << "\t\t\t Performing Blocking"<<endl;
  dmrginp.guessgenT -> start();
  // figure out if we are going forward or backwards  
  bool forward = (system.get_sites() [0] == 0);
  SpinBlock systemDot;
  int systemDotStart, systemDotEnd;
  int systemDotSize = sweepParams.get_sys_add() - 1;
  if (forward)
  {
    systemDotStart = dmrginp.spinAdapted() ? *system.get_sites().rbegin () + 1 : (*system.get_sites().rbegin ())/2 + 1 ;
    systemDotEnd = systemDotStart + systemDotSize;
  }
  else
  {
    systemDotStart = dmrginp.spinAdapted() ? system.get_sites()[0] - 1 : (system.get_sites()[0])/2 - 1 ;
    systemDotEnd = systemDotStart - systemDotSize;
  }
  vector<int> spindotsites(2); 
  spindotsites[0] = systemDotStart;
  spindotsites[1] = systemDotEnd;
  systemDot = SpinBlock(systemDotStart, systemDotEnd, stateA==stateB);

  const int nexact = forward ? sweepParams.get_forward_starting_size() : sweepParams.get_backward_starting_size();

  system.addAdditionalCompOps();
  InitBlocks::InitNewSystemBlock(system, systemDot, newSystem, stateA, stateB, sweepParams.get_sys_add(), dmrginp.direct(), DISTRIBUTED_STORAGE, dot_with_sys, true);

  pout << "\t\t\t System  Block"<<newSystem;
  if (dmrginp.outputlevel() > 0)
    newSystem.printOperatorSummary();

  std::vector<Matrix> leftrotateMatrix, rightrotateMatrix;

  LoadRotationMatrix (newSystem.get_sites(), leftrotateMatrix, stateA);
  LoadRotationMatrix (newSystem.get_sites(), rightrotateMatrix, stateB);

#ifndef SERIAL
  mpi::communicator world;
  broadcast(world, leftrotateMatrix, 0);
  broadcast(world, rightrotateMatrix, 0);
#endif

  pout <<"\t\t\t Performing Renormalization "<<endl<<endl;

  if (stateB == stateA)
    newSystem.transform_operators(leftrotateMatrix);
  else
    newSystem.transform_operators(leftrotateMatrix, rightrotateMatrix);


  if (dmrginp.outputlevel() > 0) 
    //mcheck("after rotation and transformation of block");
  if (dmrginp.outputlevel() > 0) 
    pout <<newSystem<<endl;
  if (dmrginp.outputlevel() > 0)
    newSystem.printOperatorSummary();
  //mcheck("After renorm transform");
}