In this example, we show how to perform molecular dynamics of bulk water using the popular interatomic TIP3P potential (W. L. Jorgensen et. al.) and LAMMPS.
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
from pyiron.project import Project
import ase.units as units
pr = Project("tip3p_water")
We will setup a cubic simulation box consisting of 27 water molecules density density is 1 g/cm$^3$. The target density is achieved by determining the required size of the simulation cell and repeating it in all three spatial dimensions
density = 1.0e-24 # g/A^3
n_mols = 27
mol_mass_water = 18.015 # g/mol
# Determining the supercell size size
mass = mol_mass_water * n_mols / units.mol # g
vol_h2o = mass / density # in A^3
a = vol_h2o ** (1./3.) # A
# Constructing the unitcell
n = int(round(n_mols ** (1. / 3.)))
dx = 0.7
r_O = [0, 0, 0]
r_H1 = [dx, dx, 0]
r_H2 = [-dx, dx, 0]
unit_cell = (a / n) * np.eye(3)
water = pr.create_atoms(elements=['H', 'H', 'O'],
positions=[r_H1, r_H2, r_O],
cell=unit_cell)
water.set_repeat([n, n, n])
water.plot3d()
Failed to display Jupyter Widget of type NGLWidget
.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
water.get_chemical_formula()
'H54O27'
The initial water structure is obviously a poor starting point and requires equilibration (Due to the highly artificial structure a MD simulation with a standard time step of 1fs shows poor convergence). Molecular dynamics using a time step that is two orders of magnitude smaller allows us to generate an equilibrated water structure. We use the NVT ensemble for this calculation:
job_name = "water_slow"
ham = pr.create_job("Lammps", job_name)
ham.structure = water
# Listing available potentials for this structure
ham.list_potentials()
['H2O_tip3p', 'H2O_tip3p_slab']
We will use the H2O_tip3p
potential. The H2O_tip3p_slab
should be used if you want to switch off the periodic boundary conditions along the z-axis
ham.potential = 'H2O_tip3p'
ham.calc_md(temperature=300,
n_ionic_steps=1e4,
time_step=0.01)
ham.run()
view = ham.animate_structure()
view
Failed to display Jupyter Widget of type NGLWidget
.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
At the end of this simulation, we have obtained a structure that approximately resembles water. Now we increase the time step to get a reasonably equilibrated structure
# Get the final structure from the previous simulation
struct = ham.get_structure(iteration_step=-1)
job_name = "water_fast"
ham_eq = pr.create_job("Lammps", job_name)
ham_eq.structure = struct
ham_eq.potential = 'H2O_tip3p'
ham_eq.calc_md(temperature=300,
n_ionic_steps=1e4,
n_print=10,
time_step=1)
ham_eq.run()
view = ham_eq.animate_structure()
view
Failed to display Jupyter Widget of type NGLWidget
.
If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean that the widgets JavaScript is still loading. If this message persists, it likely means that the widgets JavaScript library is either not installed or not enabled. See the Jupyter Widgets Documentation for setup instructions.
If you're reading this message in another frontend (for example, a static rendering on GitHub or NBViewer), it may mean that your frontend doesn't currently support widgets.
We can now plot the trajectories
plt.plot(ham_eq["output/generic/energy_pot"])
plt.xlabel("Steps")
plt.ylabel("Potential energy [eV]");
plt.plot(ham_eq["output/generic/temperature"])
plt.xlabel("Steps")
plt.ylabel("Temperature [K]");
We will now use the get_neighbors()
function to determine structural properties from the computed trajectories. We take advantage of the fact that the TIP3P water model is a rigid water model which means the neighbors of each molecule never change. Therefore they need to be indexed only once
final_struct = ham_eq.get_structure(iteration_step=-1)
# Get the indices based on species
O_indices = final_struct.select_index("O")
H_indices = final_struct.select_index("H")
# Getting only the first two neighbors
neighbors = final_struct.get_neighbors(num_neighbors=2)
Every O atom has two H atoms as immediate neighbors. The distribution of this bond length is obtained by:
bins = np.linspace(0.5, 1.5, 100)
plt.hist(neighbors.distances[O_indices].flatten(), bins=bins)
plt.xlim(0.5, 1.5)
plt.xlabel("d$_{OH}$ [$\AA$]")
plt.ylabel("Count");
We need to extend the analysis to go beyond nearest neighbors. We do this by controlling the cutoff distance
neighbors = final_struct.get_neighbors(cutoff=8)
neigh_indices = np.hstack(neighbors.indices[O_indices])
neigh_distances = np.hstack(neighbors.distances[O_indices])
One is often intended in an element specific pair correlation function. To obtain for example, the O-O coordination function, we do the following:
# Getting the neighboring Oxyhen indices
O_neigh_indices = np.in1d(neigh_indices, O_indices)
O_neigh_distances = neigh_distances[O_neigh_indices]
bins = np.linspace(1, 5, 100)
count = plt.hist(O_neigh_distances, bins=bins)
plt.xlim(2, 4)
plt.title("O-O pair correlation")
plt.xlabel("d$_{OO}$ [$\AA$]")
plt.ylabel("Count");