Durhamlogo                                                                                                                                        BSIlogo

 

Build a Molecule with a Forcefield


If you do not have a pdb version (or other format) of your molecule you will need to create on.
The easiest method to do this is with Scigress (there are alternatives e.g. Avogadro).

After drawing your molecule in Scigress, the bonds lengths and angles will unlikely be correct.
The best method to achieve the correct values is to select small areas, e.g. one ring, and run
Beautify. Keep doing this until the whole molecule has been covered.

After this we need to optimise the full geometry of the molecule. This is carried out by running
an experiment using AM1 geometry.

Save a .pdb file for the molecule, e.g. molecule.pdb

If this file is opened with jmol then the correct structure should be seen, but either no or
incorrect double bonds may be seen. Therefore we need to pass this molecule through mopac. Babel
will turn the pdb into a mopac input file.
babel -i pdb molecule.pdb -o mop molecule.dat

The first line of the molecule.dat file needs to be replaced with mopac keyworks
AM1 CHARGE=-2 ! where -2 should be replace with the correct overall charge of you molecule

Now run with mopac:
run_mopac7 molecule
run_mopac7 molecule > molecule.out

This file can then be converted into a mol2 file using babel:
babel -i moo molecule.out -o mol2 molecule.mol2

Opening this mol2 file with jmol should now give the correct structure and bonding.

Now we have the correct structure for the molecule we need to calculate the charges on each
atom using antechamber. Generally the best method for this is the sqm program with a bcc charge
calculation type.
antechamber -i molecule.mol2 -fi mol2 -o newmol.mol2 -fo mol2 -c bcc -at gaff -nc -2

The -nc command tells the program the real overall charge on your molecule so needs to be
changed to the correct value.

This new mol2 needs to be checked and any missing forcefield parameters found using parmchk:
parmchk -i newmol.mol2 -f mol2 -o newmol.frcmod


Now all the parameters should be ready to input to generate the amber forcefield.

To do this the Amber programs xleap or tleap are needed. xleap is a more graphical method, but
tleap is just as effective, and can be scripted if needed. This will generate the generate the
forcefield parameters and initial position file for both the protein and ligand together. tleap
is launched by:
tleap -f leaprc.ff99SB
where -f leaprc.ff99SB tells the program to use the 99SB forcefield parameters. Commands prefixed
with
> are ran in the tleap shell.

Firstly, load the general forcefield, followed by the missing values from parmchk:
> source leaprc.gaff
> loadAmberParams ligandH.frcmod

Next load the mol2 file with your ligand parameters in:
> lig = loadmol2 ligandH.mol2
The lig should be the same as the residue name in your full protein and ligand pdb file.

Now load the full protein and ligand pdb file with no Hydrogens, This is the first step needed
if you had no ligand (apart from launching tleap that is):
> protein = loadpdb protein.pdb

This will automatically add Hydrogens to the protein (and ligand if present). The Amber parameter
files can now be created for the protein, if it wants to be ran in a vacuum:
> saveAmberParm protein protein_vacuum.prmtop protein_vacuum.inpcrd