C++ SweepParams类代码示例

//Initialize stateinfo using the rotation matrices
void SpinAdapted::Sweep::InitializeStateInfo(SweepParams &sweepParams, const bool &forward, int currentstate)
{
  sweepParams.set_sweep_parameters();
  sweepParams.set_block_iter() = 0;

  std::vector<int> sites;
  int new_site, wave_site;
  if (forward) 
    new_site = 0;
  else 
    new_site = dmrginp.spinAdapted() ? dmrginp.last_site()-1 : dmrginp.last_site()/2-1;

  if (dmrginp.spinAdapted())
    sites.push_back(new_site);
  else {
    sites.push_back(2*new_site);
    sites.push_back(2*new_site+1);
    std::sort(sites.begin(), sites.end());
  }
    
  //only need statinfos
  StateInfo stateInfo1; makeStateInfo(stateInfo1, new_site);
  StateInfo::store(forward, sites, stateInfo1, currentstate);
  
  for (; sweepParams.get_block_iter() < sweepParams.get_n_iters(); ) {
      
    if (forward) 
      new_site++;
    else 
      new_site--;

    if (dmrginp.spinAdapted())
      sites.push_back(new_site);
    else {
      sites.push_back(2*new_site);
      sites.push_back(2*new_site+1);
      std::sort(sites.begin(), sites.end());
    }

    
    StateInfo siteState, newState1; 
    makeStateInfo(siteState, new_site);
    TensorProduct(stateInfo1, siteState, newState1, NO_PARTICLE_SPIN_NUMBER_CONSTRAINT);
    newState1.CollectQuanta();
  
      
    //make the newstate
    std::vector<Matrix> rotation1; 
      
    LoadRotationMatrix (sites, rotation1, currentstate);
    
    StateInfo renormState1;
    SpinAdapted::StateInfo::transform_state(rotation1, newState1, renormState1);
    StateInfo::store(forward, sites, renormState1, currentstate);
    stateInfo1 = renormState1;
    ++sweepParams.set_block_iter();
  }
  
}

void SpinAdapted::Sweep::tiny(double sweep_tol)
{
#ifndef SERIAL
  if(mpigetrank() == 0) {
#endif
    pout.precision(12);

  int nroots = dmrginp.nroots(0);
  SweepParams sweepParams;
  sweepParams.set_sweep_parameters();
  StackSpinBlock system(0,dmrginp.last_site()-1, 0, true);
  const StateInfo& sinfo = system.get_stateInfo();
  SpinQuantum hq(0,SpinSpace(0),IrrepSpace(0));
  for (int i=0; i<sinfo.totalStates; i++) {
    if (sinfo.quanta[i] == dmrginp.molecule_quantum()) {
      StackMatrix& h = system.get_op_array(HAM).get_element(0).at(0)->operator_element(i,i);
      DiagonalMatrix energies(h.Nrows()); energies = 0.0;
      diagonalise(h, energies);
      
      for (int x=0; x<nroots; x++) 
	pout << "fullci energy  "<< energies(x+1)<<endl;

      if (mpigetrank() == 0)
      {
#ifndef MOLPRO
	FILE* f = fopen("dmrg.e", "wb");
#else
	std::string efile;
	efile = str(boost::format("%s%s") % dmrginp.load_prefix() % "/dmrg.e" );
	FILE* f = fopen(efile.c_str(), "wb");
#endif
	
	for(int j=0;j<nroots;++j) {
	  double e = energies(j+1); 
	  fwrite( &e, 1, sizeof(double), f);
	}
	fclose(f);
      }

      return;
    }
  }

  pout << "The wavefunction symmetry is not possible with the orbitals supplied."<<endl;
  abort();
#ifndef SERIAL
  }
#endif
}

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];
}

void compress(double sweep_tol, int targetState, int baseState)
{
  double last_fe = 10.e6;
  double last_be = 10.e6;
  double old_fe = 0.;
  double old_be = 0.;
  SweepParams sweepParams;
  bool direction;

  sweepParams.current_root() = -1;
  //this is the warmup sweep, the baseState is used as the initial guess for the targetState
  last_fe = SweepCompress::do_one(sweepParams, true, true, false, 0, targetState, baseState);

  direction = false;
  while ( true)
    {
      old_fe = last_fe;
      old_be = last_be;
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	break;
      last_be = SweepCompress::do_one(sweepParams, false, false, false, 0, targetState, baseState);
      direction = true;
      pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
      
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	break;
      
      last_fe = SweepCompress::do_one(sweepParams, false, true, false, 0, targetState, baseState);
      direction = false;
      
      
      pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
      
    }

  //we finally canonicalize the targetState
  //one has to canonicalize the wavefunction with atleast 3 sweeps, this is a quirk of the way 
  //we transform wavefunction
  if (mpigetrank()==0) {
    Sweep::InitializeStateInfo(sweepParams, !direction, targetState);
    Sweep::InitializeStateInfo(sweepParams, direction, targetState);
    Sweep::CanonicalizeWavefunction(sweepParams, !direction, targetState);
    Sweep::CanonicalizeWavefunction(sweepParams, direction, targetState);
    Sweep::CanonicalizeWavefunction(sweepParams, !direction, targetState);
    Sweep::InitializeStateInfo(sweepParams, !direction, targetState);
    Sweep::InitializeStateInfo(sweepParams, direction, targetState);
    
  }
  
}

void dmrg(double sweep_tol)
{
  double last_fe = 10.e6;
  double last_be = 10.e6;
  double old_fe = 0.;
  double old_be = 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 dodiis = false;

  int domoreIter = 0;
  bool direction;

  //this is regular dmrg calculation
  if(!dmrginp.setStateSpecific()) {
    sweepParams.current_root() = -1;
    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()) )
    {
      old_fe = last_fe;
      old_be = last_be;
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	break;
      last_be = Sweep::do_one(sweepParams, false, false, false, 0);
      direction = true;
      pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
      
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	break;
      
      //For obtaining the extrapolated energy
      old_states=sweepParams.get_keep_states();
      new_states=sweepParams.get_keep_states_ls();
      
      last_fe = Sweep::do_one(sweepParams, false, true, false, 0);
      direction = false;
      
      new_states=sweepParams.get_keep_states();
      
      
      pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
      if (domoreIter == 2) {
	dodiis = true;
	break;
      }
      
    }
  }
  else { //this is state specific calculation  
    const int nroots = dmrginp.nroots();

    bool direction=true;
    int restartsize;
    //sweepParams.restorestate(direction, restartsize);
    //sweepParams.set_sweep_iter() = 0;
    //sweepParams.set_restart_iter() = 0;

    algorithmTypes atype;
    pout << "STARTING STATE SPECIFIC CALCULATION "<<endl;
    for (int i=0; i<nroots; i++) {
      atype = dmrginp.algorithm_method();
      dmrginp.set_algorithm_method() = ONEDOT;
      sweepParams.current_root() = i;

      p1out << "RUNNING GENERATE BLOCKS FOR STATE "<<i<<endl;

      if (mpigetrank()==0) {
	Sweep::InitializeStateInfo(sweepParams, direction, i);
	Sweep::InitializeStateInfo(sweepParams, !direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, !direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
	Sweep::InitializeStateInfo(sweepParams, direction, i);
	Sweep::InitializeStateInfo(sweepParams, !direction, i);

      }

      for (int j=0; j<i ; j++) {
	int integralIndex = 0;
	Sweep::InitializeOverlapSpinBlocks(sweepParams, direction, i, j, integralIndex);
	Sweep::InitializeOverlapSpinBlocks(sweepParams, !direction, i, j, integralIndex);
      }
      dmrginp.set_algorithm_method() = atype;

      p1out << "RUNNING GENERATE BLOCKS FOR STATE "<<i<<endl;

      SweepGenblock::do_one(sweepParams, false, !direction, false, 0, i, i);
      sweepParams.set_sweep_iter() = 0;
      sweepParams.set_restart_iter() = 0;
      sweepParams.savestate(!direction, restartsize);

      
      pout << "STATE SPECIFIC CALCULATION FOR STATE: "<<i<<endl;
//.........这里部分代码省略.........

void restart(double sweep_tol, bool reset_iter)
{
  double last_fe = 100.;
  double last_be = 100.;
  double old_fe = 0.;
  double old_be = 0.;
  bool direction;
  int restartsize;
  SweepParams sweepParams;
  bool dodiis = false;

  int domoreIter = 2;

  sweepParams.restorestate(direction, restartsize);

  if (!dmrginp.setStateSpecific()) {
    if(reset_iter) { //this is when you restart from the start of the sweep
      sweepParams.set_sweep_iter() = 0;
      sweepParams.set_restart_iter() = 0;
    }
    
    if (restartwarm)
      last_fe = Sweep::do_one(sweepParams, true, direction, true, restartsize);
    else
      last_fe = Sweep::do_one(sweepParams, false, direction, true, restartsize);
    
    
    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()) )
      {
	
	old_fe = last_fe;
	old_be = last_be;
	if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	  break;
	last_be = Sweep::do_one(sweepParams, false, !direction, false, 0);
	
	
	if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	  break;
	last_fe = Sweep::do_one(sweepParams, false, direction, false, 0);	
      }
  }
  else { //this is state specific calculation  
    const int nroots = dmrginp.nroots();

    bool direction;
    int restartsize;
    sweepParams.restorestate(direction, restartsize);

    //initialize state and canonicalize all wavefunctions
    int currentRoot = sweepParams.current_root();
    for (int i=0; i<nroots; i++) {
      sweepParams.current_root() = i;
      if (mpigetrank()==0) {
	Sweep::InitializeStateInfo(sweepParams, direction, i);
	Sweep::InitializeStateInfo(sweepParams, !direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, !direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
      }
    }

    //now generate overlaps with all the previous wavefunctions
    for (int i=0; i<currentRoot; i++) {
      sweepParams.current_root() = i;
      if (mpigetrank()==0) {
	for (int j=0; j<i; j++) {
	  int integralIndex = 0;
	  Sweep::InitializeOverlapSpinBlocks(sweepParams, !direction, i, j, integralIndex);
	  Sweep::InitializeOverlapSpinBlocks(sweepParams, direction, i, j, integralIndex);
	}
      }
    }
    sweepParams.current_root() = currentRoot;

    if (sweepParams.current_root() <0) {
      p1out << "This is most likely not a restart calculation and should be done without the restart command!!"<<endl;
      p1out << "Aborting!!"<<endl;
      exit(0);
    }
    pout << "RESTARTING STATE SPECIFIC CALCULATION OF STATE "<<sweepParams.current_root()<<" AT SWEEP ITERATION  "<<sweepParams.get_sweep_iter()<<endl;

    //this is so that the iteration is not one ahead after generate block for restart
    --sweepParams.set_sweep_iter(); sweepParams.savestate(direction, restartsize);
    for (int i=sweepParams.current_root(); i<nroots; i++) {
      sweepParams.current_root() = i;

      p1out << "RUNNING GENERATE BLOCKS FOR STATE "<<i<<endl;

      if (mpigetrank()==0) {
	Sweep::InitializeStateInfo(sweepParams, direction, i);
	Sweep::InitializeStateInfo(sweepParams, !direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, !direction, i);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, i);
	for (int j=0; j<i ; j++) {
	  int integralIndex = 0;
	  Sweep::InitializeOverlapSpinBlocks(sweepParams, direction, i, j, integralIndex);
	  Sweep::InitializeOverlapSpinBlocks(sweepParams, !direction, i, j, integralIndex);
//.........这里部分代码省略.........

int calldmrg(char* input, char* output)
{
  //sleep(15);
  streambuf *backup;
  backup = cout.rdbuf();
  ofstream file;
  if (output != 0) {
    file.open(output);
    pout.rdbuf(file.rdbuf());
  }
  license();
  ReadInput(input);
  MAX_THRD = dmrginp.thrds_per_node()[mpigetrank()];
#ifdef _OPENMP
  omp_set_num_threads(MAX_THRD);
#endif
  pout.precision (12);

   //Initializing timer calls
  dmrginp.initCumulTimer();

  double sweep_tol = 1e-7;
  sweep_tol = dmrginp.get_sweep_tol();
  bool direction;
  int restartsize;
  SweepParams sweepParams;

  SweepParams sweep_copy;
  bool direction_copy; int restartsize_copy;
  Matrix O, H;
  
  switch(dmrginp.calc_type()) {

  case (COMPRESS):
  {
    bool direction; int restartsize;
    //sweepParams.restorestate(direction, restartsize);
    //sweepParams.set_sweep_iter() = 0;
    restartsize = 0;

    int targetState, baseState, correctionVector, firstorderstate;
    {
      direction = true;

      //base state is always defined
      baseState = dmrginp.baseStates()[0];

      //if targetstate is given use it otherwise use basestate+1
      if(dmrginp.targetState() == -1)
	targetState = dmrginp.baseStates()[0]+1;
      else
	targetState = dmrginp.targetState();

      algorithmTypes atype = dmrginp.algorithm_method();
      dmrginp.set_algorithm_method() = ONEDOT;
      //initialize state info and canonicalize wavefunction is always done using onedot algorithm
      if (mpigetrank()==0) {
	Sweep::InitializeStateInfo(sweepParams, direction, baseState);
	Sweep::InitializeStateInfo(sweepParams, !direction, baseState);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, baseState);
	Sweep::CanonicalizeWavefunction(sweepParams, !direction, baseState);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, baseState);
      }
      dmrginp.set_algorithm_method() = atype;
    }

    //this genblock is required to generate all the nontranspose operators
    dmrginp.setimplicitTranspose() = false;
    SweepGenblock::do_one(sweepParams, false, false, false, restartsize, baseState, baseState);


    compress(sweep_tol, targetState, baseState);

    break;
  }
  case (RESPONSEBW):
  {
    //compressing the V|\Psi_0>, here \Psi_0 is the basestate and 
    //its product with V will have a larger bond dimension and is being compressed
    //it is called the target state
    dmrginp.setimplicitTranspose() = false;


    sweepParams.restorestate(direction, restartsize);
    algorithmTypes atype = dmrginp.algorithm_method();
    dmrginp.set_algorithm_method() = ONEDOT;
    if (mpigetrank()==0 && !RESTART && !FULLRESTART) {
      for (int l=0; l<dmrginp.projectorStates().size(); l++) {
	Sweep::InitializeStateInfo(sweepParams, direction, dmrginp.projectorStates()[l]);
	Sweep::InitializeStateInfo(sweepParams, !direction, dmrginp.projectorStates()[l]);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, dmrginp.projectorStates()[l]);
	Sweep::CanonicalizeWavefunction(sweepParams, !direction, dmrginp.projectorStates()[l]);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, dmrginp.projectorStates()[l]);
      }
      for (int l=0; l<dmrginp.baseStates().size(); l++) {
	Sweep::InitializeStateInfo(sweepParams, direction, dmrginp.baseStates()[l]);
	Sweep::InitializeStateInfo(sweepParams, !direction, dmrginp.baseStates()[l]);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, dmrginp.baseStates()[l]);
	Sweep::CanonicalizeWavefunction(sweepParams, !direction, dmrginp.baseStates()[l]);
	Sweep::CanonicalizeWavefunction(sweepParams, direction, dmrginp.baseStates()[l]);
//.........这里部分代码省略.........

void SpinAdapted::Sweep::fullci(double sweep_tol)
{
  int integralIndex = 0;
  SweepParams sweepParams;
  sweepParams.set_sweep_parameters();


  StackSpinBlock system, sysdot;
  InitBlocks::InitStartingBlock(system, true, 0, 0, sweepParams.get_forward_starting_size(),  sweepParams.get_backward_starting_size(), 0, false, true, integralIndex);
  int numsites = dmrginp.spinAdapted() ? dmrginp.last_site() : dmrginp.last_site()/2;
  int forwardsites = numsites/2+numsites%2;
  int backwardsites = numsites - forwardsites;
  SpinQuantum hq(0,SpinSpace(0),IrrepSpace(0));


  StackSpinBlock newSystem;
  for (int i=0; i<forwardsites-1; i++) {
    sysdot = StackSpinBlock(i+1, i+1, integralIndex, true);
    system.addAdditionalOps();
    newSystem.set_integralIndex() = integralIndex;
    if (i == forwardsites-2)
      newSystem.default_op_components(true, true, false, true);
    else
      newSystem.default_op_components(false, true, false, true);

    newSystem.setstoragetype(DISTRIBUTED_STORAGE);
    newSystem.BuildSumBlock (NO_PARTICLE_SPIN_NUMBER_CONSTRAINT, system, sysdot);

    long memoryToFree = newSystem.getdata() - system.getdata();
    long newsysMemory = newSystem.memoryUsed();
    if (i != forwardsites-2) {
      if (i != 0) {
	newSystem.moveToNewMemory(system.getdata());
	Stackmem[0].deallocate(newSystem.getdata()+newSystem.memoryUsed(), memoryToFree);
      }
      system.clear();
      system = newSystem;
    }
  }

  StackSpinBlock environment, newEnvironment, envdot;
  InitBlocks::InitStartingBlock(environment, false, 0, 0, sweepParams.get_forward_starting_size(),  sweepParams.get_backward_starting_size(), 0, false, true, integralIndex);
  for (int i=0;i <backwardsites-1; i++) {
    envdot = StackSpinBlock(numsites-2-i, numsites-2-i, integralIndex, true);
    environment.addAdditionalOps();
    newEnvironment.set_integralIndex() = integralIndex;
    if (i == backwardsites-2)
      newEnvironment.default_op_components(true, false, true, true);
    else
      newEnvironment.default_op_components(false, false, true, true);
    newEnvironment.setstoragetype(DISTRIBUTED_STORAGE);
    newEnvironment.BuildSumBlock (NO_PARTICLE_SPIN_NUMBER_CONSTRAINT, environment, envdot);

    if (i!=backwardsites-2) {
      if (i != 0) {
	long memoryToFree = newEnvironment.getdata() - environment.getdata();
	long newenvMemory = newEnvironment.memoryUsed();
	newEnvironment.moveToNewMemory(environment.getdata());
	Stackmem[0].deallocate(newEnvironment.getdata()+newEnvironment.memoryUsed(), memoryToFree);
      }
      environment.clear();
      environment = newEnvironment;
    }
  }

  pout <<"\t\t\t System Block :: "<< newSystem;
  pout <<"\t\t\t Environment Block :: "<< newEnvironment;
  newSystem.set_loopblock(true); newEnvironment.set_loopblock(false);
  StackSpinBlock big;
  InitBlocks::InitBigBlock(newSystem, newEnvironment, big); 


  int nroots = dmrginp.nroots(0);
  std::vector<StackWavefunction> solution(nroots);

  solution[0].initialise(dmrginp.effective_molecule_quantum_vec(), big.get_leftBlock()->get_stateInfo(), big.get_rightBlock()->get_stateInfo(), false);
  solution[0].Clear();
  if (mpigetrank() == 0) {
    for (int i=1; i<nroots; i++) {
      solution[i].initialise(dmrginp.effective_molecule_quantum_vec(), big.get_leftBlock()->get_stateInfo(), big.get_rightBlock()->get_stateInfo(), false);
      solution[i].Clear();
    }
  }


  std::vector<double> energies(nroots);
  double tol = sweepParams.get_davidson_tol();

  pout << "\t\t\t Solving the Wavefunction "<<endl;
  int currentState = 0;
  std::vector<StackWavefunction> lowerStates;
  Solver::solve_wavefunction(solution, energies, big, tol, BASIC, false, true, false, false, sweepParams.get_additional_noise(), currentState, lowerStates);

  pout << "tensormultiply "<<*dmrginp.tensormultiply<<endl;

  for (int i=0; i<nroots; i++) {
    pout << "fullci energy "<< energies[i]<<endl;
  }
  if (!mpigetrank())
  {
//.........这里部分代码省略.........

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;

}

//before you start optimizing each state you want to initalize all the overlap matrices
void Sweep::InitializeOverlapSpinBlocks(SweepParams &sweepParams, const bool &forward, int stateA, int stateB)
{
  SpinBlock system;

  sweepParams.set_sweep_parameters();
  if (forward)
    pout << "\t\t\t Starting sweep "<< sweepParams.set_sweep_iter()<<" in forwards direction"<<endl;
  else
    pout << "\t\t\t Starting sweep "<< sweepParams.set_sweep_iter()<<" in backwards direction" << endl;
  pout << "\t\t\t ============================================================================ " << endl;

  int restartSize = 0; bool restart = false, warmUp = false;
  InitBlocks::InitStartingBlock (system,forward, stateA, stateB, sweepParams.get_forward_starting_size(), sweepParams.get_backward_starting_size(), restartSize, restart, warmUp);

  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, stateA, stateB); // 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 (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
    {
      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;
      }

      SpinBlock systemDot, environmentDot;
      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;
	}
      systemDot = SpinBlock(systemDotStart, systemDotEnd, true);

      SpinBlock newSystem; // new system after blocking and decimating
      newSystem.initialise_op_array(OVERLAP, false);
      newSystem.setstoragetype(DISTRIBUTED_STORAGE);
      newSystem.BuildSumBlock (NO_PARTICLE_SPIN_NUMBER_CONSTRAINT, system, systemDot);

      std::vector<Matrix> brarotateMatrix, ketrotateMatrix;
      LoadRotationMatrix(newSystem.get_sites(), brarotateMatrix, stateA);
      LoadRotationMatrix(newSystem.get_sites(), ketrotateMatrix, stateB);
      newSystem.transform_operators(brarotateMatrix, ketrotateMatrix);

      
      system = newSystem;
      if (dmrginp.outputlevel() > 0){
	    pout << system<<endl;
      }
      
      SpinBlock::store (forward, system.get_sites(), system, stateA, stateB);	 	
      ++sweepParams.set_block_iter();
      
      sweepParams.savestate(forward, syssites.size());
      if (dmrginp.outputlevel() > 0)
         mcheck("at the end of sweep iteration");
    }

  pout << "\t\t\t ============================================================================ " << endl;

  // update the static number of iterations
  return ;
  
}

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;

  last_fe = Sweep::do_one(sweepParams, false, direction, true, restartsize);

  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() <= sweepParams.get_sweep_iter()) 
	break;

      last_be = Sweep::do_one(sweepParams, false, !direction, false, 0);
      pout << "\t\t\t 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();


      pout << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;

    }
  pout << "Converged Energy  " << sweepParams.get_lowest_energy()[0]<< std::endl;
  if(dmrginp.max_iter() <= sweepParams.get_sweep_iter()) {
    
    pout << "Maximum sweep iterations achieved " << std::endl;
  }

  //one has to canonicalize the wavefunction with atleast 3 sweeps, this is a quirk of the way 
  //we transform wavefunction
  if (mpigetrank()==0) {
    Sweep::InitializeStateInfo(sweepParams, !direction, targetState);
    Sweep::InitializeStateInfo(sweepParams, direction, targetState);
    Sweep::CanonicalizeWavefunction(sweepParams, !direction, targetState);
    Sweep::CanonicalizeWavefunction(sweepParams, direction, targetState);
    Sweep::CanonicalizeWavefunction(sweepParams, !direction, targetState);
    Sweep::InitializeStateInfo(sweepParams, !direction, targetState);
    Sweep::InitializeStateInfo(sweepParams, direction, targetState);
    
  }

}

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
//.........这里部分代码省略.........

void SpinAdapted::mps_nevpt::type1::subspace_Va(int baseState)
{
  
  double energy=0;
  double overlap=0;
  VaPerturber pb;
  MPS::siteBlocks.clear();
  int virtsize = dmrginp.spinAdapted()? dmrginp.virt_size():dmrginp.virt_size()*2;
  int virtshift = dmrginp.spinAdapted()? dmrginp.core_size()+dmrginp.act_size(): (dmrginp.core_size()+dmrginp.act_size())*2;
  //int virtsize = dmrginp.virt_size();
  //int virtshift = dmrginp.core_size()+dmrginp.act_size();
  for(int i=0; i< virtsize; i++){
    double perturberEnergy=0;
    dmrginp.calc_type() = MPS_NEVPT;
    pb.init(i+virtshift);
    pout << "Begin Va subspace with a = " << pb.orb(0)<<endl;
    SweepParams sweepParams;
    sweepParams.set_sweep_parameters();
    sweepParams.current_root() = baseState;
    //sweepParams.current_root() = -1;
    //double last_fe = Startup(sweepParams, true, true, false, 0, pb, baseState);
    Timer timer;
    Startup(sweepParams, true, pb, baseState);
    pout <<"Start up time :" << timer.elapsedwalltime();
    //sweepParams.current_root() = baseState;
    timer.start();
    while(true)
    {
      do_one(sweepParams, false, false, false, 0, pb, baseState);
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	      break;
      do_one(sweepParams, false, true, false, 0, pb, baseState);
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
      {
        cleanup(baseState, pb);
	      break;
      }
    }
    pout <<"Sweep time :" << timer.elapsedwalltime();
//    while ( true)
//      {
//        old_fe = last_fe;
//        old_be = last_be;
//        if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
//	        break;
//        last_be = do_one(sweepParams, false, false, false, 0, pb, baseState);
//        if (dmrginp.outputlevel() > 0) 
//	        pout << "Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
//        
//        if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
//	        break;
//        
//        last_fe = do_one(sweepParams, false, true, false, 0, pb, baseState);
//        
//        if (dmrginp.outputlevel() > 0)
//	        pout << "Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
//        
//      }
//

 //   if (mpigetrank()==0) {
 //     bool direction = true;
 //     Sweep::InitializeStateInfo(sweepParams, !direction, pb.wavenumber(),pb.braquanta );
 //     Sweep::InitializeStateInfo(sweepParams, direction, pb.wavenumber(),pb.braquanta );
 //     Sweep::CanonicalizeWavefunction(sweepParams, !direction, pb.wavenumber(),pb.braquanta );
 //     Sweep::CanonicalizeWavefunction(sweepParams, direction, pb.wavenumber(),pb.braquanta );
 //     Sweep::CanonicalizeWavefunction(sweepParams, !direction, pb.wavenumber(),pb.braquanta );
 //     Sweep::InitializeStateInfo(sweepParams, !direction, pb.wavenumber(),pb.braquanta );
 //     Sweep::InitializeStateInfo(sweepParams, direction, pb.wavenumber(),pb.braquanta );
 //     
 //   }
   
    MPS pbmps(pb.wavenumber());
    double o, h;
    dmrginp.calc_type() = DMRG;

    timer.start();

    calcHamiltonianAndOverlap(pbmps, h, o,pb);

    pout <<"Calculate Expectation time :" << timer.elapsedwalltime();

    if(!dmrginp.spinAdapted())
    {
      //In nonspinAdapted, alpha and beta have the results. Only one is neccessary. 
      o*=2;
      h*=2;
      i++;
    }
    if(o> NUMERICAL_ZERO){
      double fock =dmrginp.spinAdapted()? v_1[0](2*(i+virtshift),2*(i+virtshift)): v_1[0](i+virtshift,i+virtshift);
      //perturberEnergy = h/o+fock+perturber::CoreEnergy[0];
      perturberEnergy = h/o+fock;
      energy += o/(mps_nevpt::ZeroEnergy[baseState]- perturberEnergy) ;
      //overlap +=o;
      overlap += sqrt(o)/(mps_nevpt::ZeroEnergy[baseState]- perturberEnergy);
      if (dmrginp.outputlevel() > 0){
        pout << "Amplitude : " << sqrt(o)/(mps_nevpt::ZeroEnergy[baseState]- perturberEnergy) <<endl;
        pout << "Ener(only CAS part) : " << h/o<<endl;
        pout << "Energy : " << perturberEnergy<<endl;
//.........这里部分代码省略.........

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;
  
//.........这里部分代码省略.........