[gmx-users] Problem with constraints in NVT calculations.

[Justin Lemkul]

On 10/10/14 4:59 AM, Kester Wong wrote:
> Hi Justin and all,
>
>     >     > Meanwhile, is it possible to implement a self-consistent FF from scratch? One
>     >     > example I came across is from the work by Ho and Striolo
>     >     >
>     >     > titled: Polarizability effects in molecular dynamics simulations of the
>     >     > graphene-water interface
>     >     >
>     >
>     >     Of course you can implement whatever you like.  Gromacs has been able to carry
>     >     out polarizable simulations for a very long time; I've only ever cautioned
>     >     against abuse of certain models.
>     >
>     >
>     > I guess that GROMACS is capable in running polarisable sims, but for the Drude
>     > polarisable calcs, they are prone to polarisation catastrophe due to the
>     > massless shells and thermostat instability?
>
>     Polarization catastrophe is possible in any polarizable simulation.  Usually
>     very small time steps are required to avoid this, unless using an anharmonic
>     potential or a hard wall restraint.
>
>
> Using Morse = yes for the anharmonic potential option, whereas using the
> parameters below for the hard wall restraint option?
>
> pbc = xy
> nwall = 2
> wall-atomtype =    ; optional
> wall-type     = 12-6
> wall-r-linpot = 1  ; having a positive val. is esp. useful in equil. run
> wall-density  = 5 5
> wall-ewald-zfac = 3
>

No.  I'm not suggesting a Morse potential.  What I was referring to was an 
anharmonic function for the bonds, which is present in Gromacs but I'm not sure 
if it's documented.  The wall settings in Gromacs have nothing to do with this. 
  Such a function is not present in Gromacs (yet).

>
>     > In the paper mentioned above, the authors have carried out three types of cals:
>     > i) SPC/E on non-pol graphene
>     > ii) SWM4-DP on non-pol graphene: graphene in neutral or charged states
>     > iii) SWM4-DP on graphene-DP (one Drude particle per C-atom with opposite
>     > charge): graphene-DP in neutral or charged states
>     >
>     > They seemed to have simulated their systems using both additive and polarisable
>     > (0.878 angstrom^3) models?
>     > I guess this is where I got confused.
>
>     I suppose you can make any model work if you parametrize it a certain way, but
>     my point in the previous message is that you shouldn't go off trying to build a
>     force field that has SWM4-NDP water around additive CHARMM solutes.
>
>
> Yep, now I understand it.
> If I wanted to also describe graphene, is it possible to include carbon
> parameters in the SWM4-NDP force field then?
>

Well, strictly speaking, you're not introducing graphene into a SWM4-NDP force 
field, you're creating a force field that describes both.  This can certainly be 
done if you have all the parameters.

>
>     >
>     > On the side: From my previous calcs using GRAPPA force field (TIPS3P water
>     > model), graphene's polarisation (0.91 angstrom^3) resulted in spreading of water
>     > into thin layer. But that was polarisable graphene in a rigid rod model (dummy
>     > instead of shelltype particle).
>     >
>     >
>     >
>     >     >
>     >     > Pardon me if this sounds outright wrong; regarding the massless Drude particle,
>     >     > can it be replaced with an atom (assuming an induced dipole model) instead of
>     >     > the charge-on-spring model? The mass of the atom can be set to 0.4 amu with an
>     >     > opposite charge of the water oxygen atom?
>     >     >
>     >
>     >     In the Drude model with 0.4-amu particles, the Drudes are essentially just
>     >     atoms.  There's nothing conceptually special about them, we just handle them
>     >     slightly differently in the code.
>     >
>     >
>     > Well since domain decomposition will not work on shelltype calcs, I am intrigued
>     > to experiment if I can:
>     > i) replace the Drudes to atom with the same configuration - opposite charge,
>     > mass (0.4 amu), lengths, etc
>     >
>
>     The problem is that shells/Drudes have to be relaxed (SCF) or otherwise have
>     their positions integrated (extended Lagrangian) separately from "normal" atoms.
>        Conceptually, a 0.4-amu Drude is just an atom, but the integration is carried
>     out differently, so no, this sort of hacked approach probably isn't very robust.
>
> You mean the relaxation during NVT, e.g. emtol = 0.1 and niter = 30?
> As far as I know, for pure energy minimisation, the shells are treated just like
> any other particles, so what matters is the shell minimiser/integration differs
> than an "atom" during MD right?
>

Right.

> To relax a system containing Drude particles, is md=steep enough, or the more
> accurate conjugate gradient?
>
> I get unreasonable energy minimised confout.gro structure very often, if there
> are Drude-based ions included (head-scratching).
>

This suggests either unstable topology or inadequate minimizer.  Steepest 
descent should usually be good enough, if the emtol is set appropriately.

>     > OR
>     >
>     > ii) switch to the more stable SWM4-DP with the hydronium and hydroxide
>     > implementation from David van der Spoel?
>
>     I don't know how this relates to the point above about graphene, so I'm a bit
>     lost.  SWM4-NDP is a better model than SWM4-DP, FWIW.
>
>
> Absolutely, no doubt about that; SWM4-NDP describes water surf.tension better
> than SWM4-DP.
> It was just a thought, that if SWM4-NDP becomes very unstable upon the inclusion
> of polarisable ions (e.g. hydronium and hydroxide that also contain Drude
> particle), SWM4-DP could be an alternative?
>

Our group doesn't do any more work with SWM4-DP, having left it behind when 
SWM4-NDP was created.  I see no reason to use a less accurate model; I would 
instead investigate why the other species are not stable.

-Justin

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

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

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

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

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