[gmx-users] DMPC diffusion: G54A7

[Jon Kapla]

Dear users,

I'm trying out the relatively new phospholipid parameters from Poger et. 
al. 2010 (DOI: 10.1002/jcc.21396), by simulating a box of a DMPC bilayer 
and approx. 30 SPC water molecules per lipid. I'm getting very low 
lateral diffusion compared to both other force fields (like Berger) and 
experimental data.

The topology, force field (G54A7) and a preequilibrated box of DMPC and 
water are all downloaded from ATB (http://compbio.biosci.uq.edu.au/atb/) 
and corrected according to mail conversations with Poger himself (only 
some names that needs to be swapped in the itp file to make it fit the pdb).

The latest simulations I performed at 315K gave me a lateral diffusion 
coefficients between 1.2*10^-8 cm^2/s and 2.5*10^-8 cm^2/s depending on 
treatment of electrostatics and simulation length (from 50 to 350 ns). 
These numbers are around ten times smaller than numbers typically 
reported for similar systems at similar temperatures (and also raising 
the temperature to 350K increases the diffusion by a factor of 2 or 3).

Does anyone have any clue of what (if anything) I might be doing wrong?

I paste the contents of a typical input file below, I have tried to be 
as close to what Poger used in his paper as possible (using PME instead 
of RF though). I've tried some different cutoffs for rlist (0.8, 0.9 and 
1.4) and rcoulomb (0.9,1.4) but nothing of that seems to influence the 
diffusion much.

Thank you!

Kind regards,
Jon Kapla

  RUN CONTROL PARAMETERS
integrator               = md
; Start time and timestep in ps
tinit                    = 0
dt                       = 0.002
nsteps                   = 50000000
; For exact run continuation or redoing part of a run
init_step                = 0
; Part index is updated automatically on checkpointing (keeps files 
separate)
simulation_part          = 1
; mode for center of mass motion removal
comm-mode                = Linear
; number of steps for center of mass motion removal
nstcomm                  = 10
; group(s) for center of mass motion removal
comm-grps                = upper lower SOL
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                  = 2500
nstvout                  = 2500
nstfout                  = 0
; Output frequency for energies to log file and energy file
nstlog                   = 500
nstcalcenergy            = -1
nstenergy                = 500
; Output frequency and precision for .xtc file
nstxtcout                = 0
xtc-precision            = 1000
; This selects the subset of atoms for the .xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps                 =
; Selection of energy groups
energygrps               =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist                  = 5
; ns algorithm (simple or grid)
ns_type                  = grid
; Periodic boundary conditions: xyz, no, xy
pbc                      = xyz
periodic_molecules       = no
; nblist cut-off
rlist                    = 0.8
; long-range cut-off for switched potentials
rlistlong                = -1

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = PME
rcoulomb-switch          = 0
rcoulomb                 = 1.4
; Method for doing Van der Waals
vdw-type                 = Cut-off
; cut-off lengths
rvdw-switch              = 0
rvdw                     = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = no
; Extension of the potential lookup tables beyond the cut-off
table-extension          = 1
; Seperate tables between energy group pairs
energygrp_table          =
; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx               = 0
fourier_ny               = 0
fourier_nz               = 0
; EWALD/PME/PPPM parameters
pme_order                = 4
ewald_rtol               = 1e-05
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = no
;
Tcoupl                   = berendsen
tc-grps                  = DMPC SOL
tau_t                    = 0.1 0.1
ref_t                    = 315 315
ld_seed                  = -1
;
Pcoupl                   = berendsen
Pcoupltype               = Semiisotropic
tau_p                    = 1.0 1.0
compressibility          = 4.6e-5 4.6e-5
ref_p                    = 1.0 1.0
;
; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel                  = yes
gen_temp                 = 315
gen_seed                 = 173529
; OPTIONS FOR BONDS
constraints              = all-bonds
; Type of constraint algorithm
constraint_algorithm     = Lincs
; Do not constrain the start configuration
continuation             = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR                = no
; Relative tolerance of shake
shake-tol                = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order              = 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter               = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle          = 30
; Convert harmonic bonds to morse potentials
morse                    = no


-- 
_____________________________________________________

Jon Kapla
   Division of Physical Chemistry
   Dpt. of Materials and Environmental Chemistry (MMK)
   Arrhenius Laboratory
   Stockholm University
   SE-106 91 Stockholm
Pos:    PhD Student
Phone:  +46 8 16 11 79 (office)
Phone:  +46 70 304 19 89 (cell)
E-mail: jon.kapla at mmk.su.se
_____________________________________________________