Hamiltonians.hpp

#ifndef __HAMILTONIANS__
#define __HAMILTONIANS__

#include "petscsys.h"
#include "DMRGBlock.hpp"
#include <vector>
#include <iostream>

namespace Hamiltonians
{
    struct Term
    {
        PetscScalar a;      
        Op_t        Iop;    
        PetscInt    Isite;  
        Op_t        Jop;    
        PetscInt    Jsite;  
    };

    typedef enum
    {
        OpenBC=0,
        PeriodicBC=1
    } BC_t;

    class J1J2XXZModel_SquareLattice
    {

    public:

        J1J2XXZModel_SquareLattice()
        {

        }

        PetscErrorCode SetFromOptions()
        {
            PetscErrorCode ierr;
            /* Get values from command line options */
            ierr = PetscOptionsGetReal(NULL,NULL,"-J1",&_J1,NULL); CHKERRQ(ierr);
            ierr = PetscOptionsGetReal(NULL,NULL,"-J2",&_J2,NULL); CHKERRQ(ierr);
            ierr = PetscOptionsGetReal(NULL,NULL,"-Jz1",&_Jz1,NULL); CHKERRQ(ierr);
            ierr = PetscOptionsGetReal(NULL,NULL,"-Jz2",&_Jz2,NULL); CHKERRQ(ierr);
            ierr = PetscOptionsGetInt(NULL,NULL,"-Lx",&_Lx,NULL); CHKERRQ(ierr);
            ierr = PetscOptionsGetInt(NULL,NULL,"-Ly",&_Ly,NULL); CHKERRQ(ierr);
            ierr = PetscOptionsGetBool(NULL,NULL,"-verbose",&verbose,NULL); CHKERRQ(ierr);

            /* One may also specify whether to do the NN Heisenberg model and specify only the anisotropy */
            ierr = PetscOptionsGetReal(NULL,NULL,"-heisenberg",&_Jz1,&heisenberg); CHKERRQ(ierr);
            if(heisenberg)
            {
                _J1 = 0.50;
                _J2 = 0.0;
                _Jz2= 0.0;
            }

            PetscBool BCopen = PETSC_FALSE;
            ierr = PetscOptionsGetBool(NULL,NULL,"-BCopen",&BCopen,NULL); CHKERRQ(ierr);
            if(BCopen) { _BCx=OpenBC; _BCy=OpenBC; }
            PetscBool BCperiodic = PETSC_FALSE;
            ierr = PetscOptionsGetBool(NULL,NULL,"-BCperiodic",&BCperiodic,NULL); CHKERRQ(ierr);
            if(BCperiodic) { _BCx=PeriodicBC; _BCy=PeriodicBC; }
            set_from_options = PETSC_TRUE;
            return(0);
        }

        PetscErrorCode SaveAsOptions(const std::string& filename)
        {
            PetscErrorCode ierr;
            char val[4096];
            PetscBool set;
            std::vector< std::string > keys = {
                "-J1",
                "-J2",
                "-Jz1",
                "-Jz2",
                "-Lx",
                "-Ly",
                "-heisenberg",
                "-BCopen",
                "-BCperiodic"
            };

            std::ostringstream oss;
            for(std::string& key: keys) {
                ierr = PetscOptionsGetString(NULL,NULL,key.c_str(),val,4096,&set); CHKERRQ(ierr);
                if(set) {
                    std::string val_str(val);
                    if(val_str.empty()) val_str = "yes";
                    oss << key << " " << val_str << std::endl;
                }
            }

            FILE *fp;
            ierr = PetscFOpen(PETSC_COMM_SELF,filename.c_str(),"w", &fp); CHKERRQ(ierr);
            ierr = PetscFPrintf(PETSC_COMM_SELF,fp,"%s",oss.str().c_str()); CHKERRQ(ierr);
            ierr = PetscFClose(PETSC_COMM_SELF,fp); CHKERRQ(ierr);

            return(0);
        }

        PetscInt NumSites() const { return _Lx*_Ly; }

        PetscInt NumEnvSites() const { return _Ly; }

        std::vector< Term > H(
            const PetscInt& nsites  
            );

        void PrintOut() const {
            if(heisenberg)
            {
                printf( "HAMILTONIAN: HeisenbergModel_SquareLattice\n");
            }
            else
            {
                printf( "HAMILTONIAN: J1J2XXZModel_SquareLattice\n");
            }
            printf( "  Lx  : %lld\n", LLD(_Lx));
            printf( "  Ly  : %lld\n", LLD(_Ly));
            printf( "  J1  : %g\n", _J1);
            printf( "  Jz1 : %g\n", _Jz1);
            printf( "  J2  : %g\n", _J2);
            printf( "  Jz2 : %g\n", _Jz2);
            printf( "  BCx : %s\n", _BCx?"Periodic":"Open");
            printf( "  BCy : %s\n", _BCy?"Periodic":"Open");
        }

        void SaveOut(FILE *fp) const {
            fprintf(fp, "  \"Hamiltonian\": {\n");
            if(heisenberg)
            {
                fprintf(fp, "    \"label\":\"HeisenbergModel_SquareLattice\",\n");
            }
            else
            {
                fprintf(fp, "    \"label\":\"J1J2XXZModel_SquareLattice\",\n");
            }
            fprintf(fp, "    \"parameters\": {\n");
            fprintf(fp, "      \"Lx\"  : %lld,\n", LLD(_Lx));
            fprintf(fp, "      \"Ly\"  : %lld,\n", LLD(_Ly));
            fprintf(fp, "      \"J1\"  : %g,\n", _J1);
            fprintf(fp, "      \"Jz1\" : %g,\n", _Jz1);
            fprintf(fp, "      \"J2\"  : %g,\n", _J2);
            fprintf(fp, "      \"Jz2\" : %g,\n", _Jz2);
            fprintf(fp, "      \"BCx\" : \"%s\",\n", _BCx?"Periodic":"Open");
            fprintf(fp, "      \"BCy\" : \"%s\"\n", _BCy?"Periodic":"Open");
            fprintf(fp, "    }\n");
            fprintf(fp, "  }");
            fflush(fp);
        }

        PetscInt Lx() const { return _Lx; }

        PetscInt Ly() const { return _Ly; }

        PetscInt To1D(
            const PetscInt ix,
            const PetscInt jy
            ) const;

        PetscErrorCode To2D(
            const PetscInt idx,
            PetscInt& ix,
            PetscInt& jy
            ) const;

        std::vector< std::vector< PetscInt > > NeighborPairs(const PetscInt d=1) const;

    private:

        PetscBool heisenberg = PETSC_FALSE;

        PetscBool set_from_options = PETSC_FALSE;

        PetscScalar _Jz1 = 0.0;

        PetscScalar _Jz2 = 0.0;

        PetscScalar _J1 = 1.0;

        PetscScalar _J2 = 1.0;

        PetscInt _Lx = 4;

        PetscInt _Ly = 4;

        PetscBool verbose = PETSC_FALSE;

        BC_t _BCx = OpenBC;

        BC_t _BCy = PeriodicBC;

        PetscErrorCode PrintInfo() const
        {
            SETERRQ1(PETSC_COMM_SELF,1,"Function %s not implemented.", __FUNCTION__);
        }

        std::vector<PetscInt> GetNearestNeighbors(
            const PetscInt& ix, const PetscInt& jy, const PetscInt& nsites_in
            ) const;

        std::vector<PetscInt> GetNextNearestNeighbors(
            const PetscInt& ix, const PetscInt& jy, const PetscInt& nsites_in
            ) const;

        std::vector< Term > H_full;

        PetscBool H_full_filled = PETSC_FALSE;
    };
}

#endif // __HAMILTONIANS__