Pipeline for data generation and analysis

The basic workflow

The basic workflow is as follows:

compile the code
write a SLURM script to schedule the job
submit the SLURM script
copy the data to Etna
in Metis, use the Python analysis notebooks to look at the data and make relevant plots

Compiling the MC code

Once you have set the flags you want to use, you can compile the code by typing

make -f make_sigma

It can be useful to run

make -f make_sigma clean

first to clear out old .o files that might not be updated otherwise.

Once you have compiled the code, you should have an executable by the name of nonlinearsigma which you will use to run the simulation.

Running the code

Other important things to note before running the code

The simulation requires a list of parameters. We give those to the code through the text file inputs.txt. That file consists of a list of parameter keywords and a value. Be very careful not to change the format of this document, only the numbers, otherwise the simulation won’t be able to understand what values get assigned to what parameters.

The file inputs.txt should look like this

L = 10
beta = 1.6
itheta = 1.
ntherm = 1000
nMC = 1000
freq = 100

Where L gives the length of the square lattice, beta is $\beta = 1/g_{L}$ and should be set to 1.6, itheta is the imaginary value given for the topological term and is given in fractions of $\pi$ (so you should enter 0.5 if you want $i \theta = \pi/2$ ), ntherm is the number of steps you want the simulation to take for thermalization, nMC is the number of steps in the Monte Carlo loop after thermalization, and freq sets the number of steps between saved configurations.

Running the code – single job in interactive node

If you want to run this in an interactive node to test, you can request an interactive node with the following command:

salloc --cpus-per-task=1 --time=00:30:00

where you can adjust the number of cpus per task if you want to run the parallelized code and the time requested is in format hh:mm:ss

You can run the job using the command

./nonlinearsigma inputs.txt

nonlinearsigma is the name of the executable created when you compiled the code, and inputs.txt is the list of inputs mentioned above.

Running the code – single SLURM submission

Once you’re comfortable with the code and want to submit a job that’s longer than your interactive session (or one that may need to run overnight, etc), you do that by writing a SLURM script and sending that script to the scheduler.

The script is called submit_sigma.sh and looks like this:

#SBATCH --job-name=nonlinearsigma_omp_test           # Job name
#SBATCH --mail-type=ALL                              # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=cberger3@bates.edu                # Where to send mail
#SBATCH --partition=faculty                             # Which partition to use
#SBATCH --nodes=1                                    # Number of nodes
#SBATCH --cpus-per-task=1                           # Number of threads per task (OpenMP)
#SBATCH --mem=1gb                                    # Job memory request
##SBATCH --time=05:00:00                             # Time limit hrs:min:sec
#SBATCH --output=nonlinearsigma_omp_test_%j.log      # Standard output 
#SBATCH --error=err_nonlinearsigma_omp_test_%j.log   # Standard output and error log

pwd; hostname; date

export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK

echo "Running nonlinear sigma on single CPU core"

/usr/bin/time -v ./nonlinearsigma inputs.txt

date

Let’s walk through what this script does.

The first line:

#SBATCH --job-name=nonlinearsigma_omp_test           # Job name

Assigns a name to the job, which can help you keep track of what is running in the queue. I tend to name all my jobs in a certain phase of the work the same thing (e.g. “nonlinearsigma_small_L_tests” if I’m testing out the script on small lattices or “nonlinearsigma_first_production_run” if I’m starting to take data for real). This is just for your own information, so name it whatever you’d like.

The next two lines

#SBATCH --mail-type=ALL                              # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=cberger3@bates.edu                # Where to send mail

tell SLURM who to email and when to email you. I have it set to email me at my Smith email anytime a job starts, ends, or fails. I find this helpful, especially if trouble arises, but it’s your choice what to put here (just don’t leave my email address in!).

The next line

#SBATCH --partition=faculty                             # Which partition to use

sends the jobs to the partition that belongs to the physics department. We have priority on this node, but we can request time on other nodes if we really need to. We’d need to talk to CATS about that if we wanted to do it.

The next lines specify the code’s needs:

#SBATCH --nodes=1                                    # Number of nodes
#SBATCH --cpus-per-task=1                           # Number of threads per task (OpenMP)
#SBATCH --mem=1gb                                    # Job memory request

Since we are using OpenMP for parallelization, we only need one node, but we will want to change bash --cpus-per-task to something larger for OpenMP. Unless doing a very large lattice (L > 100), I tend to ask for 30 CPUs per task, which allows for 2 jobs to run at a time on each node. If you’re doing something large, you might want to consider asking for 60 CPUs per task, which will occupy an entire node.

This line

##SBATCH --time=05:00:00                             # Time limit hrs:min:sec

sets a time limit – it will cut off your code when that limit is reached, whether it is done or not. Note there is an extra # here – that means it’s commented out, so it will not have a time limit. I tend to only use the time limits when testing.

These next two lines

#SBATCH --output=nonlinearsigma_omp_test_%j.log      # Standard output 
#SBATCH --error=err_nonlinearsigma_omp_test_%j.log   # Standard output and error log

give names to the output and error files. These files are created when the code runs and can help you debug if things go wrong.

All the above were instructions for SLURM, which schedules jobs on the machine’s available resources. Now we get into the commands to actually run the code.

pwd; hostname; date

This just prints where the job is being run from and the date and time before the job starts.

export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK

this sets the number of threads in OpenMP equal to the number of CPUs per task you assigned above.

echo "Running nonlinear sigma"

/usr/bin/time -v ./nonlinearsigma inputs.txt

This prints out an annoucement that you’re starting the job, which is useful in the logfile, and then runs the simulation with the inputs file. Make sure you have the inputs file and the executable in the same folder as the SLURM script so the computer can find them.

And finally

date

we print the datetime stamp at the end of the simulation.

Running the code – batch SLURM submissions

Each inputs file is one set of parameters, and we need to get lots of data. If you’re not interested in manually setting up these SLURM scripts, input files, etc by hand, I don’t blame you. That’s why I wrote a Python code to do most of the work for you.

In this code, you specify what parameters you want to run, and the script creates all the appropriate directories and puts the correct inputs file, slurm script, and the executable in each directory. You still have to go in and submit the files yourself, but it’s much easier than writing all these scripts yourself.

The script is called create_input_files.py and here is the part you will need to modify:


#beta = 1/g = 1.6 (float)
beta = 1.6
#number of steps in thermalization (int)
ntherm = 1000
#number of monte carlo steps (int)
nMC = 10000
#number of steps between samples saved (int)
freq = 1
#list of values for lattice length L (list of ints)
L_list = [20]
#list of values for itheta (as fractions of pi) (list of floats)
itheta_list = [0.0,0.25,0.5,0.75, 1.]

#this information goes in the SLURM script -- make sure you use your own email address
script_name = "nonlinearsigma"
job_name = "nlsigma_prelim_tests"
email = "cberger3@bates.edu"
num_cpus = 60
partition = "faculty" 
time_limit = "14-00:00:00" #dd-hh:mm:ss

Beta, ntherm, nMC, and freq all take one number as input, but you can create a list of the number of lattice lengths you want, and the values for itheta (remember these are fractions of pi).

script_name is the executable, so it should be “nonlinearsigma”, but job_name is what will be put in for the job name in the SLURM script. Similarly, you can enter the email address you want included in the SLURM script and choose how many CPUs you want.

Once you’ve modified this script to have the values you want, put it in a directory for this batch. I tend to name those something like run_yyymmdd. Also in that directory should be the executable, so copy that in once you’ve compiled the code.

Then, inside the batch directory, go ahead and run the python script

python create_input_files.py

When it’s done, you should see the subdirectories created – one for each job. You need to go into each subdirectory to submit the jobs using

sbatch submit_sigma.sh

and it will submit the jobs to the scheduler.

You can run this with

sbatch submit_sigma.sh

and you can check the status of your jobs any time with the command

squeue -u your_username

Python Code for Analysis

There are a number of Jupyter Notebooks prepared to explore the data that comes out of the simulation. These are all located in the Analysis folder. There is also a Python class, LatticeData.py, which creates an analyzer object that has internal functions for analysis and visualization.

LatticeData.py - Data Analyzer Class

Importing and initializing

This is a Python class. You can create a lattice data object, which then holds all the functionality you need to analyze and visualize the data. You will need to import this into your notebook like this:

from LatticeData import *

You will want to also import the following: * Numpy * Pandas * Matplotlib * Seaborn

The simplest way to start now is by initializing the object:

analyzer = LatticeData()

You can now do a number of operations with this analyzer. The full list of functions is included later, but if you wanted to get all the data from that folder and put it into a Pandas dataframe, that would look like this:

df = analyzer.get_data()

There are a number of default settings in this class, which you can change when you initialize. The defaults are:

analyzer_default = LatticeData(datadir = "/data/", header = "nonlinearsigma_data",
                 dirheader = "nlsigma_data", Gheader = "Gij_avg_nonlinearsigma_data", 
                 tol = 0.00001, palette = "viridis")

which initializes the following internal variables:


self.path = os.getcwd()+datadir #location of data
self.header = header #set the start of the filename for the data files
self.dirheader = dirheader #set the start of the data directory name from the runs
self.Gheader = Gheader #set the start of the filename for correlation function files
self.tol = tol #set the error range for parameters -- this is for filtering
self.palette = palette #option to change seaborn palette
self.observables = ['Q_L', 'A_L', 'S_L', 'Xi_L'] #observables whose expectation values can be computed
self.parameters = ["itheta", "beta", "length","nMC", "ntherm", "freq"] #parameters read in by the simulation code

Most of these will not need to be changed, but let’s say you want to analyze a special batch of data, which you’ve stored in a directory called data_test, and you want to use a different visualization palette, you could intialize the object like this:

analyzer_special = LatticeData(datadir = "/data_test/", palette = "magma")

(You can choose any seaborn palette for this)

Then when you run the function

df = analyzer_special.get_data()

it will aggregate all the data files in the folder data_test.

Below is a complete list of functions for the LatticeData class, with a brief description. While Python doesn’t have the same public/private distinctions as C++, I’ve organized them into those same groups. Public functions are things that you may want to use. Private functions are functions that you should never need to call yourself, but are called internally.

Lattice Data Class Built-In Functions (“Public” or external)

copy_data_from_directory

copy_data_from_directory(self, src_dir, dst_path = None)

This function loops through a specified directory, finds any simulation directories (directions that begin with “nlsigma_data” or whatever you have specified under dirheader in your intialization), and copies the .csv files within those directories to some destination directory. The default destination directory is whatever you’ve specified for your data directory. The function requires a source directory to be passed as a string – this is the directory where you have all your simulation results that you want copied over – and gives you the option to specify a different destination directory using dst_path)

If some of your simulations are not complete yet (determined by testing whether the .csv has the correct number of lines), this function will not copy those files and will print out the name of the run and how many lines there are in the data output file. For example, the folder run_7_18_23_stats contains some runs that haven’t finished yet. If I try to copy the data from that directory into my data directory:

analyzer.copy_data_from_directory("run_7_18_23_stats")

here’s what appears printed out:

run L_180_beta_1.600000_itheta_0.000000_ntherm_5000_nMC_50000_freq_100 not yet complete: 378 lines
run L_180_beta_1.600000_itheta_2.356194_ntherm_5000_nMC_50000_freq_100 not yet complete: 391 lines
run L_180_beta_1.600000_itheta_3.141593_ntherm_5000_nMC_50000_freq_100 not yet complete: 237 lines
run L_180_beta_1.600000_itheta_1.570796_ntherm_5000_nMC_50000_freq_100 not yet complete: 381 lines
run L_180_beta_1.600000_itheta_0.785398_ntherm_5000_nMC_50000_freq_100 not yet complete: 240 lines

And those runs will not be in the data directory, while completed runs will have been copied in.

all_params

all_params(self)

This function collects every unique set of parameters from your default data directory and returns it as a dataframe. It does not take any inputs – if you want to know what combinations of parameters are in your directory, this function will tell you.

params = analyzer.all_params()
params.head()

>
    freq    nMC ntherm  itheta  beta    length
0   100.0   50000.0 5000.0  0.785398    1.6 20.0
1   100.0   50000.0 5000.0  0.000000    1.6 80.0
2   100.0   50000.0 5000.0  0.000000    1.6 20.0
3   100.0   50000.0 5000.0  1.570796    1.6 40.0
4   100.0   50000.0 5000.0  2.356194    1.6 10.0

It also adds a column that tells you what fraction of $\pi$ $i\theta$ is, to make it easier to determine what new jobs should be run

If you want to collect data from more than one run, you can do this by only specifying which parameters you want in your dataframe, and the function will filter the data accordingly.

param_dict = {"length": 10, "itheta":2.356194}
filtered_data = analyzer.get_data(single_run = False, suppress_output = True, **param_dict)

If you want all the data, just leave out the parameter dictionary entirely and it won’t filter anything.

all_data = analyzer.get_data(single_run = False, suppress_output = True)

NOTE: I strongly recommend suppressing output if you are collecting more than 2 or 3 runs, as it will slow the program down and produce a flood of output.

get_data

get_data(self, single_run = False, corr = False, suppress_output = True, **kwargs)

This function will collect raw data from one or more runs and return it as a Pandas dataframe.

If you want to just get data from one simulation run (e.g. to check thermalization or autocorrelation), you should set single_run to True. If you want the correlation function data from that run, you should set corr to True – otherwise it will return the observable data. If you want it to print out the parameter sets it’s putting into the dataset, set suppress_output to False

You then need to specify what the parameters are for the run you want to see. You do this by creating a dictionary. When selecting a single one, you must ensure your dictionary has all the parameter values specified. The keys for these values are: * “length”: length of the lattice in each direction * “itheta”: value of the imaginary value used for theta – actual number here, not an integer multiple of pi, but you can always use np.pi to specify it * “beta”: for our purposes this will always be 1.6, but you need to specify it anyway * “nMC”: number of steps in the Monte Carlo loop * “ntherm”: number of steps in the thermalization loop * “freq”: frequency with which the configurations were saved

If you forget one of these, the function will remind you:

param_dict = {"length": 10, "freq": 100, "itheta":2.356194, "beta":1.6, "nMC":50000}
one_run = analyzer.get_data(single_run = True, suppress_output = False, **param_dict)

Missing parameters in input: 
['ntherm']

So you know now to add in the “ntherm” you’re looking for. Now it should work:

param_dict = {"length": 10, "freq": 100, "itheta":2.356194, "beta":1.6, "nMC":50000, "ntherm": 5000}
one_run = analyzer.get_data(single_run = True, suppress_output = False, **param_dict)

freq 100
nMC 50000
ntherm 5000
itheta 2.356194
beta 1.6
length 10

one_run.head(3)

    step    |phi|   Q_L A_L S_L Xi_L    F_LRe   F_LIm   acc dt  ... Q_L_ta  A_L_ta  S_L_ta  Xi_L_ta corr_length_Re  corr_length_Im  F_Re_py F_Im_py mass_gap_Re mass_gap_Im
0   0   100.0   0.479179    -183.166830 -184.295868 45.724903   3.675149    8.145972    0.240000    0.0 ... 1   3   3   4   8.966957    -3.640801   0.916337    0.890996    0.095738    0.038872
1   100 100.0   0.159302    -195.846927 -196.222275 62.369185   3.675149    8.145972    0.203762    0.0 ... 1   3   3   4   10.472580   -4.252121   0.916337    0.890996    0.081974    0.033283
2   200 100.0   0.000000    -199.012224 -199.012224 48.119359   3.675149    8.145972    0.207164    0.0 ... 1   3   3   4   9.198745    -3.734913   0.916337    0.890996    0.093325    0.037892
3   300 100.0   0.312655    -179.261053 -179.997729 46.124447   3.675149    8.145972    0.205681    0.0 ... 1   3   3   4   9.006048    -3.656673   0.916337    0.890996    0.095322    0.038703

do_stats

do_stats(self, therm = 0., stack = False, **kwargs)

This function collects data from the directory you specified when you initialized the object. If you want to filter the data by parameters, you just need to enter a parameter dictionary just like when using get_data() above, and it will filter to include only data that matches the parameter values you’ve entered.

Once all the raw data is collected, it will perform some basic statistical analysis. All observables will have a mean and standard error calculated from data after thermalization. If you need to change the thermalization point, you can increase it by changing the therm argument in the function – give it the fractional value of the data you want it to drop from the beginning of the dataset. So for example, if you have a run with nMC = 1000 and freq = 10, you will have a total of 100 steps in your dataset. If you set therm = 0.2, it will drop the first 20 steps before computing the mean and standard error.

After calculating means and standard errors, it also computes the autocorrelation time (the step at which the observable’s autocorrelation value drops below 0.3), and then determines how long that run took to complete, saving that information in seconds, minutes, and hours.

It saves all of this along with the parameters used in that run. The result is a very comprehensive dataframe. Here’s an example of it run on a data directory with 25 runs in it:

df_stats = analyzer.do_stats()
df_stats.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 25 entries, 0 to 24
Data columns (total 41 columns):
 #   Column               Non-Null Count  Dtype  
---  ------               --------------  -----  
 0   length               25 non-null     float64
 1   itheta               25 non-null     float64
 2   beta                 25 non-null     float64
 3   nMC                  25 non-null     float64
 4   ntherm               25 non-null     float64
 5   freq                 25 non-null     float64
 6   |phi|_mean           25 non-null     float64
 7   Q_L_mean             25 non-null     float64
 8   A_L_mean             25 non-null     float64
 9   S_L_mean             25 non-null     float64
 10  Xi_L_mean            25 non-null     float64
 11  F_LRe_mean           25 non-null     float64
 12  F_LIm_mean           25 non-null     float64
 13  acc_mean             25 non-null     float64
 14  Q_L_ta               25 non-null     float64
 15  A_L_ta               25 non-null     float64
 16  S_L_ta               25 non-null     float64
 17  Xi_L_ta              25 non-null     float64
 18  corr_length_Re_mean  25 non-null     float64
 19  corr_length_Im_mean  25 non-null     float64
 20  F_Re_py_mean         25 non-null     float64
 21  F_Im_py_mean         25 non-null     float64
 22  mass_gap_Re_mean     25 non-null     float64
 23  mass_gap_Im_mean     25 non-null     float64
 24  |phi|_std            25 non-null     float64
 25  Q_L_std              25 non-null     float64
 26  A_L_std              25 non-null     float64
 27  S_L_std              25 non-null     float64
 28  Xi_L_std             25 non-null     float64
 29  F_LRe_std            25 non-null     float64
 30  F_LIm_std            25 non-null     float64
 31  acc_std              25 non-null     float64
 32  corr_length_Re_std   25 non-null     float64
 33  corr_length_Im_std   25 non-null     float64
 34  F_Re_py_std          25 non-null     float64
 35  F_Im_py_std          25 non-null     float64
 36  mass_gap_Re_std      25 non-null     float64
 37  mass_gap_Im_std      25 non-null     float64
 38  time (sec)           25 non-null     float64
 39  time (min)           25 non-null     float64
 40  time (hr)            25 non-null     float64
dtypes: float64(41)
memory usage: 8.1 KB

If I only want the analyzed data for runs where the lattice had a length of 20, I could modify this as follows:

df_stats = analyzer.do_stats(**{"length":20})
df_stats.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5 entries, 0 to 4
Data columns (total 41 columns):
 #   Column               Non-Null Count  Dtype  
---  ------               --------------  -----  
 0   length               5 non-null      float64
 1   itheta               5 non-null      float64
 2   beta                 5 non-null      float64
 3   nMC                  5 non-null      float64
 4   ntherm               5 non-null      float64
 5   freq                 5 non-null      float64
 6   |phi|_mean           5 non-null      float64
 7   Q_L_mean             5 non-null      float64
 8   A_L_mean             5 non-null      float64
 9   S_L_mean             5 non-null      float64
 10  Xi_L_mean            5 non-null      float64
 11  F_LRe_mean           5 non-null      float64
 12  F_LIm_mean           5 non-null      float64
 13  acc_mean             5 non-null      float64
 14  Q_L_ta               5 non-null      float64
 15  A_L_ta               5 non-null      float64
 16  S_L_ta               5 non-null      float64
 17  Xi_L_ta              5 non-null      float64
 18  corr_length_Re_mean  5 non-null      float64
 19  corr_length_Im_mean  5 non-null      float64
 20  F_Re_py_mean         5 non-null      float64
 21  F_Im_py_mean         5 non-null      float64
 22  mass_gap_Re_mean     5 non-null      float64
 23  mass_gap_Im_mean     5 non-null      float64
 24  |phi|_std            5 non-null      float64
 25  Q_L_std              5 non-null      float64
 26  A_L_std              5 non-null      float64
 27  S_L_std              5 non-null      float64
 28  Xi_L_std             5 non-null      float64
 29  F_LRe_std            5 non-null      float64
 30  F_LIm_std            5 non-null      float64
 31  acc_std              5 non-null      float64
 32  corr_length_Re_std   5 non-null      float64
 33  corr_length_Im_std   5 non-null      float64
 34  F_Re_py_std          5 non-null      float64
 35  F_Im_py_std          5 non-null      float64
 36  mass_gap_Re_std      5 non-null      float64
 37  mass_gap_Im_std      5 non-null      float64
 38  time (sec)           5 non-null      float64
 39  time (min)           5 non-null      float64
 40  time (hr)            5 non-null      float64
dtypes: float64(41)
memory usage: 1.7 KB

Notice we now only have 5 entries in our dataframe, not 25. We can check that this worked:

df_stats["length"].unique()

array([20.])

If you want the data returned using Pandas MultiIndex, set stack to True, but MultiIndex doesn’t always play well with seaborn and other plotting tools, so the default is False.

get_plot_data

get_plot_data(self, obs = "Q_L", L = 10, beta = 1.6, nMC = 10000, ntherm = 1000, freq = 100, stack = False)

This function allows you to get the analyzed data for one observable (you may specify which one, but the default is Q_L) as a function of itheta in order to plot it. You must specify a single value for all other parameters (L, beta, nMC, ntherm, freq) or leave them blank to use the defaults.

This pulls its data from the internally stored dataframe self.df_stats, which is created when you run the do_stats function. If you have not run that function yet, it will do it for you, with the default settings of therm = 0.0 and stack = False.

This doesn’t play nice with MultiIndex right now, so I recommend making sure to do this with the default of

stack = False

This function could be modified in order to choose your independent variable, but right now all the plots we are interested in are functions of itheta, so it’s unneccesary to plot the observable as a function of any other parameter. Seaborn can be used with the raw data to study systematic effects or other things that may be functions of nMC, ntherm or L.

This returns three items: x, y, and y_err, which can then be plotted immediately with the matplotlib errorbar function.

For example, if you wanted to plot Q_L as a function of itheta for all the different lengths you have in your data, you would do the following:

params = analyzer.all_params()
lengths = params["length"].unique()
colors = sns.color_palette("Blues", len(lengths))
observable = "Q_L"

for n,length in enumerate(lengths):
    x,y,err = analyzer.get_plot_data(obs = observable, L = length, beta = 1.6, nMC = 50000, 
                                     ntherm = 5000, freq = 100)
    plt.errorbar(x, y , yerr = err, marker = ".", ls = "none", color = colors[n], label ="L="+str(length))
plt.legend()
plt.title(observable)
plt.show()

And you would get the following output:

{width = 300}

And if you wanted to plot the magnetic susceptibility Xi_L as a function of itheta for each length, you would do:

params = analyzer.all_params()
lengths = params["length"].unique()
colors = sns.color_palette("Reds", len(lengths))
observable = "Xi_L"

for n,length in enumerate(lengths):
    x,y,err = analyzer.get_plot_data(obs = observable, L = length, beta = 1.6, nMC = 50000, 
                                     ntherm = 5000, freq = 100)
    plt.errorbar(x, y , yerr = err, marker = ".", ls = "none", color = colors[n], label ="L="+str(length))
plt.legend()
plt.title(observable)
plt.show()

And you would get the following output:

{width = 300}

get_corr_func

get_corr_func(self,suppress_output = False,**kwargs)

This function returns the average correlation function for the set of parameters specified. This requires returning one single run, so it plays by the same rules as get_data with single_run = True and it will tell you if you left out a parameter:

itheta = np.pi
beta = 1.6
length = 20
nMC = 50000
corr_params = {"itheta": itheta, "beta": beta,"length": length,"nMC": nMC}
G_ij = analyzer.get_corr_func(suppress_output = False, **corr_params)

Missing parameters in input: 
['ntherm', 'freq']

When you specify the complete set of parameters, it returns a 2D numpy array that represents the average correlation function on each lattice site. You can then plot this with imshow:

plt.imshow(G_ij, cmap = "viridis", aspect='equal')
plt.colorbar()
plt.title("Correlation function for L = "+str(corr_params["length"])+", itheta = "+str(corr_params["itheta"]))
plt.show()

and you will get something that looks like this:

{width = 300}

Currently, the simulation returns the average value of the correlation function computed after thermalization. There is no way to change this after the simulation is run, and we don’t have the code set up yet to return any error on the correlation function. That may come in future versions of the code.

Lattice Data Class Built-In Functions (“Private” or internal)

get_data_files

get_data_files(self, corr = False)

This function is referenced by all_params and get_data (which are both in turn called by other functions).

It loops through every file in the data directory (self.path) set when you initialize the object. If you have set corr to True (e.g. as done in the function get_corr_func), it uses self.Gheader to pick out all the simulation files that are correlation function data and put them in a list. If corr is False, it uses self.header to pick out all the regular data files from the simulation and create a list of those.

It returns just a list of all the relevant data filenames, which can then be used by other functions to extract the parameters (as in all_params) or make a dataframe (as in get_data).

get_file_params

get_file_params(self, file)

This function is referenced by all_params and get_data (which are both in turn called by other functions).

It takes a filename (e.g. one from the list generated by get_data_files) and extracts from that filename the parameter values used in the simulation. It returns a dictionary with the following keys:

“freq” – the number of steps between saved configurations
“nMC” – the number of steps in the Monte Carlo loop
“ntherm” – the number of thermalization steps before beginning the Monte Carlo loop
“itheta” – the value of the parameter $i \theta$ (actual value, not integer multiples of pi)
“beta” – inverse lattice coupling (currently set to 1.6 for all simulations, but could be changed)
“length” – number of lattice sites in both x and y

in_list

in_list(self,pdict,**kwargs)

This function is only called by get_data but it serves the important function of managing all the filtering we want to it. It checks each parameter fed into get_data in the keyword arguments and it checks it against the parameters in the filename for a simulation run. If all the parameter values specified in the keyword arguments match parameters in the file, it returns True, otherwise it returns False.

calc_F

calc_F(self, **kwargs)

This function calculates the correlation function at the smallest nonzero lattice momentum ( $p_{0} = 2 \pi/L$ ) . It requires the correlation function for that set of parameters ( $G$ , which it obtains by calling get_corr_func) and then performs the following calculation, summing over all lattice sites $(x,y)$ :

$\mathcal{F} = \frac{1}{2}\sum_{x,y}(e^{2 \pi i x/L} + e^{2 \pi i y/L})G(x,y)$

calc_corr_length

calc_corr_length(self,Xi,L,F_py)

Pipeline for data generation and analysis

The basic workflow

Compiling the MC code

Running the code

Other important things to note before running the code

Running the code – single job in interactive node

Running the code – single SLURM submission

Running the code – batch SLURM submissions

Python Code for Analysis

LatticeData.py - Data Analyzer Class

Importing and initializing

Lattice Data Class Built-In Functions (“Public” or external)

copy_data_from_directory

all_params

get_data

do_stats

get_plot_data

get_corr_func

Lattice Data Class Built-In Functions (“Private” or internal)

get_data_files

get_file_params

in_list

calc_F

calc_corr_length

Analysis Notebooks

AnalysisTesting.ipynb

PhiDist.ipynb

DataComparison.ipynb

SystematicsAndTiming.ipynb

CorrelationFunction.ipynb

Observables.ipynb