Using the bonded model for metal ions in CHARMM

(Needs AmberTools18 or higher)

This tutorial follows the MCPB.py tutorial. Based on the modeling files for AMBER, we can build the modeling files for CHARMM. There is a tutorial about using GAFF in CHARMM. It is recommended to finish that tutorial first before starting this one.

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.

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 (except the metal site residues) 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 CHARMM based on the prmtop and inpcrd files

Herein we need to use the amb2chm_psf_crd.py program in AmberTools18.

We can perform the command:
amb2chm_psf_crd.py -p 1OKL_solv.prmtop -c 1OKL_solv.inpcrd -f 1OKL_solv.psf -d 1OKL_solv.crd -b 1OKL_solv.crd.pdb --dict=add_res-atom_dict.txt
The generated pdb file can be used for reference. The "--dict" option was used is because of that the metal site residues can not be recognized as the amino acid residues because their specific names (Note that this file is different for different systems, you need to prepare a different file for your system). The add_res-atom_dict.txt file contains the atom information need to be changed when transfer from AMBER to CHARMM. There are four columns in the add_res-atom_dict.txt file, the first two columns are the residue name and atom name in the AMBER force field while the last two columns are the corresponding residue name and atom name in the CHARMM force field.

2. Secondly, generate the PAR file for CHARMM

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:

parmchk2 -i MNS.mol2 -o MNS_all.frcmod -f mol2 -a Y -s 1

And then we replace the atom type n2 into Y4 using the following command:

sed 's/n2/Y4/g' MNS_all.frcmod > MNS_all_new_attyp.frcmod

Note that this step was handled implicitly in MCPB.py, but we need to handle it explicitly in this tutorial.

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_1okl.prm

Please check the improper torsion terms in the parmff14SB_1okl.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. Then we need to generate RTF files for the metal site residues.

We can use the mol2rtf.py program to transfer a mol2 file to rtf file. The usage of mol2rtf.py script is as follows:

mol2rtf.py -i mol2_file -o rtf_file -r residue_name -n new_resname --ref reference

Herein, in order to generate the RTF file for metal site amino acid residues, we need the original RTF file of the ff14SB force field as a template (the "reference" term). And in the "residue_name" term we need to specify the orignial residue name of the amino acid (e.g. "HID"). And in the "new_resname" term we need to specify the new residue name of the amino acid (e.g. "HD1"). We can do the following commands to obtain the RTF files for the metal site.

mol2rtf.py -i HD1.mol2 -o HD1.rtf -r HID -n HD1 --ref parm14sb_all.rtf

mol2rtf.py -i HD2.mol2 -o HD2.rtf -r HID -n HD2 --ref parm14sb_all.rtf

mol2rtf.py -i HE1.mol2 -o HE1.rtf -r HIE -n HE1 --ref parm14sb_all.rtf

mol2rtf.py -i ZN1.mol2 -o ZN1.rtf -r ZN -n ZN1 --ref parm14sb_all.rtf

mol2rtf.py -i MS1.mol2 -o MS1.rtf -r MNS -n MS1 --ref parm14sb_all.rtf

Then we combine these RTF files into one, at the same time we need to have the atom type list from the parmff14SB_1okl.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_1okl.prm`

bdln=`awk '{if (substr($0, 1, 4)=="BOND") print NR}' parmff14SB_1okl.prm`

awk 'NR>'$atln' && NR<'$bdln'' parmff14SB_1okl.prm > parm14sb_all_part2.rtf

cat parm14sb_all_part1.rtf parm14sb_all_part2.rtf parm14sb_all_part3.rtf HD1.rtf HD2.rtf HE1.rtf MS1.rtf ZN1.rtf > parmff14SB_1okl.rtf

4. Then we can we can do a benchmark between AMBER and CHARMM.

In the benchmarks, we will perform energy and force calculations for the snapshot. Here is the input file for AMBER: bench_amber.in, in which we will calculate the total energy of the snapshot and the analytically and numerically calculated forces for certain atoms (with the atom IDs specified in the input file). Here is the input file for CHARMM: bench_charmm.in, in which we will perform calculations for the total energy and forces for each atom in the system.

Herein we can perform the AMBER calculations as follows:

sander -i bench_amber.in -o bench_amber.out -p 1OKL_solv.prmtop -c 1OKL_solv.inpcrd -r 1OKL_solv.rst -x 1OKL_solv.netcdf

And the CHARMM calculations as follows:

charmm < bench_charmm.in > bench_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 low but there are missing parameters.

My CHARMM expertise is basic, so I suggest you to check the following resources for assistances if need any:

The official CHARMM website

The CHARMM documentation webpage

CHARMM tutorials

The CHARMM Forum