Preparing input files for a peptide using AMBER

The steps given below are when you are trying to use one of the AMBER force fields to create topology file (coords.prmtop) and coordinates file (coords.inpcrd) for use with GMIN, OPTIM and PATHSAMPLE.

Step 0: Having amber on your system (nest.ch.private.cam.ac.uk)

Assuming you have AMBER in your path by having something similar to this in ~/.bashrc

export AMBERHOME=/home/crsid/amber14
export PATH=$PATH:$AMBERHOME/bin
export PATH

Then run,

source ~/.bashrc

in the command line

Now you can run tleap from anywhere on your system

Step 1: Make topology and coordinates file using tleap in AMBER.

leap.in file specifies force field, sequence and solvent model. My leap.in file has the following lines.

source leaprc.ff99SBildn
mol = sequence {ACE TYR TYR GLY GLY TYR TYR NME}
set default PBradii mbondi3
saveamberparm mol old_coords.prmtop old_coords.inpcrd 
savepdb mol old_mol.pdb                                                            
quit           

Run,

tleap -f leap.in

This gives pdb, prmtop and incprd files as output. Note that the pdb when visualised in vmd may have strange bonds. This just means you need to minimise the structure to get proper geometry later. Minimise the structure using sander in AMBER after you have checked that library files used for creating coords.prmtop file are correct.

Step 2: Check the amber library files

After running tleap, make a note of all the library files that get loaded. For example, in the above case, all_nucleic94.lib, all_amino94ildn.lib, all_aminoct94ildn.lib, all_aminont94ildn.lib, ions94.lib and solvents.lib get loaded. Check these library files with the ones given in softwarewales/AMBERTOOLS/dat/leap/lib The amber14/dat/leap/lib library files for ff99SBildn are probably correct. This checking is necessary, since the lib files for ff99SB i.e., all_aminoct94.lib have different charges for symmetrical atoms. Basically, they differ in the following lines, !entry.NHE.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x amber14/dat/leap/lib> 1 0 0 0 0 0 softwarewales/AMBERTOOLS/dat/leap/lib/all_aminoct94.lib> 1 1 0 0 0 0 !entry.NME.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x amber14/dat/leap/lib> 1 0 0 0 0 0 softwarewales/AMBERTOOLS/dat/leap/lib/all_aminoct94.lib> 1 3 0 0 0 0

In case the library files in amber are already correct, proceed to Step 3. If not, make a copy of lib files in amber14/dat/leap/lib/ somewhere and replace the lib files in the original directory with the ones in softwarewales/AMBERTOOLS Repeat Step 1 i.e., use tleap again and create new old_coords.prmtop and old_coords.inpcrd file before proceeding to Step 3.

Step 3: Minimising the structure using sander (sander does not work on nest, it will probably work on your local workstation, you need to repeat Step 0 on your local workstation)

This is to ensure that coordinates file has a physical structure without atom overlaps. For running sander, you just require, old_coords.prmtop, old_coords.inpcrd and min.in file. min.in file can have something like,

Minimization                                                                    
&cntrl                                                                          
  imin = 1,                                                                     
  ncyc = 1000,                                                                  
  maxcyc = 2000,                                                                
  igb=8, saltcon=0.1,                                                           
  ntb    = 0,                                                                   
  ntpr=100,                                                                     
  cut    = 999.0,                                                               
  rgbmax = 25.0                                                                 
  /  

Then on the command line run,

$AMBERHOME/bin/sander -O -i min.in -o min.out -p old_coords.prmtop -c old_coords.inpcrd -r min.ncrst

The min.ncrst has minimised geometry. To visualise min.ncrst and compare it with initial geometry run the following,

$AMBERHOME/bin/ambpdb -p old_coords.prmtop -c min.ncrst > minncrst.pdb 

To visualise old_coords.inpcrd, just look at old_mol.pdb or use

$AMBERHOME/bin/ambpdb -p old_coords.prmtop -c old_coords.inpcrd > old_coords_inpcrd.pdb 

Use the min.ncrst so obtained as your coords.inpcrd file. So, now you have obtained your old_coords.prmtop and coords.inpcrd file using AMBER.

Step 4: Symmetrise the topology file so obtained

Symmetrisation scripts are given in ~/softwarewales/SCRIPTS/AMBER/symmetrise_prmtop/ For the ff99SBildn force field use perm-prmtop.ff03.py script. It is written in python2, so first run.

module load anaconda/python2/5.3.0

Its usage is

perm-prmtop.ff03.py old_coords.prmtop old_symmetrised_coords.prmtop

You do want to check whether your topology file is symmetrised properly. Basically, symmetrisation means that when you permute the permutable atoms in your system the energy of the molecule does not change. The best way to check for correct symmetrisation is by first creating perm.allow file and then generating several coords.inpcrd files and calculating single point energy for each of them to get the same energy.

Step 5: Creating a perm.allow file The perm-pdb.py is a python2 script found in ~/softwarewales/SCRIPTS/make_perm.allow/ Run

perm-pdb.py old_mol.pdb AMBER

To check the perm.allow file simply read the documentation of PERMDIST in GMIN and check the atom numbers in perm.allow with the atom numbers using vmd or pymol and check yourself if those atom numbers correspond to permutable atoms.

Step 6: Check symmetrisation

Step 6a: Creation of several coords.inpcrd files with permuted atoms Step 6b: Run GMIN or A12GMIN in this case to check if their energies are the same.

Of course, you would want this process to be automated. The script I used can be found on sinister in /home/nn320/bin/symm_check.sh. If you have python2 module loaded already, you will have to remove it and load python3 before running symm_check.sh since it is a hybrid bash and python script. You would need A12GMIN executable, min.in, data, old_symmetrised_coords.prmtop (rename it as coords.prmtop for use with A12GMIN), perm.allow, submission_script for it. You will need to copy all these files to the group directories created by when you ran symm_check.sh script. The min.in file can have

Minimization                                                                    
 &cntrl                                                                          
   imin = 1,                                                                     
   ncyc = 1,                                                                     
   maxcyc = 1,                                                                   
   igb = 8, saltcon=0.1,                                                         
   ntb    = 0,                                                                   
   cut    = 999.0,                                                               
   rgbmax = 25.0                                                                 
 /   

Data file can have

STEPS 1 1.0 
DEBUG                                                                    
STEP 0.0 0.0                                                                    
SLOPPYCONV 1.0D-6                                                               
TIGHTCONV 1.0D-6                                                                
TEMPERATURE 0.5962                                                              
AMBER12  

A general submission_script can look like,

#!/bin/bash                                                                     
                                                                                 
#SBATCH --job-name=symm_check                                                   
#SBATCH --ntasks=1 # Specify the number of nodes you want to run on             
##SBATCH --requeue # Requeue job in the case of node failure                    
#SBATCH --mail-type=FAIL # Receive an email if your job fails                   
#SBATCH --time=5-00:00:00                                                       
                                                                                
echo "Time: `date`" > jobnumber                                                 
                                                                                
echo $SLURM_NTASKS > nodes.info                                                 
srun hostname >> nodes.info                                                     
echo $USER >> nodes.info                                                        
pwd >> nodes.info                                                               
                                                                                 
                                                                                 
/sharedscratch/crsid/softwarewales/GMIN/builds/gfortran_serial_nest/A12GMIN >> output
                                                                                 
echo Finished at `date` >> jobnumber 

You can then use

grep -i "Lowest minimum" group*/output

If all the energies are same, you can stop here, you have successfully created coords.prmtop and coords.inpcrd file for your peptide using ff99SBildn force field.

Step 7: Creating topology file for modified force field ff99IDPs

Since we want to use a modified force field ff99IDPs (https://github.com/chaohao2010/ADD-CMAP) Follow the steps given on the website You may have to change the permission of the python file before running them using ``chmod 755 ADD_CMAP.py`` Obtain another prmtop file which should already be symmetrised by following these steps. You will have to unload python2 and load python3 module.

python3 ADD-CMAP.py -p amber.prmtop -c ff99IDPs.para -o amber_CMAP.prmtop -s

or

python3 ADD_CMAP.py -p amber.prmtop -c ff99IDPs.para -o amber_CMAP.prmtop -s

Here, amber.prmtop represents the symmetrised prmtop (symmetrised_coords.prmtop) and amber_CMAP.prmtop is the new coords.prmtop for ff99IDPs force field. These files ADD_CMAP.py and ff99IDPs.para can be found on sinister in /home/nn320/ff99idps_files or on the github wesbite for ff99IDPs. You might like to check two things a) Symmetrisation of this new topology file by repeating the step 6 b) Whether the A12GMIN and AMBER energy for a structure agree with each other.

Step 8: Creating atomgroups file

To create atomgroups file, have a look at http://www-wales.ch.cam.ac.uk/examples/GMIN/1LE0/ and write it yourself

Note: Input files you should have for your peptide system to use with GMIN, OPTIM, PATHSAMPLE

coords.prmtop, coords.inpcrd, atomgroups, perm.allow, min.in, data min.in file can have

Minimization                                                                    
&cntrl                                                                          
  imin = 1,                                                                     
  ncyc = 1,                                                                     
  maxcyc = 1,                                                                   
  igb = 8, saltcon=0.1,                                                         
  ntb    = 0,                                                                   
  cut    = 999.0,                                                               
  rgbmax = 25.0                                                                 
/

Example data file can be

TEMPERATURE 0.5962                                                              
SLOPPYCONV 1.0D-4                                                               
TIGHTCONV 1.0D-7                                                                
MAXERISE 1.0D-4                                                                 
TRACKDATA                                                                       
ACCEPTRATIO 0.2                                                                 
DUMPINT 100                                                                     
UPDATES 1500                                                                    
MAXIT 3000 5000                                                                 
MAXBFGS 0.2D0                                                                   
STEPS 1 1.0                                                                     
STEP 0.0 0.0                                                                 
DEBUG                                                                           
RADIUS 1000.0                                                                   
ENERGY_DECOMP                                                                   
AMBER12

Miscellaneous

PLEASE PLEASE NOTE THAT SCEE values are 1.2 and SCNB values are 2.0 for AMBER. Check AMBER manual for more information. The topology files created using above method have 0.0 for improper torsions under SCEE and SCNB flags. The A12GMIN program should not try to invert these zeros. There was a bug in AMBER12 which has been corrected in AMBER20. So do not worry about SCEE and SCNB now.