Preparing input files for a peptide using AMBER

How to create input files for a peptide for use with GMIN?

The steps given below are when you are trying to use one of the AMBER force fields.

Step 0: Having amber on your system

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

export AMBERHOME=/home/nn320/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.

Run,

tleap -f leap.in

leap.in file specifies force field, sequence, 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 coords.inpcrd 
savepdb mol mol.pdb                                                            
quit           

This gives a pdb, prmtop and incprd files as output Note that 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, all_amino94ildn, all_aminoct94ildn, all_aminont94ildn, ions94, solvents.lib Check these library files with the ones given in softwarewales/AMBERTOOLS/dat/leap/lib The amber14/dat/leap/lib library files for ff99SBildn are 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 1 0 0 0 0 0 softwareWales ambertools> 1 1 0 0 0 0 !entry.NME.unit.residueconnect table int c1x int c2x int c3x int c4x int c5x int c6x 1 0 0 0 0 0 softwareWales ambertools> 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 amber lib files somewhere and replace them with the ones in softwarewales/AMBERTOOLS Repeat Step 1 i.e., use tleap again and create new coords.prmtop and coords.inpcrd file.

Step 3: Minimising the structure using sander

This is to ensure that coordinate file has a physical structure without atom overlaps. For running sander, you just require, coords.prmtop, coords.inpcrd, 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 coords.prmtop -c coords.inpcrd -r min.ncrst

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

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

To visualise older coords.inpcrd, just use the above command replacing coordinate file and the output file.

Use the min.ncrst so obtained as your new coords.inpcrd file. So, now you have obtained your 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 above force field use perm-prmtop.py script. It is written in python2. Its usage is

perm-prmtop.py old.prmtop symmetrised.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 should not change The best way to check correct symmetrisation is by first creating perm.allow file and then generating several coords.inpcrd files and calculating single point energy of each of them to get the same energy.

Step 5: Creating a perm.allow file

Run

perm-pdb.py name.pdb AMBER

The perm-pdb.py is a python2 script found in ~/softwarewales/SCRIPTS/make_perm.allow/

To check the perm.allow file simply read the documentation of PERMDIST 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

You can probably stop here, you have 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 make the python files executable chmod 755 ADD_CMAP.py Obtain another prmtop file This should already be symmetrised The command to use is, 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 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 You might like to check two things a) Symmetrisation of this new topology file by repeating the step 6 b) Whether the A12GMIN and GMIN 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. The 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.