Creating supercells with pymatgen¶

The Pymatgen python library allows to setup solid-state calculations using a flexible set of classes as well as an API to an online data base of structures. Its Structure and Lattice objects are directly supported by the DFTK load_atoms and load_lattice functions, such that DFTK may be readily used to run calculation on systems defined in pymatgen. Using the pymatgen_structure function a conversion from DFTK to pymatgen structures is also possible. In the following we use this to create a silicon supercell and find its LDA ground state using direct minimisation. To run this example Julia's PyCall package needs to be able to find an installation of pymatgen.

First we setup the silicon lattice in DFTK.

In [1]:
using DFTK

a = 10.263141334305942  # Lattice constant in Bohr
lattice = a / 2 .* [[0 1 1.]; [1 0 1.]; [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]];


Next we make a [2, 2, 2] supercell using pymatgen

In [2]:
pystruct = pymatgen_structure(lattice, atoms)
pystruct.make_supercell([2, 2, 2])
atoms = [Si => [s.frac_coords for s in pystruct.sites]];

┌ Warning: vendor() is deprecated, use BLAS.get_config() and inspect the output instead
│   caller = npyinitialize() at numpy.jl:67
└ @ PyCall /home/runner/.julia/packages/PyCall/3fwVL/src/numpy.jl:67


Setup an LDA model and discretize using a single k-point and a small Ecut of 5 Hartree.

In [3]:
model = model_LDA(lattice, atoms)
basis = PlaneWaveBasis(model; Ecut=5, kgrid=(1, 1, 1))

Out[3]:
PlaneWaveBasis discretization:
Ecut                 : 5.0 Ha
fft_size             : (32, 32, 32)
kgrid type           : Monkhorst-Pack
kgrid                : [1, 1, 1]
num. irred. kpoints  : 1

Discretized Model(lda_xc_teter93, 3D):
lattice (in Bohr)    : [0         , 10.2631   , 10.2631   ]
[10.2631   , 0         , 10.2631   ]
[10.2631   , 10.2631   , 0         ]
unit cell volume     : 2162.1 Bohr³

atoms                : Si₁₆
atom potentials      : ElementPsp(Si, psp=hgh/lda/si-q4)

num. electrons       : 64
spin polarization    : none
temperature          : 0 Ha

terms                : Kinetic()
AtomicLocal()
AtomicNonlocal()
Ewald()
PspCorrection()
Hartree()
Xc(:lda_xc_teter93)

Find the ground state using direct minimisation (always using SCF is boring ...)

In [4]:
scfres = direct_minimization(basis, tol=1e-5);

Iter     Function value   Gradient norm
0     1.120257e+02     1.630944e+00
* time: 0.6009058952331543
1     1.093380e+01     9.061809e-01
* time: 1.6579689979553223
2    -1.227509e+01     1.051355e+00
* time: 2.1529269218444824
3    -3.456244e+01     8.542085e-01
* time: 2.915565013885498
4    -4.833683e+01     5.669814e-01
* time: 3.634446859359741
5    -5.723304e+01     1.961135e-01
* time: 4.354081869125366
6    -5.996818e+01     9.536322e-02
* time: 4.851265907287598
7    -6.091622e+01     5.340587e-02
* time: 5.356735944747925
8    -6.132444e+01     6.637419e-02
* time: 5.843098878860474
9    -6.162853e+01     4.036974e-02
* time: 6.328920841217041
10    -6.183049e+01     2.926730e-02
* time: 6.824652910232544
11    -6.196628e+01     1.987717e-02
* time: 7.326385021209717
12    -6.203799e+01     1.751728e-02
* time: 7.811548948287964
13    -6.209274e+01     1.755538e-02
* time: 8.296856880187988
14    -6.212528e+01     1.441766e-02
* time: 8.78230905532837
15    -6.215178e+01     1.062990e-02
* time: 9.268942832946777
16    -6.217004e+01     1.026420e-02
* time: 9.763037919998169
17    -6.218182e+01     8.613336e-03
* time: 10.250105857849121
18    -6.219045e+01     7.332046e-03
* time: 10.73715090751648
19    -6.219641e+01     5.332251e-03
* time: 11.224604845046997
20    -6.220111e+01     5.394773e-03
* time: 11.71941590309143
21    -6.220525e+01     5.739172e-03
* time: 12.205296993255615
22    -6.220925e+01     5.741321e-03
* time: 12.69203495979309
23    -6.221301e+01     5.688421e-03
* time: 13.178782939910889
24    -6.221643e+01     5.174043e-03
* time: 13.665106058120728
25    -6.221969e+01     5.694541e-03
* time: 14.158186912536621
26    -6.222334e+01     6.432898e-03
* time: 14.64242696762085
27    -6.222793e+01     5.835084e-03
* time: 15.128541946411133
28    -6.223378e+01     7.020916e-03
* time: 15.614642858505249
29    -6.224111e+01     6.552225e-03
* time: 16.10798192024231
30    -6.224787e+01     4.763569e-03
* time: 16.594491958618164
31    -6.225253e+01     4.200039e-03
* time: 17.08444905281067
32    -6.225614e+01     3.161248e-03
* time: 17.571759939193726
33    -6.225829e+01     2.580060e-03
* time: 18.058218955993652
34    -6.225974e+01     2.135120e-03
* time: 18.552178859710693
35    -6.226061e+01     1.228211e-03
* time: 19.03783392906189
36    -6.226104e+01     1.081224e-03
* time: 19.525140047073364
37    -6.226124e+01     8.351469e-04
* time: 20.01089596748352
38    -6.226137e+01     7.456877e-04
* time: 20.50429105758667
39    -6.226146e+01     6.431108e-04
* time: 20.99058985710144
40    -6.226153e+01     5.879519e-04
* time: 21.476955890655518
41    -6.226159e+01     4.476486e-04
* time: 21.963552951812744
42    -6.226162e+01     3.205509e-04
* time: 22.457913875579834
43    -6.226164e+01     2.561542e-04
* time: 22.94336199760437
44    -6.226165e+01     2.630813e-04
* time: 23.43017601966858
45    -6.226166e+01     1.249071e-04
* time: 23.915698051452637
46    -6.226166e+01     9.041320e-05
* time: 24.402217864990234
47    -6.226166e+01     6.602975e-05
* time: 24.89622688293457
48    -6.226167e+01     4.144375e-05
* time: 25.3821439743042
49    -6.226167e+01     3.945277e-05
* time: 25.86735200881958
50    -6.226167e+01     3.160669e-05
* time: 26.353015899658203
51    -6.226167e+01     2.497509e-05
* time: 26.84541392326355
52    -6.226167e+01     2.233248e-05
* time: 27.32987904548645
53    -6.226167e+01     1.756157e-05
* time: 27.814575910568237
54    -6.226167e+01     1.340228e-05
* time: 28.299021005630493
55    -6.226167e+01     1.184454e-05
* time: 28.784889936447144
56    -6.226167e+01     6.953365e-06
* time: 29.27808690071106
57    -6.226167e+01     4.387291e-06
* time: 29.762773036956787
58    -6.226167e+01     2.944737e-06
* time: 30.24839496612549

In [5]:
scfres.energies

Out[5]:
Energy breakdown (in Ha):
Kinetic             25.7671076
AtomicLocal         -18.8557750
AtomicNonlocal      14.8522694
Ewald               -67.1831486
PspCorrection       -2.3569765
Hartree             4.8485396
Xc                  -19.3336829

total               -62.261666459050