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
Received on Fri Nov 13 2015 - 07:00:03 PST
