# C60 orbital calculation

I am currently trying to reproduce the results published by Hands et al.1 using MATLAB. They calculated the bases of the C60 wave functions of HOMO, LUMO and LUMO+1.

I did the following:

1. I implemented these functions into MATLAB syntax, e.g.  Tp1uxL5 = @(x,y,z) 35.*z.*(x.^4-6*x.^2.*y.^2+y.^4)-5*x.*(1-14*z.^2+21.*z.^4); Tp1uyL5 = @(x,y,z) 5*y.*(28.*z.^4-28*x.*z.*(x.^2-y.^2)-(1-7*z.^2).^2); Tp1uzL5 = @(x,y,z) 2*z.*(15-70*z.^2+63*z.^4)+7*x.*(5*(x.^2-y.^2).^2-4*x.^2); 

for the pentagon prone $T^p_{1ui} (L=5)$ (i=x,y,z) functions, multiplied by an exponential function $\exp(-k R)$ utilizing $k = \frac{Z_\text{eff}}{2} a_0$ with the effective charge $Z_\text{eff}=3.14$, the Bohr radius $a_0 = 52 \text{ pm}$, and $R = \sqrt{x^2+y^2+z^2}$ as proposed by the authors.

1. I normalized the functions by integrating in three dimensions using MATLAB's integral3 function, i.e.

A = integral3( FKT , -B,B , -B,B , -B,B )

with the boundary $B$ which i varied between $0.5 \text{ nm} \le B \le 10 \text{ nm}$ with no change in result. I also tried just generating the 3D "wavefunction matrix" in those boundaries and taking the sum over all values.

1. I equated each base $T$ times exponential times normalization factor $\frac{1}{\sqrt{A}}$, squared for the propability amplitude, added the corresponding functions for the probability density $\Psi^2$ of HOMO, LUMO and LUMO+1 as well as pentagon, hexagon and double-bond prone configuration.

2. I generated plots:

• A) "Constant height" STM images at $z=z_0$ by equating $\Psi^2 (x,y,z_0)$

• B) "Constant current" STM images equating $\Psi^2 (x,y,z) = const.$ and visualizing the $z$ map.

• C) 3D plots of the electron density (comparable to Fig. 4,6,8,11 in 1) by MATLABs isosurface(X,Y,Z,DOS,DOSvalue).

Unfortunately, I fail to reproduce the pictures of the report. Below you see a calculated constant current image; constant height and the 3D probability density fronts show the same structure, so the visualization part of the code I guess is fine.

I checked and checked the functions again, I am starting to get mad over this. I am not reproducing the results of the paper. What am I doing wrong? Fig. 1. Simulated constant current images. First row: LUMO, second row: HOMO, third row: LUMO+1. First column: pentagone prone, second column: hexagone prone, third column double-bond prone.

View image here.

1 "Calculation of images of oriented C60 molecules using molecular orbital theory" Phys. Rev. B 81, 205440 (2010), http://core.ac.uk/download/pdf/97962.pdf

• I think my visualization is fine, since I get the same structure in constant current, constant height as well as the 3D orbital plots. I guess that I make a systematic mistake in assuming the wave function form.. – t0xic Jul 16 '15 at 16:32
• I think that you need to describe better your question. For example, you never mention in the text that they use Hückel metho, as they write in the very first sentence of the paper. You can't just wait people to go an read the paper for you. You need to provide enough details, so they can help you. – nicoguaro Jul 22 '15 at 19:26
• I'm sorry, I am not too familiar with theoretical methods for calculating molecular orbitals or such, I assumed it is enough to start with the basis set they gave. Am I missing something fundamental? – t0xic Jul 22 '15 at 19:47
• You should start with understanding the method then. That should be done before to attempt to obtain reported results. Have you done it for simpler molecules? (H2 maybe?). – nicoguaro Jul 22 '15 at 22:03
• Than I probably should do that... – t0xic Jul 23 '15 at 20:35