C++ SweepParams::current_root方法代码示例

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 SweepTwopdm::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;
  //if (useSlater) {
  systemDot = SpinBlock(systemDotStart, systemDotEnd, system.get_integralIndex(), true);
    //SpinBlock::store(true, systemDot.get_sites(), systemDot);
    //}
    //else
    //SpinBlock::restore(true, spindotsites, systemDot);
  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.current_root(), sweepParams.current_root(), sweepParams.get_sys_add(), dmrginp.direct(), system.get_integralIndex(), DISTRIBUTED_STORAGE, true, true);
  
  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(), 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> 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
  if (sweepParams.get_block_iter() == 0)
    compute_twopdm_initial(solution, system, systemDot, newSystem, newEnvironment, big, numprocs, state);

  compute_twopdm_sweep(solution, system, systemDot, newSystem, newEnvironment, big, numprocs, state);

  if (sweepParams.get_block_iter()  == sweepParams.get_n_iters() - 1)
    compute_twopdm_final(solution, system, systemDot, newSystem, newEnvironment, big, numprocs, state);

  SaveRotationMatrix (newSystem.get_sites(), rotateMatrix, state);

  //for(int i=0;i<dmrginp.nroots();++i)
  solution[0].SaveWavefunctionInfo (big.get_stateInfo(), big.get_leftBlock()->get_sites(), state);

  newSystem.transform_operators(rotateMatrix);

}

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 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 responseSweep(double sweep_tol, int targetState, vector<int>& projectors, vector<int>& baseStates)
{
  double last_fe = 1.e6;
  double last_be = 1.e6;
  double old_fe = 0.;
  double old_be = 0.;
  SweepParams sweepParams;

  bool direction, warmUp, restart;
  int restartSize=0;
  direction = true; //forward
  warmUp = true; //startup sweep
  restart = false; //not a restart

  sweepParams.current_root() = -1;

  algorithmTypes atype = dmrginp.algorithm_method();
  dmrginp.set_algorithm_method() = ONEDOT;

  //the baseState is the initial guess for the targetState
  if (FULLRESTART) {
    sweepParams.restorestate(direction, restartSize);
    direction = !direction;
    dmrginp.setGuessState() = targetState;
    last_fe = SweepResponse::do_one(sweepParams, warmUp, direction, restart, restartSize, targetState, projectors, baseStates);
    bool tempdirection;
    sweepParams.restorestate(tempdirection, restartSize);
    sweepParams.calc_niter();
    sweepParams.set_sweep_iter() = 0;
    sweepParams.set_restart_iter() = 0;
    sweepParams.savestate(tempdirection, restartSize);
  }
  else if (RESTART) {
    dmrginp.set_algorithm_method() = atype;
    warmUp = false;
    restart = true;
    sweepParams.restorestate(direction, restartSize);
    last_fe = SweepResponse::do_one(sweepParams, warmUp, direction, restart, restartSize, targetState, projectors, baseStates);
  }
  else 
    last_fe = SweepResponse::do_one(sweepParams, warmUp, direction, restart, restartSize, targetState, projectors, baseStates);

  dmrginp.set_algorithm_method() = atype;
  restart = false;
  restartSize = 0;
  warmUp = false;
  while ( true)
    {
      old_fe = last_fe;
      old_be = last_be;
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	break;

      last_be = SweepResponse::do_one(sweepParams, warmUp, !direction, restart, restartSize, targetState, projectors, baseStates);
      p1out << "\t\t\t Finished Sweep Iteration "<<sweepParams.get_sweep_iter()<<endl;
      
      if(dmrginp.max_iter() <= sweepParams.get_sweep_iter())
	break;
      
      last_fe = SweepResponse::do_one(sweepParams, warmUp, direction, restart, restartSize, targetState, projectors, baseStates);

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

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

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