Hi Callum,
Thank you so much for the response, your example is incredibly helpful! I
will try redoing my calculations this way to see if I can get a smoother
MSD vs time curve. I do have one quick questions though - I'm not sure I
understand step 5 of the example procedure you provided, and I have read
some bilayer papers where similar language is used to describe lipid
diffusion calculations, so I hope you can help clarify for me. When you say
using different time origins and 20ns blocks, do you mean calculating MSD
from 0 to 20,000 ps, 200 to 20,200 ps, 400 to 20,400 ps, ..., 80,000 to
10,0000 ps, and then averaging over those to get an averaged MSD vs time
plot that goes from 0 to 20 ns?
Thank you again for the help,
Amy
On Wed, Jul 20, 2016 at 4:25 PM, Dickson, Callum J <
callum.dickson09.imperial.ac.uk> wrote:
> Hi Amy,
>
>
> Yes, the MSD in XY should be sqrt(X^2 + Y^2), I cannot tell why you are
> getting different results between the two codes, you could check the source
> code to see that they do the same thing.
>
>
> If in doubt, you can always use cpptraj to find the position of the COM of
> your atoms at each frame and check the result manually.
>
>
> So for a single lipid (eg. a PC lipid composed of residues 1 to 3), you
> can find the position of the COM as:
>
>
> trajin bilayer_run.nc
> vector lip1_com :1-3 :1-3 out lip1_com.out
>
> go
>
>
> The results of diffusion and stfcdiffusion (for :1-3) should track with
> movement of the COM of this lipid.
>
>
> The 'caterpillar' effect you noticed is a known issue with lipid diffusion
> calculations, see:
>
>
> http://pubs.acs.org/doi/abs/10.1021/jp0348981
>
>
> You should first create a trajectory without this COM drift before finding
> the lipid diffusion.
>
>
>
> Here is the procedure I used to get lipid diffusion:
>
>
> 1) Find the movement of the upper monolayer COM:
>
>
> trajin bilayer_run.nc
>
> vector com center out upper_movement.dat :1-192
>
> go
>
>
> [Note: this bilayer run did not use iwrap; furthermore there is no imaging
> at any point]
>
>
> 2) Create a new prmtop and trajectory of the upper monolayer only:
>
>
> parm POPC_128.prmtop
> trajin bilayer_run.nc 1 1 1
> strip !(:1-192) outprefix upper
> go
>
>
> 3) From this trajectory, for each frame, reset the monolayer by the
> corresponding movement (in X, Y, Z) to create a new trajectory without
> monolayer COM drift.
>
>
> 4) Now using the new trajectory, get the MSD of each lipid:
>
>
> trajin output_upper.traj
>
> diffusion mask :1-3 out upper_xy_1.dat com time 10 xy
>
> go
>
>
> ...and so forth for all lipids.
>
>
> 5) Take the average over all lipids in upper monolayer and using different
> time origins (I used 20 ns blocks of a 100 ns run with many time origins
> separated by 200 ps).
>
>
> 6) Repeat for the lower monolayer.
>
>
> This gave a smooth MSD vs time curve from which to calculate lipid
> diffusion.
>
>
> Best,
>
> Callum
>
> ________________________________
> From: Amy Rice <arice3.hawk.iit.edu>
> Sent: Wednesday, July 20, 2016 2:57:32 PM
> To: AMBER Mailing List
> Subject: Re: [AMBER] Calculating lipid diffusion using stfcdiffusion
>
> Sorry to reply again so quickly, but there is something else I just noticed
> about the output from stfcdiffusion. As I mentioned above, I believe MSD in
> xy should be calculated as sqrt(X^2 + Y^2) (please correct me if I am
> wrong!!). However, the xy MSD in the output from stfcdiffusion is not the
> same as what I calculate it should be, based on the x and y MSDs written to
> the output file. Attached, I've plotted MSD using the xy result and the MSD
> using sqrt(x*x+y*y). The two clearly don't match. However, the MSD using
> sqrt(x*x+y*y) with the stfcdiffusion output matches perfectly the output
> from 'diffusion'. In other words, diffusion and stfcdiffusion agree
> perfectly when they calculate the x and y MSDs, it is just that the
> stfcdiffusion xy MSD does not match sqrt(x^2+y^2). I hope I am explaining
> this in a way that makes sense!
> Below, I have included some portions of the output from stfcdiffusion,
> using the top leaflet of my LPS bilayer system. Am I misunderstanding how
> stfcdiffusion calculates the xy MSD or is there something else going on
> here?
>
>
> ####################################
> #time x y z xy
> 100.000 1.463 1.361 0.884 2.824
> 200.000 1.743 1.801 1.206 3.545
> 300.000 2.125 1.975 1.331 4.100
> 400.000 2.172 2.290 1.405 4.462
> 500.000 2.540 2.469 1.430 5.009
> 600.000 2.529 2.621 1.496 5.150
> 700.000 2.775 2.744 1.608 5.519
> [...]
> 759100.000 14.611 15.105 6.255 29.716
> 759200.000 14.636 15.450 6.124 30.086
> 759300.000 14.774 15.283 5.919 30.057
> 759400.000 14.486 15.232 6.009 29.718
> 759500.000 14.634 15.158 6.145 29.793
> 759600.000 14.554 15.131 5.664 29.686
> 759700.000 15.075 15.811 6.089 30.886
> 759800.000 14.838 15.277 5.548 30.115
> 759900.000 14.750 15.661 5.751 30.411
> 760000.000 14.519 15.309 5.383 29.828
> 760100.000 14.959 15.787 5.795 30.745
> 760200.000 14.771 15.825 6.015 30.596
> 760300.000 15.300 15.818 5.735 31.118
> 760400.000 14.885 16.014 6.068 30.899
> 760500.000 15.090 15.375 5.983 30.465
> [...]
> 1497900.000 16.380 16.944 7.148 33.324
> 1498000.000 15.848 16.386 6.894 32.234
> 1498100.000 16.048 16.246 7.261 32.295
> 1498200.000 16.097 16.529 7.385 32.626
> 1498300.000 16.226 16.498 7.298 32.724
> 1498400.000 16.553 17.061 7.175 33.614
> 1498500.000 16.119 16.694 7.823 32.813
> 1498600.000 16.538 16.925 7.840 33.464
> 1498700.000 16.112 16.480 7.824 32.591
> 1498800.000 16.378 15.839 7.585 32.217
> 1498900.000 16.629 15.982 7.637 32.611
> 1499000.000 16.595 15.833 7.688 32.428
> 1499100.000 17.022 15.900 7.549 32.921
> 1499200.000 17.061 15.925 7.964 32.987
> 1499300.000 16.638 15.984 7.585 32.622
> 1499400.000 16.829 16.049 8.463 32.878
> 1499500.000 16.473 15.913 7.610 32.386
> 1499600.000 16.968 16.213 7.663 33.182
> 1499700.000 16.372 16.421 7.631 32.793
> 1499800.000 16.211 16.946 7.495 33.157
> 1499900.000 16.485 16.430 7.782 32.915
>
>
>
>
> On Wed, Jul 20, 2016 at 1:34 PM, Amy Rice <arice3.hawk.iit.edu> wrote:
>
> > 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
> >
>
>
>
> --
> Amy Rice
> Ph.D. Student
> Physics Department
> Illinois Institute of Technology
> _______________________________________________
> 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 Mon Jul 25 2016 - 09:30:03 PDT
