Hi Callum and Dan,
Thank you both for getting back to me!
I have gone back and tried using the "diffusion" routine. For this, I use
the exact same cpptraj input as above, except instead of the
'stfcdiffusion' line, I use 'diffusion mask .1-16740 out diffusion_top.dat
time 100 noimage'. The plotted MSD in xy is calculated as: sqrt(X^2 + Y^2),
which I believe should be the correct way to determine it. As you can see,
when using the top leaflet both diffusion and stfcdiffusion result in
similar, but slightly different plots (see the attached figure
LPS_comparison.pdf); the diffusion coefficient I calculate from them is 1.3
x10^(-10) cm^2/s and 1.9 x10^(-10) cm^2/s, respectively. However, diffusion
and stfcdiffusion produce very different plots for the bottom leaflet. The
diffusion routine gives the general shape that I would expect, and very
similar to what I see when using the top leaflet. I'm not sure why
stfcdiffusion gives something wildly different for the bottom leaflet. This
diffusion/stfcdiffusion difference is consistent for both of my bilayer
systems (mLPS_comparison.pdf, attached).
As far as my thought process behind centering the trajectory after I unwrap
it- I found that this step is necessary for two reasons. First, the two
leaflets of the bilayer have a tendency to "caterpillar" across each other,
meaning they have COM motion relative to each other. I believe this is part
of the problem I'm having when I use 'stfcdiffusion' on the full bilayer,
which is why I started looking at the individual leaflets instead. The
second reason I center the trajectory is that the simulation box itself
seems to move around with respect to the origin. What I mean is that, if I
visualize the unwrapped trajectory in VMD, the bilayer patch floats around
seemingly at random in x, y, and z. I agree that it seems like I should be
able to calculate diffusion from the unwrapped system directly, but in
practice this gives MSD plots that don't really make sense (see the
attached *_nocenter.pdf for example- these were made using the exact same
commands as before, just with commenting out the 'center' step). This
behavior occurs in both of my bilayer systems, though I'm not entirely sure
why the box moves around like this.
Thank you,
Amy
On Tue, Jul 19, 2016 at 10:26 AM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> Hi,
>
> On Thu, Jul 7, 2016 at 11:49 AM, Amy Rice <arice3.hawk.iit.edu> wrote:
> >
> > This is the procedure I've been using: first, I unwrap my trajectory,
> then
> > center one of the leaflets at the origin using "center .1-16740 mass
> > origin" (here, atoms 1-16740 correspond to the top leaflet). After this,
> I
> > calculate the MSD using "stfcdiffusion mask .1-16740 out
> diffusion_top.dat
> > time 100 xy". When I plot the MSD in xy with time, it has the general
> shape
> > I would expect (see the attached figure LPS_top.pdf), and the diffusion
> > coefficient that I calculate from it matches the experimental value.
>
> I guess I'm not clear on why you are performing 'center' after you
> 'unwrap'. To me, the point of a diffusion calculation is to determine
> the rate of free motion of molecules in your system. By adding the
> 'center' command, you're artificially changing that motion. Why
> wouldn't you just calculate diffusion for the unwrapped system itself?
>
> -Dan
>
> > However, when I go back and do the exact same procedure, using now the
> > bottom leaflet as my atom mask, the plot does not behave as I would
> expect.
> > First, the xy MSD starts at ~ 1100 instead of 0, and then decreases with
> > time (see the attached figure LPS_bottom.pdf). Finally, if I use the full
> > bilayer as my atom mask, the plot starts at 0 but has large fluctuations
> > where it increases rapidly then decreases rapidly (figure LPS_all.pdf).
> I'm
> > at a bit of a loss to explain this behavior, because if I visualize the
> > unwrapped, centered trajectories, the appropriate portion of the bilayer
> is
> > centered and no strange imaging issues seem to be occurring.
> >
> > As a test, I repeated these same calculations using a bilayer simulation
> of
> > a similar lipid, and have the same type of results where the top leaflet
> > results in a nice plot and a reasonable diffusion coefficient, whereas
> > using the bottom leaflet gives a plot that starts large and decreases,
> and
> > using both leaflets has large fluctuations and is non-linear (figures
> > mLPS_top.pdf, mLPS_bot.pdf, and mLPS_all.pdf).
> >
> > Is this behavior expected? I tried following the link given in the Amber
> 16
> > manual to learn more about the diffusion routine developed by Hannes
> > Loeffler that stfcdiffusion is based off but was unable to locate any
> > pertinent documentation or publications. The data from the top leaflet
> > seems usable, but I would like to be able to average over both leaflets
> if
> > possible so that I am including more lipids and can get a better estimate
> > of the diffusion coefficients.
> >
> > Thank you,
> > - Amy
> >
> > --
> > Amy Rice
> > Ph.D. Student
> > Physics Department
> > Illinois Institute of Technology
> >
> > _______________________________________________
> > AMBER mailing list
> > AMBER.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber
> >
>
>
>
> --
> -------------------------
> Daniel R. Roe, PhD
> Department of Medicinal Chemistry
> University of Utah
> 30 South 2000 East, Room 307
> Salt Lake City, UT 84112-5820
> http://home.chpc.utah.edu/~cheatham/
> (801) 587-9652
> (801) 585-6208 (Fax)
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
--
Amy Rice
Ph.D. Student
Physics Department
Illinois Institute of Technology
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jul 20 2016 - 12:00:02 PDT
