[gmx-users] Reproducing Charmm36 POPC electron density

[Justin Lemkul]

On 2/24/16 7:25 AM, Gmx QA wrote:
> 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?
>

Why?  The results match almost perfectly.  The magnitude of error is within a 
few % and if your APL matches then I'm sure everything is fine.  We rigorously 
tested the lipid inputs, cutoff setup, etc. and everything agrees quite well 
between lipid simulations in CHARMM, GROMACS, AMBER, NAMD, and OpenMM.  There 
are minor variations, as would be expected in comparing any two software 
programs, but nothing is really aberrant.

http://dx.doi.org/10.1021/acs.jctc.5b00935

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================