Numerical project I - MD simulation (NVE)
Numerical project I - MD simulation (NVE)¶
Congratulations! You have mastered important concepts in classical and statistical mechanics, and learned about the nuts and bolts of MD simulations. You are now ready for your first project - writing your very own MD code. How exciting! The project will be written in Python, a highly popular coding language, extremely relevant for both scientific research and the current job market.
You can find the link for the first project in the course’s moodle.
If you followed correctly the Prerequisites instructions, your first assignment repository should be cloned ub your personal GitHub.
You do not need to edit the Ar_init_super.xyz, LICENSE, or README.md files. You will only need to edit the run.py and sim.py files. They both have very detailed instructions inside. For example, sim.py reads:
#sim.py
WELCOME TO YOU FIRST PROJECT! THIS BIT OF TEXT IS CALLED A DOCSTRING.
BELOW, I HAVE CREATED A CLASS CALLED "SIMULATION" FOR YOUR CONVENIENCE.
I HAVE ALSO IMPLEMENTED A CONSTRUCTOR, WHICH IS A METHOD THAT IS CALLED
EVERY TIME YOU CREATE AN OBJECT OF THE CLASS USING, FOR EXAMPLE,
>>> mysim = Simulation( dt=0.1E-15, L=11.3E-10, ftype="LJ" )
I HAVE ALSO IMPLEMENTED SEVERAL USEFUL METHODS THAT YOU CAN CALL AND USE,
BUT DO NOT EDIT THEM. THEY ARE: evalForce, dumpXYZ, dumpThermo and readXYZ.
YOU DO NOT NEED TO EDIT THE CLASS ITSELF.
YOUR JOB IS TO IMPLEMENT THE LIST OF CLASS METHODS DEFINED BELOW.
ENTER YOUR COE WHERE YOU WILL SEE THE FOLLOWING TEXT:
################################################################
####################### YOUR CODE GOES HERE ####################
################################################################
YOU ARE, HOWEVER, EXPECTED TO UNDERSTAND WHAT ARE THE MEMBERS OF THE CLASS
AND USE THEM IN YOUR IMPLEMENTATION. THEY ARE ALL EXPLAINED IN THE
DOCSTRING OF THE CONSTRUCTOR BELOW. FOR EXAMPLE, WHENEVER YOU WISH
TO USE/UPDATE THE MOMENTA OF THE PARTICLES IN ONE OF YOUR METHODS, YOU CAN
ACCESS IT BY USING self.p.
>>> self.p = np.zeros( (self.Natoms,3) )
FINALLY, YOU WILL NEED TO EDIT THE run.py FILE WHICH RUNS THE SIMULATION.
SEE MORE INSTRUCTIONS THERE.
And the instructions in run.py are:
HERE, TO RUN THE SIMULATION, YOU WILL NEED TO DO THE FOLLOWING THINGS:
1. CREATE AN OBJECT OF THE SIMULATION CLASS. A MINIMAL EXAMPLE IS
>>> mysim = Simulation( dt=0.1E-15, L=11.3E-10, ftype="LJ" )
2. DEFINE THE PARAMETERS FOR THE POTENTIAL.
USE A DICTIONARY, FOR EXAMPLE FOR LJ MODEL OF ARGON, IN SI UNITS
>>> params = { "eps": 1.656778224E-21, "sig": 3.4E-10 }
THEN, CALLING THE METHODS YOU IMPLEMENTED IN sim.py, YOU NEED TO
3. READ THE INITIAL XYZ FILE PROVIDED OR SAMPLE INITIAL COORDINATES.
4. SAMPLE INITIAL MOMENTA FROM MB DISTRIBUTION.
5. REMOVE COM MOTION.
6. RUN THE SIMULATION, INCLUDING PRINTING XYZ AND ENERGIES TO FILES.
NOTE THAT TO CALL A METHOD OF A CLASS FOR A SPECIFIC OBJECT, THE SYNTAX IS
>>> mysim.funcName( args )
FINALLY, YOU SHOULD
7. ANALYZE YOUR RESULTS. PLOT THE KINETIC AND POTENTIAL ENERGIES,
AND SEPARATE FIGURE WITH THE CHANGE IN % OF THE TOTAL ENERGY
WITH RESPECT TO t=0.
NOTE: THE INPUT FILE GIVEN HAS THE COORDINATES IN ANGSTROM.
*YOUR OUTPUT XYZ FILE SHOULD BE PRINTED IN ANGSTROM TOO*, BUT YOU CAN USE
ANY UNITS YOU WANT IN BETWEEN, I SUGGEST USING SI UNITS.
As you can see, most of your coding will be done in the sim.py file.
I have created a class called Simulation which defines the methods and data members that you can use in your implementation. Your job is twofold:
Implement the methods for class Simulation, based on the material we have covered.
Implement a several short scripts
run.pyto call the methods you have implemented and run a simulation.In methods containing the keyword
Example, write doctest examples proving your methods are correct and well implemented.
The methods you will have to implement include, for example, sampleMB() which samples initial positions from the Maxwell-Boltzmann distribution. How do I know? Because in sim.py you will see its definition, with a description of what needs to be implemented and the place to write your code:
def sampleMB( self, temp, removeCM=True ):
"""
THIS FUNCTIONS SAMPLES INITIAL MOMENTA FROM THE MB DISTRIBUTION.
IT ALSO REMOVES THE COM MOMENTA, IF REQUESTED.
Parameters
----------
temp : float
The temperature to sample from.
removeCM : bool, optional
Remove COM velocity or not. The default is True.
Returns
-------
None. Sets the value of self.p.
Tests
-----
IMPLEMENT THE FOLLOWING TESTS BELOW THE WORD Example:
1. THE MEAN MOMENTUM OF EACH COORDINATE SHOULD BE VERY CLOSE TO ZERO
2. THE STD SHOULD BE CLOSE TO SQRT(m*Kb*T)
NOTE: IN THIS CASE, WE SHOW A TEST FOR THE MEAN, WRITE THE OTHERS YOURSELF.
Example:
>>> temp = 300 # temperature in K, using default values for mass and Natoms.
>>> mysim = Simulation( dt=0.1E-15, L=11.3E-10, ftype="LJ" )
>>> mysim.sampleMB(temp)
>>> mean_momentum = np.mean(mysim.p, axis=0)
>>> np.linalg.norm(mean_momentum[0]) < 1e-20
True
>>> np.linalg.norm(mean_momentum[1]) < 1e-20
True
>>> np.linalg.norm(mean_momentum[2]) < 1e-20
True
"""
################################################################
####################### YOUR CODE GOES HERE ####################
################################################################
I highly recommend implementing one method at a time and testing that it works by running run.py that calls it, before moving on to the next method. Do not forget to commit and push along the way to create checkpoints you can go back to if something goes wrong.
The repo also has a file called ReportProjI.ipynb. This is a Jupyter Notebook which you can open by running JupyterLab or Jupyter Notebook from Anaconda. In it, you will find the definitions and details of the specific simulations you have to run. It also contains questions that you need to answer and figures you need to plot. Create a separate run script for each one of the sections in the notebook (run_harmonic.py and run_Ar.py), run the simulations, and answer the questions. Submit your report by commiting and pushing the file to your repo.