HPhi++  3.1.0
Add new lattice model into Standard mode

Overall procedure

If you want to create a new lattice file, the following procedures are needed.

  1. Copy one of lattice files such as Kagome.cpp (Probably the most similar one) and rename it.
  2. Modify lattice model file
  3. Add the function in the header file, StdFace_ModelUtil.hpp.
  4. Add entry at
    char *fname
    )
    {
    struct StdIntList *StdI;
    :
    /*>>
    Generate Hamiltonian definition files
    */
    if (strcmp(StdI->lattice, "chain") == 0
    || strcmp(StdI->lattice, "chainlattice") == 0) StdFace_Chain(StdI);
    else if (strcmp(StdI->lattice, "face-centeredorthorhombic") == 0
    || strcmp(StdI->lattice, "fcorthorhombic") == 0
    || strcmp(StdI->lattice, "fco") == 0) StdFace_FCOrtho(StdI);
    else if (strcmp(StdI->lattice, "honeycomb") == 0
    || strcmp(StdI->lattice, "honeycomblattice") == 0) StdFace_Honeycomb(StdI);
    else if (strcmp(StdI->lattice, "kagome") == 0
    || strcmp(StdI->lattice, "kagomelattice") == 0) StdFace_Kagome(StdI);
    else if (strcmp(StdI->lattice, "ladder") == 0
    || strcmp(StdI->lattice, "ladderlattice") == 0) StdFace_Ladder(StdI);
    else if (strcmp(StdI->lattice, "orthorhombic") == 0
    || strcmp(StdI->lattice, "simpleorthorhombic") == 0) StdFace_Orthorhombic(StdI);
    else if (strcmp(StdI->lattice, "pyrochlore") == 0) StdFace_Pyrochlore(StdI);
    else if (strcmp(StdI->lattice, "tetragonal") == 0
    || strcmp(StdI->lattice, "tetragonallattice") == 0
    || strcmp(StdI->lattice, "square") == 0
    || strcmp(StdI->lattice, "squarelattice") == 0) StdFace_Tetragonal(StdI);
    else if (strcmp(StdI->lattice, "triangular") == 0
    || strcmp(StdI->lattice, "triangularlattice") == 0) StdFace_Triangular(StdI);
    else if (strcmp(StdI->lattice, "wannier90") == 0) StdFace_Wannier90(StdI);
    else UnsupportedSystem(StdI->model, StdI->lattice);//<<

Modify lattice model file

To create a new lattice file, please modify the following part (Kagome.cpp as an example):

Define function as
struct StdIntList *StdI
)
{
Lattice parameters are used only in geometry.dat and lattice.gp
StdFace_PrintVal_d("a", &StdI->a, 1.0);
StdFace_PrintVal_d("Wlength", &StdI->length[0], StdI->a);
StdFace_PrintVal_d("Llength", &StdI->length[1], StdI->a);
StdFace_PrintVal_d("Wx", &StdI->direct[0][0], StdI->length[0]);
StdFace_PrintVal_d("Wy", &StdI->direct[0][1], 0.0);
StdFace_PrintVal_d("Lx", &StdI->direct[1][0], StdI->length[1] * 0.5);
StdFace_PrintVal_d("Ly", &StdI->direct[1][1], StdI->length[1] * 0.5 * sqrt(3.0));
These are unit lattice vectors.
Just call this function to initialize all lattice related parameters
StdFace_InitSite(StdI, fp, 2);
where "2" indicates 2D.
StdI->tau[0][0] = 0.0; StdI->tau[0][1] = 0.0; StdI->tau[0][2] = 0.0;
StdI->tau[1][0] = 0.5; StdI->tau[1][1] = 0.0; StdI->tau[1][2] = 0.0;
StdI->tau[2][0] = 0.0; StdI->tau[2][1] = 0.5; StdI->tau[2][2] = 0.0;
These are the fractional coordinates of internal sites. Then set parameters of Hamiltonian
StdFace_NotUsed_J("J0", StdI->J0All, StdI->J0);
StdFace_NotUsed_J("J1", StdI->J1All, StdI->J1);
StdFace_NotUsed_J("J2", StdI->J2All, StdI->J2);
StdFace_NotUsed_J("J'", StdI->JpAll, StdI->Jp);
StdFace_NotUsed_J("J0'", StdI->J0pAll, StdI->J0p);
StdFace_NotUsed_J("J1'", StdI->J1pAll, StdI->J1p);
StdFace_NotUsed_J("J2'", StdI->J2pAll, StdI->J2p);
StdFace_NotUsed_d("D", StdI->D[2][2]);
if (strcmp(StdI->model, "hubbard") == 0 ) {
StdFace_NotUsed_i("2S", StdI->S2);
StdFace_NotUsed_J("J", StdI->JAll, StdI->J);
}/*if (strcmp(StdI->model, "hubbard") == 0 )*/
else {
StdFace_PrintVal_i("2S", &StdI->S2, 1);
StdFace_InputSpin(StdI->J, StdI->JAll, "J");
}/*if (model != "hubbard")*/
}/*if (model != "spin")@@*/
to determine the default values of them and unused parameters. For more details, please see the description of each function. Then compute the upper limit of the number of Transfer & Interaction and malloc them.
if (strcmp(StdI->model, "spin") == 0 ) {//>>
ntransMax = StdI->nsite * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/);
nintrMax = StdI->NCell * (StdI->NsiteUC/*D*/ + 6/*J*/ + 6/*J'*/)
* (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1);
}
else {
ntransMax = StdI->NCell * 2/*spin*/ * (2 * StdI->NsiteUC/*mu+h+Gamma*/ + 12/*t*/ + 12/*t'*/);
nintrMax = StdI->NCell * (StdI->NsiteUC/*U*/ + 4 * (6/*V*/ + 6/*V'*/));
if (strcmp(StdI->model, "kondo") == 0) {
ntransMax += StdI->nsite / 2 * (StdI->S2 + 1/*h*/ + 2 * StdI->S2/*Gamma*/);
nintrMax += StdI->nsite / 2 * (3 * StdI->S2 + 1) * (3 * StdI->S2 + 1);
}/*if (strcmp(StdI->model, "kondo") == 0)*/
}//<<
Please estimate the number of bonds per site.
for (kCell = 0; kCell < StdI->NCell; kCell++) {
In this loop, the parameters of Hamiltonian are computed & stored. The local term is computed as follows:
/*>>
Local term
*/
isite = StdI->NsiteUC * kCell;
if (strcmp(StdI->model, "kondo") == 0 ) isite += StdI->nsite / 2;
if (strcmp(StdI->model, "spin") == 0 ) {
for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) {
StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, isite + isiteUC);
StdFace_GeneralJ(StdI, StdI->D, StdI->S2, StdI->S2, isite + isiteUC, isite + isiteUC);
}/*for (jsite = 0; jsite < StdI->NsiteUC; jsite++)*/
}/*if (strcmp(StdI->model, "spin") == 0 )*/
else {
for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)
StdFace_HubbardLocal(StdI, StdI->mu, -StdI->h, -StdI->Gamma, StdI->U, isite + isiteUC);
if (strcmp(StdI->model, "kondo") == 0 ) {
jsite = StdI->NsiteUC * kCell;
for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++) {
StdFace_GeneralJ(StdI, StdI->J, 1, StdI->S2, isite + isiteUC, jsite + isiteUC);
StdFace_MagField(StdI, StdI->S2, -StdI->h, -StdI->Gamma, jsite + isiteUC);
}/*for (isiteUC = 0; isiteUC < StdI->NsiteUC; isiteUC++)*/
}/*if (strcmp(StdI->model, "kondo") == 0 )*/
}/*if (strcmp(StdI->model, "spin") != 0 )<<*/
Probably, it is not necessary to modify this part. The non-local term is as follows:
/*>>
Nearest neighbor intra cell 0 -> 1
*/
StdFace_SetLabel(StdI, fp, iW, iL, 0, 0, 0, 1, &isite, &jsite, 1, &Cphase, dR);
if (strcmp(StdI->model, "spin") == 0 ) {
StdFace_GeneralJ(StdI, StdI->J2, StdI->S2, StdI->S2, isite, jsite);
}/*if (strcmp(StdI->model, "spin") == 0 )*/
else {
StdFace_Hopping(StdI, Cphase * StdI->t2, isite, jsite, dR);
StdFace_Coulomb(StdI, StdI->V2, isite, jsite);
}//<<
For more details, please see each function.