Hi,
There was in fact a bug in the 'stfcdiffusion' command when using the
standard calculation mode (i.e. not center of mass or with a second
mask) when the first atom mask selection did not start at one - this
is probably why the results for the second leaflet were bad with
'stfcdiffusion'. If you get a chance, could you try repeating the
calculation with the GitHub version of cpptraj
(
http://github.com/Amber-MD/cpptraj)? I will also be preparing a patch
for AmberTools.
Thanks for the report,
-Dan
PS - I want to mention that this bug is my fault (not Hannes') and is
not present in the original modified ptraj version of this command.
On Tue, Jul 26, 2016 at 1:53 PM, Daniel Roe <daniel.r.roe.gmail.com> wrote:
> Hi,
>
> It seems like there may be a bug in the 'stfcdiffusion' command. I am
> going to look into it and report back.
>
> -Dan
>
> On Wed, Jul 20, 2016 at 12:57 PM, Amy Rice <arice3.hawk.iit.edu> wrote:
>> 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
>>
>
>
>
> --
> -------------------------
> 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)
--
-------------------------
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
Received on Thu Jul 28 2016 - 13:30:02 PDT
