Potentials for the SPC/Fw water simulation

For development of tools for molecular dynamics simulations we needed to implement one of the well-known models of molecules. For that purpose I chose the SPC/Fw model of water.

The molecule structure of $H_2O$ is preserved under action of two components of intramolecular potentials. The first one is a harmonic bond potential controlling distances between atoms of hydrogen and the central atom of oxygen. The second one is also a harmonic potential preserving the desired valence angle between those atomic bonds.

$$V^{intra}=\frac{k_b}{2}\Big[(r_{OH_1}-r^0_{OH})^2+(r_{OH_2}-r^0_{OH})^2\Big] + \frac{k_a}{2}(\theta_{\angle HOH}-\theta_{\angle HOH}^0)^2,$$

where $r^0_{OH}$ and $\theta_{\angle HOH}^0$ are the equilibrium bond length and angle, $k_a$ and $k_b$ are coefficients of the bond and valence angle elasticity.

Molecules of water interact by means of the (Lennard-Jones)[https://en.wikipedia.org/wiki/Lennard-Jones_potential] and (Coulomb)[https://en.wikipedia.org/wiki/Coulomb%27s_law] intermolecular potentials:

$$V^{inter}=\sum_{i,j}^{all\,pairs}\Big\{4\epsilon_{ij}\Big[\Big(\frac{\sigma_{ij}}{r_{ij}}\Big)^{12}-\Big(\frac{\sigma_{ij}}{r_{ij}}\Big)^{6}\Big]+k\frac{q_iq_j}{r_{ij}}\Big\},$$

where $\epsilon_{ij}$ and $\sigma_{ij}$ are the LJ potential parameters, $r_{ij}$ denotes the distance between oxygen atoms of the i-th and j-th molecules, $k$ is the electrostatic constant, $q_i$ is the charge of the i-th particle.

The implemented functions for the mentioned potentials and corresponding accelerations can be used in simulations of arbitrary molecules.

Example of simulation

First of all we need to define the constants for our system parameters. Here I used the OpenMM units of measurement:

const T = 298.16 # °K
const kb = 8.3144598e-3 # kJ/(K*mol)
const ϵOO = 0.1554253*4.184 # kJ 
const σOO = 0.3165492 # nm
const ρ = 997/1.6747# Da/nm^3
const mO = 15.999 # Da
const mH = 1.00794 # Da
const mH2O = mO+2*mH
const N = 216
const L = (mH2O*N/ρ)^(1/3)
const R = 0.9 # ~3*σOO  
const Rel = 0.49*L
const v_dev = sqrt(kb * T /mH2O)
const τ = 0.5e-3 # ps
const t1 = 0τ
const t2 = 300τ # ps
const k_bond = 1059.162*4.184*1e2 # kJ/(mol*nm^2)
const k_angle = 75.90*4.184 # kJ/(mol*rad^2)
const rOH = 0.1012 # nm
const ∠HOH = 113.24*pi/180 # rad
const qH = 0.41
const qO = -0.82
const k = 138.935458 #

The next steps include building of the simulation environment:

bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
lj_parameters = LennardJonesParameters(ϵOO, σOO, R)
e_parameters = ElectrostaticParameters(k, Rel)
spc_paramters = SPCFwParameters(rOH, ∠HOH, k_bond, k_angle)
pbc = CubicPeriodicBoundaryConditions(L)
water = WaterSPCFw(bodies, mH, mO, qH, qO,  lj_parameters, e_parameters, spc_paramters);
simulation = NBodySimulation(water, (t1, t2), pbc, kb);
result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)

After the calculation is done, we can use the result variable to access the properties of the simulated system. For example, let us collect data about energy of the system:

t = t1:τ:result.solution.t[end-1]
e_kin = kinetic_energy.(result, t)/(kb*T)
e_pot = potential_energy.(result, t)/(kb*T)

During simulation of water the law of energy conservation is satisfied:

Radial distribution function

The radial distribution function (RDF) is a very important tool in molecular dynamics.

(rs, grf) = rdf(result)

Here I present results of simulations of the two well-known systems. First of all, the RDF for liquid argon:

And the RDF for the water SPC/Fw model (216 particles, from 0 till 1.5 ps by 0.5 fs time step):

The mean squared displacement may be used for estimation of particles diffusion.

(ts, dr2) = msd(result)

Visualization in VMD

In order to export results of simulations, namely trajectories, I have implemented saving positions of particles at every timestep in a PDB file.

save_to_pdb(result, path_to_the_file )

Summary

All tools, that I develop for n-body simulations, are located in DiffEqPhysics.jl project. The test and examples folders contain examples of all physical interactions already implemented. At the moment of writing this post, the following items can be used in simulations:

• gravitational interaction
• electrostatic interaction
• magnetostatic interaction
• the Lennard-Jones potential
• custom pair-wise interaction implemented via the CustomAccelerationSystem interface
• harmonic bonds between particles
• harmonic potential for valence angle between atoms
• radial distribution function
• mean squared displacement
• Andersen thermostating

Pull requests for this post: #32