Setting up
Before we can do anything, we need to set up our system.
Creating AMBER input files
Using LEaP (from scratch)
If you're going to use GMIN to find conformations to connect (for example with a small peptide) - and you don't care too much about the initial structure, you can generate it in LEaP. For the ALA-PHE peptide we're looking at here, the following LEaP input should do!
source leaprc.ff99SB
set default PBradii mbondi2
peptide = sequence {NALA CPHE}
savepdb peptide alaphe.pdb
saveamberparm peptide coords.prmtop.old coords.inpcrd
quit
This LEaP script does the following line by line:
- load the ff99SB forcefield parameters and atom types
- set the Poisson-Boltzmann radii for atom to the mbondi2 set (appropriate for use with igb=2 and igb=5 options in AMBER)
- create a peptide object with amino acid sequence ALA-PHE where ALA is at the N-terminus and PHE is at the C-terminus. The peptide is created linear with all bond lengths and angles set to their AMBER defaults.
- save a PDB file called alaphe.pdb of the peptide object so we can easily visualise the initial conformation
- save an AMBER topology file (coords.prmtop.old) and coordinate file (coords.inpcrd) for the peptide object
You can run this leap script as follows:
tleap -f tleap.in
Using PDBs from other sources
When you first get PDBs either from a database, or another source, they often will not be using the atom naming conventions AMBER will be expecting. As a result, you need to do some preparation either manually renaming atoms, or preferably - using protonate, an AMBER program.
For this tutorial, we're going to be using a folded and unfolded structure for the dipeptide ALA-PHE. Here are the initial PDBs:
alaphe_min1.pdb
REMARK NALA ATOM 1 N ALA 1 4.690 2.524 1.502 ATOM 2 1H ALA 1 5.123 1.608 1.471 ATOM 3 2H ALA 1 4.064 2.608 2.285 ATOM 4 2H ALA 1 5.478 3.159 1.619 ATOM 5 CA ALA 1 4.023 2.813 0.213 ATOM 6 HA ALA 1 3.141 2.178 0.136 ATOM 7 CB ALA 1 3.567 4.278 0.132 ATOM 8 1HB ALA 1 4.417 4.951 0.261 ATOM 9 2HB ALA 1 3.103 4.465 -0.837 ATOM 10 3HB ALA 1 2.836 4.479 0.916 ATOM 11 C ALA 1 4.949 2.387 -0.927 ATOM 12 O ALA 1 4.840 1.244 -1.353 ATOM 13 N PHE 2 5.909 3.231 -1.323 ATOM 14 H PHE 2 5.754 4.207 -1.120 ATOM 15 CA PHE 2 7.293 2.871 -0.989 ATOM 16 HA PHE 2 7.512 1.851 -1.311 ATOM 17 CB PHE 2 8.262 3.844 -1.679 ATOM 18 2HB PHE 2 8.096 3.798 -2.756 ATOM 19 3HB PHE 2 8.038 4.862 -1.355 ATOM 20 CG PHE 2 9.726 3.561 -1.399 ATOM 21 CD1 PHE 2 10.475 2.762 -2.283 ATOM 22 HD1 PHE 2 10.004 2.342 -3.160 ATOM 23 CE1 PHE 2 11.834 2.503 -2.025 ATOM 24 HE1 PHE 2 12.406 1.886 -2.704 ATOM 25 CZ PHE 2 12.447 3.041 -0.879 ATOM 26 HZ PHE 2 13.489 2.838 -0.676 ATOM 27 CE2 PHE 2 11.700 3.835 0.009 ATOM 28 HE2 PHE 2 12.165 4.240 0.896 ATOM 29 CD2 PHE 2 10.342 4.094 -0.251 ATOM 30 HD2 PHE 2 9.763 4.693 0.439 ATOM 31 C PHE 2 7.386 2.911 0.540 ATOM 32 O PHE 2 6.994 3.978 1.066 ATOM 33 OXT PHE 2 7.287 1.835 1.161 TER END
alaphe_min2.pdb
REMARK NALA ATOM 1 N ALA 1 4.769 2.964 -0.028 ATOM 2 1H ALA 1 4.610 1.986 0.154 ATOM 3 2H ALA 1 3.983 3.392 -0.491 ATOM 4 2H ALA 1 4.904 3.417 0.875 ATOM 5 CA ALA 1 5.999 3.139 -0.824 ATOM 6 HA ALA 1 6.830 2.724 -0.257 ATOM 7 CB ALA 1 5.914 2.366 -2.148 ATOM 8 1HB ALA 1 5.813 1.299 -1.945 ATOM 9 2HB ALA 1 5.058 2.705 -2.732 ATOM 10 3HB ALA 1 6.824 2.528 -2.726 ATOM 11 C ALA 1 6.286 4.634 -1.024 ATOM 12 O ALA 1 5.759 5.213 -1.967 ATOM 13 N PHE 2 7.073 5.319 -0.182 ATOM 14 H PHE 2 7.180 6.295 -0.410 ATOM 15 CA PHE 2 7.613 4.889 1.130 ATOM 16 HA PHE 2 7.722 5.785 1.741 ATOM 17 CB PHE 2 9.033 4.298 0.978 ATOM 18 2HB PHE 2 9.713 5.116 0.738 ATOM 19 3HB PHE 2 9.348 3.908 1.948 ATOM 20 CG PHE 2 9.226 3.211 -0.070 ATOM 21 CD1 PHE 2 9.116 1.852 0.281 ATOM 22 HD1 PHE 2 8.888 1.580 1.304 ATOM 23 CE1 PHE 2 9.287 0.851 -0.692 ATOM 24 HE1 PHE 2 9.205 -0.191 -0.412 ATOM 25 CZ PHE 2 9.580 1.203 -2.021 ATOM 26 HZ PHE 2 9.725 0.433 -2.766 ATOM 27 CE2 PHE 2 9.702 2.559 -2.375 ATOM 28 HE2 PHE 2 9.937 2.831 -3.395 ATOM 29 CD2 PHE 2 9.527 3.559 -1.401 ATOM 30 HD2 PHE 2 9.629 4.600 -1.674 ATOM 31 C PHE 2 6.612 4.042 1.930 ATOM 32 O PHE 2 5.489 4.562 2.089 ATOM 33 OXT PHE 2 6.759 2.807 1.993 TER END
NOTE: from now on, I will only describe the preparation process for alaphe_min1.pdb - it is IDENTICAL for min2!
protonate (with the -k flag) reads through a PDB and replaces any hydrogen atoms it does not recognise with their AMBER equivalents. If you have AMBER installed, have set the $AMBERHOME variable and have included $AMBERHOME/exe in your $PATH - you should just be able to run it like this:
protonate -k -i alaphe_min1.pdb -o alaphe_min1_amber.pdb
If not, and you have access to clust, you can use my version:
clust:~csw34/AMBER/amber9/exe/protonate
Here is the output:
alaphe_min1_amber.pdb
ATOM 1 N ALA 1 4.690 2.524 1.502 ATOM 2 H1 ALA 1 5.123 1.608 1.471 ATOM 3 H2 ALA 1 5.478 3.159 1.619 ATOM 4 H3 ALA 1 4.064 2.608 2.285 ATOM 5 CA ALA 1 4.023 2.813 0.213 ATOM 6 HA ALA 1 3.141 2.178 0.136 ATOM 7 CB ALA 1 3.567 4.278 0.132 ATOM 8 1HB ALA 1 2.836 4.479 0.916 ATOM 9 2HB ALA 1 4.417 4.951 0.261 ATOM 10 3HB ALA 1 3.103 4.465 -0.837 ATOM 11 C ALA 1 4.949 2.387 -0.927 ATOM 12 O ALA 1 4.840 1.244 -1.353 ATOM 13 N PHE 2 5.909 3.231 -1.323 ATOM 14 H PHE 2 5.754 4.207 -1.120 ATOM 15 CA PHE 2 7.293 2.871 -0.989 ATOM 16 HA PHE 2 7.512 1.851 -1.311 ATOM 17 CB PHE 2 8.262 3.844 -1.679 ATOM 18 3HB PHE 2 8.038 4.862 -1.355 ATOM 19 2HB PHE 2 8.096 3.798 -2.756 ATOM 20 CG PHE 2 9.726 3.561 -1.399 ATOM 21 CD1 PHE 2 10.475 2.762 -2.283 ATOM 22 HD1 PHE 2 10.004 2.342 -3.160 ATOM 23 CE1 PHE 2 11.834 2.503 -2.025 ATOM 24 HE1 PHE 2 12.406 1.886 -2.704 ATOM 25 CZ PHE 2 12.447 3.041 -0.879 ATOM 26 HZ PHE 2 13.489 2.838 -0.676 ATOM 27 CE2 PHE 2 11.700 3.835 0.009 ATOM 28 HE2 PHE 2 12.165 4.240 0.896 ATOM 29 CD2 PHE 2 10.342 4.094 -0.251 ATOM 30 HD2 PHE 2 9.763 4.693 0.439 ATOM 31 C PHE 2 7.386 2.911 0.540 ATOM 32 O PHE 2 6.994 3.978 1.066 ATOM 33 OXT PHE 2 7.287 1.835 1.161 END
Now we're good to create some AMBER input using LEaP.
Now that we have a PDB using the AMBER naming conventions, we can load it into LEaP to create topology and coordinate files for the forcefield we're going to be using. For this tutorial, we'll use the ff99SB forcefield and we will be aiming to use the igb=2 GB solvent model (see the AMBER9 manual!). It is important to know which solvent model we will be using in advance as it affects how we create our initial AMBER input.
Again, assuming you have AMBER installed - you can use tleap to generate the input you need. First, create a tleap.in file containing the following:
source leaprc.ff99 set default PBradii mbondi2 peptide = loadpdb alaphe_min1_amber.pdb saveamberparm peptide coords.prmtop.old coords.inpcrd quit
This LEaP script does the following line by line:
- load the ff99SB forcefield parameters and atom types
- set the Poisson-Boltzmann radii for atom to the mbondi2 set (appropriate for use with igb=2 and igb=5 options in AMBER)
- load the alaphe_min1_amber.pdb PDB into LEaP as object peptide
- save an AMBER topology file (coords.prmtop.old) and coordinate file (coords.inpcrd) for the peptide object
You can run this leap script as follows:
tleap -f tleap.in
Symmetrising the AMBER topology file
To make sure that the improper dihedral angle energies are permutationally invariant, we need to run a script to symmetrise them. Firstly though, we need to modify the coords.prmtop.old file to make sure the script is able to identify termini correctly (the improper dihedrals for terminal residues are different!). In this case, it's quite easy as there are only two residues - both of which are terminal. Open up coords.prmtop.old in your favourite editor, and find the section that looks like this:
%FLAG RESIDUE_LABEL %FORMAT(20a4) ALA PHE %FLAG RESIDUE_POINTER
The format used here is that each residue has 4 spaces in the file. To specify termini, you need to change the N-terminal ALA to NALA and the C-terminal PHE to 'CPHE', maintaining this format. Her is what it should look like:
%FLAG RESIDUE_LABEL %FORMAT(20a4) NALACPHE %FLAG RESIDUE_POINTER
If you have the SVN repository checked-out in your home directory, you can run the symmetrisation script as follows:
~/svn/SCRIPTS/AMBER/symmetrise_prmtop/perm-prmtop.ff03.py coords.prmtop.old coords.prmtop
Assuming that worked without errors, you can then open up the new coords.prmtop file and change the residue labels back from NALA to ALA and CPHE to PHE. That's it, you have a topology file ready to go. Just one more thing before we start running things - we need to make a perm.allow file for later.
Making a perm.allow file for later
The perm.allow file specifies which groups can be permuted when aligning endpoints of paths. For example, it ensures that all three H atoms of CH3 are treated as being the same, so CH3 rotation does not affect alignment. To make a perm.allow file, we need a PDB that includes the termini being specified. This means we cannot use the original PDB. We could in this case just change ALA to NALA and PHE to CPHE in the PDB manually, but you can do it automatically using the coords.prmtop.old file we made earlier which still has the termini specified. We do this using the ambpdb program:
ambpdb -p coords.prmtop.old < coords.inpcrd > alaphe_min1_amber_ter.pdb
This PDB file can then be used to generate a perm.allow file:
~/svn/SCRIPTS/make_perm.allow/perm-pdb.py alaphe_min1_amber_ter.pdb
You should now have a perm.allow file :)