C++ Lattice::addSite方法代码示例

int main(int argc, char* argv[]) {
    boost::mpi::environment env(argc,argv);
    boost::mpi::communicator comm;

    print_section("Hubbard nxn");

    int size_x, size_y, wn;
    RealType t, mu, U, beta, reduce_tol, coeff_tol;
    bool calc_gf, calc_2pgf;
    int wf_max, wb_max;
    double eta, hbw, step; // for evaluation of GF on real axis

    try { // command line parser
        TCLAP::CmdLine cmd("Hubbard nxn diag", ' ', "");
        TCLAP::ValueArg<RealType> U_arg("U","U","Value of U",true,10.0,"RealType",cmd);
        TCLAP::ValueArg<RealType> mu_arg("","mu","Global chemical potential",false,0.0,"RealType",cmd);
        TCLAP::ValueArg<RealType> t_arg("t","t","Value of t",false,1.0,"RealType",cmd);

        TCLAP::ValueArg<RealType> beta_arg("b","beta","Inverse temperature",true,100.,"RealType");
        TCLAP::ValueArg<RealType> T_arg("T","T","Temperature",true,0.01,"RealType");
        cmd.xorAdd(beta_arg,T_arg);

        TCLAP::ValueArg<size_t> x_arg("x","x","Size over x",false,2,"int",cmd);
        TCLAP::ValueArg<size_t> y_arg("y","y","Size over y",false,2,"int",cmd);

        TCLAP::ValueArg<size_t> wn_arg("","wf","Number of positive fermionic Matsubara Freqs",false,64,"int",cmd);
        TCLAP::ValueArg<size_t> wb_arg("","wb","Number of positive bosonic Matsubara Freqs",false,1,"int",cmd);
        TCLAP::SwitchArg gf_arg("","calcgf","Calculate Green's functions",cmd, false);
        TCLAP::SwitchArg twopgf_arg("","calc2pgf","Calculate 2-particle Green's functions",cmd, false);
        TCLAP::ValueArg<RealType> reduce_tol_arg("","reducetol","Energy resonance resolution in 2pgf",false,1e-5,"RealType",cmd);
        TCLAP::ValueArg<RealType> coeff_tol_arg("","coefftol","Total weight tolerance",false,1e-12,"RealType",cmd);

        TCLAP::ValueArg<RealType> eta_arg("","eta","Offset from the real axis for Green's function calculation",false,0.05,"RealType",cmd);
        TCLAP::ValueArg<RealType> hbw_arg("D","hbw","Half-bandwidth. Default = U",false,0.0,"RealType",cmd);
        TCLAP::ValueArg<RealType> step_arg("","step","Step on a real axis. Default : 0.01",false,0.01,"RealType",cmd);

        cmd.parse( argc, argv );
        U = U_arg.getValue();
        mu = (mu_arg.isSet()?mu_arg.getValue():U/2);
        boost::tie(t, beta, calc_gf, calc_2pgf, reduce_tol, coeff_tol) = boost::make_tuple( t_arg.getValue(), beta_arg.getValue(),
                                                                                            gf_arg.getValue(), twopgf_arg.getValue(), reduce_tol_arg.getValue(), coeff_tol_arg.getValue());
        boost::tie(size_x, size_y) = boost::make_tuple(x_arg.getValue(), y_arg.getValue());
        boost::tie(wf_max, wb_max) = boost::make_tuple(wn_arg.getValue(), wb_arg.getValue());
        boost::tie(eta, hbw, step) = boost::make_tuple(eta_arg.getValue(), (hbw_arg.isSet()?hbw_arg.getValue():2.*U), step_arg.getValue());
        calc_gf = calc_gf || calc_2pgf;
        calc_gf = calc_gf || calc_2pgf;
    }
    catch (TCLAP::ArgException &e)  // catch parsing exceptions
    { std::cerr << "error: " << e.error() << " for arg " << e.argId() << std::endl; exit(1);}

    int L = size_x*size_y;
    INFO("Diagonalization of " << L << "=" << size_x << "*" << size_y << " sites");
    Lattice Lat;

    /* Add sites */
    std::vector<std::string> names(L);
    for (size_t y=0; y<size_y; y++)
        for (size_t x=0; x<size_x; x++)
        {
            size_t i = SiteIndexF(size_x, x, y);
            std::stringstream s; s << i;
            names[i]="S"+s.str();
            Lat.addSite(new Lattice::Site(names[i],1,2));
        };

    INFO("Sites");
    Lat.printSites();

    /* Add interaction on each site*/
    MelemType U_complex;
    U_complex = std::complex<double> (U,0);
    MelemType mu_complex;
    mu_complex = std::complex<double> (mu,0);
    for (size_t i=0; i<L; i++) LatticePresets::addCoulombS(&Lat, names[i], U_complex, -mu_complex);

    /* Add hopping */
    for (size_t y=0; y<size_y; y++) {
        for (size_t x=0; x<size_x; x++) {
            size_t pos = SiteIndexF(size_x, x,y);
            size_t pos_right = SiteIndexF(size_x, (x+1)%size_x,y); /*if (x == size_x - 1) pos_right = SiteIndexF(0,y); */
            size_t pos_up = SiteIndexF(size_x, x,(y+1)%size_y);

            if (size_x > 1) LatticePresets::addHopping(&Lat, std::min(names[pos], names[pos_right]), std::max(names[pos], names[pos_right]), -t);
            if (size_y > 1) LatticePresets::addHopping(&Lat, std::min(names[pos], names[pos_up]), std::max(names[pos], names[pos_up]), -t);
        };
    };

    int rank = comm.rank();
    if (!rank) {
        INFO("Terms with 2 operators");
        Lat.printTerms(2);

        INFO("Terms with 4 operators");
        Lat.printTerms(4);
    };

    IndexClassification IndexInfo(Lat.getSiteMap());
    IndexInfo.prepare(false); // Create index space
    if (!rank) { print_section("Indices"); IndexInfo.printIndices(); };
    int index_size = IndexInfo.getIndexSize();
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    boost::mpi::environment env(argc,argv);
    boost::mpi::communicator world;

    
    Lattice L;
    L.addSite(new Lattice::Site("A",1,2));
    LatticePresets::addCoulombS(&L, "A", U, -mu);
    L.addSite(new Lattice::Site("B",1,2));
    LatticePresets::addCoulombS(&L, "B", U, -mu);

    LatticePresets::addHopping(&L, "A","B", -1.0);
    INFO("Sites");
    L.printSites();
    INFO("Terms");
    L.printTerms(2);
    INFO("Terms with 4 operators");
    L.printTerms(4);

    IndexClassification IndexInfo(L.getSiteMap());
    IndexInfo.prepare();
    print_section("Indices");
    IndexInfo.printIndices();
    ParticleIndex NModes = IndexInfo.getIndexSize();

    std::vector<ParticleIndex> SpinUpIndices;
    for (ParticleIndex i=0; i<NModes; i++) if (IndexInfo.getInfo(i).Spin) SpinUpIndices.push_back(i);

    print_section("Matrix element storage");
    IndexHamiltonian Storage(&L,IndexInfo);
    Storage.prepare();
    INFO("Terms");
    INFO(Storage);
    
    ParticleIndex IndexSize = IndexInfo.getIndexSize();
    OperatorPresets::Sz Sz(IndexSize, SpinUpIndices);
    INFO("Sz terms");
    Sz.printAllTerms();
    OperatorPresets::N N(IndexSize);
    INFO("N terms");
    N.printAllTerms();
    if (!(Sz.commutes(N))) return EXIT_FAILURE;

    if (!(Storage.commutes(N))) return EXIT_FAILURE;
    INFO("H commutes with N");
    if (!(Storage.commutes(Sz))) return EXIT_FAILURE;
    INFO("H commutes with Sz");
    Symmetrizer Symm(IndexInfo, Storage);
    Symm.compute();

    StatesClassification S(IndexInfo,Symm);
    S.compute();

    Hamiltonian H(IndexInfo, Storage, S);
    H.prepare();
    H.getPart(BlockNumber(4)).print_to_screen();
    H.getPart(BlockNumber(5)).print_to_screen();
    H.compute(world);
    H.getPart(BlockNumber(4)).print_to_screen();
    H.getPart(BlockNumber(5)).print_to_screen();
    INFO("The value of ground energy is " << H.getGroundEnergy());

    RealType beta = 10.0;

    DensityMatrix rho(S,H,beta);
    rho.prepare();
    rho.compute();
    for (QuantumState i=0; i<S.getNumberOfStates(); ++i) INFO(rho.getWeight(i)); 

    FieldOperatorContainer Operators(IndexInfo, S, H);
    Operators.prepareAll();
    Operators.computeAll();
    
    ParticleIndex down_index = IndexInfo.getIndex("A",0,down);

    FieldOperator::BlocksBimap c_map = Operators.getCreationOperator(0).getBlockMapping();
    for (FieldOperator::BlocksBimap::right_const_iterator c_map_it=c_map.right.begin(); c_map_it!=c_map.right.end(); c_map_it++)
        {
            INFO(c_map_it->first << "->" << c_map_it->second);
        }

    GreensFunction GF(S,H,Operators.getAnnihilationOperator(0), Operators.getCreationOperator(0), rho);

    GF.prepare();
    GF.compute();

    ComplexVectorType G_ref(10);
    G_ref << -2.53021005e-01*I,
             -4.62090702e-01*I,
             -4.32482782e-01*I,
             -3.65598615e-01*I,
             -3.07785174e-01*I,
             -2.62894141e-01*I, 
             -2.28274316e-01*I,
             -2.01170772e-01*I,
             -1.79539602e-01*I,
             -1.61950993e-01*I;
 
    bool result = true;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    boost::mpi::environment MpiEnv(argc, argv);
    world = boost::mpi::communicator();

    /* As pomerol is an ED code, it requires a finite-size lattice to be
     * provided. Here is an example of a lattice of 2 sites. */
    // First, we construct an empty lattice
    Lattice L;
    // Add a site with a name "A", that has 1 orbitals and 2 spins.
    L.addSite(new Lattice::Site("A",1,2));
    // Add one more site with a name "B". It also has 1 orbitals and 2 spins.
    L.addSite(new Lattice::Site("B",1,2));

    // Let us now connect two sites with a hopping term, with matrix element -1.
    RealType t=1.0;
    LatticePresets::addHopping(&L, "A","B", -t);
    
    /* Now add interaction. In order to provide some custom interaction, one can
     * give any custom term of 4,6 operators. This is done via 
     * Lattice.addTerm(Lattice::Term) method.
     * We will use Hubbard-type n_{up}n_{down} interaction. For this and some
     * other typical interactions, such as SzSz or SS couplings a shortcut is 
     * provided in the LatticePresets class.
     */
    RealType U = 2.0;
    RealType mu = 1.0;
    // LatticePresets::addCoulombS adds Un_{up}n_{down} for 1 orbital and 2 spins.
    LatticePresets::addCoulombS(&L, "A", U, -mu);
    LatticePresets::addCoulombS(&L, "B", U, -mu);

    // Let us now print which sites and terms are defined.
    if (!world.rank()) {
        INFO("Sites"); // equivalent to std::cout << "Sites" << std::endl;
        L.printSites();
        INFO("Terms");
        L.printTerms(2);
        INFO("Terms with 4 operators");
        L.printTerms(4);

        print_section("Indices");
    };
    /* In order to go further, we need to introduce the index space. An index
     * is a number that uniquely identifies a combination of (site,orbital,spin). 
     * The object that takes care of handling indices is called 
     * IndexClassification.
     */
    
    // Construct IndexClassification
    IndexClassification IndexInfo(L.getSiteMap());
    
    /* Important remark 1!
     * Most of the objects that are defined within Pomerol library have the 
     * following semantics. They can be constructed, prepared and computed. 
     * This means
     *   - constructed: No operations are done except from initializing links to
     *                  other objects that current class depends on.
     *   - prepared: Typically, this is when all memory allocation takes place.
     *   - computed: The actual computation. This is the most costly operation.
     *               When no heavy computation is done, it can be done during
     *               preparation stage.
     */
     
    // IndexClassification does not require much computing time, so everything
    // is calculated during prepare() method.
    IndexInfo.prepare();
    
    // Print which indices we have
    IndexInfo.printIndices();
    // Save the total number of indices.
    ParticleIndex IndexSize = IndexInfo.getIndexSize();

    print_section("Matrix element storage");
    /* We now need to construct the actual Hamiltonian. First, we need to create
     * it in the index space, e.g. write down a formula with all terms. It is
     * useful for a subsequent symmetry analysis. The corresponding class is called
     * an IndexHamiltonian. It is inherited from Operator class, which is designed
     * to represent a generic second-quantized fermionic operator and fully model
     * the algebra of fermionic operators. This means one can multiply, add,
     * subtract it, and also calculate the commutator of two operators, etc.
     */
    
    // First construct the IndexHamiltonian object.
    IndexHamiltonian Storage(&L,IndexInfo);
    // Then prepare it. As previously no long computation required, so everything
    // is done in the prepare() method.
    Storage.prepare();
    // Print out the Hamiltonian.
    INFO("Terms");
    INFO(Storage);

    // Let us make a test, that our Hamiltonian commutes with an operator, that
    // represents the total number of particles. The preset for such operator is
    // defined in OperatorPresets.h and .cpp files.
    OperatorPresets::N N(IndexSize);
    INFO("N terms");
    N.printAllTerms();
    if ((Storage.commutes(N))) INFO("H commutes with N");
    
    /* The Hamiltonian that is set now, has a set of symmetries. One can identify
//.........这里部分代码省略.........