[gmx-users] Reproducing Charmm36 POPC electron density

[Justin Lemkul]

On 2/24/16 3:36 PM, Gmx QA wrote:
> Hi Justin,
>
> Thanks, that's reassuring then. My reason for worrying was in part because
> the reported electron densities in several papers I have read seem to match
> very closely. But since as you say the APL values also agree then
> everything might be ok after all?
>

We have been rigorously testing CHARMM-GUI for years, and the GROMACS conversion 
over many months.  The input .mdp files and topologies that we provide produce 
identical forces in single-point energy evaluations, which indicate the 
simulation outcomes should be compatible.  Your concern about the electron 
density in the middle of the membrane looks to be somewhere on the order of 2-4% 
difference, which, in the absence of any sort of error/fluctuation bars, strikes 
me as a trivial difference for something as complex as a membrane.

-Justin

> Thanks
> /PK
>
>
> 2016-02-24 15:52 GMT+01:00 Justin Lemkul <jalemkul at vt.edu>:
>
>>
>>
>> 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
>>
>> ==================================================
>> --
>> Gromacs Users mailing list
>>
>> * Please search the archive at
>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
>> posting!
>>
>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>>
>> * For (un)subscribe requests visit
>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
>> send a mail to gmx-users-request at gromacs.org.
>>

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

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

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