libAtoms / matscipy

Materials science with Python at the atomic-scale
http://libatoms.github.io/matscipy/
GNU Lesser General Public License v2.1
188 stars 55 forks source link

Pair potential for the real part of the Ewald sum #245

Closed cagrikymk closed 5 months ago

cagrikymk commented 5 months ago

Hello, I am trying to understand the following code in the context of the real part of the ewald sum: from: pair_potential /calculator.py Line-396

  # Forces
  df_pc = -0.5 * de_p[_c] * r_pc / r_p[_c]

  f_nc = mabincount(j_p, df_pc, nb_atoms) - mabincount(
      i_p, df_pc, nb_atoms
  )

I do not fully understand how f_nc is computed. I expect it to produce the same results as the following code but they are different:

f_nc = np.zeros((N,3))
df_pc = -1 * de_p[_c] * r_pc / r_p[_c]
for i, j, f in zip(i_p, j_p, df_pc):
    f_nc[i] += f

For a given interaction (i,j), the neighbor list contains both (i,j) and (j,i) entries. However, since I am only using the i_p values, this should already take care of the double counting.

However, they do not produce the same results which I don't fully understand why. It would be great if someone can explain what is happening here.

pastewka commented 5 months ago

You need to add

f_nc[j] -= f

to your code for those to be equivalent. This is Newton's third law.

pastewka commented 5 months ago

Ahh, I see what you are saying. This is related to periodic boundaries.