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

void finish_tensor_trace( SpinBlock& b, SpinBlock* sysdot, SparseMatrix& sysdot_op, SparseMatrix& op, std::string& build_pattern )
{
  // Build and store new operator
  assert( ! op.get_built() );
  op.set_built() = true;
  op.set_build_pattern() = build_pattern;
//FIXME magic number 2
  op.set_deltaQuantum(1, op.get_quantum_ladder().at( build_pattern ).at(2) );
  op.allocate(b.get_stateInfo());

  SpinAdapted::operatorfunctions::TensorTrace(sysdot, sysdot_op, &b, &(b.get_stateInfo()), op);
}

double SparseMatrix::memoryUsed(const SpinBlock& b)
{
  StateInfo stateinfo = b.get_stateInfo();
  double memory = 0.0;

  for (int i=0; i < stateinfo.quanta.size(); i++)
    for (int j=0; j<stateinfo.quanta.size(); j++)
      if (allowedQuantaMatrix(i,j) ) {
	memory += 8.0*operatorMatrix(i,j).Storage();
      }
  return memory;
}

/*
void SpinAdapted::operatorfunctions::TensorMultiply(const SpinBlock *ablock, const Baseoperator<Matrix>& a, const Baseoperator<Matrix>& b, const SpinBlock *cblock, Wavefunction& c, Wavefunction& v, const SpinQuantum opQ, double scale)
{
  // can be used for situation with different bra and ket
  const int leftBraOpSz = cblock->get_leftBlock()->get_braStateInfo().quanta.size ();
  const int leftKetOpSz = cblock->get_leftBlock()->get_ketStateInfo().quanta.size ();
  const int rightBraOpSz = cblock->get_rightBlock()->get_braStateInfo().quanta.size ();
  const int rightKetOpSz = cblock->get_rightBlock()->get_ketStateInfo().quanta.size ();

  const StateInfo* lbraS = cblock->get_braStateInfo().leftStateInfo, *rbraS = cblock->get_braStateInfo().rightStateInfo;
  const StateInfo* lketS = cblock->get_ketStateInfo().leftStateInfo, *rketS = cblock->get_ketStateInfo().rightStateInfo;

  const char conjC = (cblock->get_leftBlock() == ablock) ? 'n' : 't';

  const Baseoperator<Matrix>& leftOp = (conjC == 'n') ? a : b; // an ugly hack to support the release memory optimisation
  const Baseoperator<Matrix>& rightOp = (conjC == 'n') ? b : a;
  const char leftConj = (conjC == 'n') ? a.conjugacy() : b.conjugacy();
  const char rightConj = (conjC == 'n') ? b.conjugacy() : a.conjugacy();

  Wavefunction u;
  u.resize(leftBraOpSz*leftKetOpSz, rightKetOpSz);

  int totalmem =0;

  {
    for (int lQrQPrime = 0; lQrQPrime<leftBraOpSz*rightKetOpSz; ++lQrQPrime)
    {
      int rQPrime = lQrQPrime%rightKetOpSz, lQ = lQrQPrime/rightKetOpSz;
	for (int lQPrime = 0; lQPrime < leftKetOpSz; lQPrime++)
	  if (leftOp.allowed(lQ, lQPrime) && c.allowed(lQPrime, rQPrime))
	  {
	    int lindex = lQ*leftKetOpSz+lQPrime;
	    u.allowed(lindex, rQPrime) = true;
            
	    u(lindex,rQPrime).ReSize(lbraS->getquantastates(lQ), rketS->getquantastates(rQPrime));
	    double factor = leftOp.get_scaling(lbraS->quanta[lQ], lketS->quanta[lQPrime]);
	    MatrixMultiply (leftOp.operator_element(lQ, lQPrime), leftConj, c.operator_element(lQPrime, rQPrime), 'n',
			    u.operator_element(lindex, rQPrime), factor, 0.);	      

	  }
    }
  }

  pout << "after first step in tensormultiply"<<endl;
      mcheck("before davidson but after all blocks are built");

  {
    for (int lQrQ = 0; lQrQ<leftBraOpSz*rightBraOpSz; ++lQrQ)
    {
      int rQ = lQrQ%rightBraOpSz, lQ=lQrQ/rightBraOpSz;
	if (v.allowed(lQ, rQ))
	  for (int rQPrime = 0; rQPrime < rightKetOpSz; rQPrime++)
	    if (rightOp.allowed(rQ, rQPrime))
	      for (int lQPrime = 0; lQPrime < leftKetOpSz; lQPrime++)
		if (leftOp.allowed(lQ, lQPrime) && u.allowed(lQ*leftKetOpSz+lQPrime, rQPrime))
		{
		  int lindex = lQ*leftKetOpSz+lQPrime;
		  double factor = scale;

		  factor *= dmrginp.get_ninej()(lketS->quanta[lQPrime].get_s().getirrep(), rketS->quanta[rQPrime].get_s().getirrep() , c.get_deltaQuantum(0).get_s().getirrep(), 
						leftOp.get_spin().getirrep(), rightOp.get_spin().getirrep(), opQ.get_s().getirrep(),
						lbraS->quanta[lQ].get_s().getirrep(), rbraS->quanta[rQ].get_s().getirrep() , v.get_deltaQuantum(0).get_s().getirrep());
		  factor *= Symmetry::spatial_ninej(lketS->quanta[lQPrime].get_symm().getirrep() , rketS->quanta[rQPrime].get_symm().getirrep(), c.get_symm().getirrep(), 
				       leftOp.get_symm().getirrep(), rightOp.get_symm().getirrep(), opQ.get_symm().getirrep(),
				       lbraS->quanta[lQ].get_symm().getirrep() , rbraS->quanta[rQ].get_symm().getirrep(), v.get_symm().getirrep());
		  int parity = rightOp.get_fermion() && IsFermion(lketS->quanta[lQPrime]) ? -1 : 1;
		  factor *=  rightOp.get_scaling(rbraS->quanta[rQ], rketS->quanta[rQPrime]);
		  MatrixMultiply (u.operator_element(lindex, rQPrime), 'n',
				  rightOp(rQ, rQPrime), TransposeOf(rightOp.conjugacy()), v.operator_element(lQ, rQ), factor*parity);
		}
    }
  }
	      
}
*/
void SpinAdapted::operatorfunctions::OperatorScaleAdd(double scaleV, const SpinBlock& b, const Baseoperator<Matrix>& op1, Baseoperator<Matrix>& op2)
{
  const StateInfo& s = b.get_stateInfo();
  for (int lQ = 0; lQ< op2.nrows(); lQ++)
    for (int rQ = 0; rQ<op2.ncols(); rQ++)
      if (op2.allowed(lQ, rQ) && op1.allowed(lQ,rQ))
      {
	double factor = op1.get_scaling(s.quanta[lQ], s.quanta[rQ]);
	if (op1.conjugacy() == 't')
	  MatrixScaleAdd(scaleV*factor, op1.operator_element(lQ,rQ).t(), op2.operator_element(lQ,rQ));
	else
	  MatrixScaleAdd(scaleV*factor, op1.operator_element(lQ,rQ), op2.operator_element(lQ,rQ));
      }

}

void SparseMatrix::buildUsingCsf(const SpinBlock& b, vector< vector<Csf> >& ladders, std::vector< Csf >& s) 
{
  StateInfo stateinfo = b.get_stateInfo();
  built = true;
  allocate(stateinfo);

  for (int i=0; i < stateinfo.quanta.size(); i++)
    for (int j=0; j<stateinfo.quanta.size(); j++)
      if (allowedQuantaMatrix(i,j) )
	for (int jq =stateinfo.unBlockedIndex[j]; jq < stateinfo.unBlockedIndex[j]+stateinfo.quantaStates[j]; jq++) 
	{
	  for (int iq =stateinfo.unBlockedIndex[i]; iq < stateinfo.unBlockedIndex[i]+stateinfo.quantaStates[i]; iq++) 
	    operatorMatrix(i,j)(iq-stateinfo.unBlockedIndex[i]+1, jq-stateinfo.unBlockedIndex[j]+1) = redMatrixElement(s[iq], ladders[jq], &b);
	}

}

void finish_tensor_product( SpinBlock& b, SpinBlock* sysdot, 
                            const SparseMatrix& sysdot_op1, const SparseMatrix& sysdot_op2, SparseMatrix& op, 
                            bool include_parity, std::string& build_pattern )
{
  // Build and store new operator
  assert( ! op.get_built() );
  op.set_built() = true;
  op.set_build_pattern() = build_pattern;
  op.set_deltaQuantum(1, op.get_quantum_ladder().at( build_pattern ).at(2) );
  op.allocate(b.get_stateInfo());

  // Do tensor product
  double parity = 1.0;
  if ( include_parity ) parity = getCommuteParity( sysdot_op1.get_deltaQuantum(0), sysdot_op2.get_deltaQuantum(0), op.get_deltaQuantum(0) );
  SpinAdapted::operatorfunctions::TensorProduct(sysdot, sysdot_op1, sysdot_op2, &b, &(b.get_stateInfo()), op, parity);
}

//-------------------------------------------------------------------------------------------------------------------------------------------------------------
//  (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 SparseMatrix::allocate(const SpinBlock& b)
{
  allocate(b.get_stateInfo());
}

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

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


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

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

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

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

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

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

void SweepGenblock::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");
}