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
Received on Wed Jul 20 2016 - 14:30:03 PDT
