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])
lattice = load_lattice(pystruct)
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