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

void do_4index_2_2_tensor_products( bool forwards, const opTypes& optype, const opTypes& rhsType, const opTypes& lhsType,
                                    SpinBlock& big, SpinBlock* rhsBlock, SpinBlock* lhsBlock, std::ofstream& ofs,
                                    const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  Op_component_base& rhs_array = rhsBlock->get_op_array(rhsType);
  Op_component_base& lhs_array = lhsBlock->get_op_array(lhsType);
  assert ( (rhs_array.get_size() == 1) || (lhs_array.get_size() == 1) );

  // Loop over all rhs operator indices //FIXME copy or reference?
  std::vector<boost::shared_ptr<SparseMatrix> > rhs_ops;
  for (int iidx = 0; iidx < rhs_array.get_size(); ++iidx) {
    rhs_ops = rhs_array.get_local_element(iidx);

    // Loop over all lhs operator indices
    std::vector<boost::shared_ptr<SparseMatrix> > lhs_ops;
    for (int idx = 0; idx < lhs_array.get_size(); ++idx) {
      lhs_ops = lhs_array.get_local_element(idx);
      int i = rhs_ops[0]->get_orbs()[0]; int j = rhs_ops[0]->get_orbs()[1];
      int k = lhs_ops[0]->get_orbs()[0]; int l = lhs_ops[0]->get_orbs()[1];
      // In parallel calculations not all operators are built on each proc
      if ( ! big.get_op_array(optype).has_local_index(i,j,k,l) ) continue;
//pout << "p" << mpigetrank() << "; i,j,k,l = " << i << "," << j << "," << k << "," << l << endl;
      std::vector<boost::shared_ptr<SparseMatrix> > vec = big.get_op_array(optype).get_element(i,j,k,l);

      // Loop over rhs spin-op components
      for (int jjdx=0; jjdx < rhs_ops.size(); jjdx++) {
        boost::shared_ptr<SparseMatrix>& rhs_op = rhs_ops[jjdx];
        assert( rhs_op->get_built() );
        const SpinQuantum& spin_12 = rhs_op->get_deltaQuantum(0);

        // Loop over lhs spin-op components
        for (int jdx=0; jdx < lhs_ops.size(); jdx++) {
          boost::shared_ptr<SparseMatrix>& lhs_op = lhs_ops[jdx];
          assert( lhs_op->get_built() );
          std::string build_12 = rhs_op->get_build_pattern();
          std::string build_34 = lhs_op->get_build_pattern();
          std::string build_pattern = "(" + build_12 + build_34 + ")";
          const SpinQuantum& spin_34 = lhs_op->get_deltaQuantum(0);

          // Allocate and build new operator
          for (int sx=0; sx < vec.size(); sx++) {
            boost::shared_ptr<SparseMatrix>& op = vec[sx];
            std::vector<SpinQuantum> s1 = { op->get_quantum_ladder().at(build_pattern).at(0), op->get_quantum_ladder().at(build_pattern).at(1) };
            std::vector<SpinQuantum> s2 = { spin_12, spin_34 };
            // Select relevant spin components
            if ( s1 == s2 ) {
              finish_tensor_product( big, rhsBlock, *rhs_op, *lhs_op, *op, forwards, build_pattern );
              // Renormalise operator
              op->renormalise_transform( rotateMatrix, stateinfo );
            }
          }
        }
      }
      // Store spin-batch on disk 
      if ( ! dmrginp.do_npdm_in_core() ) store_ops_on_disk( ofs, vec );
    }
  }

}

void do_4index_tensor_trace( const opTypes& optype, SpinBlock& big, SpinBlock* sysdot, std::ofstream& ofs,
                             const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  // Get pointer to sparse operator array
  Op_component_base& sysdot_array = sysdot->get_op_array(optype);
  // Open filesystem if necessary
  std::ifstream ifs;
  if ( (! dmrginp.do_npdm_in_core()) && sysdot->size() > 1 ) ifs.open( sysdot_array.get_filename().c_str(), std::ios::binary );

//FIXME need reference?  don't want to copy?
  // Loop over all operator indices
  std::vector<boost::shared_ptr<SparseMatrix> > sysdot_ops;
  for (int idx = 0; idx < sysdot_array.get_size(); ++idx) {
    if ( dmrginp.do_npdm_in_core() || sysdot->size() <= 1) 
      sysdot_ops = sysdot_array.get_local_element(idx);
    else
//FIXME size
      sysdot_ops = get_ops_from_disk( ifs, sysdot_array.get_local_element(0).size() );

    // Loop over spin-op components
    int i = sysdot_ops[0]->get_orbs()[0]; int j = sysdot_ops[0]->get_orbs()[1]; 
    int k = sysdot_ops[0]->get_orbs()[2]; int l = sysdot_ops[0]->get_orbs()[3];
    // In parallel calculations not all operators are built on each proc
    if ( ! big.get_op_array(optype).has_local_index(i,j,k,l) ) continue;
//pout << "p" << mpigetrank() << "; i,j,k,l = " << i << "," << j << "," << k << "," << l << endl;
    std::vector<boost::shared_ptr<SparseMatrix> > new_ops = big.get_op_array(optype).get_element(i,j,k,l);
    for (int jdx=0; jdx < sysdot_ops.size(); jdx++) {
      boost::shared_ptr<SparseMatrix>& sysdot_op = sysdot_ops[jdx];
      assert( sysdot_op->get_built() );
      std::string build_pattern = sysdot_op->get_build_pattern();

      // Allocate and build new operator
      for (int sx=0; sx < new_ops.size(); sx++) {
        boost::shared_ptr<SparseMatrix>& op = new_ops[sx];
        std::vector<SpinQuantum> s1 = sysdot_op->get_quantum_ladder().at(build_pattern);
        std::vector<SpinQuantum> s2 = op->get_quantum_ladder().at(build_pattern);
        // Store spin component in correct location
        if ( s1 == s2 ) {
          finish_tensor_trace( big, sysdot, *sysdot_op, *op, build_pattern );
          // Renormalise operator
          op->renormalise_transform( rotateMatrix, stateinfo );
        }
      }
    }
    // Store spin-batch on disk 
    if ( ! dmrginp.do_npdm_in_core() ) store_ops_on_disk( ofs, new_ops );
  }

  // Close filesystem if necessary
  if ( ifs.is_open() ) ifs.close();

}

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 do_4index_3_1_tensor_products( bool forwards, const opTypes& optype, const opTypes& rhsType, const opTypes& lhsType,
                                    SpinBlock& big, SpinBlock* rhsBlock, SpinBlock* lhsBlock, std::ofstream& ofs,
                                    const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  // (i,j,k | l ) partition
  //-------------------------
  Op_component_base& rhs_array = rhsBlock->get_op_array(rhsType);
  Op_component_base& lhs_array = lhsBlock->get_op_array(lhsType);
  assert ( (rhs_array.get_size() == 1) || (lhs_array.get_size() == 1) );

  // Initialize filesystem
  std::ifstream rhsifs;
  if ( (! dmrginp.do_npdm_in_core()) && rhsBlock->size() > 1 ) rhsifs.open( rhs_array.get_filename().c_str(), std::ios::binary );

  // Loop over all rhs operator indices
  for (int idx = 0; idx < rhs_array.get_size(); ++idx) {
    std::vector<boost::shared_ptr<SparseMatrix> > rhs_ops;
    if ( dmrginp.do_npdm_in_core() || rhsBlock->size() <= 1 )
      rhs_ops = rhs_array.get_local_element(idx);
    else
      rhs_ops = get_ops_from_disk( rhsifs, rhs_array.get_local_element(0).size() );

    // Loop over all lhs operator indices
    for (int iidx = 0; iidx < lhs_array.get_size(); ++iidx) {
      std::vector<boost::shared_ptr<SparseMatrix> > lhs_ops = lhs_array.get_local_element(iidx);
      int i = rhs_ops[0]->get_orbs()[0]; int j = rhs_ops[0]->get_orbs()[1]; int k = rhs_ops[0]->get_orbs()[2];
      int l = lhs_ops[0]->get_orbs()[0];
      // In parallel calculations not all operators are built on each proc
      if ( ! big.get_op_array(optype).has_local_index(i,j,k,l) ) continue;
//pout << "p" << mpigetrank() << "; i,j,k,l = " << i << "," << j << "," << k << "," << l << endl;
      std::vector<boost::shared_ptr<SparseMatrix> > vec = big.get_op_array(optype).get_element(i,j,k,l);

      // Loop over rhs spin-op components
      for (int jdx=0; jdx < rhs_ops.size(); jdx++) {
        boost::shared_ptr<SparseMatrix>& rhs_op = rhs_ops[jdx];
        assert( rhs_op->get_built() );
        std::string build_123 = rhs_op->get_build_pattern();
        std::vector<SpinQuantum> spin_123 = rhs_op->get_quantum_ladder().at(build_123);

        // Loop over lhs spin-op components //FIXME
        for (int jjdx=0; jjdx < lhs_ops.size(); jjdx++) {
          boost::shared_ptr<SparseMatrix>& lhs_op = lhs_ops[jjdx];
          assert( lhs_op->get_built() );
          std::string build_4 = lhs_op->get_build_pattern();
          std::string build_pattern = "(" + build_123 + build_4 + ")";

          // Allocate and build new operator
          for (int sx=0; sx < vec.size(); sx++) {
            boost::shared_ptr<SparseMatrix>& op = vec[sx];
            // Select relevant spin components
            std::vector<SpinQuantum> s = { op->get_quantum_ladder().at(build_pattern).at(0), op->get_quantum_ladder().at(build_pattern).at(1) };
//            if ( s == spin_123 ) finish_tensor_product( big, rhsBlock, *rhs_op, Transposeview(lhs_op), *op, forwards, build_pattern );
            if ( s == spin_123 ) {
              finish_tensor_product( big, rhsBlock, *rhs_op, *lhs_op, *op, forwards, build_pattern );
              // Renormalise operator
              op->renormalise_transform( rotateMatrix, stateinfo );
            }
          }
        }
      }
      // Store spin-batch on disk 
      if ( ! dmrginp.do_npdm_in_core() ) store_ops_on_disk( ofs, vec );
    }
  }
  if ( rhsifs.is_open() ) rhsifs.close();
}

void do_3index_1_2_tensor_products( bool forwards, const opTypes& optype, const opTypes& rhsType, const opTypes& lhsType,
                                    SpinBlock& big, SpinBlock* rhsBlock, SpinBlock* lhsBlock, std::ofstream& ofs,
                                    const std::vector<Matrix>& rotateMatrix, const StateInfo *stateinfo )
{
  // (i | j,k ) partition
  //-------------------------
  Op_component_base& rhs_array = rhsBlock->get_op_array(rhsType);
  Op_component_base& lhs_array = lhsBlock->get_op_array(lhsType);
  assert ( (rhs_array.get_size() == 1) || (lhs_array.get_size() == 1) );
//pout << "tensor_1_2: rhs_size = " << rhs_array.get_size() << "; op = " << rhs_array.get_op_string() << endl;
//pout << "tensor_1_2: lhs_size = " << lhs_array.get_size() << "; op = " << lhs_array.get_op_string() << endl;

  // Initialize filesystem
  std::ifstream lhsifs;
  if (( ! dmrginp.do_npdm_in_core()) && lhsBlock->size() > 1) lhsifs.open( lhs_array.get_filename().c_str(), std::ios::binary );

  // Loop over all lhs operator indices
  for (int idx = 0; idx < lhs_array.get_size(); ++idx) {
    std::vector<boost::shared_ptr<SparseMatrix> > lhs_ops;
    // Assume 2-index operators are available on this processor in core
    lhs_ops = lhs_array.get_local_element(idx);

    // Loop over all rhs operator indices
    for (int iidx = 0; iidx < rhs_array.get_size(); ++iidx) {
      std::vector<boost::shared_ptr<SparseMatrix> > rhs_ops = rhs_array.get_local_element(iidx);
      int i = rhs_ops[0]->get_orbs()[0];
      int j = lhs_ops[0]->get_orbs()[0]; int k = lhs_ops[0]->get_orbs()[1];
//pout << "i = " << i << endl;
//pout << "j,k = " << j << "," << k << endl;
      // In parallel calculations not all operators are built on each proc
      if ( ! big.get_op_array(optype).has_local_index(i,j,k) ) continue;
//pout << "building i,j,k = " << i << "," << j << "," << k << endl;
      std::vector<boost::shared_ptr<SparseMatrix> > vec = big.get_op_array(optype).get_element(i,j,k);
//pout << "got i,j,k\n";

      // Loop over lhs spin-op components
      for (int jdx=0; jdx < lhs_ops.size(); jdx++) {
        boost::shared_ptr<SparseMatrix>& lhs_op = lhs_ops[jdx];
        assert( lhs_op->get_built() );
        std::string build_23 = lhs_op->get_build_pattern();
//pout << build_23 << endl;
//pout << "getting spin_23\n";
        std::vector<SpinQuantum> spin_23 = lhs_op->get_quantum_ladder().at(build_23);
//pout << "got spin_23\n";

        // Loop over rhs spin-op components
        for (int jjdx=0; jjdx < rhs_ops.size(); jjdx++) {
          boost::shared_ptr<SparseMatrix>& rhs_op = rhs_ops[jjdx];
          assert( rhs_op->get_built() );
          std::string build_1 = rhs_op->get_build_pattern();
          std::string build_pattern = "(" + build_1 + build_23 + ")";

          // Allocate and build new operator
          for (int sx=0; sx < vec.size(); sx++) {
            boost::shared_ptr<SparseMatrix>& op = vec[sx];
            // Select relevant spin component
            std::vector<SpinQuantum> s = { op->get_quantum_ladder().at(build_pattern).at(0) };
            if ( s == spin_23 ) {
              finish_tensor_product( big, rhsBlock, *rhs_op, *lhs_op, *op, forwards, build_pattern );
              // Renormalise operator
              op->renormalise_transform( rotateMatrix, stateinfo );
            }
          }
        }
      }
      // Store spin-batch on disk 
      if ( ! dmrginp.do_npdm_in_core() ) store_ops_on_disk( ofs, vec );
    }
  }
  if ( lhsifs.is_open() ) lhsifs.close();
}

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