[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.