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.