another option might be to make sure the specific definitions include
something for each periodicity (up to some reasonable number that might be
in the generics, like 4). it could use amplitude 0 if that's what it wants.
might be better than leaving it undefined and assuming that means 0.
On Fri, Nov 13, 2015 at 9:54 AM, Jason Swails <jason.swails.gmail.com>
wrote:
> I *think* I see the problem here. And you're right... it can be pretty
> nasty. I'll have to look into the best way to fix this problem.
>
> You can see a chunk of code in tleap where it *explicitly* checks to see if
> an "old" torsion was assigned as a generic torsion (with wildcards) and a
> specific torsion was found, it replaces the generic with the specific:
>
> https://github.com/choderalab/ambermini/blob/master/leap/src/leap/parmSet.c#L427-L443
>
> However, this is on a per-torsion basis. What's happening here is that the
> generic torsion defines periodicity 2, while the specific torsion assigns
> periodicity 1 and 3. So there is no term for the specific torsions that
> overwrites the periodicity-2 term that is assigned by the generic torsion.
> If you change this periodicity to 1 in gaff.dat, you should see this
> problem disappear.
>
> So the bug here *seems* to be that if a generic torsion is defined *first*
> in the parameter sets (which it appears is true for gaff.dat), and a
> specific torsion does not overwrite every periodicity that the generic
> torsion specified, then that term (or those terms) will be added to the
> specific torsion terms.
>
> There are a couple fixes I can think of:
>
> 1. Move all of the generic torsions to the *end* of the torsion list in
> gaff.dat -- this still leaves you open to this problem if you override
> generic torsions with an frcmod file, don't override all generic
> periodicities, and then load the frcmod file *after* gaff.dat is loaded.
>
> 2. After all torsions have been assigned, go through and delete the ones
> that are tagged as wild-card for 4 atoms that *also* have specific torsions
> assigned
>
> 3. Same thing as 2, except just set the force constant to the wild-card
> torsions to 0 so they don't add to the potential.
>
> #1 is the easiest, but also the worst fix. It'll get the majority of
> cases, but will still exist in workflows that are fairly standard.
>
> #2 is the "best" fix, and the way I hope things can be done, but also
> involves changing the size of an array while iterating over it. And I just
> don't know if I know the leap data structures well enough to do that with
> confidence.
>
> #3 avoids the issue of changing iterator size, but could cause confusion in
> people that don't expect to see an extra zero-term.
>
> Anyway, that's my take on it. I'll open up a bugzilla issue on this once I
> get a repro.
>
> All the best,
> Jason
>
> On Thu, Nov 12, 2015 at 8:19 PM, Niel Henriksen <shireham.gmail.com>
> wrote:
>
> > Hi all,
> >
> > I've noticed there have been many email threads over the years about wild
> > card torsions and what the proper behavior should be. Here is an example
> > <http://dev-archive.ambermd.org/200506/0024.html>. They mostly focus on
> > whether/how explicit torsion definitions should overwrite the wildcard
> > torsions. But I've stumbled on an odd case where Leap is inconsistent
> > within the same molecule. It's not mission critical at the moment, so
> this
> > is more a case study.
> >
> > We have a small molecule, parameterized with GAFF, that has four
> > carboxylate groups attached to an sp3 carbon, like this: R-CH2-COO .
> >
> > The problem occurs for the hc-c3-c -o torsion parameter, for which there
> > are sixteen instances in the molecule (four per carboxylate). In
> gaff.dat,
> > the following definitions exist:
> > X -c -c3-X 6 0.000 180.000 2.000 JCC, 7,
> (1986),
> > 230
> > hc-c3-c -o 1 0.80 0.0 -1. Junmei et al,
> > 1999
> > hc-c3-c -o 1 0.08 180.0 3. Junmei et al,
> > 1999
> >
> > I'm not 100% sure what is *supposed* to happen here, but I would think
> that
> > all instances of the torsion should be treated identically. Either all
> > sixteen instances of the torsion should get two terms (with periodicity 1
> > and 3), or they should all have three terms (with periodicity 1, 2, and
> 3).
> >
> > Unfortunately, when I "printdihedrals" from the prmtop using parmed.py, I
> > find that ten of the instances get 2 terms, and six of the instance get 3
> > terms.
> >
> > Fortunately the wildcard torsion has a barrier height of zero, so it
> > doesn't really matter in this specific case. But is it possible there
> > could be other instances like this lurking in GAFF where it does matter?
> >
> > I've attached my prmtop/inpcrds for reference. And here are the parmed
> > commands to get the torsion info:
> > printdihedrals .H23 .C37 .C38 .O10
> > printdihedrals .H24 .C37 .C38 .O10
> > printdihedrals .H23 .C37 .C38 .O32
> > printdihedrals .H24 .C37 .C38 .O32
> > printdihedrals .H31 .C43 .C44 .O12
> > printdihedrals .H32 .C43 .C44 .O12
> > printdihedrals .H31 .C43 .C44 .O25
> > printdihedrals .H32 .C43 .C44 .O25
> > printdihedrals .H27 .C40 .C41 .O11
> > printdihedrals .H28 .C40 .C41 .O11
> > printdihedrals .H27 .C40 .C41 .O26
> > printdihedrals .H28 .C40 .C41 .O26
> > printdihedrals .H19 .C34 .C35 .O9
> > printdihedrals .H20 .C34 .C35 .O9
> > printdihedrals .H19 .C34 .C35 .O24
> > printdihedrals .H20 .C34 .C35 .O24
> >
> > Thanks for your thoughts,
> > --Niel
> >
> > _______________________________________________
> > AMBER-Developers mailing list
> > AMBER-Developers.ambermd.org
> > http://lists.ambermd.org/mailman/listinfo/amber-developers
> >
> >
>
>
> --
> Jason M. Swails
> BioMaPS,
> Rutgers University
> Postdoctoral Researcher
> _______________________________________________
> AMBER-Developers mailing list
> AMBER-Developers.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber-developers
>
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Fri Nov 13 2015 - 07:00:04 PST
