[gmx-users] SETTLE vs. LINCS -- different final energies of energy minimized structures

Tue Apr 19 18:07:54 CEST 2016

Hi,

On Tue, Apr 19, 2016 at 4:49 PM Rakesh Sharan <rs11chembhu at gmail.com> wrote:

> Thanks Justin and Mark!
>
> I am sorry for not being clear. Actually I am taking configurations from an
> equilibrated MD trajectory and minimizing these configuration using two
> constraints -- LINCS and SETTLE. So, basically my initial starting
> configuration is exactly same for each minimization (using SETTLE and
> LINCS), but in few cases final energy of the minimized configurations are
> different (in one case its -29146 and -29288 kj/mol) depending on the
> choice of constraint algorithm used.
>
> During energy minimization, for few configurations, in case of SETTLE I am
> getting warning "Water molecule starting at atom 1893 can not be settled"
> and in case of LINCS 4-5 warnings like "bonds that rotated more than 30
> degrees" (runs are not crashing).
>

That sounds like your actual workflow isn't what you think it is.
Equilibrated snapshots of rigid water recorded at pretty much any precision
should do EM smoothly and cleanly.

Mark, actually I am interested in the structural features of energy
> minimized (or inherent) structures. This is the reason that I want to make
> sure that I have properly energy minimized structures.
>

But "proper" doesn't have a clear meaning. If there was a single global
minimum and a smooth hypersurface then you would expect multiple
algorithms, etc. to achieve numerically the same result. But this system
has a lot of roughness, many similar minima, and there is no reason to
expect the algorithms to produce the same result, nor anything you could
deduce about the one that they found. Any single EM run is as "proper" as
you could hope for.

Mark

> Best,
> Rakesh
>
> ----------------------------------------------------
> Rakesh S. Singh, Ph.D.
> Postdoctoral Research Associate,
> Department of Chemical & Biological Engineering,
> Princeton University,
> Princeton, NJ 08544 (US).
>
>
> On Mon, Apr 18, 2016 at 12:24 PM, Mark Abraham <mark.j.abraham at gmail.com>
> wrote:
>
> > Hi,
> >
> > I would further speculate that these differences are not even
> significant.
> > You can also try rotating the whole system with gmx editconf before you
> > generate the EM input, and that will be enough to change the order of
> the
> > floating point arithmetic, so that a different final energy minimum will
> be
> > found. I don't know offhand what range of energy values you might expect,
> > but you already have some related data for that, from your multiple EM
> runs
> > from different starting points.
> >
> > I'm also sceptical about whether anything useful can be observed; the
> > potential energy surface is only somewhat relevant to the properties of
> the
> > thermalized ensemble.
> >
> > Mark
> >
> > On Mon, 18 Apr 2016 18:12 Justin Lemkul <jalemkul at vt.edu> wrote:
> >
> > >
> > >
> > > On 4/18/16 10:30 AM, Rakesh Sharan wrote:
> > > > Thanks very much Justin.
> > > >
> > > > Actually I am minimization a series of configurations from an
> > > equilibrated
> > > > trajectory of bulk TIP4P/2005 (LINCS used as constraint algorithm). I
> > > need
> > >
> > > So, to be clear, you altered the TIP4P topology to use a system of 3
> > > constraints
> > > that were held rigid by LINCS, using "constraints = all-bonds" when
> doing
> > > the
> > > simulation?  Simply setting LINCS as the constraint method in the .mdp
> > > file is
> > > insufficient to turn off SETTLE (which is always used unless you tell
> > > mdrun not to).
> > >
> > > > the final energy of the energy minimized structures, however, you can
> > see
> > > > that depending on the choice of constrained algorithm, for few
> > > > configurations, the final energy is different (e. g. for the attached
> > > plot
> > > > -29146 and -29288 kj/mol). Actually, both are following same
> > minimization
> > >
> > > The list doesn't allow attachments.
> > >
> > > Re-minimizing a configuration with different algorithms may lead to
> some
> > > differences, I would suspect.  But the internal geometry should be the
> > > same, so
> > > all that's varying will be LJ and electrostatic terms.  You can check
> by
> > > doing a
> > > further decomposition into which term is changing the most.
> > >
> > > > trajectory, however, one is finishing before other. This is bit
> > troubling
> > > > as depending on the choice of constraint algorithm average final
> energy
> > > is
> > > > different even though initial starting configuration is exactly
> same. I
> > > > rechecked topology files and they look same apart from constraint
> > > > specifications.
> > > >
> > >
> > > So let me clarify, because I'm not sure if I follow what you're saying
> > > (because
> > > here it sounds like you're talking about running different simulations,
> > > whereas
> > > above it sounds like a very different case).  Which is true:
> > >
> > > 1. You are doing two simulations from the same starting configuration
> and
> > > are
> > > trying to compare the final energy.
> > > 2. You are re-minimizing each frame of an existing trajectory and
> > comparing
> > > energies.
> > >
> > > #1 really doesn't make sense, as all MD simulations are chaotic, and
> I'm
> > > not
> > > sure what #2 would achieve.
> > >
> > > > I would very much appreciate your view over the plausible origin of
> > this
> > > > discrepancy. Also, I am wondering whether in both the cases one is
> > doing
> > > > truly constrained minimization.
> > > >
> > >
> > > If there are no bond or angle energy values, the waters were
> constrained.
> > >
> > > -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.
> > >
> > --
> > 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.
>