Amber Workshop - tutorial 5

(Note: These tutorials are meant to provide illustrative examples of how to use the AMBER software suite to carry out simulations that can be run on a simple workstation in a reasonable period of time. They do not necessarily provide the optimal choice of parameters or methods for the particular application area.)
Copyright Ross Walker 2006

Using Antechamber to Create Leap Input Files for Simulating Sustiva
using the General Amber Force Field

In this tutorial we shall use the Antechamber tools with Leap to create topology and coordinate files for the prescription drug Sustiva (efavirenz). Sustiva (http://www.sustiva.com) is a human immunodeficiency virus type 1 (HIV-1) specific, non-nucleoside, reverse transcriptase inhibitor marketed by Bristol Myers Squibb for controlling the progression of HIV infection in humans. The chemical name for Sustiva is (S)-6-chloro-(cyclopropylethynyl)-1,4-dihydro-4-(trifluoromethyl)-2H-3,1-benzoxazin-2-one. Its empirical formula is C₁₄H₉ClF₃NO₂ and it's 2D structure is:

Here is a basic 3 dimensional geometry for Sustiva from which we will start to build our topology and coordinate files.

We shall use Antechamber to assign atom types to this molecule and also calculate a set of point charges for us. Antechamber is the most important program within the set of Antechamber tools. It can perform many file conversions and can also assign atomic charges and atom types. As required by the input antechamber executes the following programs (all provide with AMBER 8): divcon, atomtype, am1bcc, bondtype, espgen, respgen and prepgen. It will also generate a series of intermediate files (all in capital letters).

Let's try using antechamber on our sustiva pdb file. To create the "prepin" file we require to define a new unit in leap we simply run the following command:

$AMBERHOME/exe/antechamber -i sustiva.pdb -fi pdb -o sustiva.prepin -fo prepi -c bcc -s 2

Here the -i sustiva.pdb specifies the name of the 3D structure file and the -fi pdb tells antechamber that this is a pdb format file (we could easily have used any number of other supported formats including Gaussian Z-Matrix [gzmat], Gaussian Output [gout], MDL [mdl], amber Restart [rst], Sybyl Mol2 [mol2]). The -o sustiva.prepin specifies the name of our output file and the -fo prepi states that we want the output file to be of amber PREP format (this is an internal format supported by Leap). The -c bcc option tells antechamber to use the BCC charge model in order to calculate the atomic point charges while the -s 2 option defines the verbosity of the status information provided by antechamber. In this case we have selected verbose output (2).

Running: /usr/local/amber8/exe/bondtype -i ANTECHAMBER_BOND_TYPE.AC0 -o ANTECHAMBER_BOND_TYPE.AC -f ac -j full
 
 
Running: /usr/local/amber8/exe/atomtype -i ANTECHAMBER_AC.AC0 -o ANTECHAMBER_AC.AC -p gaff
 
Total number of electrons: 160; net charge: 0
 
Running: /usr/local/amber8/exe/divcon
 
Running: /usr/local/amber8/exe/am1bcc -i ANTECHAMBER_AM1BCC_PRE.AC -o ANTECHAMBER_AM1BCC.AC -f ac -p /usr/local/amber8/dat/antechamber/BCCPARM.DAT -s 2 -j 1
 
Running: /usr/local/amber8/exe/atomtype -f ac -p bcc -o ANTECHAMBER_AM1BCC.AC -i ANTECHAMBER_AM1BCC_PRE.AC
 
 
Running: /usr/local/amber8/exe/atomtype -i ANTECHAMBER_PREP.AC0 -o ANTECHAMBER_PREP.AC -p gaff
 
 
Running: /usr/local/amber8/exe/prepgen -i ANTECHAMBER_PREP.AC -f int -o sustiva.prepin -rn "SUS " -rf molecule.res

The files in CAPITALS are all intermediate files used by antechamber and are not required here. You can safely delete them. These files are not deleted by default since they may be of interest if things didn't work correctly. The divcon.xxx files are input and output from the divcon quantum mechanics code used by Antechamber to calculate the atomic point charges. We are not interested in the data here except to check that the divcon calculation completed successfully:

divcon.out

 ******************************************************************************
 *                                                                            *
 *                           DIVCON 99a                                       *
 *                                                                            *
.
.
.
  CYCLE =   431         TIME =      0.400  ENERGY =    -120.030094
  GNORM =      0.537  GRDMAX =    0.1988   GRDAVR =    0.0429
 DELTAE =   -0.000057 DELTAX = 0.000029
 ENERGY TEST PASSED
 COORDINATE TEST PASSED
 GRADIENT NORM TEST PASSED
 GRADIENT COMPONENT TEST PASSED

 -- GEOMETRY OPTIMIZED --

 FINAL QUANTITIES:
 -----------------
.
.
.

The file that we are really interested in, and the reason we ran Antechamber in the first place, is the sustiva.prepin file. This contains the definition of our sustiva residue including all of the charges and atom types that we will load into Leap to when creating our prmtop and inpcrd files. Let's take a quick look at the file:

sustiva.prepin

    0    0    2

This is a remark line
molecule.res
SUS    INT  0
CORRECT     OMIT DU   BEG
  0.0000
   1  DUMM  DU    M    0  -1  -2     0.000      .0        .0      .00000
   2  DUMM  DU    M    1   0  -1     1.449      .0        .0      .00000
   3  DUMM  DU    M    2   1   0     1.522   111.1        .0      .00000
   4  C2    ca    M    3   2   1     1.540   111.208   180.000  -0.17826
   5  C1    ca    M    4   3   2     1.386    86.897  -162.979  -0.04968
   6  H1    ha    E    5   4   3     1.072   120.119    25.599   0.16937
   7  C14   ca    M    5   4   3     1.379   119.724  -154.668  -0.01913
   8  Cl1   cl    E    7   5   4     1.742   119.676  -179.857  -0.07512
   9  C13   ca    M    7   5   4     1.385   120.634     0.368  -0.06683
  10  H9    ha    E    9   7   5     1.073   120.065   179.750   0.15651
  11  C12   ca    M    9   7   5     1.380   119.705     0.192  -0.17242
  12  H8    ha    E   11   9   7     1.076   120.047   179.373   0.14949
  13  C11   ca    M   11   9   7     1.388   119.964    -0.551   0.10869
  14  N1    n     M   13  11   9     1.388   121.155   179.940  -0.47390
  15  H7    hn    E   14  13  11     0.995   120.171     3.766   0.35663
  16  C10   c     M   14  13  11     1.362   124.462  -166.513   0.84152
  17  O2    o     E   16  14  13     1.183   123.299   172.521  -0.58069
  18  O1    os    M   16  14  13     1.338   115.711    -5.631  -0.37631
  19  C3    c3    M   18  16  14     1.416   124.597   -18.703   0.31502
  20  C9    c3    3   19  18  16     1.543   105.801   -88.841   0.61999
  21  F1    f     E   20  19  18     1.321   110.093    61.264  -0.22897
  22  F2    f     E   20  19  18     1.319   110.789  -179.348  -0.23078
  23  F3    f     E   20  19  18     1.311   111.555   -58.566  -0.21487
  24  C4    c1    M   19  18  16     1.470   107.200   154.846  -0.19720
  25  C5    c1    M   24  19  18     1.186   179.487   175.816   0.01368
  26  C6    cx    M   25  24  19     1.449   179.653  -134.955  -0.07873
  27  H2    hc    E   26  25  24     1.075   114.014    53.214   0.10860
  28  C7    cx    M   26  25  24     1.507   119.531   -91.769  -0.10741
  29  H3    hc    E   28  26  25     1.075   116.902   141.614   0.08023
  30  H4    hc    E   28  26  25     1.074   117.386    -0.494   0.08282
  31  C8    cx    M   28  26  25     1.489    60.373  -109.125  -0.11247
  32  H5    hc    E   31  28  26     1.075   118.643  -106.493   0.07951
  33  H6    hc    E   31  28  26     1.075   118.172   107.331   0.08074


LOOP
  C11   C2
   C3   C2
   C8   C6

IMPROPER
   C3  C11   C2   C1
   C2  C14   C1   H1
   C1  C13  C14  Cl1
  C12  C14  C13   H9
  C11  C13  C12   H8
  C12   C2  C11   N1
  C10  C11   N1   H7
   N1   O2  C10   O1

DONE
STOP

As you can see this file contains, in internal coordinates, the 3 dimensional structure of our sustiva molecule as well as the charge on each atom, final column, the atom number (column 1), its name (column 2) and it's atom type (column 3). It also specifies loops and improper torsions. This file does not, however, contain any parameters. The GAFF parameters are all defined in $AMBERHOME/dat/leap/parm/gaff.dat. The other thing you should notice here is that all of the GAFF atom types are in lower case. This is the mechanism by which the GAFF force field is kept independent of the macromolecular AMBER force fields. All of the traditional AMBER force fields use uppercase atom types. In this way the GAFF and traditional force fields can be mixed in the same calculation.

While the most likely combinations of bond, angle and dihedral parameters are defined in this file it is possible that our molecule might contain combinations of atom types for bonds, angles or dihedrals that have not been parameterised. If this is the case then we will have to specify any missing parameters before we can create our prmtop and inpcrd files in Leap.

We can use the utility parmchk to test if all the parameters we require are available.

Run this command now and it will produce a file called sustiva.frcmod. This is a parameter file that can be loaded into Leap in order to add missing parameters. Here it will contain all of the missing parameters. If it can antechamber will fill in these missing parameters by analogy to a similar parameter. You should check these parameters carefully before running a simulation. If antechamber can't empirically calculate a value or has no analogy it will either add a default value that it thinks is reasonable or alternatively insert a place holder (with zeros everywhere) and the comment "ATTN: needs revision". In this case you will have to manually parameterise this yourself. It is hope that as GAFF is developed so the number of missing parameters will decrease. Let's look at our frcmod file:

We can see that there were a total of 4 missing angle parameters. Later on we can take a look in Xleap and see what atoms these correspond to. For the purposes of this tutorial we shall assume that the parameters Antechamber has suggested for us are acceptable. Ideally you should really test these parameters (by comparing to ab initio calculations for example) to ensure they are reasonable.

We now have everything we need to load sustiva as a unit in Leap. We just need to load Leap and ensure the GAFF force field is available. Since we can mix the traditional AMBER force fields with GAFF we could at this point load a fragment of the HIV virus and treat this using the FF99 force field while treating the sustiva molecule using the GAFF force field. For this tutorial, however, we will simply load our GAFF sustiva and embed it in a truncated octahedral box of TIP3P water. The TIP3P water parameters are loaded as part of the FF99 Leap script so we will have Xleap load this on startup:

Once Xleap is up and running we also need to ensure that it knows about the GAFF force field. There is a script in $AMBERHOME/dat/leap/cmd/ that will do this for us. We can load it into XLeap with:

If you now type list in xleap you should see a new unit called SUS. This is the same as the unit name on line 5 of our prepin file.

At this point we haven't loaded the frcmod file that parmchk gave us. Thus if we check our SUS unit we should find that there are 4 missing angle type parameters.

Let's quickly take a look at our sustiva unit and see what atoms these correspond to:

Our missing angle type parameters were ca-c3-c1, c1-c1-cx, c1-cx-hc and c1-cx-cx. These all correspond to propyl ring and the c-c triple bond. This is what we would expect since this type of system is fairly rare in organic molecules. We can now close the edit window and load our frcmod file in order to tell Xleap the parameters for these missing angle types.

If we now check out SUS unit we should find that there are no missing parameters.

This will create a truncated octahedral box of pre-equilibrated TIP3P water around the sustiva molecule with a minimum distance between any atom of our molecule and the edge of the box of 10 angstroms.

And there you have it. At this point we could then go off and run simulations in a similar fashion to how you did it in tutorials 1 to 3.

¹Wang, J., Wolf, R.M., Caldwell, J.W., Kollman, P.A., Case, D.A. "Development and Testing of a General Amber Force Field", J. Comp. Chem., 2004, 25, 1157 - 1173.

²Jakalian, A., Bush, B.L., Jack, B.D., Bayly, C.I., "Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method.", J. Comp. Chem., 2000, 21, 132-146.