Summary: TINKER vs CHARMm
Hello,
here is a summary of responses I received for my question I posted a
few weeks ago. Thank you very much for all the advice and comments
Roy
-----
S. Roy Kimura
Molecular Engineering Research Laboratory
Department of Biomedical Engineering
Boston University
email: srk # - at - # bu.edu
phone: (617)353-4697
----------------------< original question >------------------------
Dear CCLers,
Recently, I have been thinking about using the TINKER package (by Jay
Ponder) rather than CHARMm since TINKER's source looks much easier to
modify. As a check, I decided to calculate the internal energy of a 9
residue peptide in vaccuum using both the CHARMm and TINKER
programs. However, I have been getting rather different energies even
though I used the same parameter sets and identical coordinates (for
TINKER, I used the charmm.prm file supplied with the TINKER
distribution which is supposed to be equivalent to CHARMm22
parameters). Does anyone have any idea what the problem might be?
Roy
-------------------------------------------------------------------
Dear Dr. Kimura:
I recently tested TINKER on a simple MD simulation of liquid argon.
The calculated kinetic energy does not correspond to the temperature
of the simulation. There are no obvious reasons for the error in
the TINKER results).
Hubert
-------------------------------------------------------------------
Roy,
I've worked with a number of different force fields in different packages, so I
can offer a few suggestions. However, since I have not worked with either
package, I cannot comment specifically.
1) Check to see that the forcefields are defined exactly the same. The program
may say that it uses the same forcefield, but it could have been modified
slightly. This happens all the time with peptides, because everyone is
interested in getting the best result. (Ex. Improper torsional restraints may
be used in one program to hold the peptide bond flat, but this may not be the
case in the other program.)
2) Check to see that the parameters are compatible as is. In other words, check
to make sure that they don't need to be divided by 2 or something simple like
that. Also, some programs make you enter the non-bonded parameters as sigma,
epsilon, while others have you enter them a different way. You said that you
used a compatible parameter set, but is it possible that it came from a
different version of the program? Sometimes authors change these things while
they continue to improve their software.
3) Try a simple case like ethane or butane. If you are still getting different
values, then drop back to methane or similar without 1,4 interactions.
4) Check to see that dielectric constants, H-bond functions, and such are turned
on/off the same way in both programs. Alos, the non-bonded cutoffs (if used)
should be defined and set to the same value in both programs.
5) Double precision vs. Single precision. You said the coordinates were the
same, but do the packages both use double precision throughout? This would only
account for very slight differences, but it is something to keep in mind.
6) I really hate to say this, but if all else fails, trust the package that has
been around the longest.
Good luck. Trust me, I know how frustrating it came be to find these errors.
David Maxwell
Senior Scientist
Texas Biotechnology Corp.
7000 Fannin, Suite 1920
Houston, TX 77584
E-Mail: dmaxwell # - at - # tbc.com
Phone: (713) 796-8822 x102
-------------------------------------------------------------------
Dear Roy Kimura,
I just returned from some time away from St. Louis and saw your
post to CCL about "TINKER vs. CHARMm". I would be extremely
interested
in resolving (or at least rationalizing) the differences you are
seeing between energies computed with our charmm.prm parameters and
whatever CHARMm22 parameters you are using with "real" CHARMm.
There are several possibilities. First, "charmm22" means
different
things to different people, MSI put out a version of charmm22 parameters,
and Alex MacKerrell put out several versions over many years that were
all loosely called "charmm22". The charmm.prm provided with TINKER is
a
translation to TINKER format of a version of charmm22 distributed with
the related XPLOR package from Axel Brunger as of about mid-1994. At
the time we compared energies from TINKER's charmm.prm with the
corresponding XPLOR calculation, and got very close, but not
"perfect"
agreement for various polypeptide tests. The "official" paper on the
protein parameters for charmm22 has just been published (J Phys Chem,
102, 3586-3616, 1998). We certainly intend to update TINKER to this
latest version of charmm22 and validate vs. the tests reported in this
article, but we have not yet found time to do so. I would rather expect
that the differences between the new "official" parameters and our
older ones would be "small" (but not necessarily insignificant...).
If you are instead getting "large" differences between TINKER and
CHARMm, then the problem may lie elsewhere. Several people have caused
themselves trouble by mixing parameter sets from one version of TINKER
with source code and/or executables from another. This can lead to very
odd results since we have frequently made changes to both parameter
file formats and the code that are not "backward compatible".
Regardless
you should test with the most recent version of TINKER (currently 3.6).
As a relatively new package, TINKER has undergone a large number of bug
fixes as it has come into somewhat wider use. Also we have corrected
mistranslations in our charmm.prm as they have cropped up (our charmm.prm
was originally generated manually!). I know that some fairly early
versions of TINKER were floating around BU Biomed. Engineering. Of course
it is still possible that we have a "bug" in our implementation of
the
charmm function, but I tend to think we have things right at this point.
Perhaps the best way to proceed, would be for you to send me the
input and output files you are using for CHARMm (with as much detail
as possible turned on for individual energy terms in the output). I'll
take a look and see if I can't quickly find the origin of your trouble.
Obviously, I'm very interested in seeing that TINKER reproduces charmm
(and all the other force fields we support) as accurately as possible.
Jay Ponder
--------
Jay W. Ponder Phone: (314) 362-4195
Biochemistry, Box 8231 Fax: (314) 362-7183
Washington University Medical School
660 South Euclid Avenue Email: ponder # - at - # dasher.wustl.edu
St. Louis, Missouri 63110 USA WWW: http://dasher.wustl.edu/
-------------------------------------------------------------------
THIS IS MY RESPONSE TO DR. PONDER'S POST
Dear Professor Ponder,
Since I posted that message, I experimented with different "charmm22"
parameter sets, and found a set that gave charmm results very close to
those from TINKER -- that is the charmm22 parameters that I downloaded
from Axel Brunger's X-PLOR website (called top_all22_prot_b4.inp and
par_all22_prot_b4.inp).
In my initial test, I used the all atom "charmm22" parameters supplied
with MSI, which, according to some people, is known as charmm19
parameters, despite the fact that it says charmm22 in the beginning of
the file.
However, when I say very close, it means that it is much closer than
what I had been seeing when I used the MSI parameters -- but there are
some differences as you mentioned. Here are snippets from the outputs
of both programs:
......................................................................
Total Potential Energy : -40.4963 Kcal/mole
Energy Component Breakdown : Kcal/mole Interactions
Bond Stretching 18.2026 135
Angle Bending 23.3052 244
Urey-Bradley 1.8769 112
Improper Dihedral 0.1833 17
Torsional Angle 46.8472 362
1-4 van der Waals 89.1869 346
Other van der Waals -10.2909 7921
Charge-Charge -209.8075 7787
......................................................................
ENER ENR: Eval# ENERgy Delta-E GRMS
ENER INTERN: BONDs ANGLes UREY-b DIHEdrals IMPRopers
ENER EXTERN: VDWaals ELEC HBONds ASP USER
---------- --------- --------- --------- --------- ---------
ENER> 0 -36.49659 0.00000 10.32686
ENER INTERN> 16.88385 26.02777 1.56513 38.97038
0.18396
ENER EXTERN> 91.72723 -211.85492 0.00000 0.00000
0.00000
---------- --------- --------- --------- --------- ---------
......................................................................
The above are the energy calculations from a single conformation of a
9 residue peptide 2VAB (PHE ALA PRO GLY ASN TYR PRO ALA LEU). In both
cases, I used the identical coordinates (same pdb file) which I
created using charmm after minimizing the crystal structure by 50
steps. TINKER calculations were done with the latest version that I
just recently downloaded (ver 3.6). In the TINKER keyfile, I have set
the 1-4 scaling factor to 2.0 since that is what I understand charmm
does with 1-4 nonbonded interactions. It appears the largest
discrepancies lie in the torsional and Van der Waals energies (I
presume one adds your '1-4' and 'other VdW' to get the equivalent of
charmm's 'VDWaals'). Do you have any ideas as to why these
discrepancies exist?