Algorithm for calculating dielectric constant of water TIP3P

Question

I am running a NVT TIP3P simulation of water with 125 molecules of water in a 16-by-16-16 angstrom box with periodic boundary conditions on LAMMPS, with a time-step of 1 fs for 10 ps.

Once the simulation runs, I extract the positions of these particles at each time-step into a .lammpstrj file, then process the data in python to evaluate the dipole moment of my ensemble.

To evaluate dielectric constant $\epsilon$, I will make use of the following relation: $$ \epsilon = 1+\frac{4\pi}{3k_bT} \left( \langle |M|^2\rangle-\langle |M| \rangle ^2 \right)$$

From classical electrodynamics, we know that $$\mathbf{M} = \sum_{i=1}^N q_i \mathbf{r}_i$$

My question is, can I still apply this formula to evaluate dipole moment when I have periodic boundary conditions? Because of periodic boundary conditions, one portion of a molecule might be on one side of the box, while the other portion is on the other side of the box, and this I think leads to artificially large fluctuations, because the molecules are flickering on the edge of the box.

Given the position and charge of each atom in your simulation with periodic boundary conditions, what is the most effective algorithm to evaluate the dielectric constant?

Answer 1

Molecules will not be on both sides of the box at once because this is explicitly prevented by most good MD packages. You can calculate distances which take into consideration the PBC. For example, here is a code to calculate all the pairwise distances with periodic boundary conditions (x_size = [16,16,16])

This is modified from periodic boundary conditions on Wikipedia. I've essentially added a list to store all the pairwise distances, and calculated the distance as

$r = \sum_i \sqrt{(x_i-x_0)^2 + (y_i-y_0)^2 + (z_i-z_0)^2 } $ with the np.linalg.norm function.

r=[]
for i in range(0, N):
     for j in range(0, N):
        dx1 = x[j] - x[i]
        dx = np.mod(dx1, x_size * 0.5)
        r.append(np.linalg.norm(dx))

The np.mod is choosing the distance which is the smallest distance. It is the remainder of dividing the distance by x_size/2. So if the closest molecule is one image away, it is further than x_size/2, from the center of the box. Therefore, dividing by x_size/2 removes this extra amount.

Answer 2

Aren't you missing the volume in the denominator? And the order of the norm and average is probably off in the second term. The original equation should be $$ \epsilon = 1 + \frac{\langle |\mathbf{M}|^2\rangle - |\langle\mathbf{M}\rangle|^2} {3\epsilon_0 V k_B T} $$ and in units where coulombs constant is one, $k = 1 = (4\pi\epsilon_0)^{-1}$, you get $1/\epsilon_0 = 4\pi$ as in your equation. (Also I'm not sure why the term $|\langle\mathbf{M}\rangle|^2$ is there at all. The whole expression seem to be the trace of a dielectric tensor assuming rotational symmetry, and under rotational symmetry $\langle\mathbf{M}\rangle=0$. Maybe I'm misstaken or some else can comment on that. Including it could be a better approximation for finite times.)

If the molecules are broken due to wrapping the coordinates, that could be a problem! You should check that, and "heal" them it if is the case. Alternatively you can probably turn off the coordinate wrapping in Lammps.

Since the molecules are neutral, the total dipole moment is the sum of the molecular dipole moments and therefor independent of wrapping (as long a you don't break the molecules) and can be calculated with your second formula where $i$ runs over all atoms. Since the dipole of an uncharged system is independent of the origin $\mathbf{r}_i$ can just be the position vector. So if R is your Nx3 numpy-array of positions and q the numpy-array of charges, you get your system dipole vector M like:

M = 0
for i in range(len(R)):
  M += q[i] * R[i]  

or why not

M = q @ R

which should come out in units of eÅ assuming charges in units of electrons and distances in units of Ångstrom.