Copyright (C) Pengfei Li 2018
Use GAFF in CHARMM
(needs AmberTools18 or higher)
By Pengfei Li
Like AMBER, CHARMM is another widely used molecular dynamics simulation software package. AMBER and CHARMM all have their own advantages, and some specific functions one has are not available in the other. Below is a tentative tutorial for using the general AMBER force field (GAFF) in CHARMM. This can help users to take advantage of the CHARMM functions while still use the AMBER force field. AMBER ff14SB force field for protein was already supported in the CHARMM software package. The related RTF and PAR files can be found as "/toppar/non_charmm/parm14sb_all.rtf" and "/toppar/non_charmm/parm14sb_all.prm" under the CHARMM directory path. Please copy these files into your current directory for the usage of this tutorial. Note that if you are just doing simulations of system only need ff14SB force field (no need of other force fields like the general AMBER force field - GAFF, or something else), you only need to create PSF and CRD files as followings. And then use the PSF and CRD files together with the parm14sb_all.rtf and parm14sb_all.prm files from the CHARMM software package (stated above) directly. That is, you do not need to create the RTF and PRM files again as shown below.
Unlike AMBER, CHARMM has the molecular information for molecular dynamics stored in the RTF, PAR, PSF, and CRD files. Users can check the LEaP tutorial for its comparison with AMBER about this issue. This tutorial assumes users are familiar with the LEaP program in AMBER and the CHARMM program as well. It is dealing with the same protein-ligand system in the LEaP tutorial. We will build the modeling files for CHARMM based on the files we obtained for AMBER in the LEaP tutorial. If you have not finished the LEaP tutorial, please finish it first.
Before we perform the modeling, it is important to note that AMBER and CHARMM have different naming scheme for the atom names and atom types.
First, the atom names in protein residues in AMBER and CHARMM may be different. For example, in the HID residue, names of the HB2 and HB3 atoms in AMBER will be changed to HB1 and HB2 in CHARMM.
AMBER atom types usually have two characters (except for the atomic ions which can have up to four characters) and the atom types are case sensitive. For example, for protein and nucleic acid force fields, the two letters are all in upper case (e.g. "CA"). In the general amber force field (GAFF), the two letters are all in lower case (e.g. "ca"). In sugar force field, the first letter is in upper case while the second letter is in lower case (e.g. "Ca"). In lipid force field, the first letter is in lower case while the second letter is in upper case (e.g. "cA"). CHARMM atom types are usually upper cases, but it can have up to four characters for normal elements. In the following file conversions, atom types were adapted from AMBER to CHARMM if necessary. AMBER have two upper case letters for the atom types in protein and nucleic acid residues. These atom types will be kept when transfer to CHARMM except N* and C* which will be converted to NG and CG in CHARMM. The GAFF atom types in AMBER are lower case letters, and it will be converted to "GA" + its upper case form in CHARMM. The metal ions in AMBER have atom types with the style of "element + charge", they will be changed in order to be compatible with CHARMM: "Na+" will become "SOD", "K+" will become "POT", and "Cl-" will become "CLA", etc.
All these above conversions will be handled implicitly by the python codes mentioned in this tutorial. Note that this tutorial needs AmberTools18 (or higher version) and CHARMM.
1. Generate the PSF and CRD files for the protein-ligand system.
Run the following command:
amb2chm_psf_crd.py -p 1ODX.prmtop -c 1ODX.inpcrd -f 1ODX.psf -d 1ODX.crd -b 1ODX.crd.pdb
2. Generate the PAR file
ff14SB force field contains the parm10.dat, frcmod.ff14SB files in AMBER (which can be seen from the $AMBERHOME/dat/leap/cmd/leaprc.protein.ff14SB file). These files can be found under $AMBERHOME/dat/leap/parm/ directory.
Since only atom types "CU" and "FE" do not have atom van der Waals parameters in the parm10.dat file, which will cause trouble when we combine the parameter files. We will perform the following command:
awk '$1!="CU" && $1!="FE"' $AMBERHOME/dat/leap/parm/parm10.dat > parm10_noCuFe.dat
And we will copy the frcmod.ff14SB file into the current directory:
cp $AMBERHOME/dat/leap/parm/frcmod.ff14SB .
And we can use the following file for Na+, K+, and Cl- ions: NaCl.frcmod.
We need to generate a frcmod file for the ligand, which has all the force field parameters needed for this ligand. We can run the commands:
antechamber -fi prepi -fo mol2 -i 0E8_clean_H.prepi -o 0E8.mol2 -pf y
parmchk2 -i 0E8.mol2 -o 0E8_all.frcmod -f mol2 -a Y -s 1
Afterwards we will have the amb2chm_par.py program to combine and convert (in one step) these AMBER parameter files into CHARMM. We need to have an input file to tell the amb2chm_par.py program which files will be read: par_merge.in. Then we can run the command:
amb2chm_par.py -i par_merge.in -o parmff14SB_0e8.prm
Please check the improper torsion terms in the parmff14SB_0e8.prm file. The atom type sequence inside each improper torsion needs to be consistent with the improper torsions in the AMBER dat/frcmod files used. Although relevant action was taken inside the amb2chm_par.py program, however, it is guanrateed the conversion is 100% correct. Users need to be careful by themselves to check the CHARMM output file in the benchmark calculations (shown below) to see whether any parameters were missing.
3. Generate the RTF file:
Herein we need a RTF file containing both the protein residues and the ligand. We will take advantage the RTF for ff14SB force field in CHARMM, which has information for the standard protein and nucleic acid residues.
First, we will generate a RTF file for the ligand since the parm14sb_all.rtf file does not contain information for the ligand.
We can run the command:
mol2rtf.py -i 0E8.mol2 -o 0E8.rtf -r 0E8 -n 0E8
Then we combine the two RTF files into one, at the same time we need to have the atom type list from the parmff14SB_0e8.prm file in the new RTF file.
We can do the following commands:
Note that the following line numbers in the AWK commands are only for the original parm14sb_all.rtf file in the CHARMM software package. If you have modified this file by yourself, the line number may need to be adapted.
awk 'NR<272' parm14sb_all.rtf > parm14sb_all_part1.rtf
awk 'NR>341' parm14sb_all.rtf > parm14sb_all_part3.rtf
atln=`awk '{if (substr($0, 1, 4)=="ATOM") print NR}' parmff14SB_0e8.prm`
bdln=`awk '{if (substr($0, 1, 4)=="BOND") print NR}' parmff14SB_0e8.prm`
awk 'NR>'$atln' && NR<'$bdln'' parmff14SB_0e8.prm > parm14sb_all_part2.rtf
cat parm14sb_all_part1.rtf parm14sb_all_part2.rtf parm14sb_all_part3.rtf 0E8.rtf > parmff14SB_0e8.rtf
4. Benchmark calculations
Next we will have energy and force benchmark calculations for AMBER and CHARMM. The AMBER input file can be found here: amber_bench.in. The CHARMM input file can be found here: charmm_bench.in.
We can run the commands:
sander -O -i amber.in -o amber.out -p 1ODX.prmtop -c 1ODX.inpcrd -x 1ODX.netcdf -r 1ODX.rst7
charmm < charmm.in > charmm.out
In these input files we will not use periodic boundary conditions, we will use a 999 Angstrom cut off for the nonbonded interactions. After finishing the calculations, we can check the "bench_amber.out" and "bench_charmm.out" files for comparison. Note that the energy terms in AMBER and CHARMM output files are not always in one-to-one correspondence. Due to the inner algorithms differ a little (e.g. different constants were used for the electrostatics interaction), some terms may not match exactly. For example, the constant for electrostatic interaction calculations are not exactly the same in AMBER and CHARMM. Exact convergence may not be achieved but tolerence should not be large. Meanwhile, one needs to checkout about the CHARMM output file to see whether any parameters missing in the calculations. This is because CHARMM can still work when "BOMLEV" was set high but there are missing parameters.
There is another tutorial for using AMBER force field in CHARMM about Metalloprotein. One can check it if interested.
My CHARMM expertise is basic, so I suggest you to check the following resources for assistances if need any: