A Mismatching in Pair Potential Energy Measurement

Dear LAMMPS Developers and Users,

I have a 3d system of 500 charged particles interacting via Yukawa Pair potential. I am using the NVT thermostat (using ‘fix nvt’ command) to equilibrate the system at desire temperature. I am calculating the total potential energy and pair energy (both are equal for this case, as expected) using the following commands:

pair_style yukawa {kappa} {cutf}
pair_coeff * * 3.2891e-20

thermo_style custom step temp pe epair
thermo 500

Beside that, I have also dumped all the positions and velocities by following command:

dump dp1 all custom 500 dump_nvt.lammpstrj id type x y z vx vy vz

By knowing the positions of all the particle at a given time step, I can calculate total pair energy of the system at that time step by calculating

sum_(i,j) [(A./r_{ij}).*exp(-{kappa}.*r_{ij})] ,

With the condition " r_c(cutf) >= r_{ij} > 0 "

where r_i = sqrt((x_i).^2 + (y_i).^2 + (z_i).^2) and r_{ij} = abs(r_j - r_i)

I expected potential energies from these two different methods should be same or might be different by a constant normalization factor, which will be indipendent of ‘kappa’. But I am not getting this. Am I missing any conceptual issue by doing the calculation of pair potential energy?
Any help or suggestion will be appreciable.

I am using ‘lammps-17Nov16’ version of LAMMPS.

Thank you in advance.

With regards
Srimanta Maity
Institute for Plasma Research, Gandhinagar, India

Dear LAMMPS Developers and Users,

I have a 3d system of 500 charged particles interacting via Yukawa Pair
potential. I am using the NVT thermostat (using 'fix nvt' command) to
equilibrate the system at desire temperature. I am calculating the total
potential energy and pair energy (both are equal for this case, as expected)
using the following commands:

pair_style yukawa \{kappa\} {cutf}
pair_coeff * * 3.2891e-20

thermo_style custom step temp pe epair
thermo 500

Beside that, I have also dumped all the positions and velocities by
following command:

dump dp1 all custom 500 dump_nvt.lammpstrj id type x y z vx
vy vz

By knowing the positions of all the particle at a given time step, I can
calculate total pair energy of the system at that time step by calculating

sum_(i,j) [(A./r_{ij}).*exp(-{kappa}.*r_{ij})] ,

With the condition " r_c(cutf) >= r_{ij} > 0 "

where r_i = sqrt((x_i).^2 + (y_i).^2 + (z_i).^2) and r_{ij} = abs(r_j -
r_i)

1) Is this a periodic system? Did you remembering to consider
periodic images in the summation above?

I expected potential energies from these two different methods should be
same or might be different by a constant normalization factor, which will be
indipendent of 'kappa'. But I am not getting this. Am I missing any
conceptual issue by doing the calculation of pair potential energy?

2) One suggestion would be to compare differences between the energies
of two different conformations, rather than the absolute energy. You
can add a constant to the energy without effecting the physics of the
simulation.

3) Are you using pair_modify shift? (this will add a constant to the
energy of the system)

4) You can ask LAMMPS to generate a table of energy (and force) as a
function of inter-particle spacing (r) using the "pair_write" command.

5) This does not effect you, but (generally speaking) if you are using
a pair_style with a long range electrostatic contribution, then you
have to consider the contributions from periodic images carefully.

That's what comes to my mind.
I don't know if this helped at all.

Cheers

Andrew

Thank you for your suggestions. Yes, I have missed to point out that I am using periodic boundary condition in my simulation. Here I am using Yukawa pair potential, which is a short-range interaction limited up to the distance r = r_c (force cutf). So in this particular case, I wonder is it necessary to consider periodic images in the summation?

No, I am not using “pair_modify shift” command. But I have used the “pair_write” command to get the profile ® of pair energy which is exactly matching with the analytical expression [(A./r).*exp(-{kappa}.*r)]. But it does not give me the total pair potential energy of the entire system.

Thank you
With regards
Srimanta

Thank you for your suggestions. Yes, I have missed to point out that I am using periodic boundary condition in my simulation. Here I am using Yukawa pair potential, which is a short-range interaction limited up to the distance r = r_c (force cutf). So in this particular case, I wonder is it necessary to consider periodic images in the summation?

yes, ​absolutely​. please grab your MD text book and look it up.

axel.