3D03
QM(SAC-CI)/MM Calculations with Gaussian: Excited States inLiquids and Enzymes(Kyoto University) O B. Swerts, J. Hasegawa, H. Nakatsuji
When molecular systems are too large to treat on a quantum mechanical (QM) level, thecombination with molecular mechanics (MM) can be a good solution. The system is thendivided into a small region of interest and treated by QM whereas the bulk is treated by MM.
This results in the so-called QM/MM methods. Numerous implementations exist, primarily
differing in the way in which broken bonds are treated. Practically all of these methods use acombination of individual QM and MM programs to achieve the desired results. This entails
adding an interface to one of these applications and passing information through external
files.1 This is not an ideal situation as both programs have to be obtained and studied how touse. It also means that most of the time two sets of input files have to be maintained.
The aim of this study was the implementation of a QM/MM method within Gaussian, awidely used QM program. It already had the ONIOM method, which can be configured as aspecial QM/MM method, but has a wider use in multilayer and combined QM/QM
applications. For that method, a number of MM force fields were added: AMBER, CFF andDREIDING. Using these, it was possible to implement a QM/MM method all within oneprogram.
All QM methods, be it semi-empirical, Hartree-Fock, or post-Hartree-Fock, can be used todescribe the QM region. The MM region can be described either by one of the three internal
force fields, by an externally interfaced CHARMM force field or any other external program.
The basic implementation uses the same ideas as the ONIOM implementation: atoms aregiven an MM atom type and added to either the QM or MM region. During geometryoptimisations, the atoms of the MM region can be given extra optimisation cycles, so-calledmicroiterations. These are handled in a different way. Depending on whether electronic
embedding is used, ONIOM either uses the MM force field as is or it uses ESP-charges asclassical charges on the QM atoms. This leads to the use of a different set of forces during thenormal optimisation and the microiterations. This is problematic when, for example, during
the microiterations the MM region is completely optimised whereas during the normal
optimisation this is not yet the case. As such convergence can never be achieved. In theQM/MM method presented, the part of the forces on the MM atoms due to the interactions
between the QM electron density and the MM charges is saved. During the microiterations
these forces are added to the MM atoms at each iteration. This only solves the problem
exactly at the first iteration and still exhibits it at subsequent iterations. This problem,however, cannot be solved exactly without recalculating the entire wavefunction at eachmicroiteration, thereby completely eliminating the advantages of using these microiterations.
The treatment of broken bonds is based on the Linkatom method. It differs from other
implementations in three aspects: the set of bonded cross-terms described by the MM force
field, the positioning of the linkatoms during geometry optimisations and the treatment ofclassical charges close to the QM region.2 These aspects were optimised by comparing
geometries, rotational barriers, electron density maps and proton affinities of a set of smallmolecules. The bonded cross-terms are those used in the GAMESS/CHARMM interface withthe omission of the MM valence angle QM-QM-MM. The linkatom is positioned by themethod used in the DYNAMO library: it is kept at a fixed distance from the nearest QM atom
and on the broken bond. As its position depends parametrically on the neighbouring QM andMM atoms, its forces can consequently be transferred to these atoms. The charges on the MMatoms close to the QM region are scaled according to a scheme that zeroes out the charge onthe nearest MM atom and redistributes this charge to the surrounding atoms while always
recovering the total system’s charge exactly. The nature of the implementation in Gaussian
allows for easily adding different Linkatom methods for a greater diversity of available
methods.
The QM/MM method was applied to the study of the excited states of liquids and enzymes.
For highly accurate results, the QM region was described by the SAC-CI method, an efficient
Coupled-Cluster based method developed specifically for the study of excited states.
A first type of applications is the solvatochromic shift of acetone and 1-naphthol in multiplesolvents. Here the solute is described by the SAC-CI method and surrounded by a 20 Å layerof solvent molecules. Although methods exist that are suited more specifically to solvation
effects, these applications are a good test for the implementation of the QM/MM method.
For acetone, the ground state geometries were determined using the SAC/D95(d):AMBER
method. Afterwards the excitation energies were determined using SAC-CI:AMBER with thecc-pVTZ basis set augmented with one Rydberg s- and p-function on each of the Carbon andOxygen atoms. For the molecule in the gas phase, an excitation energy of 4.4871 eV wasobtained, in excellent agreement with the experimental value of 4.476 ± 0.012 eV.
1-Naphthol was first optimised using the SAC/D95(d):CHARMM method and afterwards theSAC-CI:CHARMM method was used with the cc-pVDZ basis set augmented with two
Rydberg s- and p-functions on each of the Carbon and Oxygen atoms to determine theexcitation energies. The obtained excitation energy for the molecule in the gas phase was4.0915 eV, differing slightly from the experimental value of 4.381 eV.
As in other theoretical studies of the solvatochromic shift, when water was used as a solvent,explicit water molecules had to be considered in order to correctly reproduce the type of shift:
a blue-shift instead of a red-shift. The results for these calculations will be discussed andcompared with the gas-phase structures to determine the exact shift.
Enzymes pose a more challenging situation. In this case bonds have to be broken when
dividing the system into QM and MM regions so linkatoms have to be added. The opticalproperties of the Green Fluorescent Protein (GFP) are actively studied in NakatsujiLaboratory. The chromophore is the active center and is treated as the QM region in theQM/MM description. This results in a 30 atom QM region connected via two linkatoms to the3788 atoms of the surrounding MM region. The ground state geometry was obtained using
B3LYP/6-31G*:AMBER. The geometry of the lowest excited state (HOMO-LUMO) wasdetermined using CIS/6-31G*:AMBER. The resulting geometry and fluorescence energy willbe discussed in detail and compared with experimental values.
12B. Swerts, J. Vandroogenbroeck, A. Peeters, C. Van Alsenoy, J. Phys. Chem. A 106 (2002) 424
B. Swerts, Ph. D. Thesis (2004)