Electrostatic energy calculation

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

Does anybody know if there is any way to calculate pure electrostatic between two groups of atoms? I know there is a compute group/group command that can calculate total pair energy, but I only need the electrostatic part. Moreover I have used “pair_style lj/cut/coul/cut” in my input file.

Thanks in advance.

2

Hello,

Does anybody know if there is any way to calculate pure electrostatic
between two groups of atoms? I know there is a compute group/group command
that can calculate total pair energy, but I only need the electrostatic
part. Moreover I have used "pair_style lj/cut/coul/cut" in my input file.

there is no provision for that.

you have two options:
- do your simulation and then use the rerun command to post-process
the trajectory where you replace lj/cut/coul/cut with coul/cut and
then you can use compute group/group
- modify the source code and build a custom executable. e.g. you could
modify the single() method for lj/cut/coul/cut to only return the
coulomb term. since that is independent from the regular force compute
in compute(), it will (well, it should) not affect the trajectory.

axel.

3

If you only need it once in a while you can exit the run command, switch the pair interaction to LJ/cut/coul with an epsilon of zero (pure Coulomb) and then set some charges to zero.
Disable the time integration fix

  1. Everything but group 1 to zero
  2. Everything but group 2 to zero
  3. Everything but group 1 and 2 to zero
    And do a run 0 after each step and put the potential energy into a variable
    E3-E2-E1 should be what you want.
    Then just restore all charges, the pair interaction and the time integration and continue to run.
4

If you only need it once in a while you can exit the run command, switch the
pair interaction to LJ/cut/coul with an epsilon of zero (pure Coulomb) and
then set some charges to zero.
Disable the time integration fix
1. Everything but group 1 to zero
2. Everything but group 2 to zero
3. Everything but group 1 and 2 to zero
And do a run 0 after each step and put the potential energy into a variable
E3-E2-E1 should be what you want.
Then just restore all charges, the pair interaction and the time integration
and continue to run.

To clarify...

You would run a simulation first, using the correct charges. Then (as
Axel suggested) I would analyse the resulting trajectory (dump) file
using the "rerun" command:

As Daniel suggested I would analyse the trajectory 3 times
1) setting charge of atoms to zero which do not belong to group1
2) setting charge of atoms to zero which do not belong to group1
3) setting the charge to zero of atoms which do not belong to either
group 1 or group 2

At each dump interval, I would print out the "pair" contribution to
the energy at each time (not the total energy).
(Call the pair energy: E1, E2, E3)

Then, as Daniel suggested, the Coulombic interaction energy should be:

Einteraction = E3 - (E1+E2)

I mostly just wanted to add that using "rerun" is (probably) more
convenient and faster than using ("run 0" and "read_dump" or "read
data" and "set" the coordinates.)

Cheers
Andrew

5

If you only need it once in a while you can exit the run command, switch the
pair interaction to LJ/cut/coul with an epsilon of zero (pure Coulomb) and

At each dump interval, I would print out the "pair" contribution to
the energy at each time (not the total energy).

Those should be exactly the same if you follow my recipe.
Daniel