[gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[chris . neale]

I have patched g_order version 4.5.4 to compute the order parameters  
in a different way. Despite trying, I don't understand the original  
g_order implementation. In the version that I am providing, the  
hydrogens are explicitly built and the coordinates of these hydrogens  
are output to stdout.


Compared to the original g_order: there are large differences for the  
2 carbons in the double bond and all single bonded carbons are the  
same within what appears to be rounding error.


Compared to the VMD script that Thomas Piggot mentioned, the results  
differ more, but only because that VMD script is using the real  
hydrogen positions and not some with idealized geometry. There is a  
larger difference for the double-bonded carbons, but if I use the  
idealized hydrogen positions (output by my new version of g_order)  
with the VMD .tcl script then I get the same values.


notes on the patch:
1. it is not complete for all options. Use it only with g_order -od  
and with or without -unsat.
2. There are extraneous print statement that will provide you with the  
coordinates of the atoms that are built. This is for debugging. Either  
run with g_order >/dev/null or comment out the 3 fprintf lines near  
the comments containing the word LOOK in capitals.

3. I did my best, but it could always be in error. Please use with caution.

To apply the patch (version 4.5.4)
cd src/tools
patch < g_order.patch

If it will not compile, then the spacing/hard-returns of the patch  
must have been mangled in the process of posting. Contact me off-list  
and I will provide you with a copy.


Here is the patch:

--- gmx_order.c 2011-03-04 06:10:44.0 -0500
+++ gmx_order.c 2011-06-09 17:57:11.586591000 -0400
@@ -379,6 +379,13 @@
real arcdist;
gmx_rmpbc_t  gpbc=NULL;

+// BEGIN CN MOD  


+// variables for modification by CN
+rvec ab,ac,bc,e,ce,eg,ef,g,f,cg,cf;
+float dab,def,deg,gdot,fdot,edot;
+// END CN MOD  


+
+
   /* PBC added for center-of-mass vector*/
   /* Initiate the pbc structure */
   memset(&pbc,0,sizeof(pbc));
@@ -501,98 +508,66 @@
direction[0],direction[1],direction[2]);*/
  }

-   if (bUnsat) {
- /* Using convention for unsaturated carbons */
- /* first get Sz, the vector from Cn to Cn+1 */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
- length = norm(dist);
- check_length(length, a[index[i]+j], a[index[i+1]+j]);
- svmul(1/length, dist, Sz);
-
- /* this is actually the cosine of the angle between the double bond
-and axis, because Sz is normalized and the two other components of
-the axis on the bilayer are zero */
- if (use_unitvector)
- {
-	sdbangle += acos(iprod(direction,Sz));  /*this can probably be  
optimized*/

- }
- else
-   sdbangle += acos(Sz[axis]);
-   } else {
- /* get vector dist(Cn-1,Cn+1) for tail atoms */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
- length = norm(dist);  /* determine distance between two atoms */
- check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
-
- svmul(1/length, dist, Sz);
- /* Sz is now the molecular axis Sz, normalized and all that */
-   }
-
-   /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
-  we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
-   rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
-   rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
-   cprod(tmp1, tmp2, Sx);
-   svmul(1/norm(Sx), Sx, Sx);
-
-   /* now we can get Sy from the outer product of Sx and Sz   */
-   cprod(Sz, Sx, Sy);
-   svmul(1/norm(Sy), Sy, Sy);
-
-   /* the square of cosine of the angle between dist and the axis.
-  Using the innerproduct, but two of the three elements are zero
-  Determine the sum of the orderparameter of all atoms in group
-  */
-   if (use_unitvector)
-   {
-	cossum[XX] = sqr(iprod(Sx,direction)); /* this is allowed, since Sa  
is normalized */

-   cossum[YY] = sqr(iprod(Sy,direction));
-   cossum[ZZ] = sqr(iprod(Sz,direction));
-   }
-   else
-   {
-   cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized 
*/
-   cossum[YY] = sqr(Sy[axis]);
-   cossum[ZZ] = sqr(Sz[axis]);
+// BEGIN CN MOD  


+  if (bUnsat) {
+//exactly as for the case without unsat, but here we select the  
ce vector to represent the hydrogen
+// place the 2 imaginary hydrogen atoms by enscribing a  
tetrahedron in a cube
+// a,b,c are the 3 backbone atoms, where a and b are in two cube  
corn

Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Justin A. Lemkul]

chris.ne...@utoronto.ca wrote:

Dear Thomas:

I agree with your conclusions about g_order.

My usage of the Berger POPC lipids in combination with the g_order 
command (as per the instructions in 
http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order ) yields 
order parameters of 0.193 and 0.051 for the first and second
carbons in the double bond at 300K. This is similar to but, as expected, 
slightly larger than the values described in 
http://www.sciencedirect.com/science/article/pii/S0014579305014365 where 
they simulated at 310K and also used g_order. I suspect that one could 
find tens of publications that use g_order and have published similar 
values. The Poger paper (On the Validation of Molecular Dynamics 
Simulations of Saturated and cis-Monounsaturated Phosphatidylcholine 
Lipid Bilayers: A Comparison with Experiment) gives clearly different 
values for the first carbon in the double bond than is reported by 
g_order. One probable reason that this has gone unnoticed so long is 
that there appears to be no experimental oleoyl order parameters -- as 
far as I know -- if somebody finds them, please post them.




Cited in Klauda:

dx.doi.org/10.1021/bi00338a027

I also found:

dx.doi.org/10.1007/s00249-005-0482-z

-Justin

I took a look at the source and it is clearly incorrect for the i+1 
unsaturated partner in a double bond over i->i+1. I have developed a new 
version of g_order to fix this, but I'd like to test it first. If we 
could agree on a single .xtc file to analyze and compare our values then 
that would be great. I can provide you with one of mine (Berger lipids) 
or you can provide me with one of yours (off-list).


Thanks,
Chris.




--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[chris . neale]

Dear Thomas:

I agree with your conclusions about g_order.

My usage of the Berger POPC lipids in combination with the g_order  
command (as per the instructions in  
http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order )  
yields order parameters of 0.193 and 0.051 for the first and second
carbons in the double bond at 300K. This is similar to but, as  
expected, slightly larger than the values described in  
http://www.sciencedirect.com/science/article/pii/S0014579305014365  
where they simulated at 310K and also used g_order. I suspect that one  
could find tens of publications that use g_order and have published  
similar values. The Poger paper (On the Validation of Molecular  
Dynamics Simulations of Saturated and cis-Monounsaturated  
Phosphatidylcholine Lipid Bilayers: A Comparison with Experiment)  
gives clearly different values for the first carbon in the double bond  
than is reported by g_order. One probable reason that this has gone  
unnoticed so long is that there appears to be no experimental oleoyl  
order parameters -- as far as I know -- if somebody finds them, please  
post them.


I took a look at the source and it is clearly incorrect for the i+1  
unsaturated partner in a double bond over i->i+1. I have developed a  
new version of g_order to fix this, but I'd like to test it first. If  
we could agree on a single .xtc file to analyze and compare our values  
then that would be great. I can provide you with one of mine (Berger  
lipids) or you can provide me with one of yours (off-list).


Thanks,
Chris.


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Thomas Piggot]

Further to my last message we have managed to use a tcl script in VMD to 
calculate the order parameters for the POPC chains in the CHARMM36 
simulation. I have put the values for this calculation and the values I 
obtained using g_order (over the same period in the same trajectory) 
below (this is for the sn-2 chain, for the sn-1 chains then the 2 
methods produce very similar results). I think this shows that there is 
a problem with g_order and the calculation of the order parameter 
parameters for the first carbon in the double bond of the unsaturated tail.


   carbon   g_order VMD script

   20.0967924   0.095035
   30.2081430.207425
   40.2032510.203339
   50.2142360.213071
   60.1924480.191573
   70.18472 0.182894
   80.1017350.108403
   90.1343630.061669***
  100.0444899   0.049259
  110.0703591   0.077746
  120.1166840.115886
  130.1222370.121924
  140.1264240.125857
  150.1143840.114213
  160.1038680.103317
  170.0806387   0.080544

If nobody can see that I have done anything obviously wrong then I will 
file an issue regarding this on redmine.


Cheers

Tom

Thomas Piggot wrote:

Hi,

Both the trajectories are 100 ns in length (plus I have a repeat of each 
simulation using different starting velocities). The results from a 
block analysis, and analysis of the repeat simulations, give almost 
identical results to those in my previous message. The other thing is 
that the remainder of the order parameters also give sensible numbers, 
both for the saturated carbons in the rest of the oleoyl chain and in 
the palmitoyl chain.


Regarding the mdp settings, then I used the same parameters as given in 
the papers from which I quoted the previous values. For CHARMM36, this 
includes a 1 fs timestep, the CHARMM TIP3P water (tips3p), constraints = 
hbonds and a switching of the vdW interactions between 0.8-1.2 A. For 
the GROMOS53A6L simulations then this is with a cut-off for the 
electrostatic interactions of 1.4 nm and a reaction-field correction 
applied beyond the cut-off, a vdW cut-off of 1.4 nm, rlist of 0.8, 
nstlist every 5 steps, a 2 fs time step and the SPC water model. I can 
provide full mdp's if necessary. Analysis of other membrane properties 
from the simulations (such as area and volume per lipid, area 
compressibility, electron density) all give comparable results to those 
previously published for these two force fields.


All of this leads me to suspect that there may be a problem in the way 
g_order calculates the order parameter for the unsaturated carbons, 
however I am by no means certain of this. Currently we are trying to 
calculate the order parameters from the all-atom CHARMM36 simulations 
using a different bit of code, I will report back with the results if we 
can get this to work.


Cheers

Tom

Justin A. Lemkul wrote:

Thomas Piggot wrote:

Hi Igor,

Thanks for the reply but I think you slightly missed the point I was 
trying to make. I followed the approach you mention for the double bond 
and (with the two force fields I tried) got the values I discussed in my 
last email, which when compared to the published values seem to be 
incorrect.


Just to confirm, my index file (for the unsaturated calculation) 
contains (below is just pasted from make_ndx):


For CHARMM36:

  0 C28 :72 atoms
  1 C29 :72 atoms
  2 C210:72 atoms
  3 C211:72 atoms

and for GROMOS53A6L:

  0 C1H :   128 atoms
  1 C1I :   128 atoms
  2 C1J :   128 atoms
  3 C1K :   128 atoms


My command for g_order (this is just for CHARMM36 but I use the same 
command for GROMOS53A6L on different files) is:


g_order_4.5.4 -s for-md-popc-charmm36.tpr -f md-popc-charmm36.xtc -n 
chain2_unsat.ndx -od deut_chain2_unsat.xvg -unsat



If anyone can see what might be going wrong then any help would be 
greatly appreciated.


Are your results converged?  From the g_order command, you're considering the 
whole trajectory.  I'd suggest the usual block averaging approach to see if 
you're converged.


There's nothing wrong with the way you're calculating the order parameters, so 
the only other thing I would suspect (other than lack of convergence) is that 
the .mdp settings are somehow giving unexpected results.  Membrane properties 
can be very sensitive to vdW and electrostatics methods.


-Justin





--
Dr Thomas Piggot
University of Southampton, UK.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.or

Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Thomas Piggot]

Hi,

Both the trajectories are 100 ns in length (plus I have a repeat of each 
simulation using different starting velocities). The results from a 
block analysis, and analysis of the repeat simulations, give almost 
identical results to those in my previous message. The other thing is 
that the remainder of the order parameters also give sensible numbers, 
both for the saturated carbons in the rest of the oleoyl chain and in 
the palmitoyl chain.


Regarding the mdp settings, then I used the same parameters as given in 
the papers from which I quoted the previous values. For CHARMM36, this 
includes a 1 fs timestep, the CHARMM TIP3P water (tips3p), constraints = 
hbonds and a switching of the vdW interactions between 0.8-1.2 A. For 
the GROMOS53A6L simulations then this is with a cut-off for the 
electrostatic interactions of 1.4 nm and a reaction-field correction 
applied beyond the cut-off, a vdW cut-off of 1.4 nm, rlist of 0.8, 
nstlist every 5 steps, a 2 fs time step and the SPC water model. I can 
provide full mdp's if necessary. Analysis of other membrane properties 
from the simulations (such as area and volume per lipid, area 
compressibility, electron density) all give comparable results to those 
previously published for these two force fields.


All of this leads me to suspect that there may be a problem in the way 
g_order calculates the order parameter for the unsaturated carbons, 
however I am by no means certain of this. Currently we are trying to 
calculate the order parameters from the all-atom CHARMM36 simulations 
using a different bit of code, I will report back with the results if we 
can get this to work.


Cheers

Tom

Justin A. Lemkul wrote:


Thomas Piggot wrote:

Hi Igor,

Thanks for the reply but I think you slightly missed the point I was 
trying to make. I followed the approach you mention for the double bond 
and (with the two force fields I tried) got the values I discussed in my 
last email, which when compared to the published values seem to be 
incorrect.


Just to confirm, my index file (for the unsaturated calculation) 
contains (below is just pasted from make_ndx):


For CHARMM36:

  0 C28 :72 atoms
  1 C29 :72 atoms
  2 C210:72 atoms
  3 C211:72 atoms

and for GROMOS53A6L:

  0 C1H :   128 atoms
  1 C1I :   128 atoms
  2 C1J :   128 atoms
  3 C1K :   128 atoms


My command for g_order (this is just for CHARMM36 but I use the same 
command for GROMOS53A6L on different files) is:


g_order_4.5.4 -s for-md-popc-charmm36.tpr -f md-popc-charmm36.xtc -n 
chain2_unsat.ndx -od deut_chain2_unsat.xvg -unsat



If anyone can see what might be going wrong then any help would be 
greatly appreciated.




Are your results converged?  From the g_order command, you're considering the 
whole trajectory.  I'd suggest the usual block averaging approach to see if 
you're converged.


There's nothing wrong with the way you're calculating the order parameters, so 
the only other thing I would suspect (other than lack of convergence) is that 
the .mdp settings are somehow giving unexpected results.  Membrane properties 
can be very sensitive to vdW and electrostatics methods.


-Justin



--
Dr Thomas Piggot
University of Southampton, UK.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Justin A. Lemkul]

Thomas Piggot wrote:

Hi Igor,

Thanks for the reply but I think you slightly missed the point I was 
trying to make. I followed the approach you mention for the double bond 
and (with the two force fields I tried) got the values I discussed in my 
last email, which when compared to the published values seem to be 
incorrect.


Just to confirm, my index file (for the unsaturated calculation) 
contains (below is just pasted from make_ndx):


For CHARMM36:

  0 C28 :72 atoms
  1 C29 :72 atoms
  2 C210:72 atoms
  3 C211:72 atoms

and for GROMOS53A6L:

  0 C1H :   128 atoms
  1 C1I :   128 atoms
  2 C1J :   128 atoms
  3 C1K :   128 atoms


My command for g_order (this is just for CHARMM36 but I use the same 
command for GROMOS53A6L on different files) is:


g_order_4.5.4 -s for-md-popc-charmm36.tpr -f md-popc-charmm36.xtc -n 
chain2_unsat.ndx -od deut_chain2_unsat.xvg -unsat



If anyone can see what might be going wrong then any help would be 
greatly appreciated.




Are your results converged?  From the g_order command, you're considering the 
whole trajectory.  I'd suggest the usual block averaging approach to see if 
you're converged.


There's nothing wrong with the way you're calculating the order parameters, so 
the only other thing I would suspect (other than lack of convergence) is that 
the .mdp settings are somehow giving unexpected results.  Membrane properties 
can be very sensitive to vdW and electrostatics methods.


-Justin

--


Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin


--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Thomas Piggot]

Hi Igor,

Thanks for the reply but I think you slightly missed the point I was 
trying to make. I followed the approach you mention for the double bond 
and (with the two force fields I tried) got the values I discussed in my 
last email, which when compared to the published values seem to be 
incorrect.


Just to confirm, my index file (for the unsaturated calculation) 
contains (below is just pasted from make_ndx):


For CHARMM36:

  0 C28 :72 atoms
  1 C29 :72 atoms
  2 C210:72 atoms
  3 C211:72 atoms

and for GROMOS53A6L:

  0 C1H :   128 atoms
  1 C1I :   128 atoms
  2 C1J :   128 atoms
  3 C1K :   128 atoms


My command for g_order (this is just for CHARMM36 but I use the same 
command for GROMOS53A6L on different files) is:


g_order_4.5.4 -s for-md-popc-charmm36.tpr -f md-popc-charmm36.xtc -n 
chain2_unsat.ndx -od deut_chain2_unsat.xvg -unsat



If anyone can see what might be going wrong then any help would be 
greatly appreciated.


Cheers

Tom

Igor Marques wrote:

thomas,

i've recently placed a similar question and justin asked me to show my 
index for the double bound calculation, so, i'm asking you the same.

you should have

Ci-1
Ci - the first carbon of the double bound
Ci+1 - the second carbon of the double bound
Ci+2

in that index

doing this, for C9 and C10 i've obtained the following values
 0.0738068
  -0.000533683

that you should then replace in the correct positions in the sn2 output


good luck,
and sorry if i'm missing the point :|


igor


  Igor Marques


On Wed, Jun 8, 2011 at 5:53 PM, Thomas Piggot > wrote:


Hi Everyone,

I am facing a problem when calculating the lipid deuterium order
parameters for the unsaturated carbons of the sn-2 tail of POPC
using g_order with GROMACS version 4.5.4 (although I have tried
other older versions too but they all give the same results).

Firstly, I should say the the calculation of the order parameters
for the saturated sn-1 chain (and also both chains of DPPC) behave
as I would expect, and produce order parameters that compare well to
previously published simulations and experimental values.

To calculate the order parameters of the unsaturated chain I am
following the approach as given on the GROMACS website
(http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order), so
splitting the calculation into two parts for the saturated and
unsaturated regions of the chain. The problem I am facing is that
the order parameter for carbon 9 (so the first carbon in the double
bond), calculated using the -unsat option, is much larger than
expected. By this I mean that for the two different force fields I
have tested (namely the CHARMM36 parameters of Klauda et al., and
the GROMOS 53A6L parameters of Poger et al.) the order parameter for
this carbon is much larger than the published simulation values and
also much larger than experimental values. To highlight this, I have
just put the numbers I have obtained using g_order for this carbon
below, and compared to some rough values I have estimated from
figures provided in the Klauda and Poger papers:

CHARMM36

g_order:0.133732
Klauda estimate:0.06



GROMOS53A6L:

g_order:0.199651
Poger estimate: 0.07

Myself and a colleague have tried looking into the code to determine
how the order parameters are calculated using the -unsat option,
however we couldn't quite follow the calculation. Hopefully someone
who knows something more about g_order can help with this problem.
Again I should stress that it appears that the main difference in
order parameters between what I have calculated and the published
ones is just in this one carbon, for both force fields. A similar
issue to this has also been reported previously on the list for this
carbon of POPC using the 'Berger' force field
(http://lists.gromacs.org/pipermail/gmx-users/2010-February/049072.html).

Thank you for any help anyone can give

Cheers

Tom

-- 
Dr Thomas Piggot

University of Southampton, UK.
-- 
gmx-users mailing listgmx-users@gromacs.org


http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the www
interface or send it to gmx-users-requ...@gromacs.org
.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists




--
Dr Thomas Piggot
University of Southampton, UK.
--
gmx-users mailing listgmx-users@gromacs.org

Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Igor Marques]

thomas,

i've recently placed a similar question and justin asked me to show my index
for the double bound calculation, so, i'm asking you the same.
you should have

> Ci-1
> Ci - the first carbon of the double bound
> Ci+1 - the second carbon of the double bound
> Ci+2
>
in that index

doing this, for C9 and C10 i've obtained the following values
 0.0738068
  -0.000533683

that you should then replace in the correct positions in the sn2 output


good luck,
and sorry if i'm missing the point :|


igor


  Igor Marques


On Wed, Jun 8, 2011 at 5:53 PM, Thomas Piggot  wrote:

> Hi Everyone,
>
> I am facing a problem when calculating the lipid deuterium order parameters
> for the unsaturated carbons of the sn-2 tail of POPC using g_order with
> GROMACS version 4.5.4 (although I have tried other older versions too but
> they all give the same results).
>
> Firstly, I should say the the calculation of the order parameters for the
> saturated sn-1 chain (and also both chains of DPPC) behave as I would
> expect, and produce order parameters that compare well to previously
> published simulations and experimental values.
>
> To calculate the order parameters of the unsaturated chain I am following
> the approach as given on the GROMACS website (
> http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order), so
> splitting the calculation into two parts for the saturated and unsaturated
> regions of the chain. The problem I am facing is that the order parameter
> for carbon 9 (so the first carbon in the double bond), calculated using the
> -unsat option, is much larger than expected. By this I mean that for the two
> different force fields I have tested (namely the CHARMM36 parameters of
> Klauda et al., and the GROMOS 53A6L parameters of Poger et al.) the order
> parameter for this carbon is much larger than the published simulation
> values and also much larger than experimental values. To highlight this, I
> have just put the numbers I have obtained using g_order for this carbon
> below, and compared to some rough values I have estimated from figures
> provided in the Klauda and Poger papers:
>
> CHARMM36
>
> g_order:0.133732
> Klauda estimate:0.06
>
>
> GROMOS53A6L:
>
> g_order:0.199651
> Poger estimate: 0.07
>
> Myself and a colleague have tried looking into the code to determine how
> the order parameters are calculated using the -unsat option, however we
> couldn't quite follow the calculation. Hopefully someone who knows something
> more about g_order can help with this problem. Again I should stress that it
> appears that the main difference in order parameters between what I have
> calculated and the published ones is just in this one carbon, for both force
> fields. A similar issue to this has also been reported previously on the
> list for this carbon of POPC using the 'Berger' force field (
> http://lists.gromacs.org/pipermail/gmx-users/2010-February/049072.html).
>
> Thank you for any help anyone can give
>
> Cheers
>
> Tom
>
> --
> Dr Thomas Piggot
> University of Southampton, UK.
> --
> gmx-users mailing listgmx-users@gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
> Please don't post (un)subscribe requests to the list. Use the www interface
> or send it to gmx-users-requ...@gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
-- 
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

[gmx-users] Calculation of unsaturated deuterium order parameters for POPC

[Thomas Piggot]

Hi Everyone,

I am facing a problem when calculating the lipid deuterium order 
parameters for the unsaturated carbons of the sn-2 tail of POPC using 
g_order with GROMACS version 4.5.4 (although I have tried other older 
versions too but they all give the same results).


Firstly, I should say the the calculation of the order parameters for 
the saturated sn-1 chain (and also both chains of DPPC) behave as I 
would expect, and produce order parameters that compare well to 
previously published simulations and experimental values.


To calculate the order parameters of the unsaturated chain I am 
following the approach as given on the GROMACS website 
(http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order), so 
splitting the calculation into two parts for the saturated and 
unsaturated regions of the chain. The problem I am facing is that the 
order parameter for carbon 9 (so the first carbon in the double bond), 
calculated using the -unsat option, is much larger than expected. By 
this I mean that for the two different force fields I have tested 
(namely the CHARMM36 parameters of Klauda et al., and the GROMOS 53A6L 
parameters of Poger et al.) the order parameter for this carbon is much 
larger than the published simulation values and also much larger than 
experimental values. To highlight this, I have just put the numbers I 
have obtained using g_order for this carbon below, and compared to some 
rough values I have estimated from figures provided in the Klauda and 
Poger papers:


CHARMM36

g_order:0.133732
Klauda estimate:0.06


GROMOS53A6L:

g_order:0.199651
Poger estimate: 0.07

Myself and a colleague have tried looking into the code to determine how 
the order parameters are calculated using the -unsat option, however we 
couldn't quite follow the calculation. Hopefully someone who knows 
something more about g_order can help with this problem. Again I should 
stress that it appears that the main difference in order parameters 
between what I have calculated and the published ones is just in this 
one carbon, for both force fields. A similar issue to this has also been 
reported previously on the list for this carbon of POPC using the 
'Berger' force field 
(http://lists.gromacs.org/pipermail/gmx-users/2010-February/049072.html).


Thank you for any help anyone can give

Cheers

Tom

--
Dr Thomas Piggot
University of Southampton, UK.
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-requ...@gromacs.org.

Can't post? Read http://www.gromacs.org/Support/Mailing_Lists