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

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

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

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

}