Evaluation of the Potential 
Energy Function for n-Butane

Potential Energy Function and Parameters

Terry M. Gray and Roger L. DeKock
Department of Chemistry
Calvin College
Grand Rapids, MI  49546

Introduction

Conformational analysis using molecular mechanics has become a standard tool in modern chemistry.  There are many user friendly computer programs available to students and researchers on a variety of computer platforms.  In this Mathematica  notebook the potential energy function will be examined in detail.  The potential energy function has the following general form:

     Energy    =    S bond stretch    +    S angle bend    +     S torsion angle     +
                                  all bonds                                 all  angles                             all torsions
                                         
                                             S non-bonded van der Waals   +   S non-bonded electrostatic
                                          all atoms                                                                  all atoms 
                                          

To gain a better understanding of the potential energy function we will 

 1) look at the functional form of each term and see the effect of changing the various critical 
      parameters.  
 2) create three dimensional plots of the energy as a function of different combinations of the
      terms.  
 3) use the minimization functions in Mathematica to explore the problem of local minima.
                                
This exercise uses a "united atom" model of butane.  Hydrogen atoms are not included in the this model but are implictly included with the "heavy" atom with which they are associated.  In the united atom model of butane, for example, there are only four atoms:  two -CH3 and two -CH2 united atoms.  
[The parameters are taken from the united atom forcfield developed for Amber described in Weiner et al., "A New Force Field for Molecular Mechanics simulation of Nucleic Acids and Proteins", J. Amer Chem. Soc., 106, 765 (1984)]

Bond Stretch

Harmonic Oscillator

The harmonic oscillator is the simplest model for the covalent bond.  The bond is treated as a spring in which the energy increases as the spring is either compressed or stretched from its equilibrium position.  The equation is

Ebond  = kb (b - b0)2

where  kb is the force constant for the spring (260.0 kcal/(mol Å2) for the C-C bond) and b0 is the equilibrium bond length (1.526 Å for the C-C bond)

Mathematica  input and graphs:

bcc=1.526;
kcc=260;
Clear[e1];
e1[b_] = kcc (b-bcc)^2
e1plot=Plot[e1[b],{b, 1.0, 3.0}, PlotRange->{0,80}]

Questions:

1.  What is the value of the energy at the minimum?  What is the value of b at the minimum energy?

2.  Change the value of kcc in the input cell and re-calculate the cell.  Describe what happens to the shape of the energy function as kb is made larger or smaller.  What happens to the minimum energy and the value of b at that minimum as kb is changed?

3.  What is the experimental basis for the ideal bond length (b0)?

4.  What is the most unrealistic feature of the harmonic oscillator model of the chemical bond?  [HINT:  Think about what happens to chemical bonds at large bond distances (b).]

Morse Function

The harmonic oscillator is a useful approximation for the chemical bond at bond lengths very close to the ideal bond length; it is also a very simple mathematical function and thus easy to compute.  However, more complicated functions are sometimes used to give more accurate calculations of molecular properties.   For example, the bond stretch term in the molecular mechanics force field MM3 has a cubic and a quartic term in addition to the quadratic term.  Another function which can be used to provide a more realistic bond stretching process is the Morse potential

Ebond  = Db [1 - e-a(b - b0)]2

where  Db is the depth of the potential minimum (83.0 kcal/mol for the C-C bond), b0 is the equilibrium bond length (1.526 Å for the C-C bond), and a defines the steepness of the potential  (1.5/Å).

[NOTE: This will be the only consideration of the Morse function.  All combined potential energy functions will use the harmonic oscillator for the bond stretch potential.]

Mathematica  input and graphs:

bmorse=1.526;
Dmorse=83;
amorse=1.5;
Clear[emorse];
emorse[b_] =Dmorse (1-Exp[- amorse (b-bmorse)])^2
emorseplot=Plot[emorse[b],{b, 1.0, 3.0}, PlotRange->{0,80}]
Show[e1plot, emorseplot]

Questions:

1.  Describe what happens to the shape of the energy function as a  is made larger or smaller.

2.  What is the physical meaning of the constant Db?

3.  How might the values of a, b0, and Db be obtained?

4.  Look at the plot where the Morse potential and the harmonic oscillator function are superimposed.  When b deviates only from b0 (within 0.2 Å), how does the Morse potential compare with the harmonic oscillator?  What is the bond energy associated with a 0.2 Å deviation?  

5.  Given a set of Morse potential parameters, change the value of kcc in the harmonic oscillator function to match the two functions at the minimum.  This is a "manual" version of parameterization.  In real parameterization the fit is not to some model function, but to some experimental spectroscopic data or some results of a quantum mechanic calculation.

6.  What are the advantages and disadvantages of using the Morse potential instead of the harmonic oscillator?  Why do you think that MM3 uses a polynomial expression instead of the Morse potential?

Bond Angle

The harmonic oscillator is also used to describe how the potential energy varies as the bond angle deviates from the ideal bond angle.  The bond angle is treated as a spring between two atoms bonded to a central third atom in which the energy increases as the spring is either compressed  or stretched away from its equilibrium position.  The equation is

Eangle  =  kq (q - q0)2

where  kq is the force constant for the spring (0.019 kcal / (mol degree2) for the C-C-C bond angle) and q0 is the equilibrium bond angle (112.4 deg for the C-C-C bond angle).

Mathematica  input and graphs:

thetaccc=112.4
kccc=0.019
Clear[e2]
e2[theta_]= kccc (theta-thetaccc)^2
Plot[e2[theta], {theta, 60, 160}, PlotRange->{0, 50}]

Questions:

1.   From the graph what is the value of the energy at the minimum?  What is the value of q at the minimum energy?

2.  What is the experimental basis for the ideal bond angle (q0)?

3.  Although the functional form of the bond stretch energy term and the angle bending energy term is the same, the range of energies that result from chemically meaningful deviations from the ideal bond length or bond angle is different.  What part of the mathematical function produces this effect?  Think about it from a chemical perspective.  Which is more costly energetically, increasing the C-C-C bond angle by 15 degrees or stretching the bond by 0.5 Å?

Torsional Energy

Energy barriers exist between conformations related to each other by simple bond rotations.  In our united atom representation of butane there is a torsional rotation that alters the relative position of the first and fourth atom.  The torsion energy is the term that describes this energy as a function of the torsion angle f.  An equation that is often used is the following three term Fourier series expansion:

Etor  = 0.5 V1 [1 + cos (2p f/360)] + 0.5 V2 [1 - cos (2p 2f/360)] 
+ 0.5 V3 [1 + cos (2p 3f/360)]

where  V1, V2, and V3 are the torsional force constants and f is the bond rotation angle in degrees.  For the n-butane torsional potential energy function, the coeffiecent terms are V1=0 kcal/mol, V2=0 kcal/mol, and V3=4.0 kcal/mol.  Although the minimum for this surface is not difficult to determine, we will begin to use Mathematica's built-in commands for finding the minima of functions.

Mathematica  input and graphs:

V1=0.0; V2=0.0; V3=4.0;

Clear[Term1]
Term1[phi_] = 0.5 V1 (1+Cos[2 Pi phi/360])
plotterm1 = Plot[Term1[phi],{phi,0,360}, 
PlotRange->{0, 4}]

Clear[Term2]
Term2[phi_] = 0.5 V2 (1-Cos[2 Pi 2 phi/360])
plotterm2 = Plot[Term2[phi],{phi,0,360}, 
PlotRange->{0, 4}]

Clear[Term3]
Term3[phi_] = 0.5 V3 (1+Cos[2 Pi 3 phi/360])
plotterm3 = Plot[Term3[phi],{phi,0,360}, 
PlotRange->{0, 4}]

Clear[e3]
e3[phi_] = Term1[phi] + Term2[phi] + Term3[phi]
plote3 = Plot[e3[phi],{phi, 0, 360}, 
PlotRange->{0, 6}]

FindMinimum [e3[phi], {phi, 200}]

Show[GraphicsArray[{{plotterm1,plotterm3},
            {plotterm2,plote3}}]]

Questions:

1.  In this exercise there is no first or second term in the torsional energy expression.  In many molecules or in different force fields all three terms are needed.  Let us take for example OH-C-C-OH torsional term in Amber force field.  The coefficients are V1=0.0, V2=1.0, and V3=4.0.  What is the difference from before?

For the following questions, let us take some arbitrary numbers to get an understanding of the effect of all the coefficients: V1=3.0, V2=1.0, and V1=.5.

2.  Examine the graph of the sum of the three terms, i.e. the torsion energy as a function of the rotation angle (the last graph). 
          a. At what torsional angles is the energy a minimum? 
               What is the energy at each of the minima?
          b. At what torsional angles is the energy a maximum? 
               What is the energy at each of the maxima? 

4.  Examine the graphs of each individual term of the Fourier expansion.  What is the effect of the factor 3 in front of the f in the third term?  What is the effect of the different values for V1, V2, and V3? 

5.  Notice in the FindMinimum command that the initial guess for f is 200 deg.  What minimum energy is found and what is the value f at the minimum?  Recalculate the FindMinimum cell with the initial guess at 90 deg.  What minimum energy is found and what is the value f at the minimum?  Does the minimizer find the minimum of the function?  Try this for several other initial guesses and record your guess, the final minimized energy, and the value of f at the minimum.  Use this information to explain the difference between a local minimum and the global minimum.  Why doesn't the minimization procedure always find the global minimum?  Discuss why this is a problem in molecular mechanics calculations.  Given this problem what approach might be taken to overcome this limitation?

6. The V2(f) term has a periodicity of two.  In order to examine the term separately, set V1=V3=0 and set V2=1.  Explain why the V2 parameter is important to describe the potential energy changes as a function of torsional angle for molecules containing double bonds such as ethylene.  

Make sure that you return V1=0, V2=0, and V3=4.0 and rerun the calculations with these correct coefficients for n-butane.

Non-bonded Interactions

Introduction

Non-bonded interactions are repulsive or attractive forces between atoms that are not bonded together by a chemical bond or mutually bonded to some central atom.  In principle there is a pairwise interaction between each atom and every other atom in the molecule.  For large molecules, such a proteins, this becomes a very large number.  There are two kinds of non-bonded interactions:  van der Waals and electrostatic.  We will discuss both kinds, but only use van der Waals interactions in our butane example.

van der Waals Interactions--General

Van der Waals interactions derive from slight perturbations of the electron distribution in a given atom as a result of the presence of neighboring atoms.  Classically speaking, the negatively charged electrons on one atom repel the electrons in a nearby atom causing a local positive charge in the nearby atom which then is attracted to the local negative charge of the original atom.  In other words, transient atomic dipoles are established that cause the neighboring atoms to be attracted to each other.  Van der Waals interactions are also known as induced dipoles.  

The Lennard-Jones potential is used to described van der Waals interactions.  In addition to having a distance-to-the-sixth-power attractive term, there is a distance-to-the-twelvth-power repulsive term that prevents non-bonded atoms from getting too close together.   The equation is

Evdw  = e [(r0/r)12 -  2 (r0/r)6]

where  r0 is the minimum energy separation for the two atoms (4.00 Å for a CC interaction), r is the interatomic distance, and e is the depth of the minimum (0.075 kcal/mol for a vicinal 1-4 CC interaction).

Mathematica  input and graphs:

epsilon=0.075
r0=4.00
evdw=epsilon ((r0/r)^12 - 2 (r0/r)^6)

Plot[evdw, {r, 0.01,6}, PlotRange->{-0.15,0.3}]

Questions:

1.  Change the value of e and recompute the cell to see the effect on the function.   Do the same for r0.  What is the physical meaning of e?  

2.  How large must r become before the energy term goes to zero?  How might knowledge of this be used to simplify the calculation of non-bonded terms in large molecules?

van der Waals Interactions--Butane

In the general introduction to van der Waals interactions the geometrical term, r, is simply the distance between two atoms.  In most molecular mechanics programs this distance is calculated by applying the distance formula to the xyz coordinates of the two atoms.  We have not used xyz coordinates in our model of butane, but rather internal coordinates, i.e. in terms of the bond length, bond angle, and the torsion angle (b, q, f, respectively).  We can continue this approach in our simple butane model and derive a trigonometric function that expresses the 1-4 interatomic distance also as a function of b, q, and f.  That function is

r1,4  =  b Sqrt[3 - 4 cos (2pq/360) + 2 cos2(2pq/360) - 2 sin2(2pq/360) cos(2pf/360)]

where  all the angles are given in degrees.  Three functions will be plotted:  1) the interatomic distance as a function of the torsion angle, 2) the van der Waals energy as a function of the torsion angle, and 3) the van der Waals energy as a function of the torsional and bond angle.

Mathematica  input and graphs:

Clear[r14]
r14[b_,theta_,phi_] = b Sqrt[3 
- 4 Cos[2 Pi theta/360] 
+ 2 (Cos[(2 Pi theta/360)])^2 
- 2 (Sin[(2 Pi theta/360)])^2 Cos[2 Pi phi/360]]

plotr14 = Plot[r14[bcc,thetaccc,phi],{phi, 0,360}]

Clear[e4]
e4[b_,theta_,phi_]=epsilon ((r0/r14[b,theta,phi])^12 
- 2 (r0/r14[b,theta,phi])^6)

plote4 = Plot[e4[bcc,thetaccc,phi], {phi, -180, 180}]

plote5 = Plot3D[e4[bcc,theta,phi],{theta,95,119.8},
{phi,-180,180}]

Show[GraphicsArray[{plotr14,plote4, plote5}]]
FindMinimum[e4[bcc,thetaccc,phi],{phi,180}]

Questions:

1.  What is the maximum and the minimum of the van der Waals energy as a function of the torsion angle?  At what torsion angles do these minimum and maximum occur?  Sketch the corresponding conformation.

2.  How does the van der Waals energy contribution as a function of torsion angle compare with the "intrinsic" torsion energy contribution modeled above?

Electrostatic Interactions--General

Electrostatic interactions are the energies due to the attractive force between two particles of opposite charge or to the repulsive force between two particles of the same charge.  Obviously , there are charges in ionic species, but even in neutral molecules there can be an uneven distribution of electrons among those atoms shared in chemical bonds.  This gives rise to partial charges.  In a neutral molecule all the partial charges sum to zero, but they must be considered, especially in short range interactions.  The classical equation describing these terms is simply Coulomb's Law

Eelect  = kelectq1q2/D r

where r is the interatomic distance, q1 and q2 are the charges in units of charge (not Coulombs), D is dielectric constant (78 for water, 4 for a non-polar environment, 1 for a vacuum), and kelect is the constant to convert to kcal/mol (331.6 kcal Å/(mol charge2)).  For a given set of input parameters, we will create four different plots showing two different scales for the interatomic distance and two different scales for the energy.

Mathematica  input and graphs:

kelect=331.6
Dielectric=78
q1=1
q2=1

Clear[eelect]
eelect[r_] = kelect q1 q2 / (Dielectric r)
eelect1 = Plot[eelect[r], {r, 0.1, 50}, PlotRange->{-1.0,1.0}]
eelect2 = Plot[eelect[r], {r, 0.1, 8}, PlotRange->{-1.0,1.0}]
eelect3 = Plot[eelect[r], {r, 0.1, 50}, PlotRange->{-50,50}]
eelect4 = Plot[eelect[r], {r, 0.1, 8}, PlotRange->{-50,50}]
Show[GraphicsArray[{{eelect1,eelect2},{eelect3,eelect4}}]]

Questions:

1.  Compute the energy for different values of D, different combinations of charges (+,+; -,-; and +,-), and for different magnitudes of charges (whole charges and partial charges).  Create a table showing the different terms varied and the energy for each.  

2.  How large must r become before the energy term goes to zero?  How does this compare with the van der Waals energy? 

3.  What is the effect of the larger dielectric constant (D) on the distance dependence of the electrostatic energy term?  Why does water have such a large dielectric constant?

Potential Energy Surface--Bond Term & Angle Term

So far we have been looking at the potential energy as a function of a single geometrical term: the bond length of just one bond or the bond angle between a given central atom and two atoms bonded to it.  

In this section we depart briefly from n-butane and consider a molecule having only bond stretching and angle bending terms, such as the united atom model of propane.  We will simplify the two indepedent bond stretches as two times the bond stretching term .
																																	E(b,q)= 2*E1(b) + E2(q)
where E1(b) is the bond stretching potential energy and E2(q) is the bond angle potential energy.  The constants are the same as the above sections.  

The energy will be plotted as a function of the bond angle, q, and the bond length, b.  Three different representations are given:  the 3D surface, the density plot, and a shaded contour plot.  With complex surfaces it is also more difficult to determine the minimum energy, especially if there are multiple local minima.

Mathematica  input and graphs:

Clear[e12]
e12[b_,theta_]=2*e1[b]+e2[theta]

plot3D12 = Plot3D[e12[b,theta], {b, 1.32, 1.72}, 
               {theta, 80, 140}]

plotDens12 = DensityPlot[e12[b,theta],{b, 1.32, 1.72},
      {theta, 80, 140}, PlotPoints->20, Mesh->False]

plotCont12 = ContourPlot[e12[b,theta],{b, 1.3, 1.7},
{theta, 80, 140}]

Show[GraphicsArray[{plot3D12,plotCont12,plotDens12}]]

FindMinimum [e12[b,theta], {b,1.5}, {theta, 100}]

Questions:

1.  What are the advantages and disadvantages of each kind of representation?

2.  What is the value of energy at the minimum?  What is the value of b and q at the minimum? 

3.  Notice in the FindMinimum command that an initial guess for b and q were given.  How do these initial guesses relate to the final value?  Change the initial guesses and recompute the cell.  Do this several times and record your guesses and the final minimum obtained.  Does the minimizer still find the minimum of the function?

Potential Energy Surface--Bond Term, Angle Term & Torsional Term

Now we combine three terms in our potential energy expression: bond stetching, angle bending, and torsion angle rotation.  The energy is in terms of three geometrical parameters: the bond length, b, the bond angle, q, and the torsion angle, f.  Strictly speaking, in our united atom model of butane, there are three independent bonds, two independent bond angles, and one torsion angle, giving us six independent variables altogether.  For simplicity, we will multiply the bond term by three and the angle term by two; this is a simplification since the each bond and angle can stretch or bend independently.

E(b,q,f) = 3*E1(b) + 2*E2(q) + E3(f)

where E1(b) refers to the bond stretch, E2(q) is the bond angle energy, and E3(f) is the torsional angle energy.

With many bonds and angles, as well as the other terms, the potential energy surface quickly becomes a function of three variables.  Energy versus three variables creates a four dimensional plot.  With dimensions greater than three, a plot produces a multi-dimensional hypersurface that is nearly impossible to visualize.  One way to get an idea of the four dimensional surface is to keep one of variables constant and plot the remaining variables versus energy to obtain a three dimensional surface.  In our example of n-butane, we will keep the bond length constant corresponding to the minimum energy for that term and then plot the energy as a function of the bond angle and torsion angle.
 
In the miniminization function of Mathematica,  we will allow all three variables to vary and find the minimum of the hypersurface.

Mathematica  input and graphs:

Clear[e123]
e123[b_,theta_,phi_]= 3*e1[b]+2*e2[theta]+e3[phi]
plot3D123 = Plot3D[e123[bcc,theta,phi], {theta, 103, 117},
 {phi, 0, 360}]

plotCont123 = ContourPlot[e123[bcc,theta,phi], 
                {theta, 103, 117}, {phi, 0, 360}]

Show[GraphicsArray[{plot3D123, plotCont123}]]

FindMinimum[e123[b,theta,phi], {b, 1.4}, {theta,100},
 {phi,200}]

Questions:

1.  From the graphs where is the global minimum?  Where are the local minima?  

2.  Use the FindMinimum command to determine the locations of the global and local minima in terms of b, q, and f. Record the value of the energy at each of the minima.

3.  A saddle point (called such because of its shape) is a point on the surface which is at a minimum with respect to one of the variables, but at a maximum with respect to another.  There are three saddle points in the surfaces plotted above.  Locate them approximately in terms of q and f.  What role does a saddle point play in moving from one minima to another during a conformational change?

Potential Energy Surface--Bond Term, Angle Term, Torsional Term & Non-Bonded (van der Waals) Term

Now we look at the potential energy surface with all four terms:  bond stetching, angle bending,  torsion angle rotation, and the 1-4 van der Waals interaction.  Because we have used internal coordinates to calculate the 1-4 interatomic distance, again the energy is in terms of three geometrical parameters: the bond length, b, the bond angle, q, and the torsion angle, f.  This time we will multiply the bond term by three and the angle term by two; this is still a simplification since the each bond and angle can stretch or bend independently.  With this simplification the equation is

E(b,q,f) = 3*E1(b) + 2*E2(q) +E3(f) + E4(b,q,f)

where E1(b) is the bond stretch term, E2(q) is the bond angle term, E3(f) is the torsional term and E4(b,q,f) is the non-bonding term.

Again in our visualization of the surface  the bond length will be held constant at the value corresponding to the minimum energy for that term alone, and we will plot the energy as a function of the two angles.  We will also look at a slice through the three dimensional surface at the ideal bond angle; thus, the plot will be total energy as a function of the torsion angle, f.

Mathematica  input and graphs:

Clear[e1234]
e1234[b_,theta_,phi_]=3*e1[b]+2*e2[theta]+e3[phi]+e4[b,theta,phi]

plot3D1234 = Plot3D[e1234[bcc,theta,phi], {theta, 100, 120},
 {phi, 0, 360}, PlotRange->{0,4},
 ViewPoint->{2.950, -1.880, 0.960}]

plot1234 = Plot[e1234[bcc,thetaccc,phi],{phi, 0, 360}]

plotCont1234 = ContourPlot[e1234[bcc,theta,phi], 
   {theta, 95, 125},{phi, 0, 360}, PlotRange->{0,6}, 
    Contours->12]

Show[GraphicsArray[{plot3D1234,plot1234,plotCont1234}]]

FindMinimum[e1234[b,theta,phi], {b, 1.5}, {theta, 100}, 
{phi, 130}]

Questions:

1.  From the graphs where is the global minimum?  Are there local minima?  Where are they?

2.  Use the FindMinimum command to determine the locations of the global and local minima in terms of b, q, and f. Record the value of the energy at each of the minima.  Think carefully about the results that you obtain together with any error messages that Mathematica gives you.

3.  Think carefully about the shape of the 3DPlot and the ContourPlot in terms of the structure of butane.  Why isn't the surface symmetrical with respect to bond angle as it has been in every preceding example? 

4.  Looking at the slice of the 3DPlot at the equilibrium bond angle, why does the energy plot disagree with the potential energy diagram that would be in your textbook?  The textbook has the energy barrier between the two gauche formations at about 6.9 while according to the slice this energy is 11.  Looking at the contour or 3DPlot, at what bond angle has the lowest energy for a value of f=0?  What is the energy here?

Troubleshooting problems

1.  If the program uses up to much memory for the system, change the amount of memory allotted to Mathematica  under get info function of the finder.
2.  If you still have memory problems with Mathematica, keep all the sections closed and run through the whole document before going through each section.
3.  You must go through the individual sections before you can run the last three sections or there will be an error.