Free Energy Perturbation Calculations Within Quantum Mechanical

Chapter 10

Free Energy Perturbation Calculations Within Quantum Mechanical Methodologies

Downloaded by PURDUE UNIVERSITY on March 16, 2013 | http://pubs.acs.org Publication Date: May 5, 1996 | doi: 10.1021/bk-1996-0629.ch010

Robert V. Stanton, Steven L. Dixon, and Kenneth M . Merz, Jr.

1

Department of Chemistry, 152 Davey Laboratory, Pennsylvania State University, University Park, PA 16802

N e w l y developed techniques for use i n quantum free energy and potential of mean force calculations are presented. These techniques are the first to provide a direct means of performing arbitrary perturbations in single-topology quantum mechanical systems. The theoretical and practical considerations for carrying out the associated molecular dynamics simulations are discussed. Simulations are conducted for systems that are purely quantum mechanical, and for systems that involve a combination of quantum mechanical and molecular mechanical atoms. Preliminary results demonstrate these procedures to constitute a powerful tool in free energy calculations, with the potential to significantly increase the accuracy of simulations on both large-scale and small-scale systems. Free energy perturbation (FEP)(l-3) and potential of mean force (PMF)(4, 5) calculations within molecular mechanical ( M M ) force fields have proven to be two of the most powerful computational methodologies available. These techniques, which have been developed and refined over the last two decades, allow access to important thermodynamic quantities, such as the relative binding free energies of ions within an ionophore(6), or inhibitors to an enzyme.(7) W h i l e other techniques, such as linear response theory(8), have been developed to calculate these same quantities they may require additional parameterization. F E P and P M F calculations, in contrast, are done using statistical averages of the energy differences between the various intermediate states considered i n the calculation and thus may be applied to any system for which a potential energy surface is known. M a n y elegant applications of these theoretical techniques within molecular dynamics ( M D ) calculations currently exist within the literature.(4, 5,7,9-14) Recently, quantum mechanical ( Q M ) MD(15-18), as well as coupled potential M D ( 19-22) calculations have become computationally practical. In such calculations a system is modeled using a Q M Hamiltonian, within the B o r n Oppenheimer approximation, and is propagated classically through time using Newtonian mechanics. Coupled potentials (CP), i n addition to the Q M portion of 1

Corresponding author 0097-6156/96/0629-0142$15.00A) © 1996 American Chemical Society In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.


10.

STANTON ET AL.

FEP Calculations

143

the system, also include atoms which are modeled molecular mechanically. The division of a large system into Q M and M M portions allows for a portion of the system to be modeled using the more accurate Q M potential function while still including environmental effects from the M M portion of the system. Typically, the Q M calculation involves a region of particular interest such as the solute within a solvent system or the active site within an enzyme. Ideally the entire system would be modeled using quantum mechanics, however due to the scaling of the computational expense of Q M calculations with the size o f the system (between n ^ and rP, where η is the number of basis functions) fully Q M representations of large systems remain impractical. W i t h the increased use of quantum mechanics for M D simulations, a natural extension of these simulations is their use with F E P and P M F calculations. Unfortunately, though, the extrapolation of the F E P and P M F methods used within M M calculations to Q M potentials is not as straightforward as might be hoped. In a M M F E P calculation the relative free energy difference between two states is calculated by slowly perturbing the parameters appropriate for one state into those of the other. For example, i f a ketone group were being perturbed into an alcohol, the bond, angle, torsional, and van der Waals (vdW) parameters as w e l l as the atomic charges, would be coupled to a perturbation parameter λ such that λ =1 corresponds to the ketone and λ = 0 the alcohol. T w o approaches are available for the calculation of the intermediate states. These are referred to as the single and dual topology methods. In the dual topology method the forces and energy of the fractional lambda state are generated by weighted averages of these values from independent calculations at λ = 1 and λ = 0. The single topology method, alternatively, uses weighted averages of the parameter values to calculate forces and energy. A study contrasting the single and dual topology methods showed the single topology method to converge more quickly than the dual i n many cases.(23) If single topology methods are to be used within Q M systems techniques must be developed which allow for the calculation of systems corresponding to fractional values of the perturbation parameter. Four possible system perturbations can be defined within Q M - F E P calculations: 1) changing an atom type.; 2) changing the number of atoms; 3) changing the number of electrons; and 4) changing the number of atomic orbital basis functions. Perturbations between systems can, and usually do, incorporate many of these changes simultaneously. Recently, within our lab we have developed methods allowing for each of these four types of perturbations. W e w i l l outline these methods and the results of test calculations i n the remainder of this paper. A n alternative method for the calculation of F E P s within Q M and C P systems has been demonstrated by Warshel and Wesolowski.(24) In these calculations they determine the free energy for changing the Q M starting and end point structures into their M M equivalent. Then by calculating the M M F E P they are able to complete a free energy cycle and obtain the relative free energy difference between the two Q M systems. The method presented here involves a direct Q M conversion and requires no M M intermediate.

Theoretical Approach In the next four subsections we w i l l discuss the individual techniques necessary to change the number of electrons, atoms and orbitals within a system as w e l l as the identity of certain atoms. A l l of the derivations and tests shown here are done within the semiempirical P M 3 Hamiltonian(25, 26), although they could easily be extended for use with other Q M potential functions. Number of Electrons. Here we formulate the fractional electron method to allow for non-integer numbers of electrons i n a Q M system. The approach relies on the

In Chemical Applications of Density-Functional Theory; Laird, B., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1996.

144

CHEMICAL APPLICATIONS OF DENSITY-FUNCTIONAL THEORY

use of a pseudo-closed-shell expression for the electronic energy, where fractional occupation numbers for the M O ' s is assumed at the outset. If M is the number of atomic orbitals ( A O ) , η is the total number of electrons i n the system and n k is the number of electrons i n the k molecular orbital ( M O ) then it is possible to define the electronic energy of a system using the pseudo-closed shell formula: t n

Ε^

(D

= ΣΣ(Η +β )Ρ μν

μν

μν

μ=\ v=l


where,

H„=J*>)[4 G ,

v

Vi - Σ VJirfy.W

~

(2)

(3)

= f ï \ W t e ) - k ^ ) ) P i o

λ=ΐσ=ΐΙ-

Ρμν

τ,

Σ

J

z

C vk

flk

^

H * » ) = I I U M Z ^ z l m . W τ.«* τ

2

Τη Note that i n this derivation the number, nk, of electrons in a M O may take on a noninteger value (0^nk^2). W h i l e not representative of a physical system non-integer electrons have been used previously i n , for example, the so-called half-electron method.(27) Energy gradients and thus atomic forces can be computed analytically only i f the electronic energy is variationally optimized with respect to the M O ' s . Further, the energy expression (1) implicitly assumes that the M O ' s are orthonormal. Thus we require,

clc,=5, C)SCj

(7) (8)

= Ô

t

Equation 8 is the standard ab initio Hartree-Fock formulation, while equation 7 is employed with semiempirical treatments. T o keep the discussion general we w i l l focus on the Hartree-Fock procedure and not the semiempirical methods. Extension to semiempirical approaches is straightforward. The quantity which must be minimized is the Lagrangian L , occ occ

L = E -2ê {c]Sc -S ) ekc

ij

j

ij


(9)

10.

STANTON ET AL.

145

FEP Calculations

where Qj are Lagrange multipliers. The minimization is most conveniently carried out with respect to ck while treating ck as formally independent. Thus we have, M M

0

=

§

=

é \ °^k

ΣΣΚ

ΣΣ«ι

^

+ G

~2

v=l

[μ=\

i

-δΜ

(ίο) J

j

Formal differentiation, followed by a series of algebraic manipulations yields, MM


μ=\

Xp

v=l

M

M 5fp

occ

t

j

λ=\ σ=1

V^-k

Using the Hermitian nature of G the subscripts of G can be interchanged and the expression can be rewritten as, MM

Xp

M

M Sp

° = Σ Σ Κ + G , J & + Σ Σ ^ Γ μ=1 v=l

^jt

A=l σ=1

occ σ

Α „ - 2 Σ « Λ

d ) 2

j

^k

Rearranging terms gives,

M

M

ftp

occ

MM 0

= Σ Σ Κ μ=1 v=l

occ

+ 2G,vK v* M c

e

2

occ 0 = n ( H + 2G)c - 2 ^ X S c , t

Σ */ ; ; ε

Λ

k

< > 14

(15)

A n d finally, occ

«*Fc = 2 £ e S c , . k

v

(16)

where the Fock matrix F is defined as (H+2G). For a closed-shell system, the familiar equation occ

Fc = £ e , S c , k

(17)

j

would result because nk=2 for all occupied M O ' s . In the case of the fractional electron method, however, at least one occupation number w i l l be less than 2 and equation 17 is retained. A t this point it is customary to define a unitary transformation U of the occupied M O ' s that would diagonalize the matrix of Lagrange multipliers €kj, so that the problem to solve has the general form, Fcê.SCfc


(18)

146

CHEMICAL APPLICATIONS OF DENSITY-FUNCTIONAL THEORY


This is easily justified i n the case of a closed-shell system because the wavefunction and a l l the molecular properties calculated from it are unaffected by a unitary transformation of the occupied M O ' s . Unfortunately, there may be no actual wavefunction that corresponds to the fractional electron system, and the energy w i l l be affected by any unitary transformation that mixes fully occupied M O ' s with partially occupied M O ' s . So the unitary transformation argument cannot be used to convert equation 17 to equation 18. However, as long as equation 17 is satisfied for some set of Lagrange multipliers that preserve the orthonormality of the M O ' s , the corresponding vectors ck w i l l be nontrivial solutions that coincide with an extremum of Eelec- This is all that is required i n the fractional electron method, thus a diagonal matrix of Lagrange multipliers may be assumed, and the M O ' s are solutions of.

n Fc = 2e Sc k

k

k

k

(19)

Note that equation 20 can be related to equation 19 by defining an F Ι"Α=εί8€

4

F i = ^ F 2

(20) (21)

The eigenvalues of equation 20 are related to those of equation 19 by

Thus by making nk a function of λ we can alter the number of electrons i n a system. The states corresponding to fractional λ values w i l l be not be physical, however since the transformation between states is smooth and the free energy a state function only the free energy difference between the initial and final states can be determined. Number of Atomic Orbitals. Changing the number of atomic orbitals can be done by scaling the Fock matrix elements corresponding to these orbitals with the perturbation parameter λ. This requires that the perturbations be done i n the direction i n which the number of atomic orbitals decreases. A s λ approaches zero, the energy levels of the disappearing atomic orbitals becomes very high so that they no longer contribute to the occupied molecular orbitals, and thus do not affect the energy for the system. This simple method allows the deletion of atomic orbitals without the generation of a new independent Fock matrix as would be required for a dual topology Q M approach. Identity of an Atom. If the number of valence electrons and the locations of the atomic orbitals on two molecular species are identical, the energies for these two systems w i l l exhibit the same mathematical dependence on the atomic parameters, independent of the identity of the element. A brief derivation of the P M 3 a p p r o x i m a t i o n ^ , 26, 28), which is largely based on the A M 1(29) and M N D O ( 3 0 ) method, is given here as the dependency of the parameters is essential. The total heat of formation of a molecule, within P M 3 , is given as a combination of electronic and nuclear repulsion energies E l e c t and E u c along with the experimental heat of formation of the atoms, and the calculated electronic energy for the gaseous atom e

n


10. STANTON ET AL. = E

FEP Calculations

+

tUa

-ΣΕ* A

147

+ΣΔ#;

(23)

A

Eei is taken to be a function of (1) the ground state atomic orbital ( A O ) population of i , (given by the density matrix element rî), (2) the one-electron energies for A O i of the ion resulting from removal of the valence electrons, U u , and (3) the twoelectron one-center integrals, and . (24)

EÎi = f{P«JJ {ii\jj)m)) u

The two electron one center integrals are given as =G s, =Gsp» =Gpn,

Free Energy Perturbation Calculations Within Quantum Mechanical

Recommend Documents