[gmx-users] Calculation of unsaturated deuterium order parameters for POPC
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 corners
Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC
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 t.pig...@soton.ac.uk mailto:t.pig...@soton.ac.uk 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 mailto: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 mailto: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 list
Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC
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
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
[gmx-users] Calculation of unsaturated deuterium order parameters for POPC
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
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
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
Re: [gmx-users] Calculation of unsaturated deuterium order parameters for POPC
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 t.pig...@soton.ac.uk 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