[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

            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.