Numerical Project II - MC simulation (NVT)

The project

Now that you have mastered how to sample arbitrary distribution functions using the Metropolis algorithm, you are ready to tackle the second exercise.

You will use the same repository as for the first exercise. Create a new branch named project2 from your project1 branch, since Project II reuses your Project I code: git checkout -b project2 project1 (the first time you push it: git push -u origin project2), or with the Codespaces interface as described in the README. Note that you will need your code to evaluate the LJ potential and a harmonic trap working also for the second exercise. If it did not work in the previous exercise, you have another shot at correcting it. But if you are truly stuck and it does not work, contact Barak for help as soon as possible.

Your goal is to implement a Metropolis MC simulation in the canonical ensemble. To do that, you will need to add the following two methods to your simulations class. The first, MCstep, performs a single translation move and decides whether to accept or reject it. The second, runMC, contains the instructions for performing an MC simulation, similarly to how run holds the instructions for an MD simulation. As usual, the docstring for each method contains all of the relevant information.

    def MCstep( self, **kwargs ):
        """
        THIS FUNCTIONS PERFORMS ONE METROPOLIS MC STEP IN THE NVT ENSEMBLE.
        YOU WILL NEED TO PROPOSE TRANSLATION MOVES, APPLY  
        PBC, CALCULATE THE CHANGE IN POTENTIAL ENERGY, ACCEPT OR REJECT, 
        AND CALCULATE THE ACCEPTANCE PROBABILITY.

        Returns
        -------
        None. Sets self.R, self.U and self.F.
        """
        
        ################################################################
        ####################### YOUR CODE GOES HERE ####################
        ################################################################

    def runMC( self, **kwargs ):
        """ 
        THIS FUNCTION DEFINES WHAT AN MC SIMULATION DOES, GIVEN AN INSTANCE OF 
        THE SIMULATION CLASS. YOU WILL NEED TO LOOP OVER MC STEPS, 
        PRINT THE COORDINATES AND ENERGIES EVERY PRINTFREQ TIME STEPS 
        TO THEIR RESPECTIVE FILES, SIMILARLY TO YOUR MD CODE.

        Returns
        -------
        None.

        """   
        
        ################################################################
        ####################### YOUR CODE GOES HERE ####################
        ################################################################

You will perform simulations for the two systems we have investigated before, a single Ar atom in a 1D harmonic trap and a system of 256 Ar atoms at three thermodynamic conditions. You will also need to implement a new method, evalAnharm, and perform simulations for an Ar atom in an anharmonic trap, and compare it with the harmonic example. Use \( \lambda = 7.757 \cdot 10^{23} J m^{-4} \).

    def evalAnharm(self, Lambda):
        """
        THIS FUNCTION EVALUATES THE POTENTIAL AND FORCE FOR AN ANHARMONIC TRAP.

        Parameters
        ----------
        Lambda : float
            The parameter of the trap U = 0.25 * Lambda * x**4

        Returns
        -------
        None. Sets the value of self.F and self.U.

        Tests
        -----
        CREATE AN OBJECT OF THE SIMULATION CLASS WITH FOLLOWING PARAMS:
        TIMESTEP 0.1, BOX SIDE LENGTH 10, ANHARMONIC POTENTIAL, LAMBDA 4,
        ONE ATOM WITH MASS 2.

        IMPLEMENT THE FOLLOWING TESTS BELOW THE WORD Example:
        1. CHECK THAT AT THE ORIGIN THE POTENTIAL AND THE FORCE ARE ZERO.
        2. CHECK THAT AT (0.5,0,0) THE POTENTIAL IS 0.0625 AND THE FORCE IS (-0.5,0,0).

        Example: 
        """

        ################################################################
        ####################### YOUR CODE GOES HERE ####################
        ################################################################

As before, the simulations to run and the questions to answer are in ReportProjII.ipynb. Push your work to the project2 branch by the deadline set in GitHub Classroom. Please note that while you will only need these three functions to have a working MC simulation, I assume here that you will re-evaluate the full potential after every MC step. This is OK for the single particle simulations. But, as we discussed in class, this is very wasteful for the case of 256 Ar atoms and I encourage you to come up with a solution to calculate efficiently the change in energy after randomly moving one atom at a time. This is, however, not mandatory.