Basic Workshop - Hands on Session 4

(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 2008

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 provided with AmberTools): 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 "mol2" file, required to define a new unit in leap, we simply run the following command:

$AMBERHOME/exe/antechamber -i sustiva.pdb -fi pdb -o sustiva.mol2 -fo mol2 -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.mol2 specifies the name of our output file and the -fo mol2 states that we want the output file to be of Tripos Mol2 format (this is an internal format supported by Leap). The -c bcc option tells antechamber to use the AM1-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/amber10/bin/bondtype -j full -i ANTECHAMBER_BOND_TYPE.AC0 -o ANTECHAMBER_BOND_TYPE.AC -f ac

Running: /usr/local/amber10/bin/atomtype -i ANTECHAMBER_AC.AC0 -o ANTECHAMBER_AC.AC -p gaff

Total number of electrons: 160; net charge: 0

Running: /usr/local/amber10/bin/mopac.sh

Running: /usr/local/amber10/bin/am1bcc -i ANTECHAMBER_AM1BCC_PRE.AC -o ANTECHAMBER_AM1BCC.AC -f ac
-p /usr/local/amber10/dat/antechamber/BCCPARM.DAT -s 2 -j 1

Running: /usr/local/amber10/bin/atomtype -f ac -p bcc -o ANTECHAMBER_AM1BCC.AC -i ANTECHAMBER_AM1BCC_PRE.AC

ANTECHAMBER_AC.AC      ANTECHAMBER_AM1BCC_PRE.AC  ATOMTYPE.INF  mopac.out  sustiva.mol2
ANTECHAMBER_AC.AC0     ANTECHAMBER_BOND_TYPE.AC   divcon.pdb    mopac.pdb  sustiva.pdb
ANTECHAMBER_AM1BCC.AC  ANTECHAMBER_BOND_TYPE.AC0  mopac.in      SHUTDOWN

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 mopac.xxx files are input and output from the mopac 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 mopac calculation completed successfully:

mopac.out

 *******************************************************************************
 ** FRANK J. SEILER RES. LAB., U.S. AIR FORCE ACADEMY, COLO. SPGS., CO. 80840 **
 *******************************************************************************

                                 AM1 CALCULATION RESULTS

.
.
.
 CYCLE:  27 TIME:   0.19 TIME LEFT:   3597.0 GRAD.:    17.008 HEAT:-119.3928    
 CYCLE:  28 TIME:   0.16 TIME LEFT:   3596.9 GRAD.:     8.438 HEAT:-119.4094    
 TEST ON X SATISFIED
                    HOWEVER, A COMPONENT OF GRADIENT IS LARGER THAN  0.20

 CYCLE:  29 TIME:   0.21 TIME LEFT:   3596.7 GRAD.:    11.531 HEAT:-119.4101    


          HERBERTS TEST SATISFIED - GEOMETRY OPTIMIZED

 -------------------------------------------------------------------------------
 AM1 ANALYT MMOK GEO-OK PRECISE CHARGE=0
 created by wmopcrt() for mopac

     HERBERTS TEST WAS SATISFIED IN BFGS
.
.
.

The file that we are really interested in, and the reason we ran Antechamber in the first place, is the sustiva.mol2 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.mol2

@<TRIPOS>MOLECULE
SUS
   30    32     1     0     0
SMALL
bcc


@<TRIPOS>ATOM
      1 C1          0.7280    1.4030    0.2550 ca        1 SUS     -0.046700
      2 H1         -0.1360    1.7560    0.7820 ha        1 SUS      0.166500
      3 C2          0.8030    0.0800   -0.1500 ca        1 SUS     -0.179000
      4 C3         -0.2900   -0.9320    0.1580 c3        1 SUS      0.316100
      5 C4         -1.6350   -0.3520    0.0320 c1        1 SUS     -0.189800
.
.
.

As you can see this file contains 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 6). It also specifies the bonding at the end of the file. 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 the parameter 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:

sustiva.frcmod

remark goes here
MASS

BOND

ANGLE
ca-c3-c1   64.784     110.735   Calculated with empirical approach
c1-c1-cx   56.400     177.990   same as c1-c1-c3
c1-cx-hc   48.300     109.750   same as c1-c3-hc
c1-cx-cx   64.200     111.590   same as c1-c3-c3

DIHE

IMPROPER
ca-ca-ca-ha         1.1          180.0         2.0          General improper torsional angle (2 general atom types)
n -o -c -os        10.5          180.0         2.0          General improper torsional angle (2 general atom types)
c -ca-n -hn         1.1          180.0         2.0          General improper torsional angle (2 general atom types)
ca-ca-ca-n          1.1          180.0         2.0          Using default value

NONBON

We can see that there were a total of 4 missing angle parameters and 4 missing improper dihedrals. 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. If you see any parameters listed with the comment "ATTN: NEEDS REVISION" then it means that Antechamber could not determine suitable parameters and so you must manually provide these before you can proceed with the simulation. By default Antechamber will have set the values to zero.

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 FF99SB 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 FF99SB 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:

Now we can load our sustiva unit. (sustiva.mol2), we will call this new unit SUS:

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 correspond to the 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 (use the Unit menu, DO NOT click the X to close the window since this will shut down Leap entirely) 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 in the earlier hands on sessions.

¹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.