[gmx-users] Semiisotropic pressure coupling

[Thomas Piggot]

Hi Chris/Justin,

I observed this impact upon membrane area dramatically with DPPC when 
testing both the water model (TIP3P and CHARMM TIP3P) and the switching 
mechanism (potential and force). See the table 3 in the SI of the 
following paper:

http://pubs.acs.org/doi/abs/10.1021/ct3003157

At the time, the version of GROMACS didn't have force switching and so I 
used NAMD for these simulations (as I suspected the gel phase seen with 
DPPC in GROMACS was a combination of both normal TIP3P and the switching 
off of the potential rather than force).

Cheers

Tom

On 13/02/15 14:37, Justin Lemkul wrote:
>
>
> On 2/13/15 12:18 AM, Christopher Neale wrote:
>> Dear Justin:
>>
>> If you have time, can you please provide a link / reference for your 
>> comment that potential switching leads to errors in the charmm lipid 
>> force field? I tried a google search, but couldn't come up with 
>> anything specific. You mentioned noticeable artifacts, so I'm hoping 
>> that you can point me at these (or is it just that changing anything 
>> about the LJ cutoffs changes the APL, which is pretty general across 
>> all lipid force fields in my experience).
>>
>
> No reference or link; this is just an issue that gets reported to the 
> CHARMM community/developers frequently.  I have a bunch of emails 
> related to this particular issue, but you can't really cite those :)  
> With the increase in activity on the Membrane Builder of CHARMM-GUI, 
> it seemed that people were running into issue due to 
> misinterpretations of literature (though it is clearly stated that 
> forces are switched in the parametrization).  With our introduction of 
> GROMACS 5.0-compatible inputs, the reports of problems have 
> diminished. Everyone reported the same thing - dramatic decrease (10% 
> or more) in APL. You're correct that this is not unique to CHARMM 
> lipids, it just seemed like a frequent problem for us because a lot of 
> users thought "vdwtype = switch" was correct.  Common misconception, 
> that even I screwed up in the past.
>
> -Justin
>
>> Thank you,
>> Chris.
>>
>> ________________________________________
>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se 
>> <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of 
>> Justin Lemkul <jalemkul at vt.edu>
>> Sent: 12 February 2015 20:41
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] Semiisotropic pressure coupling
>>
>> On 2/12/15 8:15 PM, shivangi nangia wrote:
>>> Hello All,
>>>
>>> I am doing all-atom simulation with charmm 36 ff (obtained from 
>>> reverse CG
>>> using the script provided by Tsjerk)
>>> The system consists of popc, 21 AA protein, water and ions.
>>>
>>> After reverse CG I did short EM and 5 ns NVT simulation.
>>> However, when I do NPT simulation, the area per lipid drops from 67 
>>> ang2 to
>>> 57 within 1 ns ( x and y dimensions decrese and z increases).
>>>
>>> I am using semiisotropic pressure coupling.
>>>
>>> I have a question regarding that.
>>> I noticed in example .mdps and the tutorial by Justin Lemkul
>>>
>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/membrane_protein/Files/npt.mdp 
>>>
>>>
>>>
>>> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/membrane_protein/Files/md.mdp 
>>>
>>>
>>>
>>> tau_p is specified, the ref_p and compresibility is same, does not that
>>> mean essentially I am doing isotropic coupling?
>>>
>>
>> No, that does not indicate isotropic pressure coupling.  The full 
>> pressure
>> tensor is used to calculate pressures in the x-y and z dimensions 
>> separately.
>> The values of tau_p and compressibility are simply used for 
>> calculating the
>> response time of the barostat.
>>
>> The problem is that your nonbonded settings are incorrect:
>>
>>> ; nblist cut-off
>>> rlist                    = 1.2
>>>
>>> ; OPTIONS FOR ELECTROSTATICS AND VDW
>>> ; Method for doing electrostatics
>>> coulombtype              = PME
>>> rcoulomb-switch          = 0
>>> rcoulomb                 = 1.2
>>> fourierspacing           = 0.16
>>> pme_order                = 4
>>> optimize-fft            = yes
>>> ; Method for doing Van der Waals
>>> vdw-type                 = switch
>>> ; cut-off lengths
>>> rvdw-switch              = 1.0
>>> rvdw                     = 1.2
>>>
>>
>> Note that potential switching leads to errors in the lipid force 
>> field.  You
>> need a buffered neighbor list with force switching.  The exact 
>> parameters that
>> you should use are listed on
>> http://www.gromacs.org/Documentation/Terminology/Force_Fields/CHARMM. 
>> For
>> lipids, you should absolutely not deviate from these settings. The 
>> protein and
>> nucleic acid force fields (and the remainder of CHARMM additive 
>> models) should
>> use these settings, though they are more forgiving (i.e. potential 
>> switching
>> does not lead to noticeable artifacts).  With lipids, the 
>> requirements are
>> pretty absolute.  People frequently report a drop in APL; every time 
>> they do,
>> it's because they're using the wrong settings.
>>
>> -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.
>>
>

-- 
Dr Thomas Piggot
University of Southampton, UK.