I've just emerged from a long debugging session, trying to understand what
happens when the HAS_10_12 option is turned on on sander and pmemd. This was
prompted by people in California who are trying to use this option for
hydrogen-bond tweaking, and noticed that sander and pmemd give divergent
results. Here is some of what I found:
1. The current (git) version of sander won't compile if you turn this
option on. Easy to fix (I'll update the repo soon), but this does point
out the dangers of having #ifdefs in the code -- options that are rarely
(or never) tested can easily become non-functional.
2. There is indeed a bug in the HAS_10_12 part of pmemd, which looks like it
has been there at least since amber10. Basically some optimization took place
in the rest of the code that led to an variable (delrinv) sometimes being
incorrectly set in the HAS_10_12 part of the code.
[Both of these problems are evils of #ifdefs, and pmemd has a large number
of these. In the past, we've had code audits in an attempt to eliminate
un-needed compiler directives. We'll have to do that again. Of course,
there is no "right" answer here, but try to leave as few #ifdefs as
possible....]
3. Beyond this, there are inconsistencies in how 1-4 nonbonds that happen
also to be 10-12 terms are handled. Peter's original decision (probably
around 1985) was that 1-4 interactions would always use 6-12 terms, even
if a 10-12 term would otherwise be used. When Tom Darden split off
the 1-4 code into separate routines (in order to handle extra points)
he didn't think about this, and essentially wrote code that would not
be correct if 1-4 10/12 terms were in the force field. In sander, the
GB code handles 1-4 interactions in the original (Kollman) way. Pmemd
follows Darden's code for both PME and GB calculations. Since for all
current Amber force fields, the only 10-12 term is OW-HW (which is never
part of a 1-4 interaction), this problem has not be noticed before.
4. So, I need to create a bugfix for at least pmemd version 11, plus updates to
the git repo. Then test everything. But I'll try to get that done in a few
days. Note that we have at least two other energy-calculating sets of code,
(in sff and mdgx) which need to be audited for compliance with these issues.
I'm half-tempted to just say "Amber doesn't support 10-12 interactions any
more" and force the routines to exit with an error if they find any. But
people outside the usual development crew (who I can browbeat) are relying on
the older behavior, so we will probably have to live with this, and hope we
are finding all the problems.
5. We could avoid some of this by being much stricter about code requirements,
but that would inhibit some of the exciting developments that go on. But
please keep complexity in mind as you write code: if there is old stuff just
hanging around, consider getting rid of it. If you had #ifdefs in the code
just for development or testing, remove them when they are no longer needed.
If you thought you needed the #ifdef for performance reasons, check again and
see if the speed gain is really worth the extra complexity.
If you have comments or experience or suggestions with regard to the
HAS_10_12 option, please speak up.
...thx...dac
_______________________________________________
AMBER-Developers mailing list
AMBER-Developers.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber-developers
Received on Wed Oct 26 2011 - 11:00:02 PDT
