Computational Spectroscopy III. Spectroscopic Hamiltonians

1. Computational SpectroscopyIII. Spectroscopic Hamiltonians (a) Introduction to MATLAB (b) Matrix operations (c) Programing (d) Diagonalizing Hamiltonian matrices Updated: March 25, 2008

2. (a) What is MATLAB? • MATLAB is the computational environment that we will use to solve spectroscopic Hamiltonians to find the resulting energy levels and wavefunctions. • MAT LAB stands for Matrix Laboratory • The basic object in MATLAB is the matrix. • Matrices are constructed, manipulated, matrix problems solved, and the results displayed. • Any problem that can be expressed in terms of matrices, that is in terms of linear algebra, can be easily solved in MATLAB. • A matrix is an m n array of numbers • A vector is a 1  n matrix • A scalar is a 1  1 matrix, that is just a number. • Complex numbers are no problem in MATLAB. • MATLAB calculations can be run in two ways: • Type in commands one at a time in the command window • Create a program, an “M-file”, that contains a series of such commands. • MATLAB comands read like ordinary algebra • Many high-level matrix operations are available. • Each command packs a big punch; you get a lot of work done without much programming.

3. Warnings about MATLAB • The hardest part is learning the user interface, and the various kinds of windows that it gives you. • You do need to work through the tutorials to learn this. • Be careful where you save your stuff. • Make a folder in your “My Documents” folder and keep all your MatLab stuff there. • Don’t leave any files in the MALAB directory or elsewhere on the computer. • The current directory is shown at the top of your main MATLAB window. • Your workspace contains all of the data and matrices that you have created. • The default name for your workspace is “matlab.mat”, but you can use other names with the “.mat” extension on the file name. • Your programs are saved separately as “M-files”, such as “myprog.m”. • These are just text files that can be opened and edited with any text editor. • You can save your matrices as text files: save ABC.txt ABC -ascii This saves the matrix ABC in the text file ABC.txt in the current directory. ABC.txt can then be opened with any text editor, with Excel, or other program. • Variable and matrix names are case-sensitive

4. How will MATLAB help us compute spectra? • To compute molecular energy levels, we need to solve the Schroedinger equation H = E . (1) • We will be concerned with the nuclear motion, therefore the Hamiltonian is H = T + Eelec(R), (2)where T is the nuclear kinetic energy operator and the effective potential energy is Eelec(R), the electronic energy calculated as a function of the nuclear coordinates R. Eelec(R) is what we could get from Gaussian or from another electronic structure program. • To solve the Schroedinger equation (1), we expand the wavefunction as a linear combination of a set of nbasis functions, • The wavefunction  then becomes a vector in the n-dimensional space of the basis functions, • the Hamiltonian H becomes an nnmatrix, E is a constant, and • H = E becomes linear algebra eigenvalue problem, that is solved in MATLAB with the single command “eig”. • Hilbert space is the infinite-dimensional vector space of wavefunctions. For practical reasons, we truncate to n dimensions. • Most of our work will be in setting up the matrix H. • We will also use matrix algebra to make the computer do most of that work. • Spectral transition frequencies are obtained by subtracting to get the difference in energy between energy levels connected by the selection rules.

5. What we didn’t get from Gaussian or Spartan • Didn’t solve the quantum mechanics of the nuclear motion. • No rotational energy levels or rotational (microwave) spectra. • No vibrational energy levels for fundamentals or for overtone and combination bands. • Got only classical harmonic frequencies • No treatment of large-amplitude motion such as • Torsion, inversion • Conformational change • Motion along a reaction coordinate • No coupling of nuclear degrees of freedom • Vibration-vibration (IVR = intramolecular vibrational redistribution) • Rotation-vibration (Coriolis coupling and centrifugal distortion) • Quantum chaos • Restricted ergodicity (Ergodicity is the core assumption of molecular mechanics calculations) • No coupling of electronic and nuclear degrees of freedom • Radiationless transitions, dynamics at conical intersections

6. (b) Matrix Operations • Live demonstration. A = [1,2,3;4,5,6;7,8,9] • Let A, B be matrices and c and d be constants. • The elements of A are Aij, i labels the rows, i=1,m, and j labels the columns, j=1,n • scalar multiplication: cA is c*A • Scalar addition: A+c is A+c • Matrix multiplication: C=AB or Cik=jAijBjk C=A*B • The number of columns in A must be equal to the number of rows in B. • Indentity matrix: I I= eye(n,n) • Always a square matrix • Matrix Inversion: I/A = inv(A)= I/A (I must have same dimension as A) • A must be square and “non-singular” to have an inverse • Matrix division: B/A= B*inv(A)= B/A • A matrix raised to a power: A2=AA; A3=AAA; etc. A^3 • A must be square • Exponentiation of a matrix: exp(A) = 1 + A + A2/2! +A3/3! + A4/4! + …= expm(A) • Many other functions of matrices are defined. Search “funm” in MATLAB help.

7. Fun with Complex Numbers • z =c+d*i with i = sqrt(-1) and c, d real • i and j are both predefined in MATLAB as the imaginary unit. • Complex conjugation: z* is conj(z) • Exponentiation of imaginary number uses Euler’s formula: • exp(i*pi) gives -1.000 + 0.000i • Functions and matrix operations work fine with complex numbers (where mathemetically allowed). • Complex numbers are user-friendly!!

8. More fun with matrices • Matrix transpose: At is A.’ • Matrix transpose conjugate: (At)* is A’. This is the more useful form. • Select a sub-matrix: A(1:2,:) selects the first two rows and all of the columns of A. • A(2,3) gives the 3rd element from the 2nd row. • Concatenate 2 matrices side-by-side [A,B] or one on top of the other [A;B]. • Inner product: A*B • Examples are matrix multiplication and “dot” product • If x is a column vector containing the variables xi and p is a column vectorthen A*x=p is a set of 3 linear equations in the 3 unkowns x1, x2, and x3.Get the solutions by typing x=A\p !! • A linear least quares fit is obtained with “\” if there are more equations (rows) than unknowns (columns). x=A\p • “\” is called “left divide”. “/” is “right divide” and can be used to solve x*A=p with x=A/p . Note that for matrices left and right divide are NOT the same.

9. (d) Programming • Priorities in programming • Documentation: Each subroutine should have extensive comments about • What it does • Definition and format of required input data • Definition and format of the output data • Description of intermediate steps and internal variables throughout the routine. • Structured programming: Each routine should be short (< 1 page including comments), and perform only one well defined activity. Refrain from using common or “global” variables. • Ease of producing error-free code: Use approaches quick to program and debug. The time and effort spent on programing is THE major limitation in producing useful computed results. • Any results saved onto disk should be in a common data format easily importable into other software. Use the capabilities of commercial software as needed; don’t try to reprogram these capabilities yourself. • Less important priorities in programming • Speed of execution: Less important for casual programmers. • Efficient use of RAM memory and disk space: Unimportant, these are now cheap and plentiful. • Object-oriented programming • Objects are data structures called classes that have particular operations or “methods” that act on that class and may have special definitions for that class. • Very popular in modern programming because it allows complicated codes to be more easily extended and maintained. • Supported by MATLAB but will not be taught in this course.

10. Kinds of programs in MATLAB:1. Scripts • A series of commands just the way you would type them into the command window. • Access the existing data in your workspace and can modify or add to that data. • May contain comments: % MYSCRIPTdoes whatever … • May contain loops. • Click on the “new M-file” icon at upper left and type in commands. • Save the file as myscriptname.m • Scripts are just a time-saving trick that saves the effort of typing the commands in by hand. • In some environments scripts are call “macros” or “batch files”. • Examples include UNIX scripts, Applescript, MS Word Macros, DOS batch files, etc. • Scripts are NOT used for structure programming, because there is no isolation of the script variables from the workspace.

11. Kinds of programs in MATLAB:2. Functions • Each function has its own workspace separate from your main MATLAB workspace. • Variables are transferred in the calling expression: C = myfunc(A, B) • Arguments in (…) are passed to the function and can be modified by the function. Arguments can be matrices, vectors, or scalars. • The “value” of the function when it is returned is C, which could be a matrix, vector, or scalar. • You can share your main (“base”) workspace with the function by using global variables of the “evalin” function but this is NOT recommended because it weakens the structured nature of your programing. • The first line of the m-file begins with function • The 2nd line is a comment line with the function name in upper case. To get help, type help MYFUNC. • The rest of the line have your executable code. • Functions can have subfunctions inside them or you can save the subfunctions as separate M-files. • Call functions from the command line or from within other functions. • MATLAB looks for your functions in the current directory. Save as myfunc.m: function C=myfunc(A,B) % MYFUNC evaluates A^2 + i*B AA=A’*A C=AA+i*B

12. MATLAB Programming Tips • Use descriptive function and variable names to make your code easier to understand. • Order subfunctions alphabetically in an M-file to make them easier to find. • Precede each subfunction with a block of help text describing what that subfunction does. This not only explains the subfunctions, but also helps to visually separate them. • Don't extend lines of code beyond the 80th column, otherwise, it will be hard to read when you print it out. • Name each function uniquely. • Be sure to document your programs well to make it easier for you or someone else to maintain them. Add comments generously, explaining each major section and any smaller segments of code that are not obvious. You can add a block of comments as shown here. %--------------------------------------------- % This function computes the ... <and so on> %--------------------------------------------- • Don't try to write the entire program all at once. Write a portion of it, and then test that piece out. When you have that part working the way you want, then write the next piece, and so on. It's much easier to find programming errors in a small piece of code than in a large program. • When making modifications to a working program, don't make widespread changes all at one time. It's better to make a few small changes, test and debug, make a few more changes, and so on. Tracking down a difficult bug in the small section that you've changed is much easier than trying to find it in a huge block of new code. • If you have a function that is called by only one other function, put it in the same M-file as the calling function, making it a subfunction.

13. Flow control:if … elseif … else … end function C = myfunc(A,B) % MYFUNC evaluates A^2 + i*B C=0; if A==B % Does A equal B? C=1; % If so, do this elseif A>B % If not, is A > B ? C=-1; % If yes, do this. else% If neither of the above is true, C=0; % then do this. end

14. Flow control:for … end and while function C = myfunc(n) % MYFUNC evaluates sums n integers C=0 for k=1:n C=C+k end function C = myfunc2(n) % MYFUNC2 evaluates sums n integers C=0 k=1 while k <= n C=C+k k=k+1 end • See alsoswitch … case … otherwise

15. Debugging your programs • Unless you end each command or program line with “;”, then the result of executing that line will appear in your command window. • Once one part is working add the “;” at the end of the relvant lines to surpress unwanted output on the screen. • Test your program on simple data and on limiting case where you know what the output should be. • Search on “debug” in the MATLAB help to find out about the sophisticated debugging tools available. • If necessary, step through the program in the MATLAB debugger while keeping a record of each line that gets executed on a printed copy of the program. Use different combinations of inputs until you have observed that every line of code is executed at least once.

16. (e) Eigenvalue problems • If Ax=cx where A is an nn matrix, x is a 1n column vector, and cis a constant, then x is an eigenvector of an matrix A and c is its eigenvalue. • The primary example for us will be H = E • Other important eigenvalue relationships include J2 = J(J+1)ℏ2 with J=0, 1/2, 1, 3/2, etc., and Jz = mℏ with m = -J, -J+1, … +J. • In quantum mechanics, the matrices that correspond to observables, such as the total energy H, the square of the total angular momentum J2, or the projection of J onto the z-axis, are Hermitian. • A Hermitian matrix is equal to its transpose conjugate. • In MATLAB notation, H is equal to H’. • The eigenvalues of a Hermitian matrix are real. • In general, there are n eigenvectors and n corresponding eigenvalues of an nn matrix • If 2 or more eigenvalues are the same, they are said to be degenerate. • The eigenvectors are orthogonal: v1v2 = iv*1iv2i = 0 v1’*v2 gives 0.00 • The eigenvectors are normalized: v1v1 = iv*1iv1i = 1 v1’*v1 gives 1.00

17. Finding eigenvalues with MATLAB • The command E = eig(H) will give a column vector E with the eigenvalues (calculated allowed energies) for the Hamiltonian H. • If you want the eigen vectors (wavefunctions) as well, then type [V,D] = eig(H) . • D is a square matrix with the eigenvalues alon the principal diagonal. Type E = diag(D) to eaxtract the eigenvalues into the column vector E as we had above. • V is a square matrix in which the column are the eigenvectors, in the same order as the eigenvalues in D and in E.

18. Diagonalizing the Hamiltonian • Means transforming to a new basis set so as to produce zeros everywhere in Hexcept along the principal diagonal. • The eigenvalue problem H=E is then solved. • The diagonal elements of H are the eigenvalues. • H is given as the matrix D in [V,D] = eig(H). • The new basis set is a linear combination of the old basis set. • The new basis functions are the eigenvectors. • Given as the volumns of the matrix V in [V,D] = eig(H) . • Diagonalizing the Hamiltonian is really applying the linear variation method using our original basis set. • See Ira N. Levine, Quantum Chemistry, 5th ed, pp 220-235.