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

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.)

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

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

- 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

- atomic_charges_read
- qmmm_add_ter
- qmmm_auto_fragmentation
- qmmm_auto_run
- qmmm_charges_around_links
- qmmm_charges_extra_info
- qmmm_core
- qmmm_cut_bonds
- qmmm_embedding
- qmmm_geometry_only
- qmmm_peptide_frag_add_pro
- qmmm_qmregion_file
- qmmm_remove_charges
- qmmm_rename_residues

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
```

```
#===============================================================================
# QMMM example 5 - Automated fragmentation
#===============================================================================
# The QM/MM interface can automatically prepare QMMM calculation of a peptide,
# adding caps to the broken backbone. The initial selection of the QM region
# should contain the region of interest which is then enlarged as needed:
# 1) Complete residues included in the selection are used
# 2) No other bonds than peptide bonds in the peptide chain can be broken
# 3) Capping residues are added to terminate the broken peptide bonds
job: energy
interface: qmmm
geometry: trpcage.pdb
# Switch on automated fragmentation
qmmm_auto_fragmentation: peptide_backbone
# Here, we define what we want in the QM region. This selection is then extended
# automatically so that the peptide backbone is properly capped capped.
qmmm_core: ':8-10'
calculation_qm:
interface: mopac
method: pm6
mopac_mozyme: yes
calculation_mm:
interface: amber
# The following leaprc packaged with cuby loads the Amber ff03 with definition
# of the cap redidues needed
amber_leaprc: "%interface/data/amberff03_pm6.leaprc"
# With this input, cuby will run the calculation immediately. If you want to check
# the generated setup and QM region geometry, set the following keyword to 'no'.
qmmm_auto_run: yes
```