<jupyter_start><jupyter_text>Using Monte Carlo on molecular systemsWe are going to use Monte Carlo (MC) simulation to solve a chemical problem. MC simulation is used in many different fields, and one of those is molecular simulation. MC methods for molecules can be quite advanced these days, but for this bootcamp, we will be using MC to simulate noble gases.In statistical mechanics, we are interested in computing averages of thermodynamic properties as a function of atom positions an momenta. A thermodynamic average depending only on configurational properties can be computed by evaluating a multivariate integral.$$ \left = \int_V Q\left(\textbf{r}^N\right)	\rho\left(\textbf{r}^N\right) d\textbf{r}^N	\label{eq.statMechAverage} $$ where $Q\left(\textbf{r}^N\right)$ is the thermodynamic quantity of interest that depends on only on the configuration, $\rho\left(\textbf{r}^N\right)$ is the probability density, and $V$ defines the volume of **configuration space** over which $\rho$ has support.This integral is very hard to compute, even for small atomic systems. For instance, a monoatomic system of 10 atoms leads to a 30 dimensional integral. We learned in the previous lesson about using Monte Carlo to evaluate integrals. We can apply this same approach to solve this integral. There are a few extra things we have to consider, however, because of the high dimensionality of the problem. We'll worry more about exactly how to implement this MC integration tomorrow. Today's lesson will focus on the model we are going to use for our thermodynamic quantity (in our case the system potential energy). We will use this in our Monte Carlo simulation in the next lesson. Modeling the systemWhen modeling atomic and molecular systems, the energy of nonbonded interactions are often approximated using the **Lennard Jones Potential**. The [Lennard Jones Potential](https://en.wikipedia.org/wiki/Lennard-Jones_potential) was proposed by <NAME> in 1924 and is commonly used to model nonbonded interactions of atoms.$$ U(r) = 4 \epsilon \left[\left(\frac{\sigma}{r}\right)^{12} -\left(\frac{\sigma}{r}\right)^{6} \right] $$This is a pairwise potential which is a function of distance between two atoms. It has two parameters, $\sigma$ and $\epsilon$, which represent the size of the particle, and the depth of the potential well (you can think of this as how strongly the particles are attracted to one another), respectively.This potential is commonly used in molecular simulation to model nonbonded interactions, not just for those of noble gases. For more complicated systems, you will have additional energy terms in addition to the LJ potential, but for our system of a Lennard Jones fluid, the is the only contribution to the potential energy. The first thing we will do is write a function which evaluates the LJ energy based on an input of distance.Commonly used parameters for simulating Argon are $\sigma = 3.4 \unicode{xC5} $ and $\epsilon/k_B = 120~K$ Dealing with unitsThroughout this project, we will be working with **reduced units**. As stated above, common values are $\sigma = 3.4 \unicode{xC5}$ and $\epsilon/k_B = 120~K$. When converted to SI units, these quickly become inconvenient to work with. Often, in simulation, we will use something called *reduced units* in order to make calculations more convenient. This approach essentially scales quantities by characteristic values to get them closer to unity.For example, when working with Argon, the distances we compute will be in units of $\sigma$ instead of angstrom. Quantity    | Expression------------|------------Length      | $$L^*=L / \sigma$$Density     | $$\rho^* = N \sigma^3 / V$$Energy      | $$U^* = U / \epsilon$$Pressure    | $$P^* = P \sigma^3 / \epsilon$$Volume      | $$V^* = V / \sigma^3 $$Temperature | $$T^* = k_{B} T / \epsilon $$Time        | $$t^* = t \sqrt{\frac{\epsilon}{ m \sigma^2}}$$Conveniently for us, Lennard Jones fluids have the surprisingly pleasant behavior of possessing **universal behavior** when expressed in terms of reduced units as (you can verify this for yourself by substituting $\sigma$ and $\epsilon$ for their reduced unit expressions):$$ 	U^*\left(r_{ij} \right) = 4 \left[\left(\frac{1}{r^*_{ij}}\right)^{12} -\left(\frac{1}{r^*_{ij}}\right)^{6} \right] $$Since we know that we will need this as a function for our MC calculation, we will write a function for it to use later.Recall that in Python, we define a function using the syntax:```pythondef function_name(function_arguments):    ** FUNCTION WHICH USES function_arguments HERE**    return return_value```<jupyter_code># write imports here
import math
# write calculate_LJ function here
def calculate_LJ(r_ij):
    r6_term = math.pow(1/r_ij, 6)
    r12_term = math.pow(1/r_ij, 12)
    
    pairwise_energy = 4 * (r12_term - r6_term)
    
    return pairwise_energy

calculate_LJ(8)<jupyter_output><empty_output><jupyter_text>Sanity checksIt is always a good idea to test the expected behavior of your function, or do some "sanity checks". We'll think about a few qualities of the LJ potential to check if our function is reasonable.By rearranging the equation for the LJ potential, we can see a few qualities that should be true if our function is implemented correctly.**Full Expression**$$ U(r) = 4 \epsilon \left[\left(\frac{\sigma}{r}\right)^{12} -\left(\frac{\sigma}{r}\right)^{6} \right] $$**Reduced Unit Expression**$$ 	U^*\left(r_{ij} \right) = 4 \left[\left(\frac{1}{r^*_{ij}}\right)^{12} -\left(\frac{1}{r^*_{ij}}\right)^{6} \right] $$1. **FOR WHAT VALUE OF r SHOULD $U(r)$ be ZERO?**.1. A minimum occurs at $r = 2^{1/6}\sigma.$ In our case, we are using reduced units and $\sigma$ is equal to one. This minimum will have a value equal to $-\epsilon$, in our case -1. For figuring out these values, you could write out the solution on paper. For 1, you would set $U(r)=0$ and solve for $r$. For 2,you would take the derivative, set it equal to zero and solve for $r$. You could also use [Wolfram Alpha](https://www.wolframalpha.com/input/?i=graph+4*%28%281%2Fr%29%5E12+-+%281%2Fr%29%5E6%29++from+r+%3D+0.9+to+r%3D3). Python Concept: Assert StatementsWe'll check these two test statements using something called an `assert` statement in Python. In an `assert` statement, you assert that something is `True`. If the statement following the word `assert` is `True`, nothing happens. However, if the statement is `False`, you will get an assertion error.<jupyter_code>calculate_LJ(1)
calculate_LJ(math.pow(2, (1/6)))
assert 1 == 1
assert calculate_LJ(1) == 0
assert calculate_LJ(math.pow(2, (1/6))) == -1<jupyter_output><empty_output><jupyter_text>Python Style - DocstringsPart of this class is cultivating good coding habits. A "best practice" when writing your code is to make sure that your code is adequately documented. This could be through code comments, but a very good way to document your code is to use something called `docstrings` in Python. In this course, we may not always write `docstrings` right away, but you should always write them.A `docstring`, or `documentation string` is a string which follows the function definition and describes a function's inputs and behavior. If you write `docstrings`, you will have a much easier time returning to your code, sharing it with others, and creating online documentation. For Python, there is a tool called `Sphinx` which can parse your code for `docstrings` and build web documentation for you (we will not discuss this in this course).We will write `Numpy` style `docstrings`. Consider our `calculate_LJ` function with a docstring. You should rewrite your `cacluate_LJ` function in the next cell or scroll up to add a docstring to your original function.<jupyter_code># write calculate_LJ function here
def calculate_LJ(r_ij):
    """
    The Lennard Jones interaction energy between two particles.
    
    Computes the pairwise Lennard Jones interaction energy based on the separation distance in reduced units.
    
    Parameters
    ----------
    r_ij : float
        The distance between the particles in reduced units.
    """
    
    r6_term = math.pow(1/r_ij, 6)
    r12_term = math.pow(1/r_ij, 12)
    
    pairwise_energy = 4 * (r12_term - r6_term)
    
    return pairwise_energy

calculate_LJ(8)
help(calculate_LJ)<jupyter_output>Help on function calculate_LJ in module __main__:

calculate_LJ(r_ij)
    The Lennard Jones interaction energy between two particles.
    
    Computes the pairwise Lennard Jones interaction energy based on the separation distance in reduced units.
    
    Parameters
    ----------
    r_ij : float
        The distance between the particles in reduced units.<jupyter_text>Lennard Jones Energy of Atomic SystemWhen applying the LJ potential to a set of atomic coordinates, the total potential energy of the system is equal to the sum of the LJ energy over all pairwise interactions:$$ U \left( \textbf{r}^N \right) = \sum_{i < j} U \left( r_{ij} \right) $$We can write some steps out for calculating the total Lennard Jones potential energy for a system.1. Calculate the distance between each particle and every other particle in the system.1. Evaluate the Lennard Jones potential for each distance1. Sum the energies to get the total potential energy.From this procedure, we can see that we will need a few more functions. 1. We will need a function that can calculate a distance between two particles based on their coordinates. 2. We will need to loop over particle pairs.So that we have a system where we know the answer, we will use some sample systems and benchmarks from the National Institutes of Standards and Technology (NIST). Some sample configurations have been included with your lesson materials.>  Exercise> Write a function called `calculate_distance` which takes in two coordinates as lists ([x, y, z]), and returns the distance between the two points.>> Use the provided docstring to write your function:<jupyter_code>def calculate_distance(coord1, coord2):
    """
    Calculate the distance between two 3D coordinates.
   
    Parameters
    ----------
    coord1, coord2: list
        The atomic coordinates
    
    Returns
    -------
    distance: float
        The distance between the two points.
    """
    distance = 0
        
    
    for i in range(3):
        distance += math.pow(coord2[i] - coord1[i], 2)
    
    
    return math.sqrt(distance)

## Check your function with assert statements
point_1 = [0,0,0]
point_2 = [1,0,0]

dist_1 = calculate_distance(point_1, point_2)
assert dist_1 == 1
## Calculate total energy function here
def calculate_total_energy(coordinates):
    
    total_energy = 0
    num_atoms = len(coordinates)
    
    for i in range(num_atoms):
        for j in range(i+1, num_atoms):
            dist_ij = calculate_distance(coordinates[i], coordinates[j])
            pairwise_energy = calculate_LJ(dist_ij)
            
            total_energy += pairwise_energy
    
    return total_energy
# Provided function

def read_xyz(filepath):
    """
    Reads coordinates from an xyz file.
    
    Parameters
    ----------
    filepath : str
       The path to the xyz file to be processed.
       
    Returns
    -------
    atomic_coordinates : list
        A two dimensional list containing atomic coordinates
    """
    
    with open(filepath) as f:
        box_length = float(f.readline().split()[0])
        num_atoms = float(f.readline())
        coordinates = f.readlines()
    
    atomic_coordinates = []
    
    for atom in coordinates:
        split_atoms = atom.split()
        
        float_coords = []
        
        # We split this way to get rid of the atom label.
        for coord in split_atoms[1:]:
            float_coords.append(float(coord))
            
        atomic_coordinates.append(float_coords)
        
    
    return atomic_coordinates, box_length

## Read in first sample configuration file here
import os
config1_file = os.path.join("lj_sample_configurations", "lj_sample_config_periodic1.txt")

sample_coords, box_length = read_xyz(config1_file)

print(box_length)
print(len (sample_coords))
## Assertions here
total_energy = calculate_total_energy(sample_coords)
assert total_energy == 4351.5
## Check total energy here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08912990296012395, 'fraction_numerical': 0.012884286063769061, 'mean_word_length': 3.8340559716200238, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 21, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '266049', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3423, 'n_tokens_neox': 3255, 'n_words': 1526, 'path': 'lj-part1.ipynb', 'repo_name': 'msse-2021-bootcamp/team2-project', 'size': '12275'}