This article illustrates the use of the Electric Permittivities Analysis library (EPA) for computing static dielectric properties from molecular dynamics trayectories.
This is not a comprensive review of the EPA capabilities. The goal here is to show how to use this Python library in some typical circumstances.
If you use this software for your projects, please cite the following articles:
To follow this article, it is necessary to have a basic understanding of
A 2 ns long MD simulation of a mixture of 511 water molecules (SPC/E) and 50 ethanol molecules (OPLS), at 298.15 K and 1 bar, will be used as an example.
The library is open source and can be download from (a single .py file) from https://gitlab.com/herchem/epa.-electric-permittivities-analysis.
The EPA library must be imported in order to be used. There are not special requirements. For example, I copied the epa.py file into the work/code/lib subdirectory of my home dir, and imported it using the following code
# Add path to the EPA library and import it
from os.path import expanduser
import sys
home = expanduser("~")
dip_dir = home+"/work/code/lib"
sys.path.append(dip_dir)
import epa
We will import Numpy for working with arrays.
import numpy as np
Data can be extracted from Gromacs trajectory files (.trr or .xtc files). The geometry file (PDB or .gro files) is always needed. The .tpr file is also requiered for extracting dipole moments.
The following example shows how to extract data from a MD trajectory. We will use results of our example simulation. The data is extracted as Numpy arrays and saved to disk.
# Get volumes from Trajectory
volumes = epa.get_volumes("trayectoria.trr", "geometria.gro")
np.save("volumes", volumes)
# Get times from Trajectory
times = epa.times_to_array("trayectoria.trr", "geometria.gro")
np.save("times", times)
# Get the total dipole moment of each stored simulation step
m = epa.dipole_moments("trayectoria.trr", "geometria.gro", tpr_file = "entrada_dinamica.tpr", charges_from = "entrada_dinamica.tpr", chunk = 500, skip = 0, stride = 1)
np.save("m", m)
# Get the water dipole moment of each stored simulation step
m_H2O = epa.dipole_moments("trayectoria.trr", "geometria.gro", tpr_file = "entrada_dinamica.tpr", charges_from = "entrada_dinamica.tpr", chunk = 500, skip = 0, stride = 1, residues="HOH")
np.save("m_H2O", m_H2O)
The code is pretty transparent. It iterates over the simulation frames in blocks of chunk frames (see the documentation of the MDTraj library for details). The skip and stride allow skip frames inside a chunk, they are rarely needed.
Let's start with the fluctuations route. We will need the dipole moments and volumes from a MD performed without applying an electric field.
Now we will load the volumes, times and dipole moments of the system and the water, corresponding to each stored simulation step. We may wish to discard many initial points (here 1000 poins = 200.0 ps) to make sure that the system reached the equilibrium state, and make a time translation in order to asign the zero value to the first time.
# skip 'skip' points due to equilibration
skip = 1000
# Load volumes
v = np.load('volumes.npy', 'r')[skip:]
# Load times and set intitial time to zero
t = np.load('times.npy', 'r')
t = t[skip:]-t[skip]
# Load total and water dipole moments, and compute for ethanol
m = np.load('m.npy', 'r')[skip:]
m_w = np.load('m_H2O.npy', 'r')[skip:]
m_et= m - m_w
The system dielectric constant satisfies
\[ \varepsilon_r = 1+ \frac{\text{Var}[\bold M]}{3\varepsilon_0 V k_b T}\]where \(\bold M \) is the total dipole moment of the system. However, the variance can only be estimated. The commonly used estimators are \( \langle \bold M ^2 \rangle - \langle \bold M \rangle ^2 \) and \( \langle \bold M ^2 \rangle \). The estimator used is chosen by the unbiased optional argument (default true) of the e_s_v2 function which is used for the calculation of the dielectric constant.
### Calculations on the entire system ###
## Fluctuation Route (FR)
e_FR = epa.e_s_v2(m[:,:], np.copy(v), 298.15, coords=[1.,1.,1.], unbiased=False)
# unbiased uses <M^2> - <M>^2 as the variance estimator (else, <M^2> is used)
# coords specifies the considered directions e.g. coords=[1.,0.,0.] compute only for dir x
# Note: Centering m has not effect in this example because it does not change this estimation of the variance
We will also compute the permittivity in each spatial direction, \(x\), \(y\) and \(z\), and also we will check that their mean equals the system estimation of the system permittivity.
e_FR_x = epa.e_s_v2(m[:,:], np.copy(v), 298.15, coords=[1.,0.,0.], unbiased=False)
e_FR_y = epa.e_s_v2(m[:,:], np.copy(v), 298.15, coords=[0.,1.,0.], unbiased=False)
e_FR_z = epa.e_s_v2(m[:,:], np.copy(v), 298.15, coords=[0.,0.,1.], unbiased=False)
print(e_FR)
print(e_FR_x, e_FR_y, e_FR_z, "Mean:", (e_FR_x+e_FR_y+e_FR_z)/3)
The output is the following
OUTPUT:
63.90882720991493
63.16907630500186 70.56741481618653 57.98999050855643 Mean: 63.90882720991493
Dielectric permittivities require many more simulation steps than other properties. Therefore, knowing the uncertainty of the calculations is particularly useful. We proposed methods for the computation of the uncertainty in the estimation of dielectric constant computed throuhg the Fluctuation and other routes [1,2]. They are more accurate than a common numerical approach obtained by subdividing the whole sumple, and can be used in EPA as shown below.
# Uncertainty
Std_e_FR = epa.std_e_r(m[:,:], np.copy(v), 298.15, coords=[1.,1.,1.])
# Std for each direction
Std_e_FR_x = epa.std_e_r(m[:,:], np.copy(v), 298.15, coords=[1.,0.,0.])
Std_e_FR_y = epa.std_e_r(m[:,:], np.copy(v), 298.15, coords=[0.,1.,0.])
Std_e_FR_z = epa.std_e_r(m[:,:], np.copy(v), 298.15, coords=[0.,0.,1.])
print("Std(e_FR): ", Std_e_FR)
print(Std_e_FR_x, Std_e_FR_y, Std_e_FR_z,
np.sqrt((Std_e_FR_x**2+Std_e_FR_y**2+Std_e_FR_z**2))/3)
OUTPUT:
Std(e_FR): 5.700936284757932
9.013063820928355 11.914013072377525 8.326286313652448 5.700936284757933
Another common method for estimating permittivities is by appling an external constant electric field (\(\bold E\)) to the MD simulation, and then using \[ \varepsilon_r = 1+ \frac{\text E[\bold M \cdot \bold u]}{\varepsilon_0 V \bold E \cdot \bold u}\] Here E is the expectation, and \(\bold u\) is an unitary vector in the considered direction. That is, we will get the scalar projection of \(\bold E\) onto \(\bold u\), and evaluate the average value of the scalar projection of \(\bold M\) onto \(\bold u\). In practice \(\bold E\) is along one coordinate axis, e.g. \(z\), and we just compute the time average value of the corresponding coordinate of \(\bold M\): \(\langle M_z \rangle\). This is not only easier but preferable, because the system is not unnecessarily perturbed in other directions.
This scheme will be called the "electric field route" (EFR). Both methods (FR and EFR) are deeply connected. In ref. [1] I proposed a mixed approach, consisting in performing a MD simulation with no electric field applied, estimate for any direction what would be the value of the mean dipole moment if a small electric field were applied, and use it in the above equation.
The code below exemplifies how to perform such calculations using EPA
## Mixed Route (OR)
# Set an small electric field strenght
E = 0.000000002 # V/m
# Estimate the hypotetical dipole moments
m_fake_x = epa.mu_est1(m[0,:], E, 298.15)
m_fake_y = epa.mu_est1(m[1,:], E, 298.15)
m_fake_z = epa.mu_est1(m[2,:], E, 298.15)
e_OR_x = epa.e_s_field(m_fake_x, E, np.copy(v))
e_OR_y = epa.e_s_field(m_fake_y, E, np.copy(v))
e_OR_z = epa.e_s_field(m_fake_z, E, np.copy(v))
print(e_OR_x, e_OR_y, e_OR_z, "Mean:", (e_OR_x+e_OR_y+e_OR_z)/3)
OUTPUT:
63.16907625060437 70.5674145736076 57.9899918700039 Mean: 63.90882756473863
Notice that the results conicide with those of the FR. This is not fortuitous, in reference 1 you can find the proofs that both methods (1-D case for the FR) coincide in their values and uncertainties.
To study the contribution of individual components, it is preferable to use partial susceptibilities instead of permittivities, as the former sum up to the system susceptibility. We will compute those properties using EPA.
As mentioned in [1], it is quite straightforward the computation if the FR is used. For a given component, we just need the values of the row (or column) of the covariance matrix in which the component is always present.
This can be easily carried out using EPA. For, example, for water
### Partial permittivities ###
## Fluctuation Route (FR)
w_numerator_x = np.var(m_w[0,:], ddof=0)+np.cov(m_w[0,:],m_et[0,:], ddof=0)[0][1]
e_w_FR_x = epa.e_s_v2(m_w[:,:], np.copy(v), 298.15, coords=[1.,0.,0.], unbiased=True, M2=w_numerator_x )
w_numerator_y = np.var(m_w[1,:], ddof=0)+np.cov(m_w[1,:],m_et[1,:], ddof=0)[0][1]
e_w_FR_y = epa.e_s_v2(m_w[:,:], np.copy(v), 298.15, coords=[0.,1.,0.], unbiased=True, M2=w_numerator_y )
w_numerator_z = np.var(m_w[2,:], ddof=0)+np.cov(m_w[2,:],m_et[2,:], ddof=0)[0][1]
e_w_FR_z = epa.e_s_v2(m_w[:,:], np.copy(v), 298.15, coords=[0.,0.,1.], unbiased=True, M2=w_numerator_z )
print(e_w_FR_x,e_w_FR_y,e_w_FR_z,(e_w_FR_x+e_w_FR_y+e_w_FR_z)/3.)
OUTPUT:
57.542253140622584 66.24266260624859 52.61626989518161 58.80039521401759The mixed route (and the EFR) make the calculations trivial. This is because of the properties of the expected value. For being accurate, the mixed route requires to shift the dipole moment of the whole system in order to nullify the estimation of the mean of the transformed values. If this were not done, the noise would be too large for the small signal (net dipole moment) due to a small electric field. If we increase too much the arbitrary electric field strenght, the theoreticaly assumed linear relationship between the dipole moment and the field would be lost (see [1] for details).
The code below shows how this calculation can be performed
## Mixed Route (OR)
# Estimate the water dipole moment
m_w_fake_x = epa.mu_est1(m_w[0,:], E, 298.15, m_ref=m[0,:])
m_w_fake_y = epa.mu_est1(m_w[1,:], E, 298.15, m_ref=m[1,:])
m_w_fake_z = epa.mu_est1(m_w[2,:], E, 298.15, m_ref=m[2,:])
# m_ref helps to set that the total in absense of field the
# total dipole momentum must be zero (and not just the water contribution)
e_w_OR_x = epa.e_s_field(m_w_fake_x, E, np.copy(v))
e_w_OR_y = epa.e_s_field(m_w_fake_y, E, np.copy(v))
e_w_OR_z = epa.e_s_field(m_w_fake_z, E, np.copy(v))
print(e_w_OR_x,e_w_OR_y,e_w_OR_z,(e_w_OR_x+e_w_OR_y+e_w_OR_z)/3.)
OUTPUT:
57.54225313629818 66.2426624932671 52.616270988318625 58.80039553929464
Again, the results coincide up to several decimal places with those of the FR.
Let's do the calculations for ethanol
## Mixed Route (OR)
# Estimate the ethanol dipole moment
m_et_fake_x = epa.mu_est1(m_et[0,:], E, 298.15, m_ref=m[0,:])
m_et_fake_y = epa.mu_est1(m_et[1,:], E, 298.15, m_ref=m[1,:])
m_et_fake_z = epa.mu_est1(m_et[2,:], E, 298.15, m_ref=m[2,:])
# m_ref helps to set that the total in absense of field the
# total dipole momentum must be zero (and not just the water contribution)
e_et_OR_x = epa.e_s_field(m_et_fake_x, E, np.copy(v))
e_et_OR_y = epa.e_s_field(m_et_fake_y, E, np.copy(v))
e_et_OR_z = epa.e_s_field(m_et_fake_z, E, np.copy(v))
print(e_et_OR_x,e_et_OR_y,e_et_OR_z,(e_et_OR_x+e_et_OR_y+e_et_OR_z)/3.)
OUTPUT:
6.626823107987436 5.324752111482954 6.373720911022359 6.108432043497583
The values for ethanol are much smaller because its mole fraction is significantly smaller Note that the sum of the values of each component does not result in the value of the whole system: \(58.8+6.1=64.9\neq 63.9\). This does not occur when working with susceptibilities. Remember that \(\chi_e = \varepsilon_r-1\), so we can obtain the susceptibilities by adding 1 to the permittivity values already obtained.
Finally, we can compute the uncertainties in partial susceptibilities (or permittivities). As shown in reference [1], this can be done for each spatial direction using the following equation
\[ \sigma(\chi_e) = \sqrt{ \frac{ \text{Var}[M_k^2] + \text{Var}[M_k M_j]+ 2\text{Cov}[M_k^2,M_k M_j] } { n_{\text{eff},K,H_1}(\varepsilon_0 V k_B T)^2 } } \]Here, \(n_{\text{eff},K,H_1}\) is the effective number of steps, which is associated with the variable \(M_k \sum M_j\).
This can be computed using EPA by following the recipe below
S_m_w_x = m_w[0,:]**2+m_et[0,:]*m_et[0,:]
X_w_std_FR_x = epa.std_er_fake_field_contrib(S_m_w_x, np.copy(v),298.15)
S_m_w_y = m_w[1,:]**2+m_et[1,:]*m_et[1,:]
X_w_std_FR_y = epa.std_er_fake_field_contrib(S_m_w_y, np.copy(v),298.15)
S_m_w_z = m_w[2,:]**2+m_et[2,:]*m_et[2,:]
X_w_std_FR_z = epa.std_er_fake_field_contrib(S_m_w_z, np.copy(v),298.15)
print("Water: ",X_w_std_FR_x,X_w_std_FR_y,X_w_std_FR_z)
S_m_et_x = m_et[0,:]**2+m_et[0,:]*m_et[0,:]
X_et_std_FR_x = epa.std_er_fake_field_contrib(S_m_et_x, np.copy(v),298.15)
S_m_et_y = m_et[1,:]**2+m_et[1,:]*m_et[1,:]
X_et_std_FR_y = epa.std_er_fake_field_contrib(S_m_et_y, np.copy(v),298.15)
S_m_et_z = m_et[2,:]**2+m_et[2,:]*m_et[2,:]
X_et_std_FR_z = epa.std_er_fake_field_contrib(S_m_et_z, np.copy(v),298.15)
print("Ethanol: ", X_et_std_FR_x,X_et_std_FR_y,X_et_std_FR_z)
OUTPUT:
Water: 6.77089441335002 9.397708025694854 8.363982796825423
Ethanol: 0.3788887688883033 0.26753525491101865 0.7649594684297192
The discussion in the section on the mixed route applies here. For this example a similar simulation was employed, but in this case an electric field of 0.01 V/nm was applied in the \(x\) direction. The new data was extracted and stored as before in the variables m, m_w, m_et and v. We will compute the permittivities and they uncertainties, for the whole system, for water and for ethanol. The method used for the computation of the uncertainties is the one proposed in [1].
E = 0.01
# Calculations on the entire system
e_x = epa.e_s_field(-m[0,:].mean(),E,np.copy(v))
std_e_x =epa.std_er_field(-m[0,:], E, np.copy(v))
print(e_x, std_e_x)
# water contribution
e_x_w = epa.e_s_field(-m_w[0,:].mean(),E,np.copy(v))
std_e_x_w =epa.std_er_field(-m_w[0,:], E, np.copy(v))
print(e_x_w, std_e_x_w)
# Ethanol contribution
e_x_et = epa.e_s_field(-m_et[0,:].mean(),E,np.copy(v))
std_e_x_et =epa.std_er_field(-m_et[0,:], E, np.copy(v))
print(e_x_et, std_e_x_et)
OUTPUT:
63.69547865899746 10.537503830726648
60.62047042608366 9.524251343436367
4.075008232913781 1.3955009700218033
Notice the - sign before the dipole moments. This is just because how Gromacs and MDTraj define the spatial directions. Results are in acordance with those of the other methods.The differences are mainly due to the short duration of the simulation and the size of the system, and secondly due to the small lost in linearity between the dipole momentum and the electric field strenght. As shown in [1], it is possible to apply the above methods (FR/mixed route) to the remaining directions.