General notes: You will have to write codes to complete this assignment. While I would prefer it if you wrote the codes using fortran 90 or 77, you are welcome to use whichever language you would like. Please send me electronic copies of all programs you write with instructions regarding how they should be compiled.

Question 1. The electronic density of states, N(E) describes the distribution of electronic energy levels as a function of energy. A plot of the density of states often exhibits bands – groups of states at similar energies that are separated by large regions of energy in which no states reside. These bands can correspond to either occupied states (e.g. valence bands) or unoccupied states (e.g. conduction bands). Integrating the density of states over the range of energies containing a band tells us how many electrons are within that band (for occupied state) or could be placed within that band (for unoccupied states).

A set of data corresponding to the density of states of Al₂O₃ is provided in Al2O3_DOS.txt. Use it to do the following:

1. Interpolate the data using (i) a piecewise sextic function and (ii) a natural cubic spline. Plot the data curves obtained in each case and use these curves to estimate the ranges of energies that contain bands.

Note you should identify three bands for this system. Two of these will be at energies below 0 eV and contain occupied states. The third band, at energies above 0 eV, contains unoccupied states.

1. Numerically integrate the DOS, and plot the integrated number of states as a function of energy. Use this plot to (i) identify the ranges of energies that contain the bands, (ii) determine the number of electrons contained in the two bands with occupied states, and (iii) the number of electrons that could be placed in the band containing unoccupied states.

Question 2. In class, we discussed linear regression in the context of fitting data to polynomials. However, this procedure can be applied to cases where the data represent forms that are not best described by a polynomial.

A set of (x,y) data corresponding to a Gaussian function ( y c= exp(−α(x x− ₀)²) is provided the file gaussian.dat. Use it to:

1. Describe how you can mathematically transform the expression for a Gaussian function into a form that can be fit with linear regression.
2. Write a linear regression program to fit the provided to the form above and determine the values of c, α, and x₀.

Question 3. n-heptane has a large number of degrees of freedom that correspond to rotations about C-C bonds. This molecule is shown below, with rotatable bonds of interest to this question labeled 1 through 4.

CH3 CH2 1 CH2 2 CH2 3 CH2 4 CH2 CH3

The energy of the system as a can be described, in kJ/mol, in terms of these torsions as:

E= 4.184^⎡⎢_⎣^∑i=⁴1 ⎡_⎣0.5cos(φ_i )+0.1cos(2(φπ_i − ))+0.5cos(3φ_i )_⎦⎤^⎤_⎦⎥

where ϕ_i is the torsion angle of bond i in radians, and is defined by the relative locations of the carbon atoms defining the torsion angle. For example, ϕ₁ = 0 when the carbon atoms in the leftmost methyl group and central methylene group are eclipsed when looking down bond 1.

The energy function described above will have several local minima and one global minimum. With this in mind:

1. Use the energy function above to determine the set of torsions at which the system is at the global minimum. You don’t need a code for this, just look at how the energy varies for each torsion and you can figure it out.
2. Write an optimization code for this system that uses the Newton-Raphson method. Use this code, and at least four properly selected points on the energy surface, to show how the NewtonRaphson method will lead you to the nearest local minimum. To illustrate this effect, provide the values of the torsions that you used at the start of the optimization, the set of torsions reached when the optimization was complete, and the energy of the optimized structure.
3. Write a simulated annealing code to locate minima for this system using the energy function above. Accept all steps that take the current structure to a lower energy structure. Accept steps that take the current structure to a higher energy structure using a Maxwell-Boltzmann factor:

P=exp⎛⎜⎝EcurrentRT−Enew ⎞⎟⎠

where R is the ideal gas constant in SI units, and T is the temperature in Kelvin. Run a series of simulations using T = 1, 10, 100, and 1000 K, with all simulations starting from some common point on the energy surface that does not correspond to the energy well containing the global minimum.

1. Report how many steps it took in each simulation for the system to reach the global minimum. Does the relationship between the number of steps and the temperature make sense to you and why?
2. Explore whether the system moves away from the global minimum once it is located. Report how many steps it takes for the system to reach a new minimum on the energy surface after it first reaches the global minimum. Does the relationship between the number of steps and the temperature make sense to you and why?
3. Modify your annealing code to start with a temperature of 1000 K and decrease it to lower temperatures over time. Explore and report how changing the temperature affects the speed with which the global minimum is located and the speed with which the system moves away from the global minimum to other minima.