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 (ASE).
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 meshUse 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); CondaPkg Found dependencies: /home/runner/.julia/packages/ASEconvert/pAQve/CondaPkg.toml
CondaPkg Found dependencies: /home/runner/.julia/packages/CondaPkg/0UqYV/CondaPkg.toml
CondaPkg Found dependencies: /home/runner/.julia/packages/PythonCall/83z4q/CondaPkg.toml
CondaPkg Resolving changes
+ ase
+ libstdcxx
+ libstdcxx-ng
+ openssl
+ python
CondaPkg Initialising pixi
│ /home/runner/.julia/artifacts/cefba4912c2b400756d043a2563ef77a0088866b/bin/pixi
│ init
│ --format pixi
└ /home/runner/work/DFTK.jl/DFTK.jl/docs/.CondaPkg
✔ Created /home/runner/work/DFTK.jl/DFTK.jl/docs/.CondaPkg/pixi.toml
CondaPkg Wrote /home/runner/work/DFTK.jl/DFTK.jl/docs/.CondaPkg/pixi.toml
│ [dependencies]
│ openssl = ">=3, <3.6"
│ libstdcxx = ">=3.4,<15.0"
│ libstdcxx-ng = ">=3.4,<15.0"
│ ase = ">=3.26,<3.27"
│
│ [dependencies.python]
│ channel = "conda-forge"
│ build = "*cp*"
│ version = ">=3.10,!=3.14.0,!=3.14.1,<4"
│
│ [project]
│ name = ".CondaPkg"
│ platforms = ["linux-64"]
│ channels = ["conda-forge"]
│ channel-priority = "strict"
└ description = "automatically generated by CondaPkg.jl"
CondaPkg Installing packages
│ /home/runner/.julia/artifacts/cefba4912c2b400756d043a2563ef77a0088866b/bin/pixi
│ install
└ --manifest-path /home/runner/work/DFTK.jl/DFTK.jl/docs/.CondaPkg/pixi.toml
✔ The default environment has been installed.
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₂
pseudopot. family : PseudoFamily("cp2k.nc.sr.pbe.v0_1.largecore.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.58786408077 -0.58 5.4 6.99s
2 -16.72503845823 -0.86 -1.01 1.0 4.56s
3 -16.73054095497 -2.26 -1.57 2.1 331ms
4 -16.73121893267 -3.17 -2.16 1.0 245ms
5 -16.73132379386 -3.98 -2.58 2.0 264ms
6 -16.73133308430 -5.03 -2.82 1.9 253ms
7 -16.73067960984 + -3.18 -2.40 2.2 278ms
8 -16.73132307738 -3.19 -2.96 2.4 264ms
9 -16.73129035651 + -4.49 -2.94 1.9 255ms
10 -16.73132088167 -4.52 -3.14 1.7 234ms
11 -16.73133648834 -4.81 -3.47 1.0 208ms
12 -16.73133861537 -5.67 -3.63 1.1 704ms
13 -16.73133955264 -6.03 -3.82 1.0 205ms
14 -16.73134016554 -6.21 -4.43 1.7 236ms
15 -16.73134014103 + -7.61 -4.39 2.1 293ms
16 -16.73134016413 -7.64 -4.49 1.0 219ms
17 -16.73134018611 -7.66 -4.67 1.0 222ms
18 -16.73134020025 -7.85 -5.27 1.4 237ms
19 -16.73134020016 + -10.09 -5.28 2.3 282ms
20 -16.73134020000 + -9.79 -5.29 1.3 241ms
21 -16.73134020042 -9.38 -5.72 1.0 227ms
22 -16.73134020039 + -10.48 -5.88 1.9 281ms
23 -16.73134020043 -10.39 -6.00 1.4 244ms
scfres.energiesEnergy breakdown (in Ha):
Kinetic 5.8593966
AtomicLocal -105.6099986
AtomicNonlocal 2.3494808
Ewald 35.5044300
PspCorrection 0.2016043
Hartree 49.5614135
Xc -4.5976634
Entropy -0.0000035
total -16.731340200431