Rayn Lakha
Abstract
Characteristics of molecular dynamics drive many notable phenomena, so I was immensely interested in gaining an understanding of their nature. Over the course of a 6 week internship, I learned the fundamental principles of molecular dynamics and endeavoured to code a simple 2D simulation of Argon atoms in Python under the guidance and teaching of Professor Sarkisov of the University of Manchester and his research group. This report chronicles the development of this simulation on a weekly basis, and explains the key concepts behind its performance, including the utilisation of Newton’s laws of motion, Lennard-Jones potentials, Maxwell-Boltzmann distributions, and periodic boundary conditions.
Introduction
The development process of the simulation is discussed, and the state of the code at each week is evaluated, with the problem-solving and iterative design process emphasised. 6 key steps in the simulation process are identified, and the code’s progress towards fulfilling these acts is regularly measured.
Ultimately, the simulation behaves successfully under certain conditions - viable courses of action are also proposed to investigate and remedy remaining flaws. A more detailed exposition of its functionality is provided in the form of annotated code (attached at the end of the report).
As well as using the lectures and advice of Professor Sarkisov and his research group to inform the code’s development, various online resources (including those noted at the end of this report) were valuable.
Discussion of simulation development
Week 1: learning the principles of molecular dynamics
I began this project by gaining some grounding in the theory of molecular dynamics: lectures on the topic delivered by Professor Sarkisov explained its core components, including the nature of microstates and the classical modelling of particle trajectories by numerically solving Newton’s laws of motion. They also expounded the practical procedure applied in simulations, as described below:
Describe the particles’ environment (including the number of particles, the box’s size, and initial configurations).
Apply the relevant rules for intermolecular interactions, using periodic bounding conditions to imitate and model bulk phase.
Determine the resultant changes in configuration using Newton’s laws of motion.
Calculate average properties.
Visualise the particles.
Repeat steps 2-5
It was also noted that checking a simulation’s conformity with the law of conservation of energy is a useful tool for gauging its accuracy.
Week 2: starting the development of a pygame-based simulation
Eager to visualise the particles’ motion, I commenced development of the simulation using pygame - a Python library used for animations and game creation, with which I had some experience. I attempted to fulfil tasks 1 and 2 of the simulation procedure described above, animating a unit cell containing Argon atoms and calculating the lennard-jones forces of interaction between them.
Unfortunately, my immediate focus on visualisation caused difficulty in the implementation of the nearest image convention: I sought to create images of each atom in 8 surrounding cells, and to calculate the Lennard-Jones (L-J) forces and potentials between the actual and image atoms when the images were closer to an actual atom than the actual atom it represents. Due to issues with mapping the true atoms to these 8 cells (specifically, copied rect objects were not mapped to their necessary locations, instead they remained superposed over the original atom from which they were copied), the code calculated the Lennard-Jones potential and forces of n lots of 8 overlapped atoms, where n is the number of atoms in the central cell (2 in the image below). This produced immense and unrealistic positive lennard-jones potentials. Moreover, scaling up the simulation by increasing the number of atoms was cumbersome in the pygame-based approach.
The image created by the pygame-based code
A review of the code also revealed unclear expository commenting practices, and Professor Sarkisov suggested using a cell-based coding system, such as IPython notebooks, to simplify debugging.
Week 3: starting the development of a non-visually focused simulation
As the visualisation-first pygame-based development strategy possessed major disadvantages, I adopted a mathematics-first approach, using numpy arrays to track particle positions and velocities in order to ease the management of multiple atoms in the simulated cell and facilitate the potential future addition of more particles.
The code randomly generates the initial positions of the particles using a continuous uniform distribution across the unit cell, and produces random initial speeds using the Maxwell-Boltzmann distribution for a specified temperature and acceptance-rejection sampling. The atoms’ initial velocities are generated using these speeds and angles taken from a continuous uniform probability distribution.
The Lennard-Jones potentials and forces of all atomic pairs were calculated, and instead of generating copies of each atom, the following mathematical method was used to apply the nearest image convention: dx is set to equal dx - L * int(2*dx/L), where dx is the distance in the x direction between particle i and particle j, L is the side length of the cell, and int() is a function which floors non-negative numbers and gives the ceiling of negative numbers. If L/2 < dx <L, L< 2*dx < 2L, so 1 < 2*dx/L < 2, so int(2*dx/L) = 1 and L*int(2*dx/L) = L. As a result, the distance between particle i and particle j is not measured, rather the distance between particle i and an ‘image’ of particle j, which is L units of distance (Angstroms) to the left of particle j (and hence less than L/2 units of distance away from particle i), is measured. The same process is applied to dy.
Interatomic force vectors were determined by multiplying their magnitude (which was found by calculating the derivative of the Lennard-Jones potential) by their corresponding unit vector (which was found by dividing the displacement vector between two particles by its magnitude).
This code was verified by setting up simple configurations of two particles, calculating the L-J forces and potentials with a calculator, and comparing the answers to the simulation’s outputs. The results were consistent with each other.
Particles’ coordinates / Å | Distance between particles / Å | U from code / kJmol^-1 (3sf) | U from calculator (3sf) | F magnitude from code / kJmol^-1Å^-1 (3sf) | F magnitude from calculator / kJmol^-1Å^-1 (3sf) |
(2,1), (6,1) | 4 | -0.937 | -0.937 | -0.554 | -0.554 |
(2, 7), (2, 3) | 4 | -0.937 | -0.937 | -0.554 | -0.554 |
(3, 6), (3, 8) | 2 | 2230 | 2230 | 13700 | 13700 |
(4, 4), (9, 4) | 5 | -0.355 | -0.355 | -0.380 | -0.380 |
The code successfully fulfils tasks 1 and 2 of the checklist from week 1, but does not advance the particles’ positions and velocities through time using Newton’s laws of motion.
Week 4: starting the development of particle advancement through time
Based on my experience with mechanics, I applied the following process to advance the particles’ positions and velocities through one timestep:
Calculate the L-J forces acting on each particle.
Divide the resultant force on each particle by its mass to find its acceleration (Newton’s 2nd law).
Change each particle’s position using ds = U*dt + 0.5*a*dt^2, where U is the particle’s initial velocity, a is its acceleration, ds is its change in position, and dt is the time step. We assume that the particles’ acceleration is constant for the duration of the time step.
Change each particle’s velocity using dv = a*dt, where dv is the change in the particle’s velocity.
Apply periodic boundary conditions - keep the particles in the unit cell.
There were two significant issues with this method and its implementation. Firstly, the Eulerian ‘SUVAT’ equations used in steps 3 and 4, although simple, are quite inaccurate (the global error is of the first order). Secondly, I used a timestep of 0.1s, which is far too large to accurately model particles moving at ~ 500 m/s in a cell that is about 100 Angstroms long - treating changes in the forces acting on the particles in that time period as negligible is highly unrealistic.
As a result, the simulation modelled the Argon atoms very inaccurately: energy was not conserved, and their speeds quickly ballooned (in one test they started at ~0.04 Å/10fs and ended at about 10 Å/10fs)
Week 5: improving particle advancement through time
To improve the simulation’s movement of particles through time, I used the Störmer-Verlet method instead of an Eulerian approach to decrease the local error to the 4th order and the global one to the 2nd order (https://www.algorithm-archive.org/contents/verlet_integration/verlet_integration.html). The timestep was also decreased to 10 fs. For each run, the total energy of the system (the sum of the particles’ kinetic and potential energy) remained constant for the first few hundred timesteps, but blew up towards the end. For example, in one test the system’s energy remained at about 5.9*10^-5 fJ for 800 iterations, and then grew to ~100 fJ by the end of the 1000 timesteps.
Week 6: creating visualisations and investigating the causes of errors
In order to gain greater insight into the lack of conservation of energy in some runs, I created a visualisation of the particles’ in the unit cell and a live-updating graph of the energy of the system using the matplotlib library.
Increases in the system’s energy coincided with one or more of two occurrences:
A particle leaves the unit cell and reappears on the opposite side of the cell. This correlates with the system’s energy spiking - rapidly increasing to about 34 fJ and then almost instantly falling back to close to its previous level. See video 1.
Two or more particles overlap. Occasionally, two particles will be close enough to each other and on such a trajectory that an advancement of one timestep using the Störmer-Verlet algorithm will cause them to overlap. This is unrealistic, and results in a very large L-J potential that causes the system’s total energy to increase significantly. See video 2.
Video 1: https://clipchamp.com/watch/9wAiuLVQz1b
The cause of phenomenon 1 is unclear, but the fact that the energy rises to the same level of about 34 fJ each time suggests that it is not caused by the distribution of the other particles in the unit cell, which varied at each observation of the spike. The fact that the energy rapidly falls back to near its previous state means that this phenomenon does not have a major negative impact on the overall accuracy of the simulation; such spikes can be ignored in any analyses of the simulation’s outputs.
If my theory of the cause of phenomenon 2 (that it is caused by two or more particles being on a trajectory that leads them to overlap) is correct, then phenomenon 2 could be prevented by decreasing the simulation’s timestep even further. If positional advancements were made near-instantaneously, particles’ van der Waal radii would overlap only slightly before they repelled each other away; particles would not be transported into a state where they overlapped enough to produce unrealistically immense L-J potentials. An issue with this approach is that the overlaps would have to be very tiny to avoid significantly increasing the system’s energy: the Lennard-Jones potential grows exceedingly quickly with minuscule decreases in r once it is positive (since the repulsive term sigma/r is taken to the twelfth power), and to facilitate sufficiently slight overlaps, the timestep would have to be infinitesimal - decreasing it from 10fs to 1fs did not have a major impact. It would then take a very long time for my machine (a typical laptop) to simulate significant particle movement.
Instead, the number of iterations and number of particles could be decreased to reduce the probability of two particles colliding and overlapping. This, however, reduces the amount of motion simulated. I thus decreased the number of iterations from 1000 to 600 and the number of particles from 10 to 6 to produce the successful simulation shown in video 1.
In regard to the visualisation of the particle movements, I first attempted to simulate the particles’ advancements through time and create an array storing their positions at each timestep, which would be used to animate the atoms’ motion later. Unfortunately, the contents of this array did not match the positions printed in each loop of the advancement function. The discrepancy may have been caused by difficulties with numpy arrays storing very many values; to rectify the situation, I updated the animation directly after each loop of the advancement function. This does, however, cause the simulation to have a longer running time - the particles have to be animated before they can advance by another timestep, instead of doing all of the advancements first and animating the whole motion later.
Proposed future action
To further improve the simulation and investigate causes of inaccuracy, the following steps could be taken.
Initialising the particles on an FCC lattice. The particles’ initial positions are currently configured using a uniform probability distribution; it is therefore possible for particles to spawn overlapping each other and prevent the simulation from functioning correctly. Using an FCC lattice would prevent this.
Calculating error accumulation of the Störmer-Verlet method and determining its growth over time to partially measure the simulation’s accuracy.
Determining the optimal box size: if it is too small, the particles will be very close together and are likely to overlap, if it is too large, a pair of particles could dominate the resultant forces on each other and just move into each other and overlap without other particles pulling them away.
Executing the simulation on a more powerful machine and using a near-instant timestep (several orders of magnitude smaller than a femtosecond). As discussed above, this could lead the simulation to more accurately reflect reality and would enable investigation of motion over longer periods of time, as more iterations could be executed without significantly increasing the chance of simulation failure due to overlapping atoms.
Rigorously and repeatedly testing the effects of increasing the number of particles and number of iterations on simulation results.
Conclusion
Under certain restrictions, including a small number of particles (n < 10 works well), a small number of iterations (N < 1000 works well), and a suitable initial configuration (with no overlapping atoms), the simulation functions successfully. The system’s energy remains roughly constant, and a useful visualisation of the Argon atoms’ trajectories is provided.
In addition to the potential improvements outlined above, more expansive future extensions could include the modelling of ions using Coulomb’s law, adding the third dimension, and modelling multi-atomic molecules, taking into account bond vibrations.
This project has substantially improved my related technical knowledge, including my understanding of Lennard-Jones potentials between electrically neutral particles, Newton’s laws of motion, and useful coding structures in python. I have also gained a strong appreciation for the importance of verifying results with manual checks, and providing clear expository comments on my code to ease collaboration.
https://drive.google.com/file/d/1GvBEV173doOmd84TdYKAiSk7wCs0daPP/view?usp=sharing
Comments
Post a Comment