Theor Chem Acc (2015) 134:132
1 3
which is parsed by the interface. Here we note that LA
calculations are performed during the preparation of a
HLSCF run to generate SLMOs in model systems.
3.2 Huzinaga equation-based local self-consistent fi eld
approach
In the HLSCF method, just as in the LA approach, the communication between AMBER and MRCC is carried out by writing and reading fi les. The major difference between the two
boundary methods as described earlier is that in the case of
HLSCF, the QM and MM regions are connected by a frozen SLMO. The coeffi cients of a frozen SLMO are taken
from a molecule containing a QMH–MMH bond preferably in a chemical environment similar to that in the system
studied. Therefore, the fi rst step in order to determine the
SLMO coeffi cients of a QMH–MMH bond is to choose a
suitable model system. Model systems are generated automatically in our implementation using selection rules (see
Sect. 3.2.2 ), but the user can also select atoms for the model
systems manually using the previously mentioned template fi le. After the model systems are set, the SLMOs are
automatically calculated without any user intervention (see
Sect. 3.2.1 ) for every QMH–MMH bond, and both the outputs of model calculations and the SLMO coeffi cients are
saved in fi les. Following these procedures, bonded terms,
which are defi ned only by QM atoms are removed, thus the
following MM bonds, angles, and dihedrals are retained:
where MM and QM are MM and QM atoms, respectively,
while QMH is a QM host atom and MMH is an MM host
atom. Note that the MMH atom, which is treated specially
in the electrostatic terms is handled as any other MM atom
in the bonded terms.
The vdW and electrostatic contributions of the nonbonded terms are calculated as in the LA approach; however, only those electrostatic point charges are zeroed
which are defi ned on QM atoms, and furthermore, the
charge of the full system is preserved by distributing the
residual charges over all remaining point charges. Following this procedure, the charge of MMH atom is redistributed according to Eq. ( 2 ): the MMH point charge is
replaced by the q MMH
core core charge of the MMH atom and
the q bond charges are automatically positioned at the midpoint of the bonds formed between the MMH and the MM
atom adjacent to the MMH atom. Interactions of the q MMH
core
Bonds = MM−MM, MM−MMH, MMH−QMH
Angles = MM−MM−MM, MM−MM−MMH,
MM−MMH−QMH, MMH−QMH−QM
Dihedrals = MM−MM−MM−MM, MM−MM−MM−MMH,
MM−MM−MMH−QMH, MM−MMH−QMH−QM,
MMH−QMH−QM−QM,
and q bond charges with MM atoms are calculated as in the
AMBER force fi eld [ 28 ]: q MMH
core and q bond interact neither
with MM charges adjacent to MMH nor with MM atoms
adjacent to these atoms (1–2, 1–3 interactions). Furthermore, interactions with three bonds away from MMH are
divided by a factor of 1.2. It is important to mention that
the interactions between q MMH
core and q bond , just as any interaction including these charges are computed only by direct
summation, and hence, periodic boundary conditions cannot be handled at this stage of development.
When the setup of the MM energy terms completed,
similarly to the LA approach, an input fi le is written for
MRCC by the interface and MRCC is executed. In the case
of HLSCF, host atoms are fl agged in the input fi le. The
aim of these fl ags is to assign the previously calculated
SLMO coeffi cients to the proper host atoms and to build
the C f matrix from the strictly localized molecular orbitals. Using these orbitals the overlap matrix of the frozen
molecular orbitals are evaluated as σ = (C f ) † SC f . The
inverse square root of the orbital overlap matrix is used
to calculate the projector R f = C f (σ f ) −1 (C f ) † and then
in each iteration cycle the Fock matrix and the Huzinaga
matrix F − SR f F − FR f S are built. The diagonalization of
the Huzinaga matrix according to Eq. ( 1 ) produces orbitals
orthogonal to the frozen orbitals. When self-consistency is
reached the calculated SCF energy is written in an output
fi le which is parsed by the interface. As to the force evaluation, note that the gradients derived in Ref. [ 13 ] for the
QM atoms are implemented and tested in MRCC ; however,
the QM/MM gradients are not available at this stage of
development.
3.2.1 Automatic generation of frozen strictly localized
orbitals
The coeffi cients of the frozen SLMOs can be taken from
model molecules that contain a chemical bond similar
to that of the SLMO. The basic idea for the generation of
the orbital coeffi cients is that after identifying a bond connecting the QM and MM subsystems, a molecular graph
containing this bond is cut from the system, and a QM calculation is performed for the model molecule generated
from the molecular graph using the link atom approach.
The wave function of the model molecule is calculated
by taking into account the electrostatic potential of the
point charges of the whole system (all atoms except those
included in the model molecule). Note that during the calculation of the model molecule bonded and vdW contributions are not calculated which is rationalized by the fact that
the only aim is to determine the polarized wavefunction of
the model molecule; nevertheless, charge neutrality is still
assured by various charge distribution schemes set by the
user. The orbitals obtained are Boys [ 29 ] or Pipek–Mezey
137
Reprinted from the journal
1 3
which is parsed by the interface. Here we note that LA
calculations are performed during the preparation of a
HLSCF run to generate SLMOs in model systems.
3.2 Huzinaga equation-based local self-consistent fi eld
approach
In the HLSCF method, just as in the LA approach, the communication between AMBER and MRCC is carried out by writing and reading fi les. The major difference between the two
boundary methods as described earlier is that in the case of
HLSCF, the QM and MM regions are connected by a frozen SLMO. The coeffi cients of a frozen SLMO are taken
from a molecule containing a QMH–MMH bond preferably in a chemical environment similar to that in the system
studied. Therefore, the fi rst step in order to determine the
SLMO coeffi cients of a QMH–MMH bond is to choose a
suitable model system. Model systems are generated automatically in our implementation using selection rules (see
Sect. 3.2.2 ), but the user can also select atoms for the model
systems manually using the previously mentioned template fi le. After the model systems are set, the SLMOs are
automatically calculated without any user intervention (see
Sect. 3.2.1 ) for every QMH–MMH bond, and both the outputs of model calculations and the SLMO coeffi cients are
saved in fi les. Following these procedures, bonded terms,
which are defi ned only by QM atoms are removed, thus the
following MM bonds, angles, and dihedrals are retained:
where MM and QM are MM and QM atoms, respectively,
while QMH is a QM host atom and MMH is an MM host
atom. Note that the MMH atom, which is treated specially
in the electrostatic terms is handled as any other MM atom
in the bonded terms.
The vdW and electrostatic contributions of the nonbonded terms are calculated as in the LA approach; however, only those electrostatic point charges are zeroed
which are defi ned on QM atoms, and furthermore, the
charge of the full system is preserved by distributing the
residual charges over all remaining point charges. Following this procedure, the charge of MMH atom is redistributed according to Eq. ( 2 ): the MMH point charge is
replaced by the q MMH
core core charge of the MMH atom and
the q bond charges are automatically positioned at the midpoint of the bonds formed between the MMH and the MM
atom adjacent to the MMH atom. Interactions of the q MMH
core
Bonds = MM−MM, MM−MMH, MMH−QMH
Angles = MM−MM−MM, MM−MM−MMH,
MM−MMH−QMH, MMH−QMH−QM
Dihedrals = MM−MM−MM−MM, MM−MM−MM−MMH,
MM−MM−MMH−QMH, MM−MMH−QMH−QM,
MMH−QMH−QM−QM,
and q bond charges with MM atoms are calculated as in the
AMBER force fi eld [ 28 ]: q MMH
core and q bond interact neither
with MM charges adjacent to MMH nor with MM atoms
adjacent to these atoms (1–2, 1–3 interactions). Furthermore, interactions with three bonds away from MMH are
divided by a factor of 1.2. It is important to mention that
the interactions between q MMH
core and q bond , just as any interaction including these charges are computed only by direct
summation, and hence, periodic boundary conditions cannot be handled at this stage of development.
When the setup of the MM energy terms completed,
similarly to the LA approach, an input fi le is written for
MRCC by the interface and MRCC is executed. In the case
of HLSCF, host atoms are fl agged in the input fi le. The
aim of these fl ags is to assign the previously calculated
SLMO coeffi cients to the proper host atoms and to build
the C f matrix from the strictly localized molecular orbitals. Using these orbitals the overlap matrix of the frozen
molecular orbitals are evaluated as σ = (C f ) † SC f . The
inverse square root of the orbital overlap matrix is used
to calculate the projector R f = C f (σ f ) −1 (C f ) † and then
in each iteration cycle the Fock matrix and the Huzinaga
matrix F − SR f F − FR f S are built. The diagonalization of
the Huzinaga matrix according to Eq. ( 1 ) produces orbitals
orthogonal to the frozen orbitals. When self-consistency is
reached the calculated SCF energy is written in an output
fi le which is parsed by the interface. As to the force evaluation, note that the gradients derived in Ref. [ 13 ] for the
QM atoms are implemented and tested in MRCC ; however,
the QM/MM gradients are not available at this stage of
development.
3.2.1 Automatic generation of frozen strictly localized
orbitals
The coeffi cients of the frozen SLMOs can be taken from
model molecules that contain a chemical bond similar
to that of the SLMO. The basic idea for the generation of
the orbital coeffi cients is that after identifying a bond connecting the QM and MM subsystems, a molecular graph
containing this bond is cut from the system, and a QM calculation is performed for the model molecule generated
from the molecular graph using the link atom approach.
The wave function of the model molecule is calculated
by taking into account the electrostatic potential of the
point charges of the whole system (all atoms except those
included in the model molecule). Note that during the calculation of the model molecule bonded and vdW contributions are not calculated which is rationalized by the fact that
the only aim is to determine the polarized wavefunction of
the model molecule; nevertheless, charge neutrality is still
assured by various charge distribution schemes set by the
user. The orbitals obtained are Boys [ 29 ] or Pipek–Mezey
137
Reprinted from the journal