Computers in Chemistry

Molecular Conformations

Eric R. Bittner

The aim of this assignment is to introduce numerical methods for conformational modelling, and for solving minimisation problems. These will be illustrated with simple examples. There is more here then you will be able to do in two hours. Do not worry about this; just work at your own rate and learn as much as you can! Do the purple parts.

$Get :: noopen : Cannot open \" BasicCalculations . nb \" .$

•Pre-requisites

Load the statistics and graphics packages:

In[13]:=

In[14]:=

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

In[15]:=

and the default style for text in graphics is set as follows.

In[16]:=

•3.0. Preliminaries

In this module we need to be able to calculate and plot in one-, two-, and three-dimensions.

•3.0.2 Plotting

The Basic Calculations palette can help you find out more about how to plot a function.

A scalar function depending on a single variable , can be plotted using

To plot a function of two variables use

If required, you can improve the resolution of a plot by increasing the number of plot points using the PlotPoints.

The more points you choose, the longer it will take to compute and render the image.

Plot the following functions for and :

(a)		O
(b)		O
(c)		O
(d)		O
(e)		O

•3.0.3 Finding Minima

To find the value and location of the minimum of a scalar function we can use the function. Here is a starting value from which the will attempt to proceed to find a local minimum. Clearly in the case of some particularly non-monotonic functions, this may or may not actually lead to the global minimum.

0. Find minima for the following functions. State whether the value found is a local or global minimum.:

(a)		O
(b)		O
(c)		O
(d)

(a)
(b)
(c)
(d)

•3.0.4 Animated Graphics

There are several ways to create an animated graphic of a static or dynamic 3D model in Mathematica. The following code generates a series of graphs showing a coil that rotates about the -axis:

$anime(θ_) := ParametricPlot3D[{r sin(x), r cos(x), (pitch x)/(2 π)}, {x, θ, θ + Φ}, SphericalRegion -> True, Axes -> False, Boxed -> False, ViewPoint -> {4, 3, 2.5}] ;$

Close the group of graphics cells (by double-clicking on the grouping cell bracket) and then double-click on one of the images to animate the sequence (or use Animate Selected Graphics under the Cell menu).

Run the above code and try varying the pitch, radius and maximum θ. [3 Points]

•3.1. Coulomb Potential

Computer modelling (or simulation) is commonly used when exact analytic solutions of mathematical models are unavailable. When used for differential equations, it produces approximate solutions only, so it is important to be aware of the potential limitations of the method.

Two point charges and , separated by distance , represent an electrical potential energy given by

where is a constant reference level corresponding to the maximally separated limit . The Coulomb force is found from by the relation

Since is independent of , we can choose our reference to be zero without affecting -- that is: set as . For a group of point charges with distances between the and charges, the total potential energy is the sum of all pairwise contributions:

Example: Consider the one dimensional case of two fixed charges and at distance apart and one movable charge at some point from . The total PE is then :

What are the practical restrictions on the value of ?

What is the minimum possible value for ?

•In order to plot we shall incorporate the unit-dependent constants and into the expression for , by defining a new potential . Then for a charge in the vicinity of a point charge we will have . ( being the valence of the charge)

0. Plot versus for the case , , and . Let , and use range for .

In[51]:=

$Vtot[r_] := (q1 q2)/R + (q1 q3)/r + (q2 q3)/(R - r) ; params = {q1 -> e, q2 -> 2 e , q3 -> e, R -> 1, e -> 1} ; Plot[Vtot[r] //. params, {r, -2, 2}]$

[Graphics:HTMLFiles/index_78.gif]

Out[53]=

0. Can you find an analytic solution for the value of between 0 and 1 that minimises ? What will the force be at this position?

In[59]:=

Out[59]=

${{x -> (q1^(1/2) R)/(q1^(1/2) - q2^(1/2))}, {x -> (q1^(1/2) R)/(q1^(1/2) + q2^(1/2))}}$

In[71]:=

$Vtot[(q1^(1/2) R)/(q1^(1/2) + q2^(1/2))] //. params // N D[Vtot[r], r] /. {r -> (q1^(1/2) R)/(q1^(1/2) + q2^(1/2))} // Simplify$

Out[71]=

Out[72]=

•

Let the two fixed charges of the previous question lie on the axis () and allow the third 'mobile' charge to have arbitrary coordinates. Display the potential energy surface using Plot3D as in §3.0.2 above. You may wish to increase the number of PlotPoints.

In[73]:=

In[122]:=

[Graphics:HTMLFiles/index_94.gif]

[Graphics:HTMLFiles/index_95.gif]

What can you now say about the value of that you found in Question 0 above that minimised for the one-dimensional case?

Repeat the 2D example with . Briefly describe what might happen if an electron or a proton were injected into the vicinity of the two fixed positive charges.

Repeat the 2D example.  This time we will make a trajectory of the electron's path.
(a) First, write down the equations of motion followed by the electron.  Do this in Mathematica by defining (x,y,vx, & vy) as the (xy) position and (vx,vy) velocity components.    Assume the electron has unit mass (m= 1) and e = 1 (we're using Atomic units in this case...so length is in Bohr and time in atomic units).
(b) Calculate the force on an electron at point (x,y) due to the charges.  You will have 2 force components, fx and fy, both functions of (x,y) position.
(c) Define a set of equations of motion and initial conditions:   Here you choose the initial coordinate point (xo,yo)

     eom = {x'[t]==vx[t], y'[t]==vy[t], vx'[t]==fx[x[t],y[t]], vy'[t]==fy[x[t],y[t]]};
     inits = {x[0]==xo, y[0]==yo, vx[0]==0,vy[0]==0};
     diffeq = Join[eom,inits];

   (d.) Solve the equations of motion using
   tmax = 10 (or so)
  traj =  NDSolve[diffeq,{x[t],y[t],vx[t],vy[t]},{t,0,tmax}];

  e.) Plot the trajectory in 2D:
  plot = ParametricPlot[Evaluate[{x[t],y[t]}/.traj],{t,0,tmax}]

  ---CAVEAT EMPTOR...you may need to play with the options in NDSolve to get the best trajectory.

  f.) Try finding both bound states and scattering states by varying the initial velocities.   Can you find "chaotic" states and "regular" or "periodic" solutions?

$fx[x_, y_] := - D[Vtot[X, Y], X] /. {X -> x, Y -> y} fy[x_, y_] := - D[Vtot[X, Y], Y] /. {X -> x, Y -> y}$

[Graphics:HTMLFiles/index_102.gif]

In[158]:=

[Graphics:HTMLFiles/index_104.gif]

In[163]:=

[Graphics:HTMLFiles/index_106.gif]

In[372]:=

[Graphics:HTMLFiles/index_108.gif]

As I mentioned in class, there is an interesting connection between the shapes of quasi-periodic orbitals (i.e. trajectories) and the corresponding quanutm wavefunction. For more information, explore Eric Heller's web-site at Harvard (http://monsoon.harvard.edu)

•3.2. Lennard-Jones Potential

When two atoms approach each other, the interaction is much more complicated than a simple Coulomb potential due to electron shielding and orbital dynamics as well as the fundamentally quantum nature of atomic systems. Solving the equivalent of Schroedinger's equation for any practical situation is extremely difficult, if not impossible. However, a level of understanding can be achieved by a semi-classical approach based on measured values of various atomic properties used to estimate parameters in a simplified model.

Empirical modeling of experimental data leads to the following approximate expression for the interaction between two atoms and : , where and are constants for each pair of atoms. We shall assume is measured in , and and take appropriate units such that is measured in Joules.


HH	0.1974	0.060
CC	1.5540	0.154
NN	1.5246	0.148
OO	1.5414	0.148
CO	1.5414	0.151
CN	1.5372	0.151
CH	0.5376	0.109
NO	1.533	0.148
HO	0.5208	0.104
NH	0.5250	0.104

Experimental values of A, B, and the Covalent Bondlength (nm), for the Lennard-Jones model.

A general expression for the Lennard-Jones potential function is then:

$V _ LJ[d_, {A_, B_}] := A/d^12 - B/d^6 ;$

Plot for the hydrogen bond (). Find the value of that minimises . Compare this value to the covalent bondlength.

•3.3. Ethane Rotational Conformation

A common method used in the modeling of large molecules involves a combination of Lennard-Jones and Coulomb potentials:

The effective charges are not the atomic charges but somewhat less due to the effect of screening and cancellation between the positive nucleus and the negative electron shells. However, in this section we restrict ourselves to the Lennard-Jones potentials of §3.2.

•3.3.1 Building the Ethane Molecule

Ethane contains two covalently bonded carbon atoms with three hydrogens connected to each carbon. X-ane is a hypothetical molecule made by replacing one of the hydrogens in each group by a large atom X. We shall choose to place one of the two carbon atoms at the origin with the other along the axis. The -C or -CX group attached to the carbon at the origin shall be assumed fixed, while the other three are free to rotate about the second carbon. Both groups are also assumed to adopt tetrahedral structure.

First we define the covalent bondlengths for the C-C bond, for the C-H bond, and for a C-X bond where X is some large atom such as bromine. and the positions of the atoms in terms of these. We introduce a rotation angle θ as a variable to assist in the minimisation of the potential energy. The positions of the three mobile atoms are therefore expressed as functions of θ. Here are the co-ordinates for ethane:

$C _ 1 = {0, 0, 0} ; C _ 2 = {0, 0, ! _ C} ;$
$H1 = {(2 2^(1/2) ! _ H)/3, 0, -! _ H/3} ;$
$H2 = {1/3 (-2^(1/2)) ! _ H, 2/3^(1/2) ! _ H, -! _ H/3} ;$
$H3 = {1/3 (-2^(1/2)) ! _ H, -2/3^(1/2) ! _ H, -! _ H/3} ;$
$H4[θ_] := {2/3 2^(1/2) ! _ H Cos[θ], 2/3 2^(1/2) ! _ H Sin[θ], ! _ C + ! _ H/3} ;$
$H5[θ_] := {2/3 2^(1/2) ! _ H Cos[(θ + (2 π)/3)], 2/3 2^(1/2) ! _ H Sin[(θ + (2 π)/3)], ! _ C + ! _ H/3} ;$
$H6[θ_] := {2/3 2^(1/2) ! _ H Cos[(θ + (4 π)/3)], 2/3 2^(1/2) ! _ H Sin[(θ + (4 π)/3)], ! _ C + ! _ H/3} ;$

In order to see the molecule in 3D as a function of θ we can use the following routine:

$spin[θ_] := Show[Graphics3D[{PointSize[0.04], RGBColor[1, 0, 0], Point[C _ 1], Point[C _ 2], RGBColor[0, 1, 0], Point[H1], Point[H2], Point[H3], Point[H4[θ]], Point[H5[θ]], Point[H6[θ]], GrayLevel[0], Thickness[0.01], Line[{C _ 1, C _ 2}], Line[{C _ 1, H _ 1}], Line[{C _ 1, H _ 2}], Line[{C _ 1, H3}], Line[{C _ 2, H4[θ]}], Line[{C _ 2, H5[θ]}], Line[{C _ 2, H6[θ]}]}], PlotRange -> ( -2 2 ), ViewCenter -> {0, 0, ! _ C/3}, Boxed -> False, ViewPoint -> {4, 3, 2.5}] ; -2 2 -1.5 2.5$

3.3.2 Exercises

0. Modify the previous method for the case where and are replaced by the large atom X as and , assuming C-X bondlength as above and the following values for LJ coefficients: H-X ; X-X ; [2 Points]

0. Plot the Lennard-Jones potential for the H-H, H-X, and X-X interactions in isolation. Find the distances for which each potential is minimised. [2 Points]

0. Assemble the expression for the total Lennard-Jones potential for the ethane and X-ane molecules as functions of the rotation angle θ. [2 Points]

[θ_]:=
[θ_]:=

0. Plot this total potential function versus over the range . Find the value(s) of for which each function is minimised. [1 Point]

Ethane: θ=
X-ane: θ=

0. What is the most likely conformation for X-ane according to this model? [1 Point]

Converted by Mathematica (September 27, 2002)