# Algorithm for calculating dielectric constant of water TIP3P

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?

• +1. Great question. I just put the actual questions in bold for people looking for what specifically you want. Aug 25, 2020 at 0:24
• You need to be simulating a much larger system and for a much longer time. You need to be simulating for greater than 10 nanoseconds. 10 picoseconds is incredibly short. You will not have even equilibrated PVT properties. Dielectric takes a MUCH longer time to equilibrate. Aug 27, 2020 at 18:12
• I completely agree with @CharlieCrown. The size of the system may affect the final result by up to 20% moreover I once have played with a 216 water molecules box and the found convergence time was on the order of 8 nano seconds. You can find an interesting article here sciencedirect.com/science/article/pii/S0009261411002740 Sep 3, 2020 at 12:44

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.

• Also known as mirror image separation :) Sep 3, 2020 at 20:42

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 D like:

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


Or why not:

D = q @ R