C7790 Introduction to Molecular Modelling -1C7790
Introduction to Molecular
Modelling
TSM Modelling Molecular Structures
Petr Kulhánek
kulhanek@chemi.muni.cz
National Centre for Biomolecular Research, Faculty of Science
Masaryk University, Kamenice 5, CZ-62500 Brno
PS/2020 Distant Form of Teaching: Rev1
Lesson 14
Potential Energy Surface II
C7790 Introduction to Molecular Modelling -2-
Context
microworldmacroworld
equilibrium (equilibrium constant)
kinetics (rate constant)
free energy
(Gibbs/Helmholtz)
partition function
phenomenological thermodynamics
statistical thermodynamics
microstates
(mechanical properties, E)
states
(thermodynamic properties, G, T,…)
microstate ≠ microworld
Description levels (model chemistry):
• quantum mechanics
• semiempirical methods
• ab initio methods
• post-HF methods
• DFT methods
• molecular mechanics
• coarse-grained mechanics
Structure
EnergyFunction
Simulations:
• molecular dynamics
• Monte Carlo simulations
• docking
• …
C7790 Introduction to Molecular Modelling -3-
PES
Potential Energy Surface
Properties
Visualization
Important Points (Stationary States)
C7790 Introduction to Molecular Modelling -4Configuration
space
)(RE R = point in 3N dimensional space (N is the number of atoms)
),,,....,,,,,,( 222111 NNN zyxzyxzyx=R
Cartesian coordinates
of the first nucleus (atom, bead).
The individual points form a configuration space.
Every point in the configuration space then
represents unique structure of the system.
individual coordinates
are degrees of freedom
of the system
C7790 Introduction to Molecular Modelling -5Potential
energy calculation
The calculation of the potential energy E(R) is possible by:
▪ approximate solution of the Schrödinger equation (quantum mechanics, QM)
▪ using empirical force fields (molecular mechanics, MM)
▪ hybrid QM/MM approach
▪ using coarse grained models
C7790 Introduction to Molecular Modelling -6Potential
energy calculation
The calculation of the potential energy E(R) is possible by:
▪ approximate solution of the Schrödinger equation (quantum mechanics, QM)
▪ using empirical force fields (molecular mechanics, MM)
▪ hybrid QM/MM approach
▪ using coarse grained models
overview of method categories
C7790 Introduction to Molecular Modelling -7Potential
energy calculation
The calculation of the potential energy E(R) is possible by:
▪ approximate solution of the Schrödinger equation (quantum mechanics, QM)
▪ HF method
▪ post HF methods (MPn, CI, CC)
▪ DFT methods (various functionals)
▪ using empirical force fields (molecular mechanics, MM)
▪ forms and parameters of force fields
▪ hybrid QM/MM approach
▪ interface, type of QM-MM interaction, link atoms, ...
▪ using coarse grained models
hundreds of methods differing in the approximations used
they do not affect general laws/properties of E(R)
)(RE
C7790 Introduction to Molecular Modelling -8Potential
energy calculation
QM (Quantum mechanics) MM (Molecular mechanics) CGM (Coarse-grained mechanics)
)(RE )(RE )(RE
R - position of atom nuclei R - position of atoms R - position of beads
approximations approximations
Potential energy surface can be calculated by various method (model chemistry)!
C7790 Introduction to Molecular Modelling -9Graphical
representation of functions
f (x, y)
x
y
f (x)
x
function of one variable functions of two variables
C7790 Introduction to Molecular Modelling -10Graphical
representation of functions
f (x, y)
x
y
f (x)
x
function of one variable
of
x
y
functions of three variables
color is a representation
of function value
f(x, y, z)
volumetric display
functions of two variables
C7790 Introduction to Molecular Modelling -11Volumetric
representation of functions
Ensing, B.; Klein, M. L. Perspective on the Reactions between F– and CH3CH2F: The Free Energy
Landscape of the E2 and SN2 Reaction Channels. PNAS 2005, 102 (19), 6755–6759.
https://doi.org/10.1073/pnas.0408094102.
This example is about FES, but the
same technique can be used for PES.
coordinates
energy as iso surfaces with
different colors
C7790 Introduction to Molecular Modelling -12Displaying
E(R)
1 atom
2 atoms
N atoms
),,( 111 zyxE
),,,,,( 222111 zyxzyxE
),,,....,,,,,,( 222111 NNN zyxzyxzyxE
cannot be displayed
only volumetrically
Example:
enzyme BsoBI has ~ 10000 atoms => 30000 degrees of freedom,
it would be necessary to use 30000 + 1 dimensional space for visualization
C7790 Introduction to Molecular Modelling -13Property
of E(R)
The potential energy is invariant to:
• displacement (translation) of entire system
• rotation (rotation) of entire system } without the action of external
force fields
(e.g., electrostatic, magnetic, etc.)
C7790 Introduction to Molecular Modelling -14Invariance
to displacement
HCHO
)( 1RE )( 2RE
)()( 21 RERE =
21 RTR =+ ,....},,,,,{ TTTTTT zyxzyxT =
translation vector
)( 1RE )( 2RE
)()( 21 RERE 
!!! Does not apply to shift in a force field !!!
C7790 Introduction to Molecular Modelling -15Invariance
to rotation
HCHO
)( 1RE
)()( 21 RERE =
21 RR =Θ
rotation matrix
)( 1RE
)()( 21 RERE 
!!! Does not apply to rotation in a force field !!!
)( 2RE
)( 2RE
C7790 Introduction to Molecular Modelling -16Diatomic
molecule
),,,,,( 222111 zyxzyxE hydrogen molecule
▪ three degrees of translational freedom
▪ two rotational degrees of freedom
(molecule is linear)
C7790 Introduction to Molecular Modelling -17Diatomic
molecule
),,,,,( 222111 zyxzyxE
)(rE
interatomic distance
hydrogen molecule
6-5 = 1
hydrogen molecule
▪ three degrees of translational freedom
▪ two rotational degrees of freedom
(molecule is linear)
C7790 Introduction to Molecular Modelling -18Thriatomic
molecule
),,,,,,,,( 333222111 zyxzyxzyxE
water molecule
▪ three degrees of translational freedom
▪ three rotational degrees of freedom
C7790 Introduction to Molecular Modelling -19Triatomic
molecule
),,,,,,,,( 333222111 zyxzyxzyxE
water molecule
▪ three degrees of translational freedom
▪ three rotational degrees of freedom
),,( 21 rrE
9-6=3
r1 r2

Internal coordinates r1, r2,
displayable only volumetrically
C7790 Introduction to Molecular Modelling -20Graphical
representation of E(R)
E (x, y)
x
y
E (x)
x
)(1 Rfx =
)(1 Rfx =
)(2 Rfy =
➢ The E(R) function is not graphically representable due to its high dimensionality.
➢ Thus, only its relevant part that best captures the problem being studied is displayed in a
two- or three-dimensional space
transformation
(projection)
functions
C7790 Introduction to Molecular Modelling -21Graphical
representation of E(R)
E (x, y)
x
y
E (x)
x
)(1 Rfx =
)(1 Rfx =
)(2 Rfy =
What about other degrees of freedom?
transformation
(projection)
functions
)(2 Rrc f= )(3 Rrc f=0
)(
=


cr
RE
value of E(R) is minimal relative to rC
➢ The E(R) function is not graphically representable due to its high dimensionality.
➢ Thus, only its relevant part that best captures the problem being studied is displayed in a
two- or three-dimensional space
C7790 Introduction to Molecular Modelling -22Illustrative
example
C7790 Introduction to Molecular Modelling -23Stationary
points
x1 x2 x3
x1, x2 and x3 are stationary points (local extrema of the function). The properties of the
quantum states of the system are derived from them.
lVRTk ExEE ,1)( +=
C7790 Introduction to Molecular Modelling -24Stationary
points
x3
tangent of the function E(x) at the
stationary point has zero slope
(a=0)
Function slope is given by gradient of function
(i.e., first derivative of the function)
x
xE


=
)(
)tan(a
a
Condition necessary for a stationary
point
0
)(
=


x
xE
C7790 Introduction to Molecular Modelling -25Stationary
points
x3
x
xE


=
)(
)tan(a
a
Condition necessary for a stationary
point
0
)(
=


x
xE
0
)(
3
=


xx
xE
tangent of the function E(x) at the
stationary point has zero slope
(a=0)
Function slope is given by gradient of function
(i.e., first derivative of the function)
C7790 Introduction to Molecular Modelling -26Types
of stationary points
local minima
local maximum
C7790 Introduction to Molecular Modelling -27Determining
stationary point type
( ) ...
)(
2
1)(
)()(
2
2
2
+


+


+=+ x
x
xE
x
x
xE
xExxE
Taylor series:
C7790 Introduction to Molecular Modelling -28Determining
stationary point type
( ) ...
)(
2
1)(
)()(
2
2
2
+


+


+=+ x
x
xE
x
x
xE
xExxE
Taylor series:
At the stationary point:
( ) ...
)(
2
1
)()(
2
2
2
+


+=+ x
x
xE
xExxE
gradient is zero
square of the
deviation is
always positive
C7790 Introduction to Molecular Modelling -29Determining
stationary point type
( ) ...
)(
2
1)(
)()(
2
2
2
+


+


+=+ x
x
xE
x
x
xE
xExxE
Taylor series:
At the stationary point:
( ) ...
)(
2
1
)()(
2
2
2
+


+=+ x
x
xE
xExxE
square of the
deviation is
always positive
The value of function increases when deviating from a stationary
point, if any second derivative at a given point has a positive
value. The stationary point is then local minimum.
The value of function decreases when deviating from a
stationary point, if any second derivative at a given point has a
negative value. The stationary point is then local maximum.
0
)(
2
2



x
xE
0
)(
2
2



x
xE
gradient is zero
C7790 Introduction to Molecular Modelling -30Stationary
points
0
)(
2
2



x
xE
Local minimum:
0
)(
=


x
xE
Local maximum:
0
)(
2
2



x
xE
0
)(
=


x
xE
!!! required condition !!!
C7790 Introduction to Molecular Modelling -31Two-dimensional
case
E (x, y)
x
y
local minima
C7790 Introduction to Molecular Modelling -32Two-dimensional
case
E (x, y)
x
y
local minima
local maxima
C7790 Introduction to Molecular Modelling -33-
saddle
points
Two-dimensional case
E (x, y)
x
y
local minima
local maxima
C7790 Introduction to Molecular Modelling -34Generalization
for E(R)
Stationary point:
0
)(
=


R
RE required condition, each component
of the gradient must be zero
gradient has 3N components
Stationary point type:
)(
)()(
)()()(
)()()(
)()()()(
2
2
1
2
2
1
2
11
2
11
2
11
2
2
1
2
11
2
1
2
11
2
11
2
2
1
2
RH=


















































NN
N
z
RE
xz
RE
z
RE
yz
RE
xz
RE
zy
RE
y
RE
xy
RE
zx
RE
zx
RE
yx
RE
x
RE





N number of atoms
Character of the stationary point is
determined by Hessian, which is a
matrix of second derivatives of
potential energy.
Not to be confused with Hamiltonian!!!
C7790 Introduction to Molecular Modelling -35Properties
of Hessian
kkk cHc =
eigenvector
eigenvalue
Nk 3,...,1=
N number of atom
• 6 (5) eigenvalues ​​are zero - it corresponds to the translation and rotation of the system
• remaining eigenvalues:
• all positive - local minimum
• one negative, other positive - first order saddle point
• two negative, other positive - saddle point of the second order
• .....
• all negative - local maximum
Diagonalization of Hessian is a
method for finding
eigenvalues ​​and eigenvectors.
C7790 Introduction to Molecular Modelling -36•
6 (5) eigenvalues ​​are zero - it corresponds to the translation and rotation of the system
• remaining eigenvalues:
• all positive - local minimum
• one negative, other positive - first order saddle point
• two negative, other positive - saddle point of the second order
• .....
• all negative - local maximum
Properties of Hessian
kkk cHc = Nk 3,...,1=
N number of atom
chemically
significant
stationary points
Diagonalization of Hessian is a
method for finding
eigenvalues ​​and eigenvectors.
eigenvector
eigenvalue
C7790 Introduction to Molecular Modelling -37-
saddle
points
Why saddle points and local minima only?
E (x, y)
x
y
local minima
local maxima
C7790 Introduction to Molecular Modelling -38Diagonalization
of Hessian
)(
)()(
)()()(
)()()(
)()()()(
2
2
1
2
2
1
2
11
2
11
2
11
2
2
1
2
11
2
1
2
11
2
11
2
2
1
2
RH=


















































NN
N
z
RE
xz
RE
z
RE
yz
RE
xz
RE
zy
RE
y
RE
xy
RE
zx
RE
zx
RE
yx
RE
x
RE





)(
)(
0
)(
00
0
)(
0
000
)(
2
2
2
3
2
2
2
2
2
1
2
Rλ=


































Nc
RE
c
RE
c
RE
c
RE





Diagonalization of Hessian is an operation,
where a rotation of the coordinate system
is searched such, that the mixed second
energy derivatives are zero.
Only the diagonal elements of the matrix
can be non-zero .
eigenvalues
Eigenvalues of Hessian then determine the
curvature of the function in direction of
axes of the new coordinate system. These
axes are determined by eigenvectors,
which are orthonormal.
1=kcijji =cc .
orthogonal normalization
condition
C7790 Introduction to Molecular Modelling -39-
Energy,
Gradient, Hessian
C7790 Introduction to Molecular Modelling -40Potential
energy calculation
The calculation of the potential energy E(R) is possible by:
▪ approximate solution of the Schrödinger equation (quantum mechanics, QM)
▪ HF method
▪ post HF methods (MPn, CI, CC)
▪ DFT methods (various functionals)
▪ using empirical force fields (molecular mechanics, MM)
▪ forms and parameters of force fields
▪ hybrid QM/MM approach
▪ interface, type of QM-MM interaction, link atoms, ...
▪ using coarse grained models
hundreds of methods differing in the approximations used
C7790 Introduction to Molecular Modelling -41Energy
gradient calculation
R
R

 )(E
Energy gradient:
)(RE 













=
Nzzyx
,...,,,
111
it is a vector; the number of components is 3N














Nz
E
z
E
y
E
x
E )(
,...,
)(
,
)(
,
)(
111
RRRR
The gradient calculation can be performed:
• analytically
• numerically
C7790 Introduction to Molecular Modelling -42Analytical/Numerical
gradient
Numerical gradient calculation is used when the analytical gradient is not available, for
example due to the complexity of the implementation of the algorithm for its calculation.
Either forward differences (FD) or central differences (CD) method can be used to
calculate the numerical gradient . In rare cases, it is also possible to use multipoint
methods.
The CD method is more accurate than the FD method and therefore the preferred
method of gradient calculation.
Analytical gradient calculation is the preferred method of calculation in cases where the
expression and subsequent calculation of energy derivatives are easy.
computationally more demanding than energy calculation
FD requires 2*3*N energy calculations
computationally more demanding than analytical gradient
calculation
C7790 Introduction to Molecular Modelling -43Calculation
of Hessian
)(
)()(
)()()(
)()()(
)()()()(
2
2
1
2
2
1
2
11
2
11
2
11
2
2
1
2
11
2
1
2
11
2
11
2
2
1
2
RH=


















































NN
N
z
RE
xz
RE
z
RE
yz
RE
xz
RE
zy
RE
y
RE
xy
RE
zx
RE
zx
RE
yx
RE
x
RE





Hessian of energy:
it is a matrix; number of
components is 3Nx3N
Calculation of Hessian can be implemented:
• analytically (memory and computationally intensive)
• numerically (by the method of central differences)
• from energies (3 x N x 3 x N x 2 energy calculations)
• from gradients (3 x N x 2 gradient calculations)
C7790 Introduction to Molecular Modelling -44-
Summary
➢ Quantum states (thermodynamic microstates) from solution of SE are characterized by
stationary points on PES.
➢ Stationary point has zero gradient.
➢ Type of stationary point can be determined from PES curvature at the stationary point
(Hessian eigenvalue analysis).
➢ The most important stationary points are local minima (stable states such as reactants,
products, intermediates) and first order saddle points (transition states).
Energy and its gradient MUST be calculated at the same
level of theory (model chemistry)!
Energy, its gradient, and Hessian MUST be calculated at
the same level of theory (model chemistry)!