[gmx-users] Reproducing Charmm36 POPC electron density

[Gmx QA]

Dear list

I have a simulation of a system of 64 POPC lipids per leaflet using the
Charmm36 lipid parameter set, and I am calculating the electron density
across the bilayer to compare to literature values. Unfortunately I cannot
quite get the same results as eg. in Klauda et al (
http://www.sciencedirect.com/science/article/pii/S0005273614002193)

I have made the starting structure using the charmm-gui, and gone through
all the equilibration steps as suggested in the gromacs tar-ball from that
website.

My production mdp file is the following:

integrator              = md
dt                      = 0.002
nsteps                  = 50000000
nstlog                  = 10000
nstxout                 = 500000
nstvout                 = 500000
nstfout                 = 0
nstcalcenergy           = 10
nstenergy               = 1000
nstxtcout               = 500000
;
cutoff-scheme           = Verlet
nstlist                 = 20
rlist                   = 1.2
coulombtype             = pme
rcoulomb                = 1.2
vdwtype                 = Cut-off
vdw-modifier            = Force-switch
rvdw_switch             = 1.0
rvdw                    = 1.2
;
tcoupl                  = Nose-Hoover
tc_grps                 = POPC SOL
tau_t                   = 1.0 1.0
ref_t                   = 303.15  303.15
;nsttcouple = 1
;
pcoupl                  = Parrinello-Rahman
pcoupltype              = semiisotropic
tau_p                   = 5.0 5.0
compressibility         = 4.5e-5 4.5e-5
ref_p                   = 1.0 1.0
;
constraints             = h-bonds
constraint_algorithm    = LINCS
continuation            = yes
;
nstcomm                 = 10
comm_mode               = linear
comm_grps               = POPC SOL
;
refcoord_scaling        = com

Simulations are run with gromacs 5.0.4, and currently 100 ns long.

Before running vmx density, I center my trajectory like this:
$ gmx trjconv -f run.trr -s em.tpr -o run_center.trr -center

selecting POPC as the group to center on.

 Then, my command for computing the electron density is the following (also
using gmx density 5.0.4):

$ gmx density -f run_center.trr -s em.tpr -o  density_traj_center.xvg -dens
electron -ei electrons.dat -center

Again, selecting POPC to center the density profile on.

The resulting profile can be found here (i hope that works)
http://www.filedropper.com/compareelectrondensity

As you can see, the depth of the central trough is not as deep as the
literature value, nor is the height of the density the same even though the
peak-peak distance appears to be similar, and values in the other regions
also match _fairly_ closely.

The average APL value over the entire trajectory is 64.0 +- 1.33, which is
in line with the reported 64.7 from Klauda, but the electron density
worries me a bit. Perhaps anyone (Justin?) can shed any light on this?

Thanks
/PK