On Mon, May 25, 2015 at 12:44 PM, Lee-Ping Wang <leeping.stanford.edu>
wrote:
> Hi everyone,
>
> I’ve been running some replica exchange simulations in Amber. I’ve got a
> simple script for cpptraj (adapted from the tutorial at
> http://ambermd.org/tutorials/advanced/tutorial7/ <
> http://ambermd.org/tutorials/advanced/tutorial7/> ) for processing the
> output trajectories, resulting in “temperature-based” trajectories where
> the temperature is constant but the conformations jump around. The
> simulations are running very nicely; thanks for writing a great software
> package.
>
> My question is how I could get potential energies to “go with” these
> trajectories? I would like to use these potential energies as input to a
> WHAM / MBAR calculation to improve the estimation of thermodynamic
> properties at each temperature.
>
The standard way is to feed those trajectories into sander, set imin=5 to
process an input trajectory, and get energies that way. But I'll provide
an alternative I think you may find more amenabl
e:
Use the sander API to compute single-point energies (and assuming you're
going to turn around and feed this directly into pymbar, you can do
everything you want in a single script). You can do something like this:
import sander
import chemistry as chem
parm = chem.amber.AmberParm('prmtop', 'inpcrd')
traj = chem.load_file('your_trajectory') # mdtraj or pytraj also works to
read trajectories
input = sander.pme_input() # or use "gas_input(NUM)" where NUM is the igb
model for implicit solvent
input.cut = 9.0 # struct members have same name as mdin variables
with sander.setup('prmtop', parm.coords, parm.box, input):
for i in range(traj.frame):
sander.set_positions(traj.coordinates(i))
ene, frc = sander.energy_forces()
# Go ahead and store your potential energy in a numpy array you set
up
You can also use OpenMM instead of the sander Python API if you're more
comfortable with that.
> My simple script is attached below.
>
> Thanks,
>
> - Lee-Ping
>
> parm ../prmtop
> ensemble ../remd.00/netcdf.00
> ensemble ../remd.01/netcdf.00
> ensemble ../remd.02/netcdf.00
> strip :WAT
>
Obviously you can't calculate the full energies if you strip waters, but I
figured I'd point this out anyway.
All the best,
Jason
--
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon May 25 2015 - 12:30:03 PDT
