[lammps-users] ewald.cpp: compute()

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps-users,

I have a quick question regarding some of the code in ewald.cpp: compute(). When computing the Coulomb part of the potential energy, the code is:

     for (k = 0; k < kcount; k++) {
       energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
       sfacim_all[k]*sfacim_all[k]);
     }
     PI = 4.0*atan(1.0);
     energy -= g_ewald*qsqsum/1.772453851 +
       0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);
     energy *= qqrd2e;

- the energy term in the for loop is the long-ranged contribution;

- the second term ( -g_ewald*qsqsum/1.772453851) is the self interaction energy;

- I am having trouble figuring out what the term

          + 0.5*PI*qsum*qsum / (g_ewald*g_ewald*volume);

corresponds to. I thought it might be the short-ranged part, but that is calculated in the pair potential code, in particular, pair_lj_cut_coul_long.cpp: compute().

Any help with this would be much appreciated!

Thanks,
Joyce

2

Hi Joyce,
That term is the contribution from the neutralizing background. It kicks in when the overall system is non-neutral. In this situation, in order for the Ewald summation method to work, one paints the simulation box with a uniform charge density, whose integral is equal to minus the total charge on the system.

best,
bala
http://www.jncasr.ac.in/bala