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: Have amber on your system
You can obtain amber14 directory from someone's home directory on sinister or from /home/nn320/amber14 on nest. You can have AMBER in your path by having something similar to this in ~/.bashrc on nest.ch.private.cam.ac.uk.
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
Please note that sander does not work on nest, it will probably work on your local workstation, you need to repeat Step 0 on your local workstation. In case, sander does not run automatically, you need to first build amber. First update your $LD_LIBRARY_PATH to include $AMBERHOME/lib by including the following in ~/.bashrc
export LD_LIBRARY_PATH="${AMBERHOME}/lib:${LD_LIBRARY_PATH}"
export PYTHONPATH="${AMBERHOME}/lib/python2.7/site-packages:${PYTHONPATH}"
Then on the command line run
source ~/.bashrc
Then run the configure script with GNU as the compiler:
./configure gnu
After configuring, run
make install
to finish building amber.
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/. Please note that, you should have a look at https://wikis.ch.cam.ac.uk/ro-walesdocs/wiki/index.php/Setting_up if you do not have your peptide capped with ACE and NME groups. 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, sbatch_GMIN_serial_nest for it. You will need to copy all these files to the group directories created when you ran symm_check.sh script. To automate this copying, you can create a directory gmin_input_files and put coords.prmtop, min.in, data, sbatch_GMIN_serial_nest files in it and create a bash script copy_gmin_input_files.sh that has
#!/bin/bash
total_groups=`head -n 1 perm.allow`
cp gmin_input_files/* all_group_one_pair/
cd all_group_one_pair
sbatch sbatch_GMIN_serial_nest
cd ../
for ((i=1;i<=$total_groups;i++));
do
cp gmin_input_files/* group$i/
cd group$i
sbatch sbatch_GMIN_serial_nest
cd ../
done
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 MAXIT 1 1 MAXBFGS 0.0D0 STEP 0.0 0.0 RADIUS 1000.0 SLOPPYCONV 1.0D-6 TIGHTCONV 1.0D-6 TEMPERATURE 0.5962 AMBER12
A general sbatch_GMIN_serial_nest 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 (old_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
The atomgroups file is required when the GROUPROTATION keyword is included in the data file. To create atomgroups file, have a look at http://www-wales.ch.cam.ac.uk/examples/GMIN/1LE0/ and write it yourself. The hint is that for an amino acid residue NH2-CH-(CH2-R)-CO- the rotation axis is defined using the central C (C alpha) and the C of side chain, the group is defined using the atoms H2 and R in the group CH2-R. You may want to visualise the molecule in pymol with atom numbers to write the atomgroups file. A pymol script label_atom_numbers.pml could have
load minncrst.pdb, mol show sticks show sphere set sphere_scale, 0.3 set label_size, 10 set label_font_id, 9 set ray_trace_fog, 0 set ray_shadows, 0 unset depth_cue bg_color white set antialias, 2 set hash_max, 300 set ray_trace_mode, 3 set ray_trace_color, black set ray_opaque_background, off
When you run the following from command line on your local workstation,
pymol label_atom_numbers.pml
A pymol window appears with the molecule loaded in. On the rightmost there is object control panel with several buttons A,S,H,L,C. From the L (label button) select atom identifiers and then select ID. If you want to save the molecule as png, on the pymol command line run,
ray 2000, 1500 png raytraced_minncrst.png, dpi=300
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 have
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.