C++ Wavefunction::Clear方法代码示例

//solves the equation (H-E)|psi_1> = -QV|\Psi_0>  where lowerState[0] contains |\Psi_0> to enforce orthogonality (Q), and lowerState[1] contains V|\Psi_0> so we can calculate its projection in the krylov space
//the algorithm is taken from wikipedia page
void SpinAdapted::Linear::ConjugateGradient(Wavefunction& xi, DiagonalMatrix& h_diag, double E0, double normtol, Davidson_functor& h_multiply, std::vector<Wavefunction> &lowerStates)
{

    pout.precision (12);
#ifndef SERIAL
    mpi::communicator world;
#endif
    int iter = 0;
    double levelshift = 0.0;

    Wavefunction b = lowerStates[1];
    //make b[0] orthogonal to lowerStates[0]
    double overlap = DotProduct(b, lowerStates[0]);
    double overlap2 = DotProduct(lowerStates[0], lowerStates[0]);
    ScaleAdd(-overlap/overlap2, lowerStates[0], b);

    //normalise all the guess roots
    if(mpigetrank() == 0) {
        Normalise(xi);
    }

    Wavefunction pi, ri;
    ri=xi;
    ri.Clear();
    h_multiply(xi, ri);
    ScaleAdd(-E0, xi, ri); // (H-E0)|psi>

    ScaleAdd(-1.0, lowerStates[1], ri);
    Scale(-1.0, ri);
    pi = ri;

    double oldError = DotProduct(ri, ri);
    while(true) {
        Wavefunction Hp = pi;
        Hp.Clear();
        h_multiply(pi, Hp);
        ScaleAdd(-E0, pi, Hp); // (H-E0)|psi>

        double alpha = oldError/DotProduct(pi, Hp);

        ScaleAdd(alpha, pi, xi);
        ScaleAdd(-alpha, Hp, ri);

        double Error = DotProduct(ri, ri);
        if (Error > normtol)
            return;
        else {
            double beta = Error/oldError;
            oldError = Error;
            ScaleAdd(1.0/beta, ri, pi);
            Scale(beta, pi);
        }
    }
}

//solves the equation (H-E)^T(H-E)|psi_1> = -(H-E)^TQV|\Psi_0>  where lowerState[1] contains |\Psi_0> to enforce orthogonality (Q), and lowerState[0] contains V|\Psi_0> so we can calculate its projection in the krylov space 
//the algorithm is taken from wikipedia page
//the algorithm is taken from wikipedia page
double SpinAdapted::Linear::ConjugateGradient(Wavefunction& xi, double normtol, Davidson_functor& h_multiply, std::vector<Wavefunction>& lowerStates)
{
  setbuf(stdout, NULL);
  p3out.precision (12);
  int iter = 0, maxIter = 100;
  double levelshift = 0.0, overlap2 = 0.0, oldError=0.0, functional=0.0, Error=0.0;

  Wavefunction& targetState = lowerStates[0];
  if (mpigetrank() == 0) {
    for (int i=1; i<lowerStates.size(); i++) {
      overlap2 = pow(DotProduct(lowerStates[i], lowerStates[i]), 0.5);
      if (fabs(overlap2) > NUMERICAL_ZERO) { 
	ScaleAdd(-DotProduct(targetState, lowerStates[i])/overlap2, 
		 lowerStates[i], targetState);
	ScaleAdd(-DotProduct(xi, lowerStates[i])/overlap2, 
		 lowerStates[i], xi);
      }
    }
  }

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

  Wavefunction pi, ri; 
  ri=xi; ri.Clear();
  h_multiply(xi, ri);  

  //Check if we should even perform CG or just exit with a zero vector.
  bool doCG = true;
  if (mpigetrank() == 0) {
    Wavefunction ricopy = ri; ricopy.Clear(); ricopy.Randomise();
    Wavefunction ricopy2 = ricopy;

    for (int i=1; i<lowerStates.size(); i++) {
      overlap2 = pow(DotProduct(lowerStates[i], lowerStates[i]), 0.5);
      if (fabs(overlap2) > NUMERICAL_ZERO) { 
	ScaleAdd(-DotProduct(ricopy, lowerStates[i])/overlap2, 
		 lowerStates[i], ricopy2);
      }
    }

    if (abs(DotProduct(ricopy2, targetState)) < NUMERICAL_ZERO) {
      pout << "The problem is ill posed or the initial guess is very bad "<<DotProduct(ricopy, targetState)<<endl;
      doCG = false;
    }
  }
#ifndef SERIAL
    mpi::broadcast(world, doCG, 0);
#endif
  if (!doCG) {
    xi.Clear();
    int success = 0;

    functional = 0.0;
    return functional;
  }

  if (mpigetrank() == 0) {
    ScaleAdd(-1.0, targetState, ri);
    Scale(-1.0, ri);
    
    for (int i=1; i<lowerStates.size(); i++) {
      overlap2 = pow(DotProduct(lowerStates[i], lowerStates[i]),0.5);
      if (fabs(overlap2) > NUMERICAL_ZERO) { 
	ScaleAdd(-DotProduct(ri, lowerStates[i])/overlap2, 
		 lowerStates[i], ri);
      }
    }
    
    pi = ri;

    oldError = DotProduct(ri, ri);
    printf("\t\t\t %15s  %15s  %15s\n", "iter", "Functional", "Error");
  }

#ifndef SERIAL
    mpi::broadcast(world, Error, 0);
    mpi::broadcast(world, oldError, 0);
    mpi::broadcast(world, functional, 0);
#endif

    if (oldError < normtol) {
      if (mpigetrank() == 0) {
	functional = -DotProduct(xi, ri) - DotProduct(xi, targetState);
	printf("\t\t\t %15i  %15.8e  %15.8e\n", 0, functional, oldError);
      }
#ifndef SERIAL
    mpi::broadcast(world, functional, 0);
#endif
      return functional;
    }

  while(true) {
#ifndef SERIAL
    mpi::communicator world;
//.........这里部分代码省略.........

void SpinAdapted::mps_nevpt::type1::calcHamiltonianAndOverlap(const MPS& statea, double& h, double& o, perturber& pb) {

#ifndef SERIAL
  mpi::communicator world;
#endif


  SpinBlock system, siteblock;
  bool forward = true, restart=false, warmUp = false;
  int leftState=0, rightState=0, forward_starting_size=1, backward_starting_size=1, restartSize =0;
  InitBlocks::InitStartingBlock(system, forward, leftState, rightState, forward_starting_size, backward_starting_size, restartSize, restart, warmUp, 0,statea.getw().get_deltaQuantum(), statea.getw().get_deltaQuantum()); 

  if (dmrginp.outputlevel() > 0)
    pout << system<<endl;
  system.transform_operators(const_cast<std::vector<Matrix>&>(statea.getSiteTensors(0)), 
      		       const_cast<std::vector<Matrix>&>(statea.getSiteTensors(0)), false, false );

  int sys_add = true; bool direct = true; 

  for (int i=0; i<mps_nevpt::sweepIters-1; i++) {
    SpinBlock newSystem;
    system.addAdditionalCompOps();

    InitBlocks::InitNewSystemBlock(system, siteBlocks_noDES[i+1], newSystem, leftState, rightState, sys_add, direct, 0, DISTRIBUTED_STORAGE, false, true,NO_PARTICLE_SPIN_NUMBER_CONSTRAINT,statea.getw().get_deltaQuantum(),statea.getw().get_deltaQuantum());

    newSystem.transform_operators(const_cast<std::vector<Matrix>&>(statea.getSiteTensors(i+1)), 
      			    const_cast<std::vector<Matrix>&>(statea.getSiteTensors(i+1)), false );

    system = newSystem;
  }

  SpinBlock newSystem, big;
  system.addAdditionalCompOps();
  //system.printOperatorSummary();
  //To set implicit_transpose to true.
  //The last site spinblock should have implicit_transpose true.
  InitBlocks::InitNewSystemBlock(system, siteBlocks_noDES[mps_nevpt::sweepIters], newSystem,  leftState, rightState, sys_add, direct, 0, DISTRIBUTED_STORAGE, false, true,NO_PARTICLE_SPIN_NUMBER_CONSTRAINT,statea.getw().get_deltaQuantum(),statea.getw().get_deltaQuantum());
  
  newSystem.set_loopblock(false); system.set_loopblock(false);
  newSystem.addAdditionalCompOps();
  //siteBlocks_noDES[mps_nevpt::sweepIters+1].set_loopblock(false);
  InitBlocks::InitBigBlock(newSystem, siteBlocks_noDES[mps_nevpt::sweepIters+1], big,statea.getw().get_deltaQuantum(),statea.getw().get_deltaQuantum()); 
  

  //FIXME
  //Assume statea.getw() has only one deltaquantum.
  //Spin Embeding for zero order wavefunction is needed.

  Wavefunction temp = statea.getw();
  temp.Clear();
  big.multiplyH(const_cast<Wavefunction&>(statea.getw()), &temp, 1);

  if (mpigetrank() == 0)
    h = DotProduct(statea.getw(), temp);

  temp.Clear();
  big.multiplyOverlap(const_cast<Wavefunction&>(statea.getw()), &temp, 1);
  if (mpigetrank() == 0)
    o = DotProduct(statea.getw(), temp);
  if(dmrginp.spinAdapted())
  {

    double cg= 0.0;
    //TODO
    //Assume the zero order wavefunction has a spin zero. 
    //Spin Embeding must be used.
    SpinQuantum wQ= statea.getw().get_deltaQuantum(0);
    //cg*= dmrginp.get_ninej()(wQ.get_s().getirrep(), 1, dmrginp.effective_molecule_quantum().get_s().getirrep(), 
    //      					   pb.delta.get_s().getirrep(), pb.delta.get_s().getirrep(), 0,
    //                          dmrginp.effective_molecule_quantum().get_s().getirrep(), 0, dmrginp.effective_molecule_quantum().get_s().getirrep());
    //cg*= Symmetry::spatial_ninej(wQ.get_symm().getirrep(), -pb.delta.get_symm().getirrep(), dmrginp.effective_molecule_quantum().get_symm().getirrep(), 
    //      					   pb.delta.get_symm().getirrep(), -pb.delta.get_symm().getirrep(), 0,
    //                          dmrginp.effective_molecule_quantum().get_symm().getirrep(), 0, dmrginp.effective_molecule_quantum().get_symm().getirrep());
    cg += pow(clebsch(wQ.get_s().getirrep(),-1,1,1,0,0),2);
    cg += pow(clebsch(wQ.get_s().getirrep(),1,1,-1,0,0),2);
    //cout << "cg coefficient: " <<cg<<endl;
    h*= cg*cg;
    o*= cg*cg;
  }


#ifndef SERIAL
  mpi::broadcast(world, h, 0);
  mpi::broadcast(world, o, 0);
#endif

  return;
}