Interface qmmm

This interface implements a hybrid QM/MM scheme. An subtractive approach similar to ONIOM[1] is used, this makes it possible to interface any QM and MM code without modifications.

Dapprich S., Komaromi I., Suzie Byun K., Morokuma K., Frisch M.J., J. Mol. Struct. Teochem 1999

Important note

Two implementations of the electrostatic embedding are available: 'electrostatic_2way' and 'electrostatic'.

The electrostatic_2way embedding uses point charges on the MM atoms to polarize the QM part. Reciprocally, the gradient on the point charges is calculated in the QM code. This is the correct scheme that yields exact gradient of the QM/MM energy. However, only some QM codes can produce the gradient on point charges (following interfaces are available in cuby: DFTB+, deMon, turbomole). orca).

The electrostatic embedding yields the same energy as electrostatic_2way but the gradient on the MM atoms is only approximate, based just on the MM electrostatics. Therefore, this interface is recommended only for single-point claculations, for optimizations and MD, electrostatic_2way should be used. This was the only option for electrostatic embedding in cuby3.

(The mechanical embedding can be used with any code, support for external charges is not needed.)

Methods and capabilities

The interface implements a single unnamed method; 'method' keyword not necessary

Input structure

The interface requires following blocks in the input:

calculation_qm - Setup for the QM calculation
calculation_mm - Setup for the MM calculation

Optionally, following blocks can be defined in the input:

calculation_qmregion_mm - MM setup applied only to the QM region
calculation_system_mm - MM setup applied only to the calulation of the whole system

Keywords used

Keywords specific for this interface:

Other keywords used by this interface:

Examples

The following examples, along with all other files needed to run them, can be found in the directory cuby4/interfaces/qmmm/examples

#===============================================================================
# QMMM example 1 - simple QM/MM calculation
#===============================================================================

# Very little setup is needed when the QM region is a separate molecule. Here,
# a water dimer is calculated, the first water is QM and the second is MM.

job: optimize

geometry: water_dimer.pdb

interface: qmmm

# Electrostatic embedding:
qmmm_embedding: electrostatic_2way

# Definition of the QM region - the first residue is QM
qmmm_core: ":1"

# Setup of the QM calculation
calculation_qm:
  interface: dftb
  method: scc-dftb
  charge: 0

# Setup for the MM calculation, shared by both the calculation of the whole system
# and the MM calculation of the QM region
calculation_mm:
  interface: amber
  amber_leaprc: "%amberhome/dat/leap/cmd/leaprc.ff03.r1"

#===============================================================================
# QMMM example 2 - QMMM boundary across a covalent bond
#===============================================================================

# When a covalent bond is cut, it is capped by a hydrogen link atom. In the input,
# it is necessary to define which bond is cut and what is the position and name
# of the link atom, it is then built automatically.

# The MM forcefield has to be able to describe the QM region, in this case,
# it was necessary to crate a new residue for the alanine sidechain.

job: optimize
optimizer: rfo
opt_quality: 0.5
maxcycles: 500

geometry: ace-ala-nme.pdb

interface: qmmm

# Electrostatic embedding:
qmmm_embedding: electrostatic_2way
# Remove charges on QM atoms close to the link atom
qmmm_remove_charges: "16,17"
# The same can be achieved automatically by removing charges on atoms
# separated by N bonds (N = 1):
#qmmm_charges_around_links: 1

# Definition of the QM region
qmmm_core: 1-14
# Bonds top cut (selection, distance and PDB name of the link atom created)
qmmm_cut_bonds:
  - {bond: 9-15, link_ratio: 0.729, link_type: HL}
# Rename the QM part of the ALA residue to "ala"
qmmm_rename_residues:
  - ":2 ala"

# Setup of the QM calculation
calculation_qm:
  interface: turbomole
  method: dft
  basisset: SVP
  functional: b-lyp
  charge: 0

# Setup for the MM calculation, shared by both the calculation of the whole system
# and the MM calculation of the QM region
calculation_mm:
  interface: amber
  amber_leaprc: leaprc.ala # Customized forcefield with "ala" residue defined

#===============================================================================
# QMMM example 3 - QM/QM calculation
#===============================================================================

# The interface can be used also for QM/QM calculations, in this case DFT for
# the inner part and semiempirical QM for the rest of the system.

job: energy

geometry: ace-ala-nme.xyz

interface: qmmm

# Mechanical embedding must be used
qmmm_embedding: mechanical

# Definition of the QM region
qmmm_core: 1-14
# Bonds top cut (selection, distance and PDB name of the link atom created)
qmmm_cut_bonds:
  - {bond: 9-15, link_ratio: 0.729, link_type: HL}

# Setup of the QM calculation
calculation_qm:
  interface: turbomole
  method: dft
  basisset: SVP
  functional: b-lyp
  charge: 0

# Setup for the MM calculation, shared by both the calculation of the whole system
# and the MM calculation of the QM region
calculation_mm:
  interface: mopac
  method: pm6
  charge: 0

#===============================================================================
# QMMM example 4 - Advanced MM settings
#===============================================================================

# This example is identical to Example 2 but additional setup ui used for
# the MM calculations

# This is the same as in example 2:
job: optimize
optimizer: rfo
opt_quality: 0.5
maxcycles: 500

geometry: ace-ala-nme.pdb

interface: qmmm

qmmm_embedding: electrostatic_2way
qmmm_remove_charges: "16,17"

qmmm_core: 1-14
qmmm_cut_bonds:
  - {bond: 9-15, link_ratio: 0.729, link_type: HL}
qmmm_rename_residues:
  - ":2 ala"

calculation_qm:
  interface: turbomole
  method: dft
  basisset: SVP
  functional: b-lyp
  charge: 0

# Common MM setup used for the calculations of both the QM region and
# the whole system
calculation_mm:
  interface: amber
  # Here, we use the unmodified ff03 forcefield
  amber_leaprc: "%amberhome/dat/leap/cmd/leaprc.ff03.r1"

# Setup for the MM calculation of the QM region
calculation_qmregion_mm:
  # Here, we load the customized forcefield for the QM region,
  # it overrides the settings in block 'calculation_mm'
  amber_leaprc: leaprc.ala

# Setup for the MM calculation of the whole system
calculation_system_mm:
  # Use implicit solvent. This is a crude approximation as the
  # solvent does not affect the QM calculation in other way than mechanically
  # but it is an useful tool to keep the native structure of larger biomolecules
  solvent_model: GBM