Code for MC simulation

Simulation code for Monte Carlo approach

This approach uses a conventional Monte Carlo simulation with a Metropolis step to simulate the $2D$ $\mathcal{O}(3)$ nonlinear sigma model where the topological term $\theta$ is assigned an imaginary value.

Code documentation: nonlinearsigma.cpp

** Coming soon **

Code documentation: lattice.cpp

** last updated November 8, 2023 **

Go through and figure out which - if any - of these functions are never referenced. Move them to a testing class or remove them entirely

Constructor UPDATE THIS WITH CHANGES FROM 11/29

    Lattice::Lattice(int length, double beta, double itheta){
        Lattice::setLength(length); //set length
        Lattice::setBeta(beta); //set beta
        Lattice::setiTheta(itheta); //set itheta
        Lattice::setnTherm(1000); //set therm steps to default number
        Lattice::setnMC(1000); //set Monte Carlo steps to default number
        Lattice::setFreq(100); //set frequency between saved configs to default number
        Lattice::generateFilename_();
        fixedr_ = false; //this should only be set to true when testing
        use_arcsin_ = true;
        maxExc_ = 0;
    }

The constructor takes in the length of the lattice (we are using a square lattice, so length is used for x and y lengths), the value of beta, and the value of itheta. It then sets those three parameters, sets a default number of thermalization steps, Monte Carlo steps, and frequency for saving the lattice configurations and observables.

The functions setLength, setBeta, setiTheta, setnTherm, setnMC, and setFreq exist to prevent a segmentation fault when running this code. They provide a way into the memory where these parameters are stored, in order to update and save them properly.

Something to consider: should you take in the thermalization steps, MC steps, and freq at this point? What purpose does waiting on those have, or is this a holdover from early testing activities?

The next thing the constructor does is generate the filename, which is described in more detail in this section. This way, as soon as the lattice object is constructed, a filename exists for the output.

Then it sets fixedr_ to false – this means that we will use the random number generator as usual instead of forcing a pre-set value for the two random numbers (we would only want to fix those numbers in very rare testing conditions). This is currently used in fixRNG, freeRNG, and exceptionalConfig.

What is this doing in exceptionalConfig? Do we need it there?

Next, it sets our preference for whether to use arccos or arcsin to calculate QL on the triangles – the default is to use arcsin, but we can change that setting with the function setTrig

And finally, it sets the counter for max number of exceptional configurations removed (maxExc_) to 0. This number is updated every time we use the removeExceptional function and is output with the configurations in saveConfig so we can track how many times the code has to iterate to get a clean configuration.

various set param functions

Each of the functions below allows us to change one of the parameters given for the simulation. These functions are useful when testing the code, as it allows us to update the default values of each of the simulation parameters.

Each of these functions has the same basic steps, which are to update the stored attribute in the object and then generate a new filename (as the parameters are listed in the filename). In addition, setLength also initializes the lattice, since the length of the lattice determines how many spins are on the lattice and the whole lattice must be reset when the length changes.

void Lattice::setLength(int length){
        //tested 6/1/2023
        length_ = length;
        Lattice::generateFilename_();
        Lattice::initialize();
    }

    void Lattice::setBeta(double beta){
        //tested 6/1/2023
        beta_ = beta;
        Lattice::generateFilename_();
    }

    void Lattice::setiTheta(double itheta){
        //tested 6/1/2023
        itheta_ = itheta;
        Lattice::generateFilename_();
    }

    void Lattice::setnTherm(int ntherm){
        nTherm_ = ntherm;
        Lattice::generateFilename_();
    }

    void Lattice::setnMC(int nMC){
        nMC_ = nMC;
        Lattice::generateFilename_();
    }

    void Lattice::setFreq(int freq){
        freq_ = freq;
        Lattice::generateFilename_();
    }

setTrig

This function allows us to change our default use of arcsin to arccos by passing use_arcsin = false into this function. We use this when performing tests on QL.

    void Lattice::setTrig(bool use_arcsin){
        use_arcsin_ = use_arcsin;
    }

fixRNG

This function allows us to set our random numbers to pre-specified values. This is only used when testing to remove the randomness inherent in the model.

    void Lattice::fixRNG(double r1, double r2){
        //tested 6/1/2023
        fixedr_ = true;
        r1_ = r1;
        r2_ = r2;
    }

freeRNG

This allows us to undo fixing the random numbers, should we wish to perform a test afterwards that requires a random number.

    void Lattice::freeRNG(){
        //tested 6/1/2023
        fixedr_ = false;
    }

various get parameter functions

The below functions all perform the same operation: fetching the internal attribute.

Investigate whether this is the proper practice – I recall it helps prevent seg faults, but also am able to use length_ directly in the for loops later? Is that an internal v external access issue?

    int Lattice::getLength(){
        //tested 6/1/2023
        return length_;
    }
    
    double Lattice::getBeta(){
        //tested 6/1/2023
        return beta_;
    }
    
    double Lattice::getiTheta(){
        //tested 6/1/2023
        return itheta_;
    }
    
    int Lattice::getnTherm(){
       return nTherm_; 
    }
    
    int Lattice::getnMC(){
        return nMC_;
    }
    
    std::string Lattice::getFilename(){
        return filename_;
    }

getFilename

This fetches the internal class attribute of the filename, which is generated in generateFilename. See this section for more on why we use getFilename instead of referencing the attribute directly.

std::string Lattice::getFilename(){
        return filename_;
    }

getPhi

This fetches the three-component vector phi at site (i,j). See this section for more on why we use getPhi instead of referencing the attribute directly.

   
    Lattice::field Lattice::getPhi(int i, int j){
        //tested 5/30/2023
        return grid_[i][j];
    }

getRandNums

This fetches the two random numbers most recently used to generate our field vectors.

Consider removing this or moving it to a testing class if you make one

    double* Lattice::getRandNums(){
        //tested 6/1/2023
        static double r[2] = {r1_, r2_};
        return r;
    }

getPhiTot

This calculates the total magnitude of phi on the lattice.

Currently this is one of our outputs for our observables, but is used nowhere else inside the code. Consider removing this or moving it to a testing class if you make one

    double Lattice::getPhiTot(){
        //tested 6/1/2023
        //double phi_tot = 0.;
        double phi_tot(0.);//optimization 7/4/23
        #pragma omp parallel for collapse(2) default(none) shared(length_) reduction(+:phi_tot)
        for(int i = 0; i < length_; i++){
            for(int j = 0; j < length_; j++){
                field phi(Lattice::getPhi(i,j)); 
                phi_tot += dot(phi,phi);
            }
        }
        return phi_tot;
    }

getAvgG

This fetches the internal attribute which is the average correlation function at site (i,j).

    double Lattice::getAvgG(int i, int j){
        //double Gij = Gij_[i][j];
        double Gij(Gij_[i][j]);//optimization 7/4/23
        return Gij;
    }

removeExceptional UPDATE WITH CHANGES FROM 11/29

This function starts with the condition that the lattice is not clean (exceptional_config = true), and begins a count of exceptional configurations.

It starts a while loop, which only exits when exceptional_config becomes false or we reach a certain number of attempts (set by exc_lim, which is passed into the function from outside). This loop checks if the lattice site we are on violates the conditions required for regularization. If either condition is violated, the configuration is exceptional, and the boolean is updated to reflect that, the count of exceptional configirations is incremented, and a new field is generated at our location. If the configuration is valid, the boolean is updated to false, thereby exiting our while loop.

At the end, it checks to see if the number of exceptional configurations the code had to throw away is larger than the previous amount saved. If so, it updates that number.

void Lattice::removeExceptional(int i, int j, int exc_lim){
        bool exceptional_config = true;
        int exc_count = 0;
        while (exceptional_config){
            if (Lattice::exceptionalConfig(i,j,0) or Lattice::exceptionalConfig(i,j,1)){
                exceptional_config = true;
                exc_count++;
                //update lattice
                field phi_new = Lattice::makePhi_();
                grid_[i][j][0] = phi_new[0];
                grid_[i][j][1] = phi_new[1];
                grid_[i][j][2] = phi_new[2]; 
            }
            else{
                exceptional_config = false;
            }
            if (exc_count >= exc_lim){
                break;
            }
        }//while still exceptional at i,j
        //std::cout << "Num attempts at non-exceptional at site (i,j) = ";
        //std::cout << i << "," << j << " was " << exc_count << std::endl;
        if (exc_count > maxExc_){
            maxExc_ = exc_count;//update number of exceptional configs removed
        }
    }

clean

This function creates a shuffled array of all the lattice sites, so it moves randomly on the lattice.

It sets a limit on number of attempts at fixing exceptional configurations, and then loops over the entire lattice (using the randomized array). It also resets our max number of exceptional configurations to zero, so we track it for each configuration and not over all iterations.

At each site, it checks whether the triangles violate the conditions that allow us to generate integer values for QL.

If either condition is violated, we mark that configuration as exceptional and generate a new field at the site.

We then check again, continuing until we find a non-exceptional configuration or we reach our limit.

We then print to the console how many attempts were required to find a non-exceptional configuration. This is just in hopes that we might be able to figure out where things are going particularly poorly.

This function should not be used during the Metropolis step, but only during the initialization step, as the Metropolis step has additional conditions for accepting or rejecting a change to the configuration.

To do: make a removeExceptional function that modularizes that while loop.

     void Lattice::clean(){
        //cleaning the lattice means removing exceptional configurations
        int nsites(length_*length_); 
        std::vector<int> site_arr(nsites);
        std::iota(site_arr.begin(), site_arr.end(), 0);     
        shuffle(site_arr.begin(), site_arr.end(), std::default_random_engine(1232));
        
        int exc_lim(1000);
        maxExc_ = 0;

        for(unsigned int n = 0; n < site_arr.size(); n++){
            int i(site_arr[n]/length_);
            int j(site_arr[n]%length_);
            Lattice::removeExceptional(i, j, exc_lim);
        }//loop over sites
    }

initialize

This function first initializes a 2D grid (using C++ vectors), where each entry in our grid is of the custom field type (a 3D array of doubles, representing our 3 components of $\phi$ ).

It then initializes a similar 2D grid, but each entry is just a double. This will store the value of the correlation function at that point.

Next, it loops over i and j (x and y), and at each site (i,j), it generates a 3-component field phi using the makePhi function and saves that field into the 2D grid.

It also initializes the correlation function to 0 at each site (i,j), as the correlation function cannot be calculated yet.

After the loops over i and j are complete, the lattice grid of phi vectors and the lattice grid of correlation function, G, are saved to their respective internal attributes (grid_ and Gij_).

Next, the triangles that we use to calculate $Q_L$ are generated, and their coordinates saved. This happens in the makeTriangles function.

This is the current version of this function:

void Lattice::initialize(){
        //tested 5/30/2023
        std::vector < std::vector < Lattice::field > > grid;
        std::vector < std::vector < double > > Gij;
        for(int i = 0; i < length_; i++){
            std::vector <double> Gj;
            std::vector < Lattice::field > gridj;
            for (int j = 0; j<length_; j++){
                field phi = Lattice::makePhi_();
                gridj.push_back(phi);
                Gj.push_back(0.);
            }
            Gij.push_back(Gj);
            grid.push_back(gridj);
        }
        grid_ = grid;
        Gij_ = Gij;
        Lattice::makeTriangles_();
        Lattice::zeroCount();
    }

The function clean() follows this, which goes through the entire lattice and checks for exceptional configurations and tries to remove them. We should consider adding that to initialization itself, and writing a modular function that checks for exceptional configurations that we can use in both the clean() function and in the metropolis step.

Consider also calculating Gij here? You initialize everything to 0, but it’s not actually 0.

makeTriangles

generateFilename

exceptionalConfig

Code documentation: mathlib.cpp

** coming soon **

Running the simulation code

Makefile flags

The code is written in C++ and OpenMP and can be compiled with numerous flags.

USE_OMP ?= TRUE
USE_GPROF ?= FALSE
USE_TEST_PRINT_STATEMENTS ?= FALSE
USE_EXTREME_TEST_CONDITION ?= FALSE
USE_CHECK_QL_COS ?= FALSE
USE_CONST_RN ?= FALSE

The first flag “USE_OMP” toggles whether to implement the parallelization in the code. It should be set to “TRUE” if you want to run the simulation in parallel. This is highly recommended for large lattices as the scaling is very poor in series.

If you wish to profile the code, set the second flag “USE_GPROF” to “TRUE”. This sets the correct compiler flags so that you can generate the profiling output. To view the output after the code has run, go to the directory in which you have the executable and run the command

gprof -l nonlinearsigma gmon.out > profiling_results.txt

You will then be able to see the profiling report. In general, you should set “USE_GPROF” to false, unless you are looking to optimize the code or troubleshoot it.

The flag “USE_TEST_PRINT_STATEMENTS” activates print statements throughout the code. This is useful for debugging, but should generally be set to FALSE as it slows down the code.

If you run into major problems, set “USE_EXTREME_TEST_CONDITION” to TRUE. This will run a testing suite built into the code, but will not run the usual simulation. This can help you identify problems in the code, and the testing suite is a function inside the main function, which can be modified as needed to add more tests.

To switch from using arcsin to calculate $Q_{L}$ to using arccos, set “USE_CHECK_QL_COS” to TRUE. In general, this should be set to FALSE< as we want to use arcsin due to its useful symmetry.

Finally, if you want to remove the random number generation and use a constant value for the random numbers, set “USE_CONST_RN” to true. This should only be done when testing the code.