Theor Chem Acc (2015) 134:132
1 3
The main steps of the algorithm are demonstrated in
Fig. 2 and can be summarized as follows.
1. Repeat the model molecule search below for every
atom pairs defi ning single bonds between the QM and
MM subsystems
2. Find the base atoms and the border atoms connected to
them.
3. Extend the model molecule using the following rules
(a) If more than one border atom is connected to a
base atom, then all of them are attached to the
model molecule, unless all but one base atom–
border atom connections are terminal (the border atoms have no further bonds to heavy atoms
not included in the model system). In this case,
only the terminal atoms are added to the model
system, and the remaining one bond is cut and it
is not used in the subsequent base atom–border
atom search.
(b) If only one border atom is connected to a base
atom and
(i) this bond is a breakable one, then terminate this extension path
(ii) this bond is a non-breakable one, then add its border
atom to the model system.
(c) If a border atom is connected to two different
base atoms (crossing), then the common border
atom is also attached to the model system.
4. If every base atom has at the most one non-terminal
bond to the border atoms and this bond is breakable,
and furthermore no crossings occur, then the model
molecule is defi ned, otherwise go to step 2.
5. Add H atoms to terminate dangling bonds.
Examples indicating the selection of model molecules
for several systems are presented as Supplementary Material. The described protocol was applied to generate
model molecules in all calculation presented in the current
contribution.
4 Sample calculations
In this section results obtained both with the widely used
link atom approach and with the newly implemented Huzinaga equation-based local self-consistent fi eld method are
presented. Both QM/MM boundary methods were tested
to model typical applications of hybrid methods. However,
the systems selected are suffi ciently small so that reference
full QM calculations are feasible.
In all the calculations presented below, the wave function of the QM region was calculated with the Hartree–
Fock–Roothaan equation for the LA approach and with the
Huzinaga equation for the HLSCF approach. Pople-type
Cartesian 6-31G* basis set [ 32 ] was used. The parameters
defi ned in the ff14SB [ 28 ] and in the GAFF [ 33 ] force
fi elds were used for peptide and hexanoic acid calculations, respectively, to evaluate the MM energy function.
The MM parameters of the hexanoic acid were generated
by the ANTECHAMBER program [ 34 ] using the AM1-BCC
charge scheme [ 35 , 36 ]. Input fi les were generated by using
the Maestro [ 37 ] program and the TELEAP program which
is part of the AMBERTOOLS package. Figures were prepared
using Gnuplot [ 38 ] and Marvin [ 39 ].
As a fi rst step, the full system was optimized at the
HF/6-31G* level in all test cases, and then the obtained
geometries were modifi ed to perform single-point calculations along the defi ned reaction coordinates. The choice
of the computational level for geometry calculations was
motivated by generating reasonable geometries that are
consistent with the calculations performed to test the QM/
MM schemes. The tests were calculated in vacuum, without periodic boundary conditions and the cutoff for electrostatic and vdW terms was set to 999 Å. In the case of the
link atom approach, the capping hydrogen atom was positioned 1.09 Å away from the QMH atom; furthermore, the
residual charges, which were originated from the zeroing
of the QM and MMH atoms, were distributed among all
MM atoms. Those bonded terms which contained at least
one MM or MMH atoms were retained. The SLMOs of the
HLSCF method were only calculated for a single confi guration, and these SLMO coeffi cients were used for all other
confi gurations. This procedure was found to better reproduce reference results than the recalculation of the SLMOs
for each individual confi guration. We selected the confi guration of the global potential minimum for calculating the
SLMO coeffi cients. In this way, the SLMO coeffi cients are
treated on an equal footing with the MM parameters in the
sense that they both provide a constant environment for
the rest of the QM subsystem. The SLMOs were generated
with automatically selected model systems, using Boys
localization and with the setup described for the link atom
approach. The model systems utilized for the SLMO calculations are appended in the Supplementary information.
In all examples presented below, the effect of the QM
subsystem size is investigated by varying the subsystem
boundary. The boundary is defi ned by the bond cut by the
subsystem separation. Note that “cut” in the LA approach
indeed means bond cutting between the QM and MM
subsystems and saturating the dangling bond of the QM
139
Reprinted from the journal
1 3
The main steps of the algorithm are demonstrated in
Fig. 2 and can be summarized as follows.
1. Repeat the model molecule search below for every
atom pairs defi ning single bonds between the QM and
MM subsystems
2. Find the base atoms and the border atoms connected to
them.
3. Extend the model molecule using the following rules
(a) If more than one border atom is connected to a
base atom, then all of them are attached to the
model molecule, unless all but one base atom–
border atom connections are terminal (the border atoms have no further bonds to heavy atoms
not included in the model system). In this case,
only the terminal atoms are added to the model
system, and the remaining one bond is cut and it
is not used in the subsequent base atom–border
atom search.
(b) If only one border atom is connected to a base
atom and
(i) this bond is a breakable one, then terminate this extension path
(ii) this bond is a non-breakable one, then add its border
atom to the model system.
(c) If a border atom is connected to two different
base atoms (crossing), then the common border
atom is also attached to the model system.
4. If every base atom has at the most one non-terminal
bond to the border atoms and this bond is breakable,
and furthermore no crossings occur, then the model
molecule is defi ned, otherwise go to step 2.
5. Add H atoms to terminate dangling bonds.
Examples indicating the selection of model molecules
for several systems are presented as Supplementary Material. The described protocol was applied to generate
model molecules in all calculation presented in the current
contribution.
4 Sample calculations
In this section results obtained both with the widely used
link atom approach and with the newly implemented Huzinaga equation-based local self-consistent fi eld method are
presented. Both QM/MM boundary methods were tested
to model typical applications of hybrid methods. However,
the systems selected are suffi ciently small so that reference
full QM calculations are feasible.
In all the calculations presented below, the wave function of the QM region was calculated with the Hartree–
Fock–Roothaan equation for the LA approach and with the
Huzinaga equation for the HLSCF approach. Pople-type
Cartesian 6-31G* basis set [ 32 ] was used. The parameters
defi ned in the ff14SB [ 28 ] and in the GAFF [ 33 ] force
fi elds were used for peptide and hexanoic acid calculations, respectively, to evaluate the MM energy function.
The MM parameters of the hexanoic acid were generated
by the ANTECHAMBER program [ 34 ] using the AM1-BCC
charge scheme [ 35 , 36 ]. Input fi les were generated by using
the Maestro [ 37 ] program and the TELEAP program which
is part of the AMBERTOOLS package. Figures were prepared
using Gnuplot [ 38 ] and Marvin [ 39 ].
As a fi rst step, the full system was optimized at the
HF/6-31G* level in all test cases, and then the obtained
geometries were modifi ed to perform single-point calculations along the defi ned reaction coordinates. The choice
of the computational level for geometry calculations was
motivated by generating reasonable geometries that are
consistent with the calculations performed to test the QM/
MM schemes. The tests were calculated in vacuum, without periodic boundary conditions and the cutoff for electrostatic and vdW terms was set to 999 Å. In the case of the
link atom approach, the capping hydrogen atom was positioned 1.09 Å away from the QMH atom; furthermore, the
residual charges, which were originated from the zeroing
of the QM and MMH atoms, were distributed among all
MM atoms. Those bonded terms which contained at least
one MM or MMH atoms were retained. The SLMOs of the
HLSCF method were only calculated for a single confi guration, and these SLMO coeffi cients were used for all other
confi gurations. This procedure was found to better reproduce reference results than the recalculation of the SLMOs
for each individual confi guration. We selected the confi guration of the global potential minimum for calculating the
SLMO coeffi cients. In this way, the SLMO coeffi cients are
treated on an equal footing with the MM parameters in the
sense that they both provide a constant environment for
the rest of the QM subsystem. The SLMOs were generated
with automatically selected model systems, using Boys
localization and with the setup described for the link atom
approach. The model systems utilized for the SLMO calculations are appended in the Supplementary information.
In all examples presented below, the effect of the QM
subsystem size is investigated by varying the subsystem
boundary. The boundary is defi ned by the bond cut by the
subsystem separation. Note that “cut” in the LA approach
indeed means bond cutting between the QM and MM
subsystems and saturating the dangling bond of the QM
139
Reprinted from the journal