17. STM Simulations - Bardeen’s

This tutorial walks through a Bardeen-formalism STM current calculation on a simple graphene-sample / tungsten-tip system. For the underlying theory, see the Scanning Tunneling Microscopy section.

We now turn to a simple sample tip system to introduce the basic keywords relevant to Bardeen current calculations. As customary, we begin with the calculation of the ground state density. Let’s create input file c2_lcao_scf.input as follows

LCAO.status = 1
info.calculationType = 'self-consistent'
info.savepath = 'results/c2_lcao_scf'
info.outfile = 'fc2_lcao_scf.out'
la = 4.648725932
atom.element = [1 1 2]
atom.fracxyz = 'c2.xyz'
domain.latvec = [[1/2  sqrt(3)/2     0]*5
[1/2 -sqrt(3)/2     0]*5
[  0          0    27]/la]*la
domain.lowres = 0.3
element.species = 'C'
element.path = './C_LDA_TM.mat'
element(2).species = 'W'
element(2).path = './W_LDA_TM.mat'
domain.boundary = [1,1,2]
mixing.type = 'density'

The atomic coordinates can be generated using the following code

da = 3.0425/27;
xyz = [2,1,(0.5-da)*3;1,2,(0.5-da)*3]/3;
repxyz = [];
for ii = 0:4
for jj = 0:4
repxyz = [repxyz; bsxfun(@plus, xyz, [ii,jj,0])];
end
end
repxyz = repxyz/diag([5,5,1]);
fid = fopen('c2.xyz','w');
fprintf(fid,'%d\n\n',size(repxyz,1)+1);
fprintf(fid,'C %1.7e %1.7e %1.7e\n', repxyz');
fprintf(fid,'W %1.7e %1.7e %1.7e\n',[0,0,0.5+da]);
fclose(fid);

The atomic coordinates can also be generated using the following Python script:

import numpy as np

a = 3.0425
da = a / 27
xyz = np.array([
    [2, 1, (0.5 - da) * 3],
    [1, 2, (0.5 - da) * 3]
]) / 3

repxyz = []
for i in range(5):
    for j in range(5):
        shifted = xyz + [i, j, 0]
        repxyz.append(shifted)

repxyz = np.vstack(repxyz)
repxyz = repxyz @ np.diag([1/5, 1/5, 1])

with open('c2.xyz', 'w') as f:
    f.write(f"{len(repxyz) + 1}\\n\\n")
    for coord in repxyz:
        f.write(f"C {coord[0]:.7e} {coord[1]:.7e} {coord[2]:.7e}\\n")
    f.write(f"W {0:.7e} {0:.7e} {0.5 + da:.7e}\\n")

You can visualize the system to make sure everything looks reasonable and then execute RESCU typing

rescu -s c2_lcao_scf.input
rescu -i c2_lcao_scf.input

Next, the input file of the stm-current calculation is defined as c2_lcao_stm.input as follows

info.calculationType = 'stm-current'
info.savepath = 'results/c2_lcao_stm'
info.outfile = 'fc2_lcao_stm.out'
rho.in = 'results/c2_lcao_scf'
stm.z0 = 14
stm.voltage = 1.2
units.energy = 'eV'

The important keywords of c2_lcao_stm.input are:

• info.calculationType When set to stm-current, RESCU computes the current passing from the sample to the tip using (9.3).

• rho.in Provides the path the to files containing the ground state density.

• stm.z0 If stm.z0 is set to \(z_0\), then the tunneling matrix elements are computed on the plane \(z = z_0\) using Eq. (9.2). Moreover, the system is decomposed into orbitals below (sample) and above (tip) the \(z_0\)-plane. Thence, the PDOS of the sample and tip (\(n^S\) and \(n^T\)), and the sample and tip wavefunctions (\(\psi^S\) and \(\psi^T\)) depend on the value of stm.z0.

• stm.voltage Bias between the sample and tip.

Execute RESCU typing

rescu -i c2_lcao_stm.input

Upon completion of the calculation, the results may be read from c2_lcao_stm.mat. RESCU may also save temporary results to c2_lcao_stm_stm_tmp.mat, or more generally to whereever the savepath keyword points plusstm_tmp appended. In this way, computationally expensive results can be reused to compute the current for different bias voltages; note that you should begin by calculating the highest bias voltage for that to work however. In particular, one will find the structure stm which contains

• I The current, in ampere (amp).

• M The tunneling matrix elements.

• voltage The bias between the sample and tip, converted to atomic units.

• z0 The separation surface \(z\) value, converted to atomic units.

In the above example, you should get an answer around \(-1.46\times 10^{-5}\) amp.