When running QM/MM calculations, Q-Chem prints in the output three types of energies:
- E(QM + QM/MM)
- E(MM)
- E(Tot)
What is the physical meaning of these energies and how are they computed when using either ONIOM or Janus embedding?
Here is what we learned from Yihan Shao:
ONIOM (mechanical embedding):
E(Tot) = E(QM, model system in vacuum) + E(MM, total system) - E(MM, model system in vacuum)
E(QM+QM/MM) = E(QM, model system in vacuum)
That means that E(QM+QM/MM) does not have any contribution from MM. It is pure QM calculation.
E(MM) = E(MM, total system) - E(MM, model system in vacuum), where the model system refers to the selected QM region.
Electrostatic embedding (with JANUS and without Ewald)
E(QM+QM/MM) = E(QM, model system in MM point charges), that is QM energy + QM/MM electrostatics
E(MM) = E(MM, mm region) + E(QM/MM Lennard-Jones between QM and MM regions)
E(Tot) = E(QM+QM/MM) + E(MM), that is E(Tot) includes QM energy, MM energy, QM/MM electrostatics, and QM/MM Lennard-Jones energies.
1 Like
As Yihan said, what this means is different for ONIOM vs. Janus, but you can follow it pretty easily in the output. The following ONIOM sample job is provided with Q-Chem:
$comment
ONIOM (B3LYP//6-31G* and CHARMM) calculation
$end
$molecule
0 1
H 2.718389 -.562114 -.017129 1
C 1.897351 .162470 .022304 27
H 1.990860 .719423 .962720 1
H 2.042772 .873439 -.800288 1
C .537586 -.536469 -.081363 26
H .435743 -1.266310 .736293 1
H .471743 -1.105916 -1.016946 1
C -.639342 .438825 -.024882 26
H -.571490 1.156042 -.850898 1
H -.600233 1.022020 .911235 1
O -1.902602 -.197282 -.181728 76
H -2.023054 -.799333 .569479 8
$end
$rem
jobtype opt
!ideriv 0
method b3lyp
basis 6-31G*
qm_mm_interface oniom
force_field charmm27
!XC_GRID 1
!forceman_print 2
$end
$qm_atoms
8:12
$end
In the output, each step of the geometry optimization has 3 clearly labeled steps (omitting some output and just showing important parts):
**** Step 1: MM calculation on the entire system ****
Etot: 2.995247675365 Kcal/mol ( 0.004773231025 hartree)
**** Step 2: MM calculation on the model system ****
Etot: 16.488274367349 Kcal/mol ( 0.026275737864 hartree)
**** Step 3: QM calculation on the model system ****
E(QM + QM/MM): -115.7139710440 E(MM): -0.0215025 E(Tot): -115.735473550853
In the final step, the E(QM+QM/MM) is just the SCF energy, as is clear from two lines above in the output. E(MM) is the sum of the two MM energies, and E(Tot) = E(QM + QM/MM) + E(MM). This means that E(Tot) is the energy that is being optimized, and this formula is consistent with general ONIOM energy formula (see Q-Chem manual), which is
E(Tot) = E(MM,full) - E(MM,model) + E(QM,model).
Steps 1,2,3 compute these three terms, in order.
For Janus, a sample job provided with Q-Chem is the following (for the same system):
$comment
JANUS (B3LYP//6-31G* and CHARMM) calculation, CH3-CH2-CH2-OH
$end
$molecule
0 1
H 2.718389 -.562114 -.017129 1 2 0 0 0
C 1.897351 .162470 .022304 27 1 3 4 5
H 1.990860 .719423 .962720 1 2 0 0 0
H 2.042772 .873439 -.800288 1 2 0 0 0
C .537586 -.536469 -.081363 26 2 6 7 8
H .435743 -1.266310 .736293 1 5 0 0 0
H .471743 -1.105916 -1.016946 1 5 0 0 0
C -.639342 .438825 -.024882 26 5 9 10 11
H -.571490 1.156042 -.850898 1 8 0 0 0
H -.600233 1.022020 .911235 1 8 0 0 0
O -1.902602 -.197282 -.181728 76 8 12 0 0
H -2.023054 -.799333 .569479 8 11 0 0 0
$end
$rem
jobtype opt
!ideriv 0
qm_mm_interface janus
method b3lyp
basis 6-31G*
!XC_GRID 1
force_field charmm27
user_connect true
$end
$qm_atoms
8:12
$end
$nuclear_charges
5 0.820
$end
In this case there are two steps at each geometry optimization cycle:
**** Step 1: MM force contributions ****
Etot: 125.659942266614 Kcal/mol ( 0.200251865627 hartree)
**** Step 2: QM calculation on the model system ****
E(QM + QM/MM): -115.4934320940 E(MM): 0.2002519 E(Tot): -115.293180228417
As in the ONIOM job, E(QM + QM/MM) is the SCF energy. Here E(MM) is the MM energy and E(Tot) is just the sum of these two terms. As above, E(Tot) is the energy surface on which the optimization is performed.
2 Likes