Modelling a gallium arsenide surface

This example shows how to use the atomistic simulation environment or ASE for short, to set up and run a particular calculation of a gallium arsenide surface. ASE is a Python package to simplify the process of setting up, running and analysing results from atomistic simulations across different simulation codes. For more details on the integration DFTK provides with ASE, see Atomistic simulation environment.

In this example we will consider modelling the (1, 1, 0) GaAs surface separated by vacuum.

Parameters of the calculation. Since this surface is far from easy to converge, we made the problem simpler by choosing a smaller Ecut and smaller values for n_GaAs and n_vacuum. More interesting settings are Ecut = 15 and n_GaAs = n_vacuum = 20.

miller = (1, 1, 0)   # Surface Miller indices
n_GaAs = 2           # Number of GaAs layers
n_vacuum = 4         # Number of vacuum layers
Ecut = 5             # Hartree
kgrid = (4, 4, 1);   # Monkhorst-Pack mesh

Use ASE to build the structure:

using ASEconvert
using PythonCall

a = 5.6537  # GaAs lattice parameter in Ångström (because ASE uses Å as length unit)
gaas = ase.build.bulk("GaAs", "zincblende"; a)
surface = ase.build.surface(gaas, miller, n_GaAs, 0, periodic=true);

Get the amount of vacuum in Ångström we need to add

d_vacuum = maximum(maximum, surface.cell) / n_GaAs * n_vacuum
surface = ase.build.surface(gaas, miller, n_GaAs, d_vacuum, periodic=true);

Write an image of the surface and embed it as a nice illustration:

ase.io.write("surface.png", surface * pytuple((3, 3, 1)), rotation="-90x, 30y, -75z")

Python: None

Use the pyconvert function from PythonCall to convert the ASE atoms to an AtomsBase-compatible system. This can then be used in the same way as other AtomsBase systems (see AtomsBase integration for details) to construct a DFTK model:

using DFTK
using PseudoPotentialData

pseudopotentials = PseudoFamily("cp2k.nc.sr.pbe.v0_1.largecore.gth")
model = model_DFT(pyconvert(AbstractSystem, surface);
                  functionals=PBE(),
                  temperature=1e-3,
                  smearing=DFTK.Smearing.Gaussian(),
                  pseudopotentials)

Model(gga_x_pbe+gga_c_pbe, 3D):
    lattice (in Bohr)    : [7.55469   , 0         , 0         ]
                           [0         , 7.55469   , 0         ]
                           [0         , 0         , 40.0648   ]
    unit cell volume     : 2286.6 Bohr³

    atoms                : As₂Ga₂
    atom potentials      : ElementPsp(Ga, "/home/runner/.julia/artifacts/9a2a5dc89d1b33bff2ad61eaf2d000191050d15c/Ga.gth")
                           ElementPsp(As, "/home/runner/.julia/artifacts/9a2a5dc89d1b33bff2ad61eaf2d000191050d15c/As.gth")
                           ElementPsp(Ga, "/home/runner/.julia/artifacts/9a2a5dc89d1b33bff2ad61eaf2d000191050d15c/Ga.gth")
                           ElementPsp(As, "/home/runner/.julia/artifacts/9a2a5dc89d1b33bff2ad61eaf2d000191050d15c/As.gth")

    num. electrons       : 16
    spin polarization    : none
    temperature          : 0.001 Ha
    smearing             : DFTK.Smearing.Gaussian()

    terms                : Kinetic()
                           AtomicLocal()
                           AtomicNonlocal()
                           Ewald(nothing)
                           PspCorrection()
                           Hartree()
                           Xc(gga_x_pbe, gga_c_pbe)
                           Entropy()

In the above we use the pseudopotential keyword argument to assign the respective pseudopotentials to the imported model.atoms. Try lowering the SCF convergence tolerance (tol) or try mixing=KerkerMixing() to see the full challenge of this system.

basis  = PlaneWaveBasis(model; Ecut, kgrid)
scfres = self_consistent_field(basis; tol=1e-6, mixing=LdosMixing());

n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -16.58731591661                   -0.58    5.3    456ms
  2   -16.72506460759       -0.86       -1.01    1.0    237ms
  3   -16.73052439263       -2.26       -1.57    2.2    290ms
  4   -16.73121756464       -3.16       -2.16    1.0    242ms
  5   -16.73131973274       -3.99       -2.58    1.8    280ms
  6   -16.73133162805       -4.92       -2.82    2.0    278ms
  7   -16.73087948461   +   -3.34       -2.48    2.1    282ms
  8   -16.73118267309       -3.52       -2.64    2.4    285ms
  9   -16.73102793209   +   -3.81       -2.55    2.1    264ms
 10   -16.73131582679       -3.54       -3.07    1.1    219ms
 11   -16.73133869759       -4.64       -3.53    2.0    256ms
 12   -16.73133788603   +   -6.09       -3.54    1.3    228ms
 13   -16.73133952188       -5.79       -3.79    1.9    281ms
 14   -16.73134005216       -6.28       -4.13    1.0    220ms
 15   -16.73134019533       -6.84       -4.72    2.0    261ms
 16   -16.73134019748       -8.67       -4.87    2.0    289ms
 17   -16.73134019997       -8.60       -5.18    1.0    217ms
 18   -16.73134020040       -9.37       -5.55    1.2    240ms
 19   -16.73134020034   +  -10.22       -5.76    1.2    232ms
 20   -16.73134020041      -10.17       -5.98    2.0    259ms
 21   -16.73134020043      -10.64       -6.40    1.0    212ms

scfres.energies

Energy breakdown (in Ha):
    Kinetic             5.8593979 
    AtomicLocal         -105.6100207
    AtomicNonlocal      2.3494813 
    Ewald               35.5044300
    PspCorrection       0.2016043 
    Hartree             49.5614343
    Xc                  -4.5976640
    Entropy             -0.0000035

    total               -16.731340200433