Re: [AMBER] Calculating electrostatic potential in cpptraj

From: Jason Swails <jason.swails.gmail.com>
Date: Fri, 15 Aug 2014 20:00:46 -0400

On Fri, Aug 15, 2014 at 4:21 PM, Covington, Cody Lance <
cody.l.covington.vanderbilt.edu> wrote:

> Hello all
>
> Is there a way to calculate the electrostatic potential at a given atom
> position or point in space in cpptraj?


​cpptraj cannot (yet) do PME electrostatics. There are a couple actions
that do electrostatic calculations, but they do direct space pair-pair
interactions only using minimum image convention.


> ​​
> Or any other analysis program?
>

​None in Amber that I know about.


> ​​
> I can calculate it using my own programs, but I want to use PME to get
> long range electrostatics that would go
> ​​
> beyond where periodic boundary conditions end. This might seem silly, but
> I have a surfactant system that has a
> ​​
> high ionic strength.
> ​​
>

​​
>
> ​​
> I have also tried to use
> ​​
> image xoffset 0 yoffset 0 zoffset 0
> ​​
> image xoffset 0 yoffset 0 zoffset 1
> ​​
> ... etc
> ​​
> To make the system more like and 'infinite' one, but this is quite
> cumbersome.
>

​It will be _very_ hard to make the system large enough that you can
compute an electrostatic potential using just a direct sum-over-pairs (and
it will get _very_ slow very quickly). Your best bet really is to use an
Ewald-based method. There are a couple options:

1) Write your own PME or Ewald code. The electrostatic potential at any
given point is given by the sum of the direct space ESP and the reciprocal
space ESP. But I believe that the convolution of the charge grid ​with the
electrostatic potential is done in reciprocal space and the result is
back-transformed to real space. So I think you would want to
back-transform the reciprocal space pair potential before convoluting it
with the charge density (if what you're after is the ESP, that is).

2) Use an existing PME or Ewald code to extract what you need. The
electric field is easy enough to get by dividing the electrostatic force by
the charge of that particle. The ESP seems to me to be a little bit more
challenging. Since I don't think you can pull out the reciprocal space
pair potential from any of the existing PME/Ewald codes, you would
effectively need to find the electrostatic energy of putting a unitary
point charge at a particular location. So compute the full electrostatic
energy (and only the electrostatic energy), put a unit charge at the
location you want the ESP at, and compute another electrostatic energy. I
believe the difference is the electrostatic energy of that point charge in
the field generated by all of the particles in all periodic boxes. Divide
by that charge and you get the ESP.

For #2, the easiest way I can think of doing this (and the way I would use
if I had to do it) would be to use the OpenMM capabilities in ParmEd and
write a Python script to compute the ESPs. The general steps here would be:

1) Use createSystem to create an OpenMM system with all of the relevant
forces

2) Create two new, separate System objects with only a single
NonbondedForce object. Populate the new NonbondedForce objects from the
NonbondedForce object in the System you created in step 1, making sure to
set all of the Lennard-Jones parameters to 0 so the only energies computed
come from electrostatics. Add a new particle to the second system (and
second NonbondedForce) that will be the unitary charge at the point in
space you want to compute the ESP.

3) Create Context objects for each of these systems and compute their total
potential energies.

4) Take the difference between the energies computed in (3) and divide by
the charge you assigned to the particle you added.

The advantage of (1) is that you'll almost certainly learn a lot more about
PME and Ewald. The advantage of (2) is that ParmEd provides Python classes
for most of Amber's input and output files (restarts, trajectories, NetCDF
restarts and trajectories, prmtop files, etc.) which should make your job
easier, and you'll add a powerful tool to your repertoire (and you won't
have to program and validate your own PME implementation, of course).

This was quite lengthy, but hopefully it was helpful,
Jason

-- 
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Aug 15 2014 - 17:30:02 PDT
Custom Search