Contents
Methods
Build the system
We calculate the number of molecules from the experimental liquid density. To build the initial configuration, we use gmx insert-molecules to randomly insert 1167 copies of a single molecule (here, dichloromethane), into a 5x5x5 nm simulation box:
gmx insert-molecules -ci dichloromethane.pdb -nmol 1167 -box 5. 5. 5. -o system.gro
The molecule topology files (.itp file), using the OPLS-AA force field, were taken from the database at virtualchemistry.org.
Fig. 2. Simulation box, with 1167 dichloromethane molecules.
Equilibrate and run
MD run parameters are typical for all-atom simulations; we use particle mesh Ewald (PME) to compute electrostatic interactions, and 1.2 nm Lennard-Jones cutoffs. The LINCS algorithm is used to fix the length of all hydrogen bonds, which allows a 2 fs timestep.
To equilibrate the system, we perform energy minimization, followed by short (100 ps) NVT and NPT runs at 300K. The production run is 10 ns, with frames saved every 10 ps (by setting nstxout = 5000), for 1000 total frames in the trajectory file. The simulations run at about 33 ns/day, on 4 cores on a single node.
Fluctuation method
From a zero-field simulation trajectory md.trr, we obtain a time series for the total system polarization M(t) using gmx dipoles:
gmx dipoles -f md.trr -s md.tpr -o dipole.xvg
This produces a file dipole.xvg, that includes the system polarization in x, y, z direction (Mx, My, Mz) versus time. Then we calculate the polarizability \(\chi\) by Eqn. 3, and \(\epsilon\) from \(\chi\) by Eqn. 1.
Linear response method
The linear response method evaluates \(\chi\) directly, by applying a constant electric field and measuring the total polarization \(M = V \langle P \rangle\), and using the definition \(P = \epsilon_0 \chi E\).
We must stay in the linear response regime, such that \(P\) is proportional to \(E\); if we apply too large a field, the polarization will begin to saturate. Roughly speaking, the energy difference \(\Delta U = 2 d E \) between an aligned and antialigned dipole \(d\) in the field \(E\) should be at most kBT, the thermal energy. To explore linear response, we apply electric fields such that the dipole alignment energy is kBT, 0.5 kBT and 0.25 kBT.
To estimate the maximum field to apply, we need a value for \(d\). These can be computed using a quantum chemistry package; we used Gaussian to perform these calculations for our molecules, using DFT with the B3LYP functional and a modest 6-31G basis set. For example, the calculated dipole moment of dichloromethane is 2.23 Debye, and the corresponding maximum field is 0.2784 V/nm.
To apply an electric field in Gromacs, we add a single line to the .mdp file. Here we apply a field of 0.2784 V/nm along x:
E-x = 1 0.2784 0.
We compute the polarization response using gmx dipoles,
gmx dipoles -f field.trr -s field.tpr -o dipole.xvg
followed by gmx analyze to extract the average polarization and its error bar (estimated from the variance and number of independent samples as (σ2/N)1/2 ):
gmx analyze -f dipole.xvg
The polarizability is the slope of polarization density versus electric field E, and the dielectric constant follows from Eqn. 1.
Adjusting partial charges
The partial charges on a molecule can be adjusted by editing its .itp file. An example of a modified .itp file for dichloromethane, in which the default OPLS charges are multiplied by 1.43, looks like this:
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 opls_151 1 DCM CL1 1 -0.143 35.453
2 opls_152 1 DCM C 2 -0.0086 12.011
3 opls_153 1 DCM H1 3 0.1473 1.008
4 opls_153 1 DCM H2 4 0.1473 1.008
5 opls_151 1 DCM CL2 5 -0.143 35.453