# How to get electrostatic potential at same point in pppm calculations?

Hello, I want to get electrostatic potential at same point in pppm
calculations of solution but have no idea how to do it even after
search the mailing list a lot. Could Lammps dump such quantities by a
simple command or should I calculate it manually? As far as I know,
the quantity is the sum of the interactions between all the charged
particles with the chosen point, am I right? If so, I think I’d like
to do it myself than other inconvenient tools.

Thank you very much.

Shijun Zhao

LAMMPS doesn't compute such a quantity since it isn't needed
for dynamics. Perhaps Stan can comment on how hard it would
be to add it as a diagnostic, i.e. a compute.

Steve

I worked on this a few years back, and will have to recommence in the
near future. The code has changed quite a bit since back then (2008),
but the general approach should still be valid:

Look at pppm.cpp in void PPPM::poisson_ik()
The electrical field in all spatial directions is calculated using a
fourier backtransform of -ik*V(k) (this corresponds to a derivative of
the potential in the real space direction corresponding to k)

You should be able to obtain the potential by simply performing a
fourier backtransform of V(k) at this point! Be careful with the units
and remember that you may have to add the short range contributions
for the potential.

Daniel

Alternately, why not add a test charge (or dipole) of unit magnitude
and measure how the potential energy changes in response?
You can do this using the "rerun" command.

Details below:

---- For the electrostatic potential ----
Suppose you have a system with "n" atoms.
Suppose you wish to calculate the eletrostatic potential at N
locations in space during your simulation.

Create a LAMMPS DATA file which has N+n atoms.
N of these "atoms" could be located at the positions you want to
calculate the electrostatic potential.
These N particles would be test/ghost particles which do not interact
with the system during the simulation. (You could give them zero
charge and zero Lennard-Jones epsilon parameters. To immobilize these
particles, you could either set their initial velocities to zero, or
simply exclude them from the group of atoms which are passed to the
"fix nvt" command, for example.)

After running the simulation, you can turn on the charge of one of the
ghost particles. (Assign it's charge to 1.0 in the "Atoms" section of
the DATA file.)
Then recalculate the potential energy of the entire system using the
"rerun" command.

http://lammps.sandia.gov/doc/rerun.html

The difference between the energies (at every frame in the trajectory)
is the electrostatic potential. You could repeat this process for
each of the N ghost particles you added to the simulation (each time
setting one of the charges to 1.0).

---- For the electric-field ----

You could do the equivalent thing using a "test dipole":

0) place an electric point dipole at the location(s) where you want to
query the electric field, but set their dipole moments to 0.
Run the simulation.
1) align the unit dipole along the x axis, calculate the energy using "rerun"
2) align the unit dipole along the y axis, calculate the energy using "rerun"
3) align the unit dipole along the z axis, calculate the energy using "rerun"

Subtract 0) from 1) and 2), and 3).
The three energy differences should be proportional to the x,y,z
components of the electric field at that location. If you have N
"test dipoles" in your simulation, you can repeat steps 1),2),3) for
all N of these particles. Again, when you run the simulation during
step 0, you would turn off dipole moments so they don't interact with
the other particles in the simulation (and you would immobilize them).

This way, users can even move the locations in space where the
electric field is queried by allowing the ghost particles to move, or
by writing a script to modify the corresponding lines of the "dump"
file.

Cheers
Andrew

Thank all your enthusiasm!
Now I realize that it is not a trivial task if I want to draw an
electrostatic potential distribution map. I think I will spend much
more time on it if I following your kindly advices. It is really a big
challenge!

2012/12/4, Andrew Jewett <[email protected]>:

Shouldn't be too hard, but you would have to write your own compute. You just need to interpolate from the PPPM mesh to the point you care about. For example, this is already done at the location of the atoms to get the per-atom energy.

Stan