QMMM tutorial

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.

Input file for the same calculation using CHARMM is listed at the end of this page, without detailed information of preparation of the forcefield.

Files

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

AMBER

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 constant 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)
  - distance;1.10;680;%atomname(HL);%atomname(CB)

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

QMMM input

We are now able to run all the partial calculations, following input for QM/MM (qmmm.yml) just combines them together:

interface: qmmm

# Charge for the QM part
qmmm_qm_charge: 0

# Selection of atoms in the QM part
# Here as atom numbers (in the full system)
qmmm_core: 130-147

# List of bonds replaced by link atoms
# Link_type should be the name of the atom in
# the modified residue
qmmm_cut_bonds:
  - {bond: 128-130, link_ratio: 0.721, link_type: HL}

# List of residues to be renamed when PDB for the
# QM part is constructed
qmmm_rename_residues:
  - ":TRP trp"

# QM method specification
qmmm_method_hi:
  interface: dftb

# MM method specifiactaion
qmmm_method_lo:
  interface: amber
  leaprc: leaprc.qmmm

# Additional setup for MM calculation
# of the QM part only
qmmm_method_lo_qmpart:
  modifiers: restraints
  restraints:
    - distance;1.09;-680;%atomname(HL);%atomname(CB)
    - distance;1.10;680;%atomname(HL);%atomname(CB)

# Optimization setup
optimize_region: "%within(4;130-147)"
maxcycles: 400

An important thing is the "link_ratio" value in the definition of link atoms. It is the ratio between equilibrium bond length of the link atom and the bond it replaces. The C-C bond length is the taken from the forcefield parameter file.

This input also specified partial optimization of a region of 4 A around the QM part. It can be run by typing

cuby_optimize qmmm.yaml trpcage.pdb

Electrostatic embedding

To run QM/MM calculation with electrostatic embedding (QM calculation includes point charges for all MM atoms), it is necessary to prepare list of these charges first. MM calculation with following input (charges.yaml) should build the list:

interface: amber
leaprc: leaprc.qmmm
atomic_charges_write: charges.txt

by running

cuby_charges charges.yaml trpcage.pdb

In the QMMM input, electrostatic embedding is then specified by keywords:

polarized: yes
atomic_charges_read: charges.txt

By default, all MM atoms except the ones that are replaced by link atoms are used as point charges. More charges can be excluded using keyword qmmm_remove_charges

CHARMM version

Following input runs the same calculation with CHARMM. The PSF files have to be provided by the user.

interface: qmmm
qmmm_core: 130-147
charge: 0
qmmm_qm_charge: 0

qmmm_cut_bonds:
  - {bond: 128-130, link_ratio: 0.741, link_type: HL}
qmmm_rename_residues:
  - ":9 trp"

qmmm_method_hi:
  interface: dftb

qmmm_method_lo:
  interface: charmm
  charmm_ff_top: /home/rezac/bin/charmm35/c35b3/toppar/top_all27_prot_na.rtf
  charmm_ff_par: /home/rezac/bin/charmm35/c35b3/toppar/par_all27_prot_na.prm

qmmm_method_lo_qmpart:
  charmm_psf_file: qmmm_cluster_ok.psf
  # Sidechain terminated with link atom defined here:
  charmm_stream_files:
    - link_atoms.str
  modifiers: restraints
  restraints:
    - distance;1.111;-618;%atomname(HL);%atomname(CB)
    - distance;1.100; 618;%atomname(HL);%atomname(CB)


qmmm_method_lo_mmpart:
  charmm_psf_file: trpcage_ok.psf

optimize_region: "%within(4;130-147)"
maxcycles: 400