VMD-L Mailing List

From: saam (saam_at_ks.uiuc.edu)
Date: Mon May 04 2009 - 14:31:31 CDT

Dear Richard,

I think what you are seeing is a boundary effect.
I assume you simulated the water box with periodic boundary conditions
in the NPT ensemble. Then you computed the free energy map over the
whole box. However, since you were not telling "volmap ligand" to take
periodic boundaries into account you will get wrong results for all
gridpoints that are within 12 A of the system boundaries because atoms
in that region won't see their interaction partner that lie in
neighboring periodic cells.
The easiest way for you to fix this would be to average the occupancy
only over the central 16x16x16 region of your box where the map is
correct. I don't know if the 40x40x40 box Jordi reported in his paper
refers to the actual simulation box or the smaller box he averaged
over (at that time there was no option to use periodic boundaries
yet). I suppose his simulation was bigger and he averaged over
40x40x40.

In VMD 1.8.6 there is are undocumented options for aligning the
molecule and for taking PBC into account.
I listed the documentation below.

However, we have a much improved version of ILS in our CVS tree that
is able to use GPUs giving you a speedup of up to 20. Also you can now
use probes with more than two atoms. If you are interested in getting
the new version let me know and I'll update the documentation for
this bleeding edge version and send it to you.

I hope this helps,
Jan

PS. there is one other thing in your ILS setup that makes me
suspicious: You use "segname WAT" instead of "all" as a selection.
What are the atoms you omit? Really you should include all atoms and
use the atomselections only to select spatial regions. this is why in
the new version you provide just the minmax of the box defining the
grid and the selection is assumed to be all atoms in that box and all
atoms within cutoff distance of the box.

Additional options for volmap ligand:

-alignsel <sel>
Use the provided selection to align all trajectory frames to the first frame
Note that the results are slightly different from what you get when you align
all frames manually prior to the computation. The differences are
quite noticeable when comparing the files as ascii texts but actually
they are of numerical nature and when you display the maps they
actualy look the same.

-transform <matrix>
Suppose you want to align your trajectory to a reference frame from a
different molecule then you can specify the alignment matrix (as
returned by "measure fit") that would align the first frame to the
reference.

-pbccenter <vector>|<list of vectors>
This option signals that you want a PBC aware PMF calculation. Depending
on the desired target grid size image atoms from neighboring PBC cells
are taken into account for the computation. The atoms used for the
calculation are chosen from a box that exceeds the target grid size by
the interaction cutoff in each direction. I.e. there's no need to
choose a smaller grid.
Since the PBC cell origin is not stored by VMD you have to specify it.
In many cases it's just the same value for all frames, e.g. {0 0 0},
then you can just specify a single coordinate vector. If your origin
varies over the frames you must specify a list with a vector for each
frame. (It would be nice if we could add an origin field to the timestep
data structure so that one doesn't have to fiddle with very long TCL
lists. Then I could write a simple TCL script that reads the origins
from an xst file and fills them in timestep.)

-pbccutoff <vector>
If this option is omitted or {0 0 0} then the target grid will be the
rectangular box that encloses the PBC cell unless the target grid size
was explicitely given with -minmax. You can shrink or increase the
target grid in each dimension by specifying negative or positive values
respectively. For example {-5 0 8} shrinks the box by 5 A in each
direction of the x-dimension and increases it by 8 A in z.

On Mon, May 4, 2009 at 12:35 PM, Klaus Schulten <kschulte_at_ks.uiuc.edu> wrote:
> Dear Jan:
> Can you look into this please?
> Thanks, Klaus
>
>
> Begin forwarded message:
>
> From: Richard Daigle <richard.daigle.3_at_ulaval.ca>
> Date: May 4, 2009 12:21:41 PM CDT
> To: kschulte_at_ks.uiuc.edu
> Cc: vmd-l_at_ks.uiuc.edu
> Subject: Solvation energies using implicit ligand sampling
> Dear all,
>
> I am writing you because I am unable to reproduce NO/O2 solvation energies
> by
> implicit ligand sampling.
>
> I am getting 1.34 kcal/mol and 1.72 kcal/mol for NO and O2, respectively.
> In,
> the orignal paper, solvations energies of 1.60 kcal/mol and 1.97 kcal/mol
> were reported for NO and O2.
>
> For this calculation, I believe that I did correctly the same procedure as
> published in the original paper of Cohen and al. (Biophys. J. 2006, 91,
> 1844-1857). Briefly, I performed implicit ligand sampling (NO and O2) of a 5
> ns trajectory of a water box (40Å x 40Å x 40Å). PMF of each grid point was
> converted to its associated occupancy probability ( e^(-PMF) ) using
> volutil.
> Then, occupation probabilities were averaged over all grid points and the
> solvation free energy was calculated according to the formula :
> deltaG = -RT ln(computed average)
>
> Here is the volmap command I used :
>
> volmap ligand [atomselect $molid "segname WAT"] -allframes -combine
> pmf -minmax {{ -8 -8 -8 } { 8 8 8 }} -dx "./box_$job.dx" -dihetero -conf
> 50 -probe1 -0.20 1.85 -probe2 -0.12 1.70 -bond 1.15 -cutoff 12 -subres 2
>
> Help will be very appreciated !
>
> Thank you !
>
> Richard Daigle
> Étudiant au doctorat en biochimie
> Laboratoire du professeur Patrick Lagüe
> Université Laval, PROTEO
> Québec, Canada
>
> 

-- 
****************************************
Jan Saam
Theoretical and Computational Biophysics Group 3061 Beckman Institute
University of Illinois
405 N. Mathews Ave
Urbana, IL 61801
Phone: (217) 244-1928
Fax:   (217) 244-6078
saam_at_ks.uiuc.edu