# 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. – Nike Dattani Aug 25 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. – Charlie Crown Aug 27 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 – Pierpy Sep 3 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 :) – Charlie Crown Sep 3 at 20:42