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

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

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

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

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

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

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

}

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

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

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

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

}

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

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


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

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

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

  DiagonalMatrix e;


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

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

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

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

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

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

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

void compute_one_pdm_1_1(Wavefunction& wave1, Wavefunction& wave2, const SpinBlock& big, Matrix& onepdm)
{
  SpinBlock* leftBlock = big.get_leftBlock();
  SpinBlock* rightBlock = big.get_rightBlock();

  int ketS = dmrginp.total_spin_number().getirrep();
  int braS = dmrginp.bra_spin_number().getirrep();


  for (int j = 0; j < rightBlock->get_op_array(CRE).get_size(); ++j)	
   {
    boost::shared_ptr<SparseMatrix> op2 = rightBlock->get_op_array(CRE).get_local_element(j)[0]->getworkingrepresentation(rightBlock);
    int jx = op2->get_orbs(0);
    for (int i = 0; i < leftBlock->get_op_array(DES).get_size(); ++i)
    {
      boost::shared_ptr<SparseMatrix> op1 = leftBlock->get_op_array(DES).get_local_element(i)[0]->getworkingrepresentation(leftBlock);
      int ix = op1->get_orbs(0);

      vector<SpinQuantum> opQ = op2->get_deltaQuantum(0)+op1->get_deltaQuantum(0);
      Wavefunction opw2;
      vector<SpinQuantum> dQ = wave1.get_deltaQuantum();
      opw2.initialisebra(dQ, &big, true);

      operatorfunctions::TensorMultiply(rightBlock, *op2, *op1, &big, wave2, opw2, opQ[0], 1.0);
      double sum = DotProduct(wave1, opw2);

    pout << " right  CRE  left DES "  <<endl;

      onepdm(2*jx+1, 2*ix+1) = sum/sqrt(2.0);
      onepdm(2*jx+2, 2*ix+2) = sum/sqrt(2.0);
    pout << "onepdm(2*jx+1, 2*ix+1)   "<<  "ix   "<< ix<<"  jx    "<< jx  <<"   "<< 2*jx+1 <<"   "<< 2*ix+1<<"  "<<onepdm(2*jx+1, 2*ix+1)<<endl;
    pout << "onepdm(2*jx+2, 2*ix+2)   "<<  "ix   "<< ix<<"  jx    "<< jx  <<"   "<< 2*jx+2 <<"   "<< 2*ix+2<<"  "<<onepdm(2*jx+2, 2*ix+2)<<endl;
    }
  }

//-----
  for (int j = 0; j < rightBlock->get_op_array(DES).get_size(); ++j)
  {
    boost::shared_ptr<SparseMatrix> op2 = rightBlock->get_op_array(DES).get_local_element(j)[0]->getworkingrepresentation(rightBlock);
    int jx = op2->get_orbs(0);
    for (int i = 0; i < leftBlock->get_op_array(CRE).get_size(); ++i)
    {
      boost::shared_ptr<SparseMatrix> op1 = leftBlock->get_op_array(CRE).get_local_element(i)[0]->getworkingrepresentation(leftBlock);
      int ix = op1->get_orbs(0);

      vector<SpinQuantum> opQ = op1->get_deltaQuantum(0)+op2->get_deltaQuantum(0);
      Wavefunction opw2;
      vector<SpinQuantum> dQ = wave1.get_deltaQuantum();
      opw2.initialisebra(dQ, &big, true);
      operatorfunctions::TensorMultiply(leftBlock, *op1, *op2, &big, wave2, opw2, opQ[0], 1.0);
      double sum = DotProduct(wave1, opw2);

    pout << " left  CRE  right DES "  <<endl;

      onepdm(2*ix+1, 2*jx+1) = sum/sqrt(2.0);
      onepdm(2*ix+2, 2*jx+2) = sum/sqrt(2.0);
    pout << "onepdm(2*ix+1, 2*jx+1)   "<<  "ix   "<< ix<<"  jx    "<< jx  <<"   "<< 2*ix+1 <<"   "<< 2*jx+1<<"  "<<onepdm(2*ix+1, 2*jx+1)<<endl;
    pout << "onepdm(2*ix+2, 2*jx+2)   "<<  "ix   "<< ix<<"  jx    "<< jx  <<"   "<< 2*ix+2 <<"   "<< 2*jx+2<<"  "<<onepdm(2*ix+2, 2*jx+2)<<endl;
   }
  }
}

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

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

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


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

  std::vector<Matrix> rotateMatrix;


  if (!dmrginp.get_fullrestart()) {
    //this should be done when we actually have wavefunctions stored, otherwise not!!
    SpinBlock environment, environmentDot, newEnvironment;
    int environmentDotStart, environmentDotEnd, environmentStart, environmentEnd;
    InitBlocks::InitNewEnvironmentBlock(environment, systemDot, newEnvironment, system, systemDot,
					sweepParams.get_sys_add(), sweepParams.get_env_add(), forward, dmrginp.direct(),
					sweepParams.get_onedot(), nexact, useSlater, true, true, true);
    SpinBlock big;
    InitBlocks::InitBigBlock(newSystem, newEnvironment, big); 
    DiagonalMatrix e;
    std::vector<Wavefunction> solution(1);
    
    GuessWave::guess_wavefunctions(solution[0], e, big, sweepParams.get_guesstype(), true, state, true, 0.0); 
    solution[0].SaveWavefunctionInfo (big.get_stateInfo(), big.get_leftBlock()->get_sites(), state);


    DensityMatrix tracedMatrix;
    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 = newSystem.makeRotateMatrix(tracedMatrix, rotateMatrix, sweepParams.get_keep_states(), sweepParams.get_keep_qstates());
    
  }
  else
    LoadRotationMatrix (newSystem.get_sites(), rotateMatrix, state);

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

  if (!dmrginp.get_fullrestart())
    SaveRotationMatrix (newSystem.get_sites(), rotateMatrix, state);

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