QMMM tutorial

Show page Syntax

In this tutorial we will prepare [QM/MM|QMMM interface] 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|http://cuby.molecular.cz/?page=QMMM+tutorial#charmm], without detailed information of preparation of the forcefield.

!! Files
Tryptofan cage miniprotein is used. PDB ready for AMBER calculation is available [here|http://cuby.molecular.cz/downloads/tutorial_qmmm/trpcage.pdb].

!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::

{html}<a name="charmm"></a>{/html}
!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}}

Password Edit summary