Numerical Project II - MC simulation (NVT)
Contents
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. 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.
"""
################################################################
####################### YOUR CODE GOES HERE ####################
################################################################
def runMC( self, **kwargs ):
"""
THIS FUNCTION DEFINES 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.
"""
if( self.debug ):
print( "Called evalAnharm with Lambda = " + str(Lambda) )
################################################################
####################### YOUR CODE GOES HERE ####################
################################################################
As before, report your results in a file called ReportProjII.ipynb. 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.