Subject: Summary: TINKER vs CHARMm

From: "S. Roy Kimura" <[email-masked]>
Date: Thu, 6 Aug 1998 10:58:19 -0400 (EDT)

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: [email-masked]
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: [email-masked]
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:	[email-masked]
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?