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
2 0.0967924 0.095035
3 0.208143 0.207425
4 0.203251 0.203339
5 0.214236 0.213071
6 0.192448 0.191573
7 0.18472 0.182894
8 0.101735 0.108403
9 0.134363 0.061669 ***
10 0.0444899 0.049259
11 0.0703591 0.077746
12 0.116684 0.115886
13 0.122237 0.121924
14 0.126424 0.125857
15 0.114384 0.114213
16 0.103868 0.103317
17 0.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 list gmx-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
