MC_heisenbergEquilibriumQMC.hpp

//
//  MC_heisenbergEquilibriumQMC.hpp
//  
//
//  Created by Marton Kanasz-Nagy on 10/28/16.
//
//

#ifndef MC_heisenbergEquilibriumQMC_hpp
#define MC_heisenbergEquilibriumQMC_hpp

#include <stdio.h>
#include <iostream>
#include <vector>
#include "MC_random.hpp"
#include "MC_finiteTemperatureEquilibriumQMC.hpp"

using namespace std;

// --------------------------------------------------------------------------------------------
// Child class of squareLatticeFiniteTemperatureEquilibriumQMC implementing an equilibrium
// quantum Monte Carlo calculation of the SU(2) symmetric AFM quantum Heisenberg model, without
// magnetic field. Lattice size is set in MC_coord.hpp by the global variable LATTICE_SIZE.
// --------------------------------------------------------------------------------------------
// This algorithm was written using the tutorials on Anders Sandwik's website and his papers.
// http://physics.bu.edu/~sandvik/vietri/sse.pdf
// Since most details are on the website, I will only write short comments to the code.
// ------------------------------------------------------------------------------------------
class heisenbergEquilibriumQMC : public squareLatticeFiniteTemperatureEquilibriumQMC {
protected:
    double T;                       // temperature
    double J;                       // spin coupling (positive since we assume antiferromagnetic ordering)
    
private:
    int hamiltonianPower;           // power of the Hamiltonian in its Taylor expansion
    int hamiltonianPowerMax;        // maximal Hamiltonian power used so far (to update hamiltonianPowerCutoff)
    int hamiltonianPowerCutoff;     // maximal allowed power of the Hamiltonian
    vector<int> oparray;            // array containing which operator we have at time p = 0...L-1
                                    // oparray[p] == 0: no operator
                                    // oparray[p] >  0: oparray[p] = a(p) + 2*bondindex(p) + 1
                                    // denoted s(p) in Sandwik's notes
                                    // Sandwik uses s(p) = a(p) + 2*bondindex(p) - 1, but in our case, the sign
                                    // shoud be +1, since we number the bonds from 0, whereas he numbers them
                                    // starting from 1
    vector<double> pInsertDiagonal; // pInsertDiagonal[n] = probability of accepting a diagonal
                                    // insertion update when the Hamiltonian is on the power n
    vector<double> pDeleteDiagonal; // pDeleteDiagonal[n] = probability of accepting a diagonal
                                    // deletion update when the Hamiltonian is on the power n
    vector<int> linkedVertexList;   // auxiliary array containing the connectivity of operator vertex legs.
                                    // Has as many elements as many vertex legs we have.
                                    // Note: v == linkedVertexList[linkedVertexList[v]]
    vector<int> vfirst;             // auxiliary array: vfirst[i] == the first operator vertex leg at site i
    vector<int> vlast;              // auxiliary array: vlast[i] == the last operator vertex leg at site i
    
public:
    // constructor and destructor
    heisenbergEquilibriumQMC(double Jval=1., double Tval=100.);
    ~heisenbergEquilibriumQMC();
    
    // reset data and reinitialize simulation
    void reset(double Jval, double Tval);
    
    // QMC routines
    void thermalize(int nBurnInSteps);
    void run(int nMonteCarloSteps);
    
    // print out all data, using the operator << of the parent class as well
    friend ostream& operator<<(ostream &, const heisenbergEquilibriumQMC &);
    
private:
    // initialize the private variables
    void initialize(void);
    
    // insertion or deletion of an S_z(i)*S_z(j) diagonal operator, see Sandwik's notes
    void diagonalUpdate(void);
    
    // flip all spins and operators in each operator loop with probability 1/2, see Sandwik's notes
    void operatorLoopUpdate(void);
    
    // auxiliary function to update the values of pInsertDiagonal and pDeleteDiagonal
    // needed when the maximal cutoff of powers of the Hamiltonian is changed
    void calculateDiagonalUpdateProbabilities(void);
    
    // auxiliary routine to update L during the thermalization process, if nmax becomes too large
    // According to Sandwik, the ideal cutoff is L = nmax + nmax/3
    void updatePowerCutoff(void);
    
    // update average energy and magnetic susceptibility
    void updateAvgEnergy(void);
    void updateMagneticSusceptibility(void);
    void updateSpecificHeat(void);
};


#endif /* MC_heisenbergEquilibriumQMC_hpp */