So I'm trying to use the LCAO method to calculate molecular orbitals on a square lattice of hydrogen atoms. To this end, I need to compute the Hamiltonian matrix. With the nearest neighbor interaction approximation, the matrix should be pentadiagonal, analogous to the Laplacian matrix.

Now the problem is that I don't want to assume orthogonality of nearest neighbor atomic orbitals so there will be an overlap matrix. I believe it's also pentadiagonal but when I solved for molecular orbitals and their energies, the result is wrong. So I wonder if anyone has done similar things and could provide me some insight? Thank you!

If needed, I can post my Matlab codes and the computed Hamiltonian and overlap matrix here.

PS: The LCAO approach I described here can be found in Simon’s Oxford Solid State Basics textbook, Chapter 6, exercise (6.5)

PPS: I'm not constructing Bloch functions because I don't want to impose boundary periodicity constraints. I'm trying to calculate band structures of finite-sized nanoparticles but not bulk material in my project.

My code:

N = 10; %Size of the lattice is N^2
S = eye(N) + diag(zeros(N-1,1)+s12,1) + diag(zeros(N-1,1)+s12,-1); % Initialize the overlap block
Hop = 4*(eye(N)) + diag(zeros(N-1,1)+v12,1) + diag(zeros(N-1,1)+v12,-1); % Initialize the hopping block 

B0 = S*e_atom + Hop;   % block matrix for interactions between same row of atoms
B1 = eye(N)*(e_atom*s12+v12);  % block matrix for interactions between adjacent row of atoms
Stmp = kron(eye(N),S)+ kron(diag(1+zeros(N-1,1),1),eye(N)*s12) + kron(diag(1+zeros(N-1,1),-1),eye(N)*s12); % Pentadiagonal Overlap matrix 

H = kron(eye(N),B0) + kron(diag(1+zeros(N-1,1),1),B1) + kron(diag(1+zeros(N-1,1),-1),B1); %Hamiltonian matrix, justification see my notebook

[EVec, EVal] = eigs(H,Stmp, N^2,'sm'); %EVec contains eigenvectors, EVal contains eigenenergy on diagonal
E = real(diag(EVal));
[E, ind] = sort(E);
tmp = EVec(:,ind);
EVec = tmp;

For two adjacent hydrogen atoms, denote their $1s$ orbital by $|1\rangle$ and $|2\rangle$. $s_{12}=\langle1|2\rangle$, $v_{12}=\langle1|V_1|2\rangle$, $v_{22}=\langle2|V_1|2\rangle$ where $V_1$ is the Coulombic potential centered on the first hydrogen nuclei, $e_{atom}=-13.6eV$ is the energy of electron in an isolated $1s$ orbital.

  • $\begingroup$ It probably would help to see your code. At this point, what you are doing sounds right in principle, so it is likely some type of coding error. If that's the case, the question maybe better on the scicomp SE $\endgroup$ – Tyberius Dec 1 '19 at 16:23
  • $\begingroup$ @Tyberius Hi Tyberius, I've attached relevant codes to my problem. Perhaps I'm making wrong assumptions somewhere? $\endgroup$ – Macrophage Dec 1 '19 at 16:31
  • 1
    $\begingroup$ What are you comparing your calculation against? Or I guess phrased another way, how do you know your result is wrong? $\endgroup$ – Tyberius Dec 2 '19 at 2:07
  • $\begingroup$ @Tyberius Hi Tiberius, sorry for the late reply! I'm checking my density of state distribution with that in the textbook Orbital Interactions in Chemistry by Albright. I'm not sure if the DOS plot presented in that book is for wavevectors in the full first Brillouin zone or just along the specific path through high symmetry points. $\endgroup$ – Macrophage Dec 5 '19 at 21:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.