{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import additional_code as code\n",
"import warnings\n",
"from pylj import md, sample, comp, examples, util\n",
"warnings.filterwarnings('ignore')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ISIS Neutron Training Course\n",
"## Computational Simulation for Neutron Scattering\n",
"\n",
"#### Andrew R. McCluskey \n",
"##### University of Bath/Diamond Light Source - arm61@bath.ac.uk\n",
"\n",
"2018-03-08\n",
"\n",
">The speed and development of computers is now so rapid, and the advances in modelling and informatics are so dramatic that is 15 years' time no chemist will be doing an experiments at the bench without trying to model them first.\n",
"\n",
"> **Prof. Dominic Tildesley, Former RSC President, 2013**\n",
"\n",
"### Introduction \n",
"\n",
"Molecular dynamics simulation is a common and powerful technique in computational chemistry/statistical physics for the determination of equilibrium & transport properties of many-body systems. These systems consist of particles interacting through force-fields, whose motions are described by classical Newtonian mechanics, for molecular dynamics this comes down to two equations,\n",
"\n",
"$$ \\mathbf{r}(t+\\Delta t) = \\mathbf{r}(t) + \\mathbf{v}(t)\\Delta t + \\dfrac{1}{2} \\mathbf{a}(t)\\Delta t^2, $$\n",
"\n",
"where, $\\mathbf{r}(t)$ is the position, $\\mathbf{v}(t)$ is the velocity, and $\\mathbf{a}(t)$ is the acceleration of the particle at a given time $t$, while $\\mathbf{r}(t+\\Delta t)$ is the position after the timestep, $\\Delta t$, and, \n",
"\n",
"$$ \\mathbf{f}(t) = m\\mathbf{a}(t), $$\n",
"\n",
"where $\\mathbf{f}(t)$ is the force on some particle at time $t$, and $m$ is the particle mass.\n",
"\n",
"Generally speaking the use of classical mechanics compared to quantum mechanics is a very good approximation for many molecular materials, where quantum mechanical effects are weak -- these only become significant for very light atoms, e.g. $\\text{H}_2$ or $\\text{He}_2$ molecules. \n",
"\n",
"A molecular dynamics simulation can be considered in a similar way to a real experiment. A general recipe for a real experiment would be: \n",
"\n",
"1. Prepare the sample -- buy the constituent chemicals and mix them approporiately to get your single crystal of a photovoltaic material, grow and purify the protein in deuterared medium, etc. \n",
"2. Allow the sample to interact with some measuring instrument -- fly to the ILL and put the sample on D9, bring it to ISIS and put it on SANS2D.\n",
"3. Finally, measure some property of the sample for a given time -- collect the diffraction, or SANS pattern until the counting is good. \n",
"\n",
"A similar recipe can be applied to a molecular dynamics simulation:\n",
"\n",
"1. Prepare the sample -- build the particle coordinates, determine a suitable force-field and equilibriate the sample by solving Newton's equations of motion until some property of the system no longer changes beyond random fluctuations.\n",
"2. Allow the sample to interact with some measuring instrument -- within molecular dynamics this involves determining a quantity that can be measured as a function of the positions, velocities, and accelerations of the particles. \n",
"3. Finally, measure that property of the sample for a given time -- this means calculating Newton's equations of motion iteratively, for a long enough time that good statistics can be obtained of the measurement quantity. \n",
"\n",
"It is important to be aware that this means that the problems present in real experiments are possible in molecular dynamics simulations; preparing the sample incorrectly by using an inaccurate force field, or measuring for too short a time. \n",
"\n",
"For using molecular dynamics simulation to understand neutron scattering, the property that is being measured in the simulation should be as similar as possible to that measured in the neutron experiment. For example, small angle neutron scattering (SANS) measured the position of clusters of atoms with respect to each other, therefore to determine a SANS profile from molecular simulation we must sample the atomic coordinates during the molecular dynamics simulation. \n",
"\n",
"### Force-fields \n",
"\n",
"In molecular dynamics, the electrons are not explicitly simulated, instead the system is modelled as a series of point particles that interact based on some *force-field*. The force-field is a function that describes the potential energy, $u$, for a given configuration of the system, typically for a molecular system this consists of bonded and non-bonded components, \n",
"\n",
"$$ u = u_{\\text{bonded}} + u_{\\text{non-bonded}} $$\n",
"\n",
"The potential energy depends on the functions used, their parameterisation and the configuration of the system at a given time. The functions and their parameterisation depend on the underlying chemistry of the system, the particular functions that are described herein are those associated with the OPLS2005 force field.[1]\n",
"\n",
"#### Non-bonded interactions \n",
"\n",
"In the OPLS2005 force-field, two forms of non-bonded interaction are present, van der Waals and Coulombic. The van der Waals interaction is modelled with a Lennard-Jones potential, and has the following form, \n",
"\n",
"$$ u_{\\text{vdW}} = 4\\epsilon_{ij}\\Bigg[\\bigg(\\frac{\\sigma_{ij}}{\\mathbf{r}_{ij}}\\bigg)^{12} - \\bigg(\\frac{\\sigma_{ij}}{\\mathbf{r}_{ij}}\\bigg)^6\\Bigg], $$\n",
"\n",
"where, $\\epsilon_{ij}$ and $\\sigma_{ij}$ are parameters specific to the interaction, hence the atom types, such that the van der Waals interaction vaires with distance between atoms, $\\mathbf{r}_{ij}$. The Coulombic interaction is that between ions, modelled with a Coulombic potential, \n",
"\n",
"$$ u_{\\text{Coulombic}} = \\frac{q_iq_j}{4\\pi\\epsilon\\mathbf{r}_{ij}} $$\n",
"\n",
"where, $q_i$ and $q_j$ are the charges on the interacting particles -- it should be noted that the charge does not need to be an integer essentially allowing for the presense of polarity, and $\\epsilon$ is the permittivity of the solvent. The Lennard-Jones, Coulombic, and total non-bonded potentials for a sodium-chloride ion pair are determined an shown below. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"sigma_ij = 3.7552\n",
"epsilon_ij = 0.0756\n",
"\n",
"r = np.arange(0.01, 10.01, 0.01)\n",
"\n",
"u_vdW = 4. * epsilon_ij * ((sigma_ij / r) ** 12 - (sigma_ij / r) ** 6)\n",
"u_Coulombic = (-1 * 1) / (4. * np.pi * r) / (1.036E-2) # The division is to correct the units\n",
"u = u_vdW + u_Coulombic\n",
"\n",
"\n",
"fig, ax = code.nonbond_plots()\n",
"plt.plot(r, u_vdW, label='$u_{vdW}$')\n",
"plt.plot(r, u_Coulombic, label='$u_{Coulombic}$')\n",
"plt.plot(r, u, label='$u_{Total}$')\n",
"plt.legend() \n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Bonded interactions\n",
"\n",
"The bonded interactions describe the nature of the models used to form chemical bond between particles. In the OPLS2005 force-field, these generally consist of bond stretchs, angle bends, and dihedral flexes; with the functional forms, \n",
"\n",
"$$ \\begin{array} \n",
"Uu_{\\text{bonded}} & = \\sum_{\\text{bonds}}K_b(b-b_0)^2 + \\sum_{\\text{angles}} K_\\theta(\\theta - \\theta_0)^2 \\\\\n",
" & = \\sum_{\\text{dihedrals}}\\dfrac{1}{2}\\Big\\{A_1\\big(1 + \\cos{[\\phi]}\\big)+A_2\\big(1-\\cos{[2\\phi]}\\big) + A_3\\big(1+\\cos{[3\\phi]}\\big)\\Big\\}\\end{array}$$\n",
" \n",
"where $K_b$, $b_0$, $K_\\theta$, $\\theta_0$, $A_1$, $A_2$, and $A_3$ are potential model dependent parameters for the bonds, angles, and dihedrals respectively, while $b$, $\\theta$, and $\\phi$ are the bond length and siize of the angles and dihedrals in the system. It can be seen that both the bond and angle interactions are modeled with harmonic functions, while the dihedral angle interaction consists of a series of cosines. The form of each of these potentials for a chain of carbon atoms are shown below. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"fig, ax = code.bond_plots()\n",
"\n",
"K_b = 1121.3\n",
"b_0 = 1.529\n",
"\n",
"b = np.arange(1, 2, 0.001)\n",
"u_b = K_b * (b - b_0) ** 2\n",
"\n",
"ax[0].plot(b, u_b)\n",
"\n",
"K_theta = 244.1\n",
"theta_0 = 112.7\n",
"\n",
"theta = np.arange(110, 120, 0.001)\n",
"u_theta = K_theta * (theta - theta_0) ** 2\n",
"\n",
"ax[1].plot(theta, u_theta)\n",
"\n",
"A1 = 1.135\n",
"A2 = -0.151\n",
"A3 = 0.400\n",
"\n",
"phi = np.arange(0, 360, 0.01)\n",
"u_phi = 0.5 * (A1 * np.cos(np.deg2rad(phi)) + A2 * np.cos(np.deg2rad(2*phi)) + A3* np.cos(np.deg2rad(3*phi)))\n",
"\n",
"ax[2].plot(phi, u_phi)\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Off-the-shelf force-fields\n",
"\n",
"Many \"off-the-shelf\"force-fields exist, which have been parametrised, from experimental data or quantum calculations. These should be used with caution, as the application of an off-the-shelf force-field can lead to poor agreement with experimental data, if the new system is significantly different from those used in the parameterisation process -- generally the off-the-shelf force-field should be used to reproduce some simple experiment parameters before use, e.g. density. Common off-the-shelf force-fields include: \n",
"- OPLS2005/OPLS-AA -- this is the Optimised Potential for Liquid Simulations, a general purpose force-field for the simulation of liquid systems (OPLS2005 is an update to OPLS-AA). \n",
"- CHARMM -- Chemistry at HARvard Molecular Mechanics, widely used for small molecules and proteins. \n",
"- GROMOS -- the GROningen MOlecular Simulation comes with the GROMOS/GROMACS software, particularly useful for biomolecules. \n",
"- AMBER -- the Assisted Model Building and Energy Refinement, again popular for proteins and DNA. \n",
"- MARTINI -- this is a coarse-grained force field where each 4 heavy atoms ($Z > 1$) is represented with one coarse-grained \"bead\", parameterisation to reproduce lipid and protein interactions. \n",
"- Many, many others for many types of application. \n",
"\n",
"### Periodic boundary conditions\n",
"\n",
"The purpose of molecular simulation for neutron scattering is to reproduce the properties of some experimental system. However, current computational systems cannot efficiently model such large systems -- limited to about 10s of millions of particles ($1\\times10^{-17}$ moles).[2] Therefore simplications must be made; the most routinely applied is that of the periodic box. This is where a boundary condition is applied to the edges of the simulation box to mimic an infinite system -- by creating the situation that the simulation box is surrounded by infinite identical images of itself.[3] This is achieved by having a particle that reaches the edge of the simulation box appear on the other side of the box; as though it came in from the adjecent periodic image. Additionally, when the force is calculated on each particle, the periodic nature is considered and if particles are more than half of the box apart, it is recognised that they are closer by replacing one of the particles with that from the adjacent periodic image. This shorter distance is then used in the force calculation.\n",
"\n",
"Below is an interactive example that shows a single particle that has been given some initial velocity in a random direction. It will travel across the periodic boundary of the simulation box appearing on the other side. The first number (initially 10000) is the number of steps in the simulation, while the second number (initially 50) is the temperature of the particle, in Kelvin. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"code.periodic_boundary(10000, 50)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A molecular dynamics program\n",
"\n",
"One of the best ways to understand how a molecular dynamics simulation works is to consider how a molecular dynamics program would work. Below is an example of a molecular dynamics simulation, of Lennard-Jones particles interacting in a 2D simulation box. The property we are sampling is the radial distribution function of each particle. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"# Assigning some parameters\n",
"number_of_particles = 100\n",
"temperature = 273.15\n",
"number_of_steps = 10000\n",
"sample_frequency = 10\n",
"box_length = 80\n",
"\n",
"# Initialise the system\n",
"system = md.initialise(number_of_particles, temperature, box_length, 'square')\n",
"# This sets the sampling object\n",
"sample_system = sample.RDF(system)\n",
"# Start at time 0\n",
"system.time = 0\n",
"# Begin the molecular dynamics loop\n",
"for i in range(0, number_of_steps):\n",
" # At each step, calculate the forces on each particle and get acceleration\n",
" system.particles, system.distances, system.forces, system.energies = comp.compute_forces(system.particles,\n",
" system.box_length, \n",
" system.cut_off)\n",
" # Run the equations of motion integrator algorithm\n",
" system.particles = md.velocity_verlet(system.particles, system.timestep_length, system.box_length, \n",
" system.cut_off)\n",
" # Sample the parameters of interest for the simulation\n",
" system = md.sample(system.particles, system.box_length, system.initial_particles, system)\n",
" # Iterate the time and step\n",
" system.time += system.timestep_length\n",
" system.step += 1\n",
" # At a given frequency sample the positions and plot the RDF\n",
" if system.step % sample_frequency == 0:\n",
" sample_system.update(system)\n",
"sample_system.average()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Initialisation\n",
"\n",
"The initialisation process involves getting the starting positions and velocities for the particles in the system. Generally speaking the initial positions are system dependent and can be obtained from a source or generated using appropriate software. For example, the simulations of proteins often involves taking a starting structure from the Protein Data Bank [4] crystallographic information, this requires a crystal structure of the protein to have been resolved. Other methods would include generating energy minimised structures *in-silico*, e.g. with Jmol, [5] and orienting the structures appropriately using a program such as Packmol. [6]\n",
"In the above 2D example the initial positions were generated by placing the particles on a square lattice, however the particles could also be placed in random positions. However, this may be problematic as some particles may be too close to others, resulting in the simulation \"exploding\". Therefore it is necessary to preform some energy minimsation which usually takes the form of a steepest descent algorithm. This is not discussed in detail herein, however Frenkel and Smith [3] has detail for the interested reader. \n",
"\n",
"The generation of initial velocities is a more general process than the positions. The velocities are randomly selected from a uniform distribution, $\\mathbf{v}_i(0) \\in [-0.5, 0.5]$, it is necessary to then shift them such as to obtain a total momentum of zero. The velocites are also scaled by a factor, that includes the simulation temperature, $T$.\n",
"\n",
"$$ \\mathbf{v}_i(0) = \\sqrt{2 T}\\times\\cos\\big[2\\pi \\mathbf{v}_i(0)\\big]. $$\n",
"\n",
"There is no need to intitialise the velocites with a Maxwell-Boltzmann distribution as this will happen naturally within the simulation. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"v_i = np.random.uniform(-0.5, 0.5, 10000)\n",
"temperature = 1.\n",
"fig, ax = code.velocity_plot()\n",
"v_i = (v_i - np.average(v_i)) * np.sqrt(3 * temperature / (np.average(v_i) ** 2.))\n",
"ax.hist(v_i, histtype='step', lw=2)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The final step of the initialisation process is to determine the particle positions before the first timestep, e.g. at $t = -\\Delta t$. This can be achieved using the following relation, \n",
"\n",
"$$ \\mathbf{r}_i(-\\Delta t) = \\mathbf{r}_i(0) - \\mathbf{v}_i(0)\\Delta t. $$\n",
"\n",
"Having initialised both the positions and velocites we can now start the iterative procedure of the molecular dynamics simulation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Force calculations\n",
"\n",
"The calculation of the forces on each of the atoms is the most computationally time consuming part of a molecular dynamics simulation. This is due to the fact that we must consider the pairwise interactions over all of the particles in our system -- this means calculating the force for $\\frac{(N^2 - N)}{2}$ pairs. This results in the speed of the algorithm associated with the force calculation scale as order $N^2$ meaning that quite quickly the simulation can become size limited. \n",
"\n",
"The force on a particle is equal to the negative of the gradient of the potential energy for a given system configuration, this means that the force is found from the first derivative of the potential energy function, \n",
"\n",
"$$ \\mathbf{f} = -\\frac{\\delta u}{\\delta \\chi} $$\n",
"\n",
"where, $\\mathbf{f}$ is the force on the particle, $u$ is the particle potential energy -- from the force-field, and $\\chi$ is the current configuration of the particle. This calculation is achieved by finding the distance in the $x$, $y$, and $z$ dimensions for every pair of particles, $i$ and $j$. If the particles are further apart in a particular direction than half the length of the simulation box in that dimension the periodic boundary correction is applied; where one of the particles is swapped with that from an adjacent image, such that the particles are closer together. At this stage, the system cut-off is also applied, this is where if the particles are beyond a certain distance apart, the force that they apply on each other is considered to be effectively zero so the energy and force are not calculated -- this is a method of improving the speed of the force calculation. \n",
"\n",
"#### Equations of motion integration \n",
"\n",
"Due to the fact that the potential energy is a function of the entire system configuration ($3N$ in 3-dimensions), Newton's equations of motion are no longer analytically solvable -- instead an approximate Newtonian integrator must be used. Many different integrators exist, the most common come from the Velocity family, the simplest of which is the Verlet algorithm, \n",
"\n",
"$$ \\mathbf{r}_i(t+\\Delta t)\\approx 2\\mathbf{r}_i(t) - \\mathbf{r}_i(t-\\Delta t) + \\dfrac{\\mathbf{f}_i(t)}{m_i}\\Delta t^2 $$\n",
"\n",
"where, $\\mathbf{r}(t+\\Delta t)$ is the new particle position, $\\mathbf{r}(t)$ is the current particle position, $\\mathbf{r}(t-\\Delta t)$ is the previous particle position, $\\mathbf{f}(t)$ is the force calculated in the current step, and $m$ is the particle mass. \n",
"\n",
"This algorithm does not compete the velocity explicitly; calculation of the kinetic energy requires the velocity, $\\mathbf{v}(t)$, therefore it can be found from knowledge of the trajectory, \n",
"\n",
"$$ \\mathbf{v}_i(t) = \\dfrac{\\mathbf{r}_i(t+\\Delta t) - \\mathbf{r}_i(t - \\Delta t)}{2 \\Delta t}. $$\n",
"\n",
"The kinetic energy -- therefore velocity, is used to find the instantaneous temperature, $T_{\\text{inst}}$, of the system, this is important in the constant temperature ensembles, \n",
"\n",
"$$ T_{\\text{inst}} = \\dfrac{\\sum_{i=1}^{N} \\mathbf{v}_i(t)}{3N} $$\n",
"\n",
"Following the integration, the system configuration is sampled for some desired information at a given frequency. \n",
"\n",
"Within molecular dynamics, one of the significant limitations is the size of the time-step, $\\Delta t$. Generally speaking the time-step should be less than the highest frequency vibration in the system to ensure that the force-field is modelling the system accurately. For simulations containing explicitly modelled hydrogen atoms, the time-step should be 1 fs as a maximum. Larger time-steps can be accessed using coarse-grained force-fields such as MARTINI, or by fixing the positions of the hydrogen atoms with respect to the heavy atom to which it is bound.\n",
"\n",
"#### Ensembles \n",
"\n",
"The molecular dynamics simulation that is discussed above was in the NVE ensemble, where the number of particles (N), volume of the simulation box (V), and energy of the system (E) are held constant. Other ensembles exist, such as the NVT ensemble, where T indicates constant temperature, or NPT, constant pressure (P). \n",
"\n",
"In an NVT simulation, the simulation is coupled with a heat bath at a given temperature, for the Anderson thermostat, the simulation algorithm is modified as follows:\n",
"\n",
"1. The initial set of positions and velocities of the particles are generated as normal.\n",
"2. After the first step, $\\Delta t$, a number of particles are selected to undergo a collision with a heat bath. The probability that a particle is selected depends on some coupling frequency.\n",
"3. If particle $i$ is selected to undergo a collision, a new velocity is selected from a Maxwell-Boltzmann distribution corresponding to the desired temperature, $T$. \n",
"\n",
"This allows the system to stay approximately equal to the defined temperature, with random fluctuations in the instanteous pressure. The below simulation shows the effect of the heat bath within the Velocity Verlet algorithm. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"# Assigning some parameters\n",
"number_of_particles = 100\n",
"temperature = 273.15\n",
"number_of_steps = 10000\n",
"sample_frequency = 10\n",
"box_length = 80\n",
"\n",
"# Initialise the system\n",
"system = md.initialise(number_of_particles, temperature, box_length, 'square')\n",
"# This sets the sampling object\n",
"sample_system = sample.Interactions(system)\n",
"# Start at time 0\n",
"system.time = 0\n",
"# Begin the molecular dynamics loop\n",
"for i in range(0, number_of_steps):\n",
" # At each step, calculate the forces on each particle and get acceleration\n",
" system.particles, system.distances, system.forces, system.energies = comp.compute_forces(system.particles,\n",
" system.box_length, \n",
" system.cut_off)\n",
" # Run the equations of motion integrator algorithm\n",
" system.particles = md.velocity_verlet(system.particles, system.timestep_length, system.box_length, \n",
" system.cut_off)\n",
" # Sample the parameters of interest for the simulation\n",
" system = md.sample(system.particles, system.box_length, system.initial_particles, system)\n",
" # Allow the system to interate with a heat bath\n",
" system.particles = comp.heat_bath(system.particles, system.temperature_sample, temperature)\n",
" # Iterate the time and step\n",
" system.time += system.timestep_length\n",
" system.step += 1\n",
" # At a given frequency sample the positions and plot the RDF\n",
" if system.step % sample_frequency == 0:\n",
" sample_system.update(system)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Molecular dynamics simulation packages\n",
"\n",
"Similar to the number of force-fields, there are also a wide range of molecular dynamics software packages available -- the choice of software can depend on personal preference and what is being simulated. Common molecular dynamics software include:\n",
"\n",
"- LAMMPS -- possible the most common molecular dynamics package due to the flexibility offered, popular amoung hardcore theoreticians.\n",
"- NAMD -- popular package from the USA for the simulation of biological systems in particular.\n",
"- GROMACS -- effectively the European challenger to NAMD, popular for coarse-grained simulation due to the inclusion of the MARTINI force-field. \n",
"- DL_POLY -- developed through CCP5-EPSRC funding, popular for the development of force-fields. \n",
"\n",
"### Application to neutron scattering\n",
"\n",
"Shown below is a molecular dynamics simulation, from which the neutron scattering/diffraction profile is being determined. The scattering is found from a Fourier transform of the radial distribution function, however could also be found by the Debye relation. \n",
"\n",
"$$ I(q) = \\sum_j^N \\sum_k^N b_jb_k \\frac{sin(q|\\mathbf{r}_j - \\mathbf{r}_k|)}{q|\\mathbf{r}_j - \\mathbf{r}_k|}, $$\n",
"\n",
"where $b_j$ and $b_k$ are the scattering length densities of particles $j$ and $k$, $\\mathbf{r}_j$ and $\\mathbf{r}_k$ are the positions of particles $j$ and $k$, $q$ is the scattering vector, and $I(q)$ is the determined intensity. It is possible to use this to observe phase transitions in simulation by varying the temperature, pressure or volume and considering the sharpness of the peaks in the scattering profile. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"# Assigning some parameters\n",
"number_of_particles = 100\n",
"temperature = 600\n",
"number_of_steps = 1000\n",
"sample_frequency = 10\n",
"box_length = 30\n",
"\n",
"# Initialise the system\n",
"system = md.initialise(number_of_particles, temperature, box_length, 'square')\n",
"# This sets the sampling object\n",
"sample_system = sample.Scattering(system)\n",
"# Start at time 0\n",
"system.time = 0\n",
"# Begin the molecular dynamics loop\n",
"for i in range(0, number_of_steps):\n",
" # At each step, calculate the forces on each particle and get acceleration\n",
" system.particles, system.distances, system.forces, system.energies = comp.compute_forces(system.particles,\n",
" system.box_length, \n",
" system.cut_off)\n",
" # Run the equations of motion integrator algorithm\n",
" system.particles = md.velocity_verlet(system.particles, system.timestep_length, system.box_length, \n",
" system.cut_off)\n",
" system = md.sample(system.particles, system.box_length, system.initial_particles, system)\n",
" # Allow the system to interate with a heat bath\n",
" system.particles = comp.heat_bath(system.particles, system.temperature_sample, temperature)\n",
" # Iterate the time and step\n",
" system.time += system.timestep_length\n",
" system.step += 1\n",
" # At a given frequency sample the positions and plot the RDF\n",
" if system.step % sample_frequency == 0:\n",
" sample_system.update(system)\n",
"sample_system.average()"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"1. J. L. Banks, *et al.* , *J. Comput. Chem.*, 2005, **26**, 1752–1780.\n",
"2. K. Y. Sanbonmatsu and C. S. Tung, *J. Struct. Biol*, 2007, **157**, 470-480.\n",
"3. D. Frenkel and B. Smit, *Understanding Molecular Simulation: From Algorithms to Applications*, Academic Press, San Diego, USA, 2nd Edn., 2002. \n",
"4. *RCSB Protein Data Bank -- RCSB PDB*, http://www.rscb.org/pdb/home/home.do, Accessed 2017-03-02.\n",
"5. *Jmol: an open-source Java viewer for chemical in 3D*, http://www.jmol.org/, Accessed: 2017-03-02.\n",
"6. L. Martínez, *et al.*, *J. Comput. Chem*, 2009, **30**, 2157-2164."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:pylj]",
"language": "python",
"name": "conda-env-pylj-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": false,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}