Theor Chem Acc (2015) 134:132
1 3
As in the previous test, both QM/MM boundary methods well approximate the reference; moreover, the quality
of the energy curves also improves as the size of the QM
subsystem increases. In the case of cut1, the mean absolute error and the error of barrier height at 1.5 Å are both
1.1 kcal/mol for the LA method, whereas these errors are
0.6 and 0.8 kcal/mol, respectively, for the HLSCF method.
The positions of the two minima are also very well reproduced. Applying a larger QM region (cut2), both methods
approximate the reference even better: For both methods,
both the mean error and the error at the potential barrier are
0.1 kcal/mol (well below 1 %). Summarizing the results, it
was found that the both methods perform well even with
the smaller QM subsystem, and they reproduce both the
positions and the relative heights of the reference minima
extremely well with the larger QM subsystem.
5 Conclusions
A mixed QM/MM program is developed and tested
with sample calculations. The coupling of the molecular
mechanics program AMBER and the quantum mechanics
program MRCC is based on the interface of Walker et al.
[ 23 ]. Besides the traditional link atom approach, the Huzinaga equation-based local self-consistent fi eld method was
also implemented in the AMBER – MRCC code. In addition
to the implementation of the Huzinaga equation, an automatic model system selector algorithm was also developed,
which enables the automatic generation of the required
strictly localized molecular orbital coeffi cients without any
user intervention. This is an important step toward an effective and user-friendly QM/MM code since the necessity of
the generation of system specifi c frozen SLMOs was considered as an obstacle that is eliminated by the proposed
general method to derive the SLMOs.
The current implementation makes it possible to extensively test the HLSCF method, to analyze the effect of its
parameters and to compare its results with those of the link
atom and of the full QM calculations.
It was found from the deprotonation energy and from
the rotational energy profi le of the hexanoic acid that both
boundary methods monotonously converge to the full QM
calculations with increasing QM subsystem size. Calculations performed for peptides with various QM subsystem sizes also confi rm this conclusion. The deprotonation
energy was reproduced within a few tenth of kcal/mol
by both boundary methods when the subsystem boundary is separated from the site of deprotonation by fi ve
bonds. The HLSCF results were found to be sensitive to
the charge parameters at the boundary, and this sensitivity
is attributed to the calculation of the energy difference of
differently charged systems. The potential energy curves
calculated by rotating the carboxyl group of the hexanoic
acid in one example, and the imidazole group in the Ace–
His–Nme peptide in the other example showed that the
reference QM results are well reproduced by both methods with a maximal error not exceeding 5 % (1 kcal/mol
and less then 0.1 kcal/mol in the two examples) when the
rotating bond is separated from the QM–MM boundary by
three bonds. The proton transfer energy curve between an
Asp and His residues was very well reproduced considering both the positions and the energies of the minima. Even
with a subsystem boundary separated by two bonds from
the proton acceptors, the error is within 2 kcal/mol (7 %)
and it reduces signifi cantly, to less than 0.1 kcal/mol when
the QM subsystem is further extended by three bonds.
The results obtained are in line with former studies [ 10 ]
that compared the link atom and the LSCF approaches at
semiempirical level and concluded that they give comparable results when applied cautiously. However, the LSCF
approach is conceptually more appealing and its current Huzinaga equation-based implementation allows for
the fi rst time to perform calculations at the ab initio level
without orthogonalizing the basis set to the frozen orbitals.
Further developments, such as gradient calculations, inclusion of electron correlation effects from the QM side and
the treatment of periodic boundary conditions and to perform molecular dynamics simulations on the MM side will
be performed. These developments will allow a more thorough exploration of the optimal parameters of the HLSCF
method with the aim of providing a versatile QM/MM tool
to the community.
Acknowledgments This paper is dedicated to Professor Péter R.
Surján on the happy occasion of his sixtieth birthday. M.K. expresses
his gratitude to Professor Surján for mentoring him at the early stages
of his career and for continuous support. G.G.F. is grateful to Professor Peter R. Surjan for his contribution to creating a highly motivating
research atmosphere and for his friendly support the author enjoyed
as a PhD student. The computing time granted on the Hungarian
HPC Infrastructure at NIIF Institute is gratefully acknowledged. The
research work has been accomplished in the framework of the “BME
R+D+I project,” supported by the grant TÁMOP 4.2.1/B-09/1/KMR2010-0002. The authors are grateful for the fi nancial support from the
Hungarian Scientifi c Research Fund (OTKA, Grant No. K111862).
References
1. Náray-Szabó G, Surján P (1983) Chem Phys Lett 96:499–501
2. Ferenczy GG, Rivail JL, Surján PR, Náray-Szabó G (1992) J
Comput Chem 13:830–837
3. Náray-Szabó G, Surján P (1985) Theochem 123:85–95
4. Warshel A, Levitt M (1976) J Mol Biol 103:227–249
5. Théry V, Rinaldi D, Rivail JL, Maigret B, Ferenczy GG (1994) J
Comput Chem 14:269–282
6. Philipp DM, Friesner RA (1999) J Comput Chem 20:1468–1494
7. Murphy RB, Philipp DM, Friesner RA (2000) J Comput Chem
21:1442–1457
143
Reprinted from the journal
1 3
As in the previous test, both QM/MM boundary methods well approximate the reference; moreover, the quality
of the energy curves also improves as the size of the QM
subsystem increases. In the case of cut1, the mean absolute error and the error of barrier height at 1.5 Å are both
1.1 kcal/mol for the LA method, whereas these errors are
0.6 and 0.8 kcal/mol, respectively, for the HLSCF method.
The positions of the two minima are also very well reproduced. Applying a larger QM region (cut2), both methods
approximate the reference even better: For both methods,
both the mean error and the error at the potential barrier are
0.1 kcal/mol (well below 1 %). Summarizing the results, it
was found that the both methods perform well even with
the smaller QM subsystem, and they reproduce both the
positions and the relative heights of the reference minima
extremely well with the larger QM subsystem.
5 Conclusions
A mixed QM/MM program is developed and tested
with sample calculations. The coupling of the molecular
mechanics program AMBER and the quantum mechanics
program MRCC is based on the interface of Walker et al.
[ 23 ]. Besides the traditional link atom approach, the Huzinaga equation-based local self-consistent fi eld method was
also implemented in the AMBER – MRCC code. In addition
to the implementation of the Huzinaga equation, an automatic model system selector algorithm was also developed,
which enables the automatic generation of the required
strictly localized molecular orbital coeffi cients without any
user intervention. This is an important step toward an effective and user-friendly QM/MM code since the necessity of
the generation of system specifi c frozen SLMOs was considered as an obstacle that is eliminated by the proposed
general method to derive the SLMOs.
The current implementation makes it possible to extensively test the HLSCF method, to analyze the effect of its
parameters and to compare its results with those of the link
atom and of the full QM calculations.
It was found from the deprotonation energy and from
the rotational energy profi le of the hexanoic acid that both
boundary methods monotonously converge to the full QM
calculations with increasing QM subsystem size. Calculations performed for peptides with various QM subsystem sizes also confi rm this conclusion. The deprotonation
energy was reproduced within a few tenth of kcal/mol
by both boundary methods when the subsystem boundary is separated from the site of deprotonation by fi ve
bonds. The HLSCF results were found to be sensitive to
the charge parameters at the boundary, and this sensitivity
is attributed to the calculation of the energy difference of
differently charged systems. The potential energy curves
calculated by rotating the carboxyl group of the hexanoic
acid in one example, and the imidazole group in the Ace–
His–Nme peptide in the other example showed that the
reference QM results are well reproduced by both methods with a maximal error not exceeding 5 % (1 kcal/mol
and less then 0.1 kcal/mol in the two examples) when the
rotating bond is separated from the QM–MM boundary by
three bonds. The proton transfer energy curve between an
Asp and His residues was very well reproduced considering both the positions and the energies of the minima. Even
with a subsystem boundary separated by two bonds from
the proton acceptors, the error is within 2 kcal/mol (7 %)
and it reduces signifi cantly, to less than 0.1 kcal/mol when
the QM subsystem is further extended by three bonds.
The results obtained are in line with former studies [ 10 ]
that compared the link atom and the LSCF approaches at
semiempirical level and concluded that they give comparable results when applied cautiously. However, the LSCF
approach is conceptually more appealing and its current Huzinaga equation-based implementation allows for
the fi rst time to perform calculations at the ab initio level
without orthogonalizing the basis set to the frozen orbitals.
Further developments, such as gradient calculations, inclusion of electron correlation effects from the QM side and
the treatment of periodic boundary conditions and to perform molecular dynamics simulations on the MM side will
be performed. These developments will allow a more thorough exploration of the optimal parameters of the HLSCF
method with the aim of providing a versatile QM/MM tool
to the community.
Acknowledgments This paper is dedicated to Professor Péter R.
Surján on the happy occasion of his sixtieth birthday. M.K. expresses
his gratitude to Professor Surján for mentoring him at the early stages
of his career and for continuous support. G.G.F. is grateful to Professor Peter R. Surjan for his contribution to creating a highly motivating
research atmosphere and for his friendly support the author enjoyed
as a PhD student. The computing time granted on the Hungarian
HPC Infrastructure at NIIF Institute is gratefully acknowledged. The
research work has been accomplished in the framework of the “BME
R+D+I project,” supported by the grant TÁMOP 4.2.1/B-09/1/KMR2010-0002. The authors are grateful for the fi nancial support from the
Hungarian Scientifi c Research Fund (OTKA, Grant No. K111862).
References
1. Náray-Szabó G, Surján P (1983) Chem Phys Lett 96:499–501
2. Ferenczy GG, Rivail JL, Surján PR, Náray-Szabó G (1992) J
Comput Chem 13:830–837
3. Náray-Szabó G, Surján P (1985) Theochem 123:85–95
4. Warshel A, Levitt M (1976) J Mol Biol 103:227–249
5. Théry V, Rinaldi D, Rivail JL, Maigret B, Ferenczy GG (1994) J
Comput Chem 14:269–282
6. Philipp DM, Friesner RA (1999) J Comput Chem 20:1468–1494
7. Murphy RB, Philipp DM, Friesner RA (2000) J Comput Chem
21:1442–1457
143
Reprinted from the journal