From srk-!at!-engc.bu.edu Thu Aug 6 10:58:19 1998 Received: from engc.bu.edu (xmw39ACBuQdV2acrwV895n5Z5WLBuzYS-: at :-ENGC.BU.EDU [128.197.185.241]) by www.ccl.net (8.8.3/8.8.6/OSC/CCL 1.0) with ESMTP id KAA19245 Thu, 6 Aug 1998 10:58:19 -0400 (EDT) Received: (from srk |-at-| localhost) by engc.bu.edu ((8.8.8.buoit.b4)/8.8.8/(BU-W-06/24/98-b7)) id KAA10963; Thu, 6 Aug 1998 10:58:19 -0400 (EDT) From: "S. Roy Kimura" Message-Id: <199808061458.KAA10963 ^at^ engc.bu.edu> Subject: Summary: TINKER vs CHARMm To: chemistry |-at-| www.ccl.net Date: Thu, 6 Aug 1998 10:58:19 -0400 (EDT) Cc: srk _-at-_)engc.bu.edu X-Mailer: ELM [version 2.4ME+ PL19 (25)] MIME-Version: 1.0 Content-Type: text/plain; charset=US-ASCII Content-Transfer-Encoding: 7bit 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?