QMMM tutorial

This revision is from 2010/12/02 11:07. You can Restore it.

WORK IN PROGRESS

In this tutorial we will prepare QM/MM calculation of a peptide where one amino acid (tryptofan) sidechain is calculated at QM level. AMBER is used for the MM calculation, DFTB for QM. The preparation of the MM calculation is specific for AMBER, but DFTB can be replaced with any other QM method.

Files

Tryptofan cage miniprotein is used. PDB ready for AMBER calculation is available here.

Forcefield preparation

The MM calculation of the whole peptide will run as it is. To perform the calculation on the QM region, we have to prepare residue for the sidechain.

Firstly, wee need to determine the bond length for the link atom in the QM method we will use. The sidechain can be selected and saved using the geometry script:

geometry -s ':TRP&%not(%atomname(N,H,HA,C,O))' trpcage.pdb > selection.pdb

In the file, we replace the C alpha:

ATOM    128  CA  TRP     9       0.273   2.606  -0.574

with hydrogen link atom

ATOM    128  HL  TRP     9       0.273   2.606  -0.574

For further use, we also rename the residue "TRP" in this file to "trp".

To optimize the bond length, we prepare following input (opt_link.yaml):

interface: dftb
charge: 0
freeze_atoms: "%not(1)" # Only the first atom in the geometry is optimized
restart_file: optimized_link.xyz
# default setup is used for the optimization

and run the optimization:

cuby_optimize opt_link.yaml selection.pdb

FInally, we measure the optimized bond length on the resulting structure (optimized_link.xyz); it is 1.10 A.

Forcefield preparation - AMBER specific

Now, we will prepare new residue for the sidechain using AMBER's xleap. We'll start by copying TRP to trp. In the editor, we'll delete atoms N,H,HA,C,O. Atom CA will be changed to hydrogen, named HL, atom type HC (The same as the other hydrogens on CB atom.) Finally, we'll save the new residue in OFF format.

trp = copy TRP
edit trp
saveoff trp trp.off

In the file trp.off, the residue is still named TRP. It should be changed to "trp" in text editor.

New leaprc for leap that loads both the original forcefield and the trp residue should be created, we'll name it leaprc.qmmm:

source leaprc.ff03.r1
loadoff /home/rezac/cuby3/temp/qmmm/tutorial/trp.off"

With this leaprc file, we should be able to load the file selection.pdb we prepared earlier. When it works in leap, we can try to do the MM calculation in cuby, using following input (amber.yaml):

interface: amber
leaprc: leaprc.qmmm

We can run both full system and the QM part with this:

cuby_energy amber.yaml trpcage.pdb
cuby_energy amber.yaml selection.pdb

There is still one catch. The bond length for the link atom in the forcefield is not equal to the one from our QM calculation; this has to be fixed. It is possible to modify the forcefield, but there is an easier workaround in cuby. Using the Restraints interface, we can subtract the original harmonic potential used in the forcefield and replace it with new one. The bond length and force constant can be found in the AMBER file defining the forcefield. However, the force consstant has to be multiplied by 2 for the use in cuby. Here is the input for MM calculation of the QM part:

interface: amber
leaprc: leaprc.qmmm
modifiers: restraints
restraints:
  - distance;1.09;-680;%atomname(HL);%atomname(CB) # Negative force constant, potential subtracted
  - distance;1.10;680;%atomname(HL);%atomname(CB)

# And replaced by new one

Optimization run with this input should yield bond length of 1.10 A.

Not covered here

Multiple residues in pdb - adding TERs