Re: [AMBER] Plotting heat maps in AMBER aMD simualtions

From: Jason Swails <jason.swails.gmail.com>
Date: Tue, 16 Apr 2013 11:10:05 -0400

On Tue, Apr 16, 2013 at 3:19 AM, Neha Gandhi <n.gandhiau.gmail.com> wrote:

> Dear Amber Developers and Users,
>
> How do I extract free energies from aMD simulations? Is there a way to plot
> free energy contour map/heat map (phi/psi and coloring according to free
> energies)?
>

You need to re-weight each point in order to get its corresponding
'contribution' to the unbiased canonical ensemble. The prescription for
doing that is in the original Accelerated MD paper: J. Chem. Phys. 120,
11919 (2004); http://dx.doi.org/10.1063/1.1755656. So what I would do is
use cpptraj to calculate whatever properties you want from the trajectories
(phi/psi maps, for instance), then histogram the data, taking each point
and adding its weight into the bin in which it resides. You can calculate
the weight at whatever step you want based on the information written to
the AMD log.

The reweighting formula may be in other papers as well, but it's written
pretty clearly in the Hamelberg paper I showed. Admittedly the existing
documentation is not super clear about what terms mean what in the AMD log
file, (nor do I see a tutorial anywhere describing it), but I believe the
format is:

<traj writing frequency>
<current step #>
<unbiased total potential energy>
<total unbiased dihedral energy>
<force weighting factor>
<force weighting factor from dihedrals>
<total boost energy>
<dihedral boost energy>

Hopefully Romelia will correct me here if I'm wrong. In any case, once you
have the weighted value of each point along the time series of your
dihedrals, you can create a histogram using these weights---numpy can do
this easily provided a 'weight' array to numpy.histogram2d. Then simply
generate a free energy-ified heat map from that data (i.e., -kT ln(bin /
max) -- where max is the most-populated bin), to give a relative free
energy surface with '0' as the most stable state).

Ben Roberts wrote a nice little Python script that can generate heat maps (
http://www.merzgroup.org/heatmap.html) -- you can modify it to also read in
a set of weights and give that to the numpy.histogram2d call, and analyze
the source to see how he does it.

HTH,
Jason

-- 
Jason M. Swails
Quantum Theory Project,
University of Florida
Ph.D. Candidate
352-392-4032
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Apr 16 2013 - 08:30:03 PDT
Custom Search