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
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.
Run,
tleap -f leap.in
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 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.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
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. Note that you will have to run this in a separate directory with old_coords.prmtop and old_coords.inpcrd files renamed as coords.prmtop and coords.inpcrd respectively. Softwares like AMBER and GMIN are very particular about file names and they accept coords.prmtop and coords.inpcrd as input files. 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 min.ncrst 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, i.e.
$AMBERHOME/bin/ambpdb -p coords.prmtop -c coords.inpcrd > old_inpcrd.pdb
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 ff99SBildn force field use perm-prmtop.py script. It is written in python2. Its usage is
perm-prmtop.py old_coords.prmtop 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 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
Run
perm-pdb.py mol.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 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
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 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.