[gmx-users] Gromos54a7 electrostatics interactions ?
Sim gmx
simgmx at gmail.com
Mon Sep 25 11:06:26 CEST 2017
Thanks again for your reply and all the info that I've read with great
interest!
I would like to simulate lipid bilayers that contain my molecules and have
a look at my molecules' location and orientation as well as how they can
affect lipids properties. My plan would be to try two different
forcefields: Gromos54a7 and Berger lipids. For gromos54a7, I would use a
mdp file (for run production) like this one:
"integrator = md
nsteps = 125000000
dt = 0.002
nstxout = 500000
nstvout = 500000
nstxtcout = 10000
nstenergy = 1000
nstlog = 1000
continuation = yes
constraint_algorithm = lincs
constraints = all-bonds
lincs_iter = 1
lincs_order = 4
ns_type = grid
nstlist = 5
nstcalclr = 1
vdw-type = cut-off
rlist = 0.8
rcoulomb = 1.4
rvdw = 1.4
coulombtype = reaction-field
epsilon-rf = 62
tcoupl = Nose-Hoover
tc-grps = MOLECULE_LIPID SOL
tau_t = 0.5 0.5
ref_t = 298 298
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 2.0
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
pbc = xyz
DispCorr = No
gen_vel = no
nstcomm = 1
comm-mode = Linear
comm-grps = MOLECULE_LIPID SOL
cutoff-scheme = group"
Berger lipids simulations would be performed with similar mdp options,
except for a few parameters, namely:
"integrator = md
nsteps = 125000000
dt = 0.002
nstxout = 500000
nstvout = 500000
nstxtcout = 10000
nstenergy = 1000
nstlog = 1000
continuation = yes
constraint_algorithm = lincs
constraints = all-bonds
lincs_iter = 1
lincs_order = 4
ns_type = grid
nstlist = 5
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
coulombtype = PME
pme_order = 4
fourierspacing = 0.12
tcoupl = Nose-Hoover
tc-grps = MOLECULE_LIPID SOL
tau_t = 0.5 0.5
ref_t = 298 298
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 2.0
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
pbc = xyz
DispCorr = EnerPres
gen_vel = no
nstcomm = 1
comm-mode = Linear
comm-grps = MOLECULE_LIPID SOL
cutoff-scheme = Verlet"
I would do one repetition of each system with each forcefield. Then I would
compare the results between each other and with experimental results.
Hopefully, both forcefields would give similar results and would be in good
agreement with experimental results, I would then perform two more
repetitions with the faster setup. If both forcefields don't give similar
results, then I would select the one that is in better agreement with
experimental results (hopefully there would be at least one huhu).
The forcefield choice is also guided by the possibility of adding some
lipids (sterols...) if useful (depending on experimental results). With
gromos54a7 I would use sterol topology available on ATB and for Berger
lipids I would start from sterol topology from Höltje et al..
I have for sure to make my own opinion about the different forcefields.
However, I am lucid enough to understand that some people are much more
experienced and qualified than I am. I hope I will have the opportunity to
progress, and I would not blindly follow anyone's opinion though. But I
don't think Mark Abraham's opinion and yours are 'anyone's opinion', so I
carefully read your posts :)
Thanks!
2017-09-22 20:37 GMT+02:00 Thomas Piggot <t.piggot at soton.ac.uk>:
>
>
> On 22/09/17 16:59, Sim gmx wrote:
>
>> First of all, thanks for your replies!
>>
>> 1) Indeed, both epsilon-r and epsilon-rf are such constants. However, in
>> the paper I guess that they refer to the relative dielectric constant of
>> the reaction field, i.e. epsilon-rf. Or do I miss something ?
>>
> Yes, they are referring to epsilon-rf.
>
>> 2) I will for sure contact them to collect all these info if I go further
>> with this forcefield!
>>
> I've tested things pretty extensively with this force field so I don't
> think you need to contact them. If you post your mdp it will make things
> much easier as we can see exactly what you're are doing and can let you
> know if there is anything that doesn't look correct.
>
>>
>> 3) I don't think I understand everything there. What I think I get is:
>> using a verlet scheme with parameters that were not initially designed for
>> this scheme is risky and should be carefully assessed.
>>
> This is what I was saying I have already assessed and, for DPPC at least,
> the membrane does not behave well with Verlet and reaction-field. If you
> use PME, Verlet and group give the same results (but the membrane is more
> ordered with PME compared to RF, irrespective of the cutoff scheme).
>
>> Plus, parameters
>> used by Poger are very likely to be incompatible with verlet scheme.
>> From this, I conclude it's not even worth giving it a try. So my options
>> are:
>> - use a group scheme, that will be slower
>>
> Well, this depends upon your computer. On one HPC I use that is CPU only,
> I get nearly identical performance with Verlet and a straight 1.4 nm cutoff
> with RF (i.e the settings that don't work well) and the fastest group
> twin-range setup that gives appropriate membrane properties. If you are
> wanting to use a GPU this will obviously be different.
>
>> - use a different force field, but which one since I do have the
>> topologies
>> of my small molecules for gromos 53a6/54a7 ? Would Berger Lipids be a good
>> option here ?
>>
> This really depends upon what you are wanting to do. There are other
> GROMOS compatible lipid force fields out there (including Berger) but, as
> with pretty much everything, there are limitations/considerations with them
> all. In addition, you say you already have small molecule topologies.
> Depending upon what these structures are, it may be fairly trivial to make
> ones for other force fields too. This would open up lots of other possible
> lipid force fields too. Again depending upon what you are wanting to look
> at, these could be better/more appropriate for your study. Personally, I'd
> suggest looking at more than once force field, that way you can be more
> confident in any results you obtain. This, obviously, takes more time and
> effort though.
>
> You should probably do a bit of reading about the advantages/disadvantages
> of different (lipid) force fields, that way you could make your own mind up
> rather than just relying on someones opinion on a mailing list! Although if
> you do want my opinion, I'm happy to share it if you can give more details
> of what you are trying to do.
>
>>
>> About Tom's reply:
>>
>> 1. That's good to know, especially as it does not appear anywhere in the
>> original papers.
>>
> The original simulations in the papers were performed with a really old
> version of GROMACS (3.2.1 IIRC) and so this was not an issue. The paper (by
> the same authors) I linked to in my previous message provides all the
> details about this.
>
>> 2. Is using a setup with all cutoffs = 1.4 more "adequate" than the
>> twin-range setup then ? Or is it just a "good to know" stuff ? I would
>> naturally use the fastest correct setup which seems to be the twin-range
>> one (with group-scheme), am I wrong ?
>>
> With the group cutoff scheme, simulations using both the twin-range
> 0.8/1.4 and the straight 1.4 cutoffs result in the same properties of a
> DPPC membrane (again see the paper I linked to; I have also tested this
> too). So yes, you should use the faster setup.
>
>
>> Again, thanks to both of you.
>>
>> 2017-09-22 14:29 GMT+02:00 Piggot T. <T.Piggot at soton.ac.uk>:
>>
>> Hi,
>>>
>>> In addition to Mark's comments, there are a few other points to be aware
>>> of:
>>>
>>> 1. If you use the twin-range scheme (i.e. 0.8/1.4) with nstlist 5, you
>>> need to set nstcalclr to 1 to get results matching those as reported by
>>> Poger and co-workers (see http://pubs.acs.org/doi/abs/
>>> 10.1021/acs.jctc.7b00178?src=recsys). This is the fastest setup that
>>> will
>>> give you the 'correct' membrane properties with the group scheme.
>>>
>>> 2. The combination of reaction field with Verlet and all cutoffs set to
>>> 1.4 does not give the correct membrane properties for DPPC (these are my
>>> own results). Using the group scheme with this setup is fine but slower
>>> than the twin-range setup mentioned above. You should therefore use the
>>> group setup if using reaction-field with this force field. If you wish to
>>> use PME the membrane will be more ordered. However, Verlet and group
>>> schemes give the same results for PME (again my results),
>>>
>>> Regarding the second point, can anyone shed some more light onto why
>>> using
>>> the Verlet scheme would result in substantially different membrane
>>> properties to the group scheme with a single 1.4 nm cutoff and
>>> reaction-field (PME gives the same membrane properties using these two
>>> setups)?
>>>
>>> Cheers
>>>
>>> Tom
>>> ________________________________________
>>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [
>>> gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of Mark
>>> Abraham [mark.j.abraham at gmail.com]
>>> Sent: 22 September 2017 12:41
>>> To: gmx-users at gromacs.org; gromacs.org_gmx-users at maillist.sys.kth.se
>>> Subject: Re: [gmx-users] Gromos54a7 electrostatics interactions ?
>>>
>>> Hi,
>>>
>>> On Fri, Sep 22, 2017 at 10:21 AM Sim gmx <simgmx at gmail.com> wrote:
>>>
>>> Hi!
>>>>
>>>> I would like to do simulations of lipids bilayers with gromos54a7
>>>> parameters. To do so, I want to use the same mdp parameters as Poger et
>>>>
>>> al.
>>>
>>>> (the authors).
>>>>
>>>> In their papers, they write : "Nonbonded interactions were evaluated
>>>>
>>> using
>>>
>>>> a twin-range cutoff scheme: interactions falling within the 0.8-nm
>>>> short-range cutoff were calculated every step, whereas interactions
>>>>
>>> falling
>>>
>>>> within the 1.4-nm long-range cutoff were updated every 10 fs, together
>>>>
>>> with
>>>
>>>> the pair list. A reaction-field correction was applied to account for
>>>> the
>>>> truncation of the electrostatic interactions beyond the long-range
>>>> cutoff
>>>> using a relative dielectric permittivity constant of 62, as appropriate
>>>>
>>> for
>>>
>>>> SPC water"
>>>>
>>>> Note that such twin-range multiple stepping schemes have only ever been
>>> supported in GROMACS in the group cutoff scheme, and not at all in the
>>> most
>>> recent versions. See (older versions of) the reference manual for
>>> discussion of known issues with the twin-range scheme. However,
>>> parameters
>>> derived with this irreversible (and thus unphysical) scheme have usually
>>> been shown to be usable outside that scheme, but obviously there's a
>>> higher
>>> burden of validation bourne by a the user in that case.
>>>
>>> From this, I conclude that my mdp file should contain something like
>>> (among
>>>
>>>> other parameters):
>>>>
>>>> dt=0.002
>>>> vdw-type = cut-off
>>>> rlist = 0.8
>>>> nstlist = 5
>>>> rvdw = 1.4
>>>> coulombtype = reaction-field
>>>> epsilon-rf = 62
>>>>
>>>> My questions are:
>>>>
>>>> 1) Does "relative dielectric permittivity constant" indeed mean
>>>> "epsilon-rf" here (not epsilon-r) ?
>>>>
>>>> That's not the right question, both epsilon and epsilon-rf are such
>>> constants. See
>>> http://manual.gromacs.org/documentation/2016.4/user-
>>> guide/mdp-options.html#electrostatics
>>> for
>>> the different parts of the coulomb contributions to which they refer.
>>>
>>>
>>> 2) What about the coulomb cut-off (rcoulomb) ? I think that in
>>>>
>>> gromos53a6,
>>>
>>>> rcoulomb is often set to 0.9 (here I guess it would be 0.8 to match the
>>>> rlist). However, if the reaction-field correction is applied "to account
>>>> for the truncation of the electrostatic interactions beyond *the
>>>>
>>> long-range
>>>
>>>> cutoff*", I would set rcoulomb to 1.4. As no explicit reference to
>>>>
>>> rcoulomb
>>>
>>>> is done in the papers, I fear to miss something here.
>>>>
>>>> They refer to non-bonded interactions, that means both Coulomb and VDW.
>>> Your interpretation is correct for the group scheme. Their text implies
>>> rlist=0.8 and the others 1.4. You should also ask them if they used
>>> GROMACS
>>> and to supply a sample .mdp file so you can reproduce their reported
>>> results.
>>>
>>> 3) If rcoulomb has indeed to be set to 1.4 as rvdw, it becomes possible
>>> to
>>>
>>>> use a verlet cut-off scheme to benefit from GPU acceleration. However,
>>>>
>>> when
>>>
>>>> verlet cut-off scheme is used, mdrun shifts the nstlist and rlist to
>>>>
>>> higher
>>>
>>>> values. Isn't it likely to impact the simulation reliability in such
>>>>
>>> cases
>>>
>>>> ?
>>>>
>>>> Be careful with terms here. *Reliability* is whether the combination of
>>> code and force field correctly reproduces the results for which it was
>>> parameterized. Since you cannot use the twin-range implementation of the
>>> group scheme in the Verlet scheme, you will be using a different setup.
>>> You
>>> should absolutely evaluate reliability yourself - can you reproduce the
>>> experimental observables that the parameterization claims, and do they
>>> accompany plausible values for the observables in which you are
>>> interested?
>>> The Verlet scheme uses and interprets rlist differently (and if we had
>>> our
>>> decisions over again, we should not have been re-using rlist and nstlist
>>> for schemes with rather different behaviour).
>>>
>>> The question of which setup is more likely to produce a useful physical
>>> model for given computational cost is open, but consensus in the field
>>> suggests that the features of
>>>
>>> a) symplectic integration (which implies reversibility),
>>> b) suitable conservation of a quantity (thus adequate buffering),
>>> c) explicit handling of non-isotropic long-range effects (including at
>>> least electrostatic, but preferably also dispersion PME)
>>>
>>> are highly desirable (particulary if the system is not isotropic).
>>> Unfortunately the scheme that Poger used in their parameterization has
>>> none
>>> of those properties.
>>>
>>> Mark
>>>
>>>
>>> Thank you in advance for your replies!
>>>> --
>>>> 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.
>>>>
>>>> --
>>> 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.
>>> --
>>> 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
> Visiting Fellow
> University of Southampton, UK.
>
>
> --
> 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.
>
More information about the gromacs.org_gmx-users
mailing list