Theor Chem Acc (2015) 134:132
1 3
where q MMH
core
is the MMH atom core charge, −3 accounts
for the two core electrons and the single valence electron
included in the SLMO, and Q
MMH
MM is the MM charge of
the MMH atom. The magnitude of q MMH
core is an adjustable
parameter that determines q bond via Eq. ( 2 ). Moreover, the
position of the bond charges is, in principle, an additional
parameter. These charges are placed at the midpoint of the
bond formed by the MMH and the MM atom adjacent to
the boundary to provide a sensible representation of the
MMH valence electrons which are not accounted for quantum mechanically.
The MM parameters at the boundary are used similarly as
they are implemented in AMBER for the link atom approach.
Internal coordinate force fi eld terms of bonds, angles,
and dihedral angles are used if any MM or MMH atom is
involved. Electrostatic and van der Waals (vdW) interactions are also calculated in accordance with the AMBER
force fi eld scheme [ 28 ]. Detailed description of the retained
bonded and nonbonded terms, in particular the interactions
with bond charges, is discussed in the following section.
3 Implementation
In this section, specifi c features of the QM/MM code
obtained with the joint implementation of the molecular mechanics program AMBER and the quantum mechanics program MRCC are presented. The coupling of the two
codes is mainly based on the interface of Walker et al. [ 23 ];
therefore, relevant parts of their work will be reviewed
with special emphasis on the similarities and differences
of the link atom and the frozen strictly localized molecular orbital approaches. In the AMBER – MRCC interface, as in
most implementations of MM and QM program couplings
[ 24 ], the necessary data exchange is realized by writing
and reading of fi les. All tasks, except the solution of the
QM orbital equations, are driven by SANDER , the molecular
dynamics engine program of AMBER .
In order to preform a QM/MM calculation, the fi rst step
for the user is to prepare a parameter and topology fi le, a
coordinate fi le, and a control fi le which are required for a
regular molecular dynamics (MD) simulation. This procedure may not be straightforward in the case of QM/MM
runs, because the QM subsystem may not be represented
by classical chemical structures and corresponding force
fi eld parameters. Nevertheless, the construction of the
topology fi le facilitates the QM/MM run since it allows the
use of the MM program infrastructure with minimal modifi cations; moreover, it may support system preparations,
for example, MM energy minimization or equilibration.
It is worth also noting that the only consequence of the
topology defi nition to the QM subsystem is the assignment
of vdW parameters. Owing to the short-range nature of the
vdW interactions, an appropriate choice of the QM subsystem can minimize the effect of potentially inconsistent
parameters.
Changes in the control fi le with respect to a standard
MM setup consist of marking the QM atoms, specifying
the charge and multiplicity of the QM region as well as
defi ning the theory and basis set. For compatibility reasons,
only the most general MRCC keywords can be denominated
in the control fi le; however, the user can also set up the
full range of options by referring in the control fi le to an
MRCC template fi le. Having done these arrangements, all
other procedures required for a simulation are automatically executed by SANDER for the link atom (LA) approach
and for the Huzinaga equation-based local self-consistent
fi eld (HLSCF) method, as well.
3.1 Link atom approach
In the case of link atom approach, the procedure for an
AMBER-MRCC calculation closely follows that formerly implemented in SANDER . Namely, it automatically
determines the place where to cut bonds and also positions the hydrogen atoms as a function of the QMH and
MMH positions to saturate the dangling bonds of the QM
subsystem. (In the context of the link atom approach the
MMH and QMH atoms are those whose bond is cut.) It is
important to note that the LA method produces a goodquality wave function when an apolar bond (e.g., bond
between two sp3 carbon atoms) is cut; moreover, link
atoms connected to a common host atom is disadvantageous owing to their short spatial separation The latter
criterion is checked by the program, and hence the user
cannot mark arbitrary QM subsystems; furthermore, both
criteria are applied in the rules of automatic selection
of model molecules for the HLSCF method (see later).
After the link atoms are defi ned, bonded terms in the QM
region are removed, and in the case of electronic embedding, the point charges are zeroed on the QM atoms as
well as, to avoid overpolarization, on the MMH atoms.
In order to preserve the total charge of the system, residual charges are distributed among the nonzero charges of
the MM region; however, there are several schemes for
both charge distribution and bonded term removal which
can be specifi ed by the user. If the MM energy terms are
defi ned, an input fi le will be written for MRCC by the
AMBER – MRCC interface in every simulation step, which
consists of QM atom coordinates, atomic numbers, MM
point charges and their coordinates, and the controlling
keywords. MRCC is then executed via a system call, and
after a successful run, the energy of the QM region and
forces acting on atoms are written into an output fi le
136
Reprinted from the journal
1 3
where q MMH
core
is the MMH atom core charge, −3 accounts
for the two core electrons and the single valence electron
included in the SLMO, and Q
MMH
MM is the MM charge of
the MMH atom. The magnitude of q MMH
core is an adjustable
parameter that determines q bond via Eq. ( 2 ). Moreover, the
position of the bond charges is, in principle, an additional
parameter. These charges are placed at the midpoint of the
bond formed by the MMH and the MM atom adjacent to
the boundary to provide a sensible representation of the
MMH valence electrons which are not accounted for quantum mechanically.
The MM parameters at the boundary are used similarly as
they are implemented in AMBER for the link atom approach.
Internal coordinate force fi eld terms of bonds, angles,
and dihedral angles are used if any MM or MMH atom is
involved. Electrostatic and van der Waals (vdW) interactions are also calculated in accordance with the AMBER
force fi eld scheme [ 28 ]. Detailed description of the retained
bonded and nonbonded terms, in particular the interactions
with bond charges, is discussed in the following section.
3 Implementation
In this section, specifi c features of the QM/MM code
obtained with the joint implementation of the molecular mechanics program AMBER and the quantum mechanics program MRCC are presented. The coupling of the two
codes is mainly based on the interface of Walker et al. [ 23 ];
therefore, relevant parts of their work will be reviewed
with special emphasis on the similarities and differences
of the link atom and the frozen strictly localized molecular orbital approaches. In the AMBER – MRCC interface, as in
most implementations of MM and QM program couplings
[ 24 ], the necessary data exchange is realized by writing
and reading of fi les. All tasks, except the solution of the
QM orbital equations, are driven by SANDER , the molecular
dynamics engine program of AMBER .
In order to preform a QM/MM calculation, the fi rst step
for the user is to prepare a parameter and topology fi le, a
coordinate fi le, and a control fi le which are required for a
regular molecular dynamics (MD) simulation. This procedure may not be straightforward in the case of QM/MM
runs, because the QM subsystem may not be represented
by classical chemical structures and corresponding force
fi eld parameters. Nevertheless, the construction of the
topology fi le facilitates the QM/MM run since it allows the
use of the MM program infrastructure with minimal modifi cations; moreover, it may support system preparations,
for example, MM energy minimization or equilibration.
It is worth also noting that the only consequence of the
topology defi nition to the QM subsystem is the assignment
of vdW parameters. Owing to the short-range nature of the
vdW interactions, an appropriate choice of the QM subsystem can minimize the effect of potentially inconsistent
parameters.
Changes in the control fi le with respect to a standard
MM setup consist of marking the QM atoms, specifying
the charge and multiplicity of the QM region as well as
defi ning the theory and basis set. For compatibility reasons,
only the most general MRCC keywords can be denominated
in the control fi le; however, the user can also set up the
full range of options by referring in the control fi le to an
MRCC template fi le. Having done these arrangements, all
other procedures required for a simulation are automatically executed by SANDER for the link atom (LA) approach
and for the Huzinaga equation-based local self-consistent
fi eld (HLSCF) method, as well.
3.1 Link atom approach
In the case of link atom approach, the procedure for an
AMBER-MRCC calculation closely follows that formerly implemented in SANDER . Namely, it automatically
determines the place where to cut bonds and also positions the hydrogen atoms as a function of the QMH and
MMH positions to saturate the dangling bonds of the QM
subsystem. (In the context of the link atom approach the
MMH and QMH atoms are those whose bond is cut.) It is
important to note that the LA method produces a goodquality wave function when an apolar bond (e.g., bond
between two sp3 carbon atoms) is cut; moreover, link
atoms connected to a common host atom is disadvantageous owing to their short spatial separation The latter
criterion is checked by the program, and hence the user
cannot mark arbitrary QM subsystems; furthermore, both
criteria are applied in the rules of automatic selection
of model molecules for the HLSCF method (see later).
After the link atoms are defi ned, bonded terms in the QM
region are removed, and in the case of electronic embedding, the point charges are zeroed on the QM atoms as
well as, to avoid overpolarization, on the MMH atoms.
In order to preserve the total charge of the system, residual charges are distributed among the nonzero charges of
the MM region; however, there are several schemes for
both charge distribution and bonded term removal which
can be specifi ed by the user. If the MM energy terms are
defi ned, an input fi le will be written for MRCC by the
AMBER – MRCC interface in every simulation step, which
consists of QM atom coordinates, atomic numbers, MM
point charges and their coordinates, and the controlling
keywords. MRCC is then executed via a system call, and
after a successful run, the energy of the QM region and
forces acting on atoms are written into an output fi le
136
Reprinted from the journal