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]];

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(1)
                               AtomicLocal()
                               AtomicNonlocal()
                               Ewald()
                               PspCorrection()
                               Hartree(1)
                               Xc([:lda_xc_teter93], 1, nothing)

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.119745e+02     1.534611e+00
 * time: 0.6858079433441162
     1     1.001989e+01     9.067000e-01
 * time: 1.9077880382537842
     2    -1.231277e+01     9.831490e-01
 * time: 2.4594340324401855
     3    -3.443302e+01     7.422740e-01
 * time: 3.296710968017578
     4    -4.795090e+01     5.380520e-01
 * time: 4.098114013671875
     5    -5.707327e+01     2.564140e-01
 * time: 4.9100000858306885
     6    -5.976870e+01     1.840169e-01
 * time: 5.471502065658569
     7    -6.080879e+01     5.363202e-02
 * time: 6.012686014175415
     8    -6.123132e+01     7.463058e-02
 * time: 6.555546045303345
     9    -6.150772e+01     6.634163e-02
 * time: 7.09138298034668
    10    -6.173787e+01     2.861154e-02
 * time: 7.646059036254883
    11    -6.192799e+01     2.190902e-02
 * time: 8.187505006790161
    12    -6.201083e+01     1.840368e-02
 * time: 8.728402137756348
    13    -6.208871e+01     1.830214e-02
 * time: 9.27440094947815
    14    -6.213202e+01     1.338024e-02
 * time: 9.814677000045776
    15    -6.216817e+01     1.165805e-02
 * time: 10.351267099380493
    16    -6.218615e+01     9.073019e-03
 * time: 10.891649007797241
    17    -6.219941e+01     9.727073e-03
 * time: 11.439162015914917
    18    -6.220877e+01     8.768039e-03
 * time: 11.978865146636963
    19    -6.221671e+01     7.544395e-03
 * time: 12.518681049346924
    20    -6.222363e+01     5.484427e-03
 * time: 13.057928085327148
    21    -6.223126e+01     5.991926e-03
 * time: 13.601854085922241
    22    -6.223941e+01     5.173361e-03
 * time: 14.139701128005981
    23    -6.224732e+01     4.641030e-03
 * time: 14.678477048873901
    24    -6.225340e+01     4.111441e-03
 * time: 15.219855070114136
    25    -6.225719e+01     3.559551e-03
 * time: 15.766541004180908
    26    -6.225913e+01     2.216527e-03
 * time: 16.306998014450073
    27    -6.226013e+01     1.443151e-03
 * time: 16.849133014678955
    28    -6.226064e+01     1.136548e-03
 * time: 17.39190101623535
    29    -6.226095e+01     1.126307e-03
 * time: 17.93880295753479
    30    -6.226117e+01     1.032155e-03
 * time: 18.477463006973267
    31    -6.226133e+01     9.955401e-04
 * time: 19.019612073898315
    32    -6.226144e+01     7.110819e-04
 * time: 19.56872797012329
    33    -6.226154e+01     5.018042e-04
 * time: 20.11106014251709
    34    -6.226160e+01     3.773442e-04
 * time: 20.65448307991028
    35    -6.226163e+01     2.736419e-04
 * time: 21.19935703277588
    36    -6.226165e+01     1.873975e-04
 * time: 21.75040292739868
    37    -6.226166e+01     1.508332e-04
 * time: 22.29786205291748
    38    -6.226166e+01     9.836398e-05
 * time: 22.84190607070923
    39    -6.226166e+01     7.053136e-05
 * time: 23.38502812385559
    40    -6.226167e+01     5.482125e-05
 * time: 23.93679714202881
    41    -6.226167e+01     4.691477e-05
 * time: 24.480637073516846
    42    -6.226167e+01     3.542361e-05
 * time: 25.02352213859558
    43    -6.226167e+01     2.930704e-05
 * time: 25.567615032196045
    44    -6.226167e+01     2.349085e-05
 * time: 26.116332054138184
    45    -6.226167e+01     1.926572e-05
 * time: 26.660706043243408
    46    -6.226167e+01     1.565035e-05
 * time: 27.204591035842896
    47    -6.226167e+01     1.035712e-05
 * time: 27.850672960281372
    48    -6.226167e+01     8.508926e-06
 * time: 28.509856939315796
    49    -6.226167e+01     5.721727e-06
 * time: 29.13272500038147
    50    -6.226167e+01     3.787734e-06
 * time: 29.77399206161499
    51    -6.226167e+01     2.304309e-06
 * time: 30.459808111190796
In [5]:
scfres.energies
Out[5]:
Energy breakdown (in Ha):
    Kinetic             25.7671062
    AtomicLocal         -18.8557671
    AtomicNonlocal      14.8522649
    Ewald               -67.1831486
    PspCorrection       -2.3569765
    Hartree             4.8485363 
    Xc                  -19.3336817

    total               -62.261666460927