Crystal structure

Modules: pbc.gto

Initializing a crystal

Initializing a crystal unit cell is very similar to the initialization of a molecule. Instead of the gto.Mole class, one uses the pbc.gto.Cell class to define a cell:

>>> from pyscf.pbc import gto
>>> cell = gto.Cell()
>>> cell.atom = '''H  0 0 0; H 1 1 1'''
>>> cell.basis = 'gth-dzvp'
>>> cell.pseudo = 'gth-pade'
>>> cell.a = numpy.eye(3) * 2
>>> cell.build()

The other two equivalent ways to initialize a molecule introduced in Molecular structure also apply here:

>>> from pyscf.pbc import gto
>>> cell = gto.Cell()
>>> cell.build(
...   atom = '''H  0 0 0; H 1 1 1''',
...   basis = 'gth-dzvp',
...   pseudo = 'gth-pade',
...   a = numpy.eye(3) * 2)

>>> import pyscf
>>> cell = pyscf.M(
...   atom = '''H  0 0 0; H 1 1 1''',
...   basis = 'gth-dzvp',
...   pseudo = 'gth-pade',
...   a = numpy.eye(3) * 2)

>>> from pyscf.pbc import gto
>>> cell = gto.M(
...   atom = '''H  0 0 0; H 1 1 1''',
...   basis = 'gth-dzvp',
...   pseudo = 'gth-pade',
...   a = numpy.eye(3) * 2)

Lattice vectors

The crystal initialization requires an additional parameter a, which represents the lattice vectors. The format of a is array-like:

cell.a = numpy.eye(3) * 2
cell.a = [[2,0,0],[0,2,0],[0,0,2]]

Each row of the 3-by-3 matrix of a represents a lattice vector in Cartesian coordinate, with the same unit as the input atom parameter.

Note

The input lattice vectors (row by row) should form a right-handed coordinate system, as otherwise some integrals may be computed incorrectly in PySCF.

Basis set and pseudopotential

PySCF uses crystalline Gaussian-type orbitals as basis functions for solid calculations. The predefined basis sets and ECPs for molecular calculations can be used in solid calculations as well. In addition, the predefined basis sets include valence basis sets that are optimized for the GTH pseudopotentials (a whole list can be found in pyscf/pbc/gto/basis and pyscf/pbc/gto/pseudo). The input format of basis sets for the Cell object is the same as that for the Mole object.

#!/usr/bin/env python

'''
Input pseudo potential using functions pbc.gto.pseudo.parse and pbc.gto.pseudo.load

It is allowed to mix the Quantum chemistry effective core potentail (ECP) with
crystal pseudo potential (PP).  Input ECP with .ecp attribute and PP with
.pseudo attribute.

See also
pyscf/pbc/gto/pseudo/GTH_POTENTIALS for the GTH-potential format
pyscf/examples/gto/05-input_ecp.py for quantum chemistry ECP format
'''

import numpy
from pyscf.pbc import gto

cell = gto.M(atom='''
Si1 0 0 0
Si2 1 1 1''',
             a = '''3    0    0
                    0    3    0
                    0    0    3''',
             basis = {'Si1': 'gth-szv',  # Goedecker, Teter and Hutter single zeta basis
                      'Si2': 'lanl2dz'},
             pseudo = {'Si1': gto.pseudo.parse('''
Si
    2    2
     0.44000000    1    -6.25958674
    2
     0.44465247    2     8.31460936    -2.33277947
                                        3.01160535
     0.50279207    1     2.33241791
''')},
             ecp = {'Si2': 'lanl2dz'},  # ECP for second Si atom
            )

#
# Some elements have multiple PP definitions in GTH database.  Add suffix in
# the basis name to load the specific PP.
#
cell = gto.M(
    a = numpy.eye(3)*5,
    atom = 'Mg1 0 0 0; Mg2 0 0 1',
    pseudo = {'Mg1': 'gth-lda-q2', 'Mg2': 'gth-lda-q10'})

#
# Allow mixing quantum chemistry ECP (or BFD PP) and crystal PP in the same calculation.
#
cell = gto.M(
    a = '''4    0    0
           0    4    0
           0    0    4''',
    atom = 'Cl 0 0 1; Na 0 1 0',
    basis = {'na': 'gth-szv', 'Cl': 'bfd-vdz'},
    ecp = {'Cl': 'bfd-pp'},
    pseudo = {'Na': 'gthbp'})

#
# ECP can be input in the attribute .pseudo
#
cell = gto.M(
    a = '''4    0    0
           0    4    0
           0    0    4''',
    atom = 'Cl 0 0 1; Na 0 1 0',
    basis = {'na': 'gth-szv', 'Cl': 'bfd-vdz'},
    pseudo = {'Na': 'gthbp', 'Cl': 'bfd-pp'})

Finally, custom basis sets can be defined just like in molecular calculations

#!/usr/bin/env python

'''
Basis can be input the same way as the finite-size system.
'''

#
# Note pbc.gto.parse does not support NWChem format.  To parse NWChem format
# basis string, you need the molecule gto.parse function.
#

import numpy
from pyscf import gto
from pyscf.pbc import gto as pgto
cell = pgto.M(
    atom = '''C     0.      0.      0.    
              C     0.8917  0.8917  0.8917
              C     1.7834  1.7834  0.    
              C     2.6751  2.6751  0.8917
              C     1.7834  0.      1.7834
              C     2.6751  0.8917  2.6751
              C     0.      1.7834  1.7834
              C     0.8917  2.6751  2.6751''',
    basis = {'C': gto.parse('''
# Parse NWChem format basis string (see https://bse.pnl.gov/bse/portal).
# Comment lines are ignored
#BASIS SET: (6s,3p) -> [2s,1p]
O    S
    130.7093200              0.15432897       
     23.8088610              0.53532814       
      6.4436083              0.44463454       
O    SP
      5.0331513             -0.09996723             0.15591627       
      1.1695961              0.39951283             0.60768372       
      0.3803890              0.70011547             0.39195739       
                                ''')},
    pseudo = 'gth-pade',
    a = numpy.eye(3)*3.5668)

Low-dimensional systems

The PySCF pbc module also supports low-dimensional periodic systems. You can initialize the attribute Cell.dimension to specify the dimension of the system:

>>> cell.dimension = 2
>>> cell.a = numpy.eye(3) * 2
>>> cell.build()

When dimension is smaller than 3, a vacuum of infinite size will be applied in certain direction(s). For example, when dimension is 2, the z-direction will be treated as infinitely large and the xy-plane constitutes the periodic surface. When dimension is 1, the y and z axes are treated as vacuum and thus the system is a wire in the x direction. When dimension is 0, all three directions are treated as vacuum, and this is equivalent to a molecular system.

K-points

The k-points used in solid calculations can be obtained through the Cell.make_kpts() method. The minimal input is to specify the number of k-points in each lattice vector direction:

>>> kpts = cell.make_kpts([1,2,2])
>>> print(kpts.shape)
(4,3)

By default, this will return the shifted Monkhorst-Pack mesh which is centered at the \(Gamma\)-point. To get the non-shifted Monkhorst-Pack mesh, one can call:

>>> kpts = cell.make_kpts([1,2,2], with_gamma_point=False)

To get a shifted Monkhorst-pack mesh centered at a specific point, one can call:

>>> kpts = cell.make_kpts([1,2,2], scaled_center=[0.,0.25,0.25])

where scaled_center is defined in the units of lattice vectors.

The obtained k-points are used as input for crystalline electronic structure calculations. For example, k-point sampled RHF can be invoked as follows:

>>> from pyscf.pbc import scf
>>> kpts = cell.make_kpts([2,2,2])
>>> kmf = scf.KRHF(cell, kpts = kpts)
>>> kmf.kernel()

More details about k-point sampling for each method can be found in the corresponding chapters.

Other parameters

The attribute precision of the Cell object determines the integral accuracy. The default value is 1e-8 hartree. Some other parameters are set automatically based on the value of precision, which include:

  • mesh – length-3 list or 1x3 array of int

    • The numbers of grid points in the FFT-mesh in each direction.

  • ke_cutoff – float

    • Kinetic energy cutoff of the plane waves in FFT-DF

  • rcut – float

    • Cutoff radius (in Bohr) of the lattice summation in the integral evaluation

The parameters above can also be set manually by the user.

The other attributes of the Mole class such as verbose, max_memory, spin, etc., have the same meanings in the Cell class.

Note

Currently, point group symmetries for crystals are not supported.

Accessing AO integrals

The Mole.intor() method can only be used to evaluate integrals with open boundary conditions. When the periodic boundary conditions of crystalline systems are used, one needs to use the pbc.Cell.pbc_intor() function to evaluate the integrals for short-range operators, such as the overlap and kinetic matrix:

overlap = cell.pbc_intor('int1e_ovlp')
kin = cell.pbc_intor('int1e_kin')

By default, the pbc.Cell.pbc_intor() function only returns the \(\Gamma\)-point integrals. If k-points are specified, it will evaluate the integrals at each k-point:

kpts = cell.make_kpts([2,2,2])  # 2x2x2 Monkhorst-pack mesh
overlap = cell.pbc_intor('int1e_ovlp', kpts=kpts)

Note

pbc.Cell.pbc_intor() can only be used to evaluate the short-range integrals. PBC density fitting methods have to be used to compute the integrals for long-range operators such as nuclear attraction and Coulomb repulsion.

The two-electron Coulomb integrals can be evaluated with the PBC density fitting methods:

from pyscf.pbc import df
eri = df.DF(cell).get_eri()

See Density fitting for crystalline calculations for more details of the PBC density fitting methods.