# 21.1. pbc.gto — Crystal cell structure¶

This module provides functions to setup the basic information of a PBC calculation. The pyscf.pbc.gto module is analogous to the basic molecular pyscf.gto module. The Cell class for crystal structure unit cells is defined in this module and is analogous to the basic molecular Mole class. Among other details, the basis set and pseudopotentials are parsed in this module.

## 21.1.1. Cell class¶

The Cell class is defined as an extension of the molecular pyscf.gto.mole.Mole class. The Cell object offers much of the same functionality as the Mole object. For example, one can use the Cell object to access the atomic structure, basis functions, pseudopotentials, and certain analytical periodic integrals.

Similar to the input in a molecular calculation, one first creates a Cell object. After assigning the crystal parameters, one calls build() to fully initialize the Cell object. A shortcut function M() is available at the module level to simplify the input.

#!/usr/bin/env python

import numpy
import pyscf.lib
from pyscf.pbc import gto

#
# Simliar to the initialization of "Mole" object, here we need create a "Cell"
# object for periodic boundary systems.
#
cell = gto.Cell()
cell.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'''
cell.basis = 'gth-szv'
#
# Note the extra attribute ".a" in the "cell" initialization.
# .a is a matrix for lattice vectors.  Each row of .a is a primitive vector.
#
cell.a = numpy.eye(3)*3.5668
cell.build()

#
# pbc.gto module provided a shortcut initialization function "gto.M", like the
# one of finite size problem
#
cell = gto.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 = 'gth-szv',
a = numpy.eye(3)*3.5668)



Beyond the basic parameters atom and basis, one needs to set the unit cell lattice vectors a (a 3x3 array, where each row is a real-space primitive vector) and the numbers of grid points in the FFT-mesh in each positive direction gs (a length-3 list or 1x3 array); the total number of grid points is 2 gs +1.

In certain cases, it is convenient to choose the FFT-mesh density based on the kinetic energy cutoff. The Cell class offers an alternative attribute ke_cutoff that can be used to set the FFT-mesh. If ke_cutoff is set and gs is None, the Cell initialization function will convert the ke_cutoff to the equivalent FFT-mesh according to the relation $$\mathbf{g} = \frac{\sqrt{2E_{\mathrm{cut}}}}{2\pi}\mathbf{a}^T$$ and will overwrite the gs attribute.

Many PBC calculations are best performed using pseudopotentials, which are set via the pseudo attribute. Pseudopotentials alleviate the need for impractically dense FFT-meshes, although they represent a potentially uncontrolled source of error. See Pseudo potential for further details and a list of available pseudopotentials.

The input parameters .a and .pseudo are immutable in the Cell object. We emphasize that the input format might be different from the internal format used by PySCF. Similar to the convention in Mole, an internal Python data layer is created to hold the formatted .a and .pseudo parameters used as input.

_pseudo
The internal format to hold PBC pseudo potential parameters. It is represented with nested Python lists only.

Nuclear-nuclear interaction energies are evaluated by means of Ewald summation, which depends on three parameters: the truncation radius for real-space lattice sums rcut, the Gaussian model charge ew_eta, and the energy cutoff ew_cut. Although they can be set manually, these parameters are by default chosen automatically according to the attribute precision, which likewise can be set manually or left to its default value.

Besides the methods and parameters provided by Mole class (see Chapter gto — Molecular structure and GTO basis), there are some parameters frequently used in the code to access the information of the crystal.

kpts

The scaled or absolute k-points (nkpts x 3 array). This variable is not held as an attribute in Cell object; instead, the Cell object provides functions to generate the k-points and convert the k-points between the scaled (fractional) value and absolute value:

# Generate k-points
n_kpts_each_direction = [2,2,2]
abs_kpts = cell.make_kpts(n_kpts_each_direction)

# Convert k-points between two convention, the scaled and the absoulte values
scaled_kpts = cell.get_scaled_kpts(abs_kpts)
abs_kpts = cell.get_abs_kpts(scaled_kpts)

Gv
The (N x 3) array of plane waves associated to gs. gs defines the number of FFT grids in each direction. Cell.Gv() or get_Gv() convert the FFT-mesh to the plane waves. Gv are the the plane wave bases of 3D-FFT transformation. Given gs = [nx,ny,nz], the number of vectors in Gv is (2*nx+1)*(2*ny+1)*(2*nz+1).
vol
Cell.vol gives the volume of the unit cell (in atomic unit).
reciprocal_vectors
A 3x3 array. Each row is a reciprocal space primitive vector.
energy_nuc
Similar to the energy_nuc() provided by Mole class, this function also return the energy associated to the nuclear repulsion. The nuclear repulsion energy is computed with Ewald summation technique. The background contribution is removed from the nuclear repulsion energy otherwise this term is divergent.
pbc_intor
PBC analytic integral driver. It allows user to compute the PBC integral array in bulk, for given integral descriptor intor (see also Mole.intor() function moleintor). In the Cell object, we didn’t overload the intor() method. So one can access both the periodic integrals and free-boundary integrals within the Cell object. It allows you to input the cell object into the molecule program to run the free-boundary calculation (see Connection to Mole class).

Note

pbc_intor() does not support Coulomb type integrals. Calling pbc_intor with Coulomb type integral descriptor such as cint1e_nuc_sph leads to divergent integrals. The Coulomb type PBC integrals should be evaluated with density fitting technique (see Chapter pbc.df — PBC denisty fitting).

### 21.1.1.1. Attributes and methods¶

class pyscf.pbc.gto.Cell(**kwargs)[source]

A Cell object holds the basic information of a crystal.

Attributes:
a
: (3,3) ndarray
Lattice primitive vectors. Each row represents a lattice vector Reciprocal lattice vectors are given by b1,b2,b3 = 2 pi inv(a).T
mesh
: (3,) list of ints
The number G-vectors along each direction. The default value is estimated based on precision
pseudo
: dict or str
To define pseudopotential.
precision
: float
To control Ewald sums and lattice sums accuracy
rcut
: float
Cutoff radius (unit Bohr) in lattice summation. The default value is estimated based on the required precision.
ke_cutoff
: float
If set, defines a spherical cutoff of planewaves, with .5 * G**2 < ke_cutoff The default value is estimated based on precision
dimension
: int
Default is 3

** Following attributes (for experts) are automatically generated. **

ew_eta, ew_cut
: float
The Ewald ‘eta’ and ‘cut’ parameters. See get_ewald_params()

(See other attributes in Mole)

Examples:

>>> mol = Mole(atom='H^2 0 0 0; H 0 0 1.1', basis='sto3g')
>>> cl = Cell()
>>> cl.build(a='3 0 0; 0 3 0; 0 0 3', atom='C 1 1 1', basis='sto3g')
>>> print(cl.atom_symbol(0))
C

bas_rcut(cell, bas_id, precision=1e-08)

Estimate the largest distance between the function and its image to reach the precision in overlap

precision ~ int g(r-0) g(r-R)

build(dump_input=True, parse_arg=True, a=None, mesh=None, ke_cutoff=None, precision=None, nimgs=None, ew_eta=None, ew_cut=None, pseudo=None, basis=None, h=None, dimension=None, rcut=None, ecp=None, low_dim_ft_type=None, *args, **kwargs)[source]

Setup Mole molecule and Cell and initialize some control parameters. Whenever you change the value of the attributes of Cell, you need call this function to refresh the internal data of Cell.

Kwargs:
a
: (3,3) ndarray
The real-space unit cell lattice vectors. Each row represents a lattice vector.
mesh
: (3,) ndarray of ints
The number of positive G-vectors along each direction.
pseudo
: dict or str
To define pseudopotential. If given, overwrite Cell.pseudo
dumps(cell)

Serialize Cell object to a JSON formatted str.

energy_nuc(cell, ew_eta=None, ew_cut=None)

Perform real (R) and reciprocal (G) space Ewald sum for the energy.

Formulation of Martin, App. F2.

Returns:
float
The Ewald energy consisting of overlap, self, and G-space sum.
pyscf.pbc.gto.get_ewald_params
eval_gto(eval_name, coords, comp=None, kpts=None, kpt=None, shls_slice=None, non0tab=None, ao_loc=None, out=None)[source]

Evaluate PBC-AO function value on the given grids,

Args:

eval_name : str

Function Expression
“GTOval_sph” sum_T exp(ik*T) |AO>
“GTOval_ip_sph” nabla sum_T exp(ik*T) |AO>
“GTOval_cart” sum_T exp(ik*T) |AO>
“GTOval_ip_cart” nabla sum_T exp(ik*T) |AO>
atm
: int32 ndarray
libcint integral function argument
bas
: int32 ndarray
libcint integral function argument
env
: float64 ndarray
libcint integral function argument
coords
: 2D array, shape (N,3)
The coordinates of the grids.
Kwargs:
shls_slice
: 2-element list
(shl_start, shl_end). If given, only part of AOs (shl_start <= shell_id < shl_end) are evaluated. By default, all shells defined in cell will be evaluated.
non0tab
: 2D bool array
mask array to indicate whether the AO values are zero. The mask array can be obtained by calling dft.gen_grid.make_mask()
out
: ndarray
If provided, results are written into this array.
Returns:
A list of 2D (or 3D) arrays to hold the AO values on grids. Each element of the list corresponds to a k-point and it has the shape (N,nao) Or shape (*,N,nao).

Examples:

>>> cell = pbc.gto.M(a=numpy.eye(3)*4, atom='He 1 1 1', basis='6-31g')
>>> coords = cell.get_uniform_grids([20,20,20])
>>> kpts = cell.make_kpts([3,3,3])
>>> ao_value = cell.pbc_eval_gto("GTOval_sph", coords, kpts)
>>> len(ao_value)
27
>>> ao_value[0].shape
(100, 2)
>>> ao_value = cell.pbc_eval_gto("GTOval_ig_sph", coords, kpts, comp=3)
>>> print(ao_value.shape)
>>> len(ao_value)
27
>>> ao_value[0].shape
(3, 100, 2)

ewald(cell, ew_eta=None, ew_cut=None)

Perform real (R) and reciprocal (G) space Ewald sum for the energy.

Formulation of Martin, App. F2.

Returns:
float
The Ewald energy consisting of overlap, self, and G-space sum.
pyscf.pbc.gto.get_ewald_params
format_basis(basis_tab)[source]

Convert the input Cell.basis to the internal data format:

{ atom: (l, kappa, ((-exp, c_1, c_2, ..), nprim, nctr, ptr-exps, ptr-contraction-coeff)), ... }

Args:
basis_tab
: dict
Similar to Cell.basis, it cannot be a str
Returns:
Formated basis

Examples:

>>> pbc.format_basis({'H':'gth-szv'})
{'H': [[0,
(8.3744350009, -0.0283380461),
(1.8058681460, -0.1333810052),
(0.4852528328, -0.3995676063),
(0.1658236932, -0.5531027541)]]}

format_pseudo(pseudo_tab)[source]

Convert the input Cell.pseudo (dict) to the internal data format:

{ atom: ( (nelec_s, nele_p, nelec_d, ...),
rloc, nexp, (cexp_1, cexp_2, ..., cexp_nexp),
nproj_types,
(r1, nproj1, ( (hproj1[1,1], hproj1[1,2], ..., hproj1[1,nproj1]),
(hproj1[2,1], hproj1[2,2], ..., hproj1[2,nproj1]),
...
(hproj1[nproj1,1], hproj1[nproj1,2], ...        ) )),
(r2, nproj2, ( (hproj2[1,1], hproj2[1,2], ..., hproj2[1,nproj1]),
... ) )
)
... }

Args:
pseudo_tab
: dict
Similar to Cell.pseudo (a dict), it cannot be a str
Returns:
Formatted pseudo

Examples:

>>> pbc.format_pseudo({'H':'gth-blyp', 'He': 'gth-pade'})
{'H': [[1],
0.2, 2, [-4.19596147, 0.73049821], 0],
'He': [[2],
0.2, 2, [-9.1120234, 1.69836797], 0]}

from_ase(ase_atom)[source]

Update cell based on given ase atom object

Examples:

>>> from ase.lattice import bulk
>>> cell.from_ase(bulk('C', 'diamond', a=LATTICE_CONST))

gen_uniform_grids(cell, mesh=None, **kwargs)

Generate a uniform real-space grid consistent w/ samp thm; see MH (3.19).

Args:
cell : instance of Cell
Returns:
coords
: (ngx*ngy*ngz, 3) ndarray
The real-space grid point coordinates.
get_Gv(cell, mesh=None, **kwargs)

Calculate three-dimensional G-vectors for the cell; see MH (3.8).

Indices along each direction go as [0...N-1, -N...-1] to follow FFT convention.

Args:
cell : instance of Cell
Returns:
Gv
: (ngrids, 3) ndarray of floats
The array of G-vectors.
get_Gv_weights(cell, mesh=None, **kwargs)

Calculate G-vectors and weights.

Returns:
Gv
: (ngris, 3) ndarray of floats
The array of G-vectors.
get_SI(cell, Gv=None)

Calculate the structure factor (0D, 1D, 2D, 3D) for all atoms; see MH (3.34).

Args:

cell : instance of Cell

Gv
: (N,3) array
G vectors
Returns:
SI
: (natm, ngrids) ndarray, dtype=np.complex128
The structure factor for each atom at each G-vector.
get_abs_kpts(scaled_kpts)[source]

Get absolute k-points (in 1/Bohr), given “scaled” k-points in fractions of lattice vectors.

Args:
scaled_kpts : (nkpts, 3) ndarray of floats
Returns:
abs_kpts : (nkpts, 3) ndarray of floats
get_bounding_sphere(cell, rcut)

Finds all the lattice points within a sphere of radius rcut.

Defines a parallelipiped given by -N_x <= n_x <= N_x, with x in [1,3] See Martin p. 85

Args:
rcut
: number
real space cut-off for interaction
Returns:
cut : ndarray of 3 ints defining N_x
get_ewald_params(cell, precision=1e-08, mesh=None)

Choose a reasonable value of Ewald ‘eta’ and ‘cut’ parameters. eta^2 is the exponent coefficient of the model Gaussian charge for nucleus at R: frac{eta^3}{pi^1.5} e^{-eta^2 (r-R)^2}

Choice is based on largest G vector and desired relative precision.

The relative error in the G-space sum is given by

precision ~ 4pi Gmax^2 e^{(-Gmax^2)/(4 eta^2)}

which determines eta. Then, real-space cutoff is determined by (exp. factors only)

precision ~ erfc(eta*rcut) / rcut ~ e^{(-eta**2 rcut*2)}
Returns:
ew_eta, ew_cut
: float
The Ewald ‘eta’ and ‘cut’ parameters.
get_kpts(cell, nks, wrap_around=False, with_gamma_point=True, scaled_center=None)

Given number of kpoints along x,y,z , generate kpoints

Args:
nks : (3,) ndarray
Kwargs:
wrap_around
: bool
To ensure all kpts are in first Brillouin zone.
with_gamma_point
: bool
Whether to shift Monkhorst-pack grid to include gamma-point.
scaled_center
: (3,) array
Shift all points in the Monkhorst-pack grid to be centered on scaled_center, given as the zeroth index of the returned kpts. Scaled meaning that the k-points are scaled to a grid from [-1,1] x [-1,1] x [-1,1]
Returns:
kpts in absolute value (unit 1/Bohr). Gamma point is placed at the first place in the k-points list

Examples:

>>> cell.make_kpts((4,4,4))

get_lattice_Ls(cell, nimgs=None, rcut=None, dimension=None)

Get the (Cartesian, unitful) lattice translation vectors for nearby images. The translation vectors can be used for the lattice summation.

get_nimgs(cell, precision=None)

Choose number of basis function images in lattice sums to include for given precision in overlap, using

precision ~ int r^l e^{-alpha r^2} (r-rcut)^l e^{-alpha (r-rcut)^2} ~ (rcut^2/(2alpha))^l e^{alpha/2 rcut^2}

where alpha is the smallest exponent in the basis. Note that assumes an isolated exponent in the middle of the box, so it adds one additional lattice vector to be safe.

get_scaled_kpts(abs_kpts)[source]

Get scaled k-points, given absolute k-points in 1/Bohr.

Args:
abs_kpts : (nkpts, 3) ndarray of floats
Returns:
scaled_kpts : (nkpts, 3) ndarray of floats
get_uniform_grids(cell, mesh=None, **kwargs)

Generate a uniform real-space grid consistent w/ samp thm; see MH (3.19).

Args:
cell : instance of Cell
Returns:
coords
: (ngx*ngy*ngz, 3) ndarray
The real-space grid point coordinates.
has_ecp()[source]

Whether pseudo potential is used in the system.

kernel(dump_input=True, parse_arg=True, a=None, mesh=None, ke_cutoff=None, precision=None, nimgs=None, ew_eta=None, ew_cut=None, pseudo=None, basis=None, h=None, dimension=None, rcut=None, ecp=None, low_dim_ft_type=None, *args, **kwargs)

Setup Mole molecule and Cell and initialize some control parameters. Whenever you change the value of the attributes of Cell, you need call this function to refresh the internal data of Cell.

Kwargs:
a
: (3,3) ndarray
The real-space unit cell lattice vectors. Each row represents a lattice vector.
mesh
: (3,) ndarray of ints
The number of positive G-vectors along each direction.
pseudo
: dict or str
To define pseudopotential. If given, overwrite Cell.pseudo
lattice_vectors()[source]

Convert the primitive lattice vectors.

Return 3x3 array in which each row represents one direction of the lattice vectors (unit in Bohr)

loads(molstr)[source]

Deserialize a str containing a JSON document to a Cell object.

make_kpts(cell, nks, wrap_around=False, with_gamma_point=True, scaled_center=None)

Given number of kpoints along x,y,z , generate kpoints

Args:
nks : (3,) ndarray
Kwargs:
wrap_around
: bool
To ensure all kpts are in first Brillouin zone.
with_gamma_point
: bool
Whether to shift Monkhorst-pack grid to include gamma-point.
scaled_center
: (3,) array
Shift all points in the Monkhorst-pack grid to be centered on scaled_center, given as the zeroth index of the returned kpts. Scaled meaning that the k-points are scaled to a grid from [-1,1] x [-1,1] x [-1,1]
Returns:
kpts in absolute value (unit 1/Bohr). Gamma point is placed at the first place in the k-points list

Examples:

>>> cell.make_kpts((4,4,4))

pack(cell)

Pack the input args of Cell to a dict, which can be serialized with pickle

pbc_eval_gto(eval_name, coords, comp=None, kpts=None, kpt=None, shls_slice=None, non0tab=None, ao_loc=None, out=None)[source]

Evaluate PBC-AO function value on the given grids,

Args:

eval_name : str

Function Expression
“GTOval_sph” sum_T exp(ik*T) |AO>
“GTOval_ip_sph” nabla sum_T exp(ik*T) |AO>
“GTOval_cart” sum_T exp(ik*T) |AO>
“GTOval_ip_cart” nabla sum_T exp(ik*T) |AO>
atm
: int32 ndarray
libcint integral function argument
bas
: int32 ndarray
libcint integral function argument
env
: float64 ndarray
libcint integral function argument
coords
: 2D array, shape (N,3)
The coordinates of the grids.
Kwargs:
shls_slice
: 2-element list
(shl_start, shl_end). If given, only part of AOs (shl_start <= shell_id < shl_end) are evaluated. By default, all shells defined in cell will be evaluated.
non0tab
: 2D bool array
mask array to indicate whether the AO values are zero. The mask array can be obtained by calling dft.gen_grid.make_mask()
out
: ndarray
If provided, results are written into this array.
Returns:
A list of 2D (or 3D) arrays to hold the AO values on grids. Each element of the list corresponds to a k-point and it has the shape (N,nao) Or shape (*,N,nao).

Examples:

>>> cell = pbc.gto.M(a=numpy.eye(3)*4, atom='He 1 1 1', basis='6-31g')
>>> coords = cell.get_uniform_grids([20,20,20])
>>> kpts = cell.make_kpts([3,3,3])
>>> ao_value = cell.pbc_eval_gto("GTOval_sph", coords, kpts)
>>> len(ao_value)
27
>>> ao_value[0].shape
(100, 2)
>>> ao_value = cell.pbc_eval_gto("GTOval_ig_sph", coords, kpts, comp=3)
>>> print(ao_value.shape)
>>> len(ao_value)
27
>>> ao_value[0].shape
(3, 100, 2)

pbc_intor(intor, comp=None, hermi=0, kpts=None, kpt=None, shls_slice=None, **kwargs)[source]

One-electron integrals with PBC.

$\sum_T \int \mu(r) * [intor] * \nu (r-T) dr$

reciprocal_vectors(norm_to=6.283185307179586)[source]
\begin{split}\begin{align} \mathbf{b_1} &= 2\pi \frac{\mathbf{a_2} \times \mathbf{a_3}}{\mathbf{a_1} \cdot (\mathbf{a_2} \times \mathbf{a_3})} \\ \mathbf{b_2} &= 2\pi \frac{\mathbf{a_3} \times \mathbf{a_1}}{\mathbf{a_2} \cdot (\mathbf{a_3} \times \mathbf{a_1})} \\ \mathbf{b_3} &= 2\pi \frac{\mathbf{a_1} \times \mathbf{a_2}}{\mathbf{a_3} \cdot (\mathbf{a_1} \times \mathbf{a_2})} \end{align}\end{split}
to_mol()[source]

Return a Mole object using the same atoms and basis functions as the Cell object.

tot_electrons(cell, nkpts=1)

Total number of electrons

unpack(moldic)[source]

Convert the packed dict to a Cell object, to generate the input arguments for Cell object.

### 21.1.1.2. Connection to Mole class¶

Cell class is compatible with the molecule pyscf.gto.mole.Mole class. They shared most data structure and methods. It gives the freedom to mix the finite size calculation and the PBC calculation. If you feed the cell object to molecule module/functions, the molecule program will not check whether the given Mole object is the true Mole or not. It simply treats the Cell object as the Mole object and run the finite size calculations. Because the same module names were used in PBC program and molecule program, you should be careful with the imported modules since no error message will be raised if you by mistake input the Cell object into the molecule program.

Although we reserve the flexibility to mix the Cell and Mole objects in the same code, it should be noted that the serialization methods of the two objects are not completely compatible. When you dumps/loads the cell object in the molecule program, informations of the Cell object or the faked Mole object may be lost.

### 21.1.1.3. Serialization¶

Cell class has two set of functions to serialize Cell object in different formats.

• JSON format is the default serialization format used by pyscf.lib.chkfile module. It can be serialized by Cell.dumps() function and deserialized by Cell.loads() function.
• In the old version, Mole.pack() and Mole.unpack() functions are used to convert the Mole object to and from Python dict. The Python dict is then serialized by pickle module. This serialization method is not used anymore in the new PySCF code. To keep the backward compatibility, the two methods are defined in Cell class.

## 21.1.2. Basis set¶

The pbc module supports all-electron calculation. The all-electron basis sets developed by quantum chemistry community can be directly used in the pbc calculation. The Cell class supports to mix the QC all-electron basis and PBC basis in the same calculation.

#!/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
''')},
a = numpy.eye(3)*3.5668)


Note

The default PBC Coulomb type integrals are computed using FFT transformation. If the all-electron basis are used, you might need very high energy cutoff to converge the integrals. It is recommended to use mixed density fitting technique (pbc.df — PBC denisty fitting) to handle the all-electron calculations.

## 21.1.3. Pseudo potential¶

Quantum chemistry community developed a wide range of pseudo potentials (which are called ECP, effective core potential) for heavy elements. ECP works quite successful in finite system. It has high flexibility to choose different core size and relevant basis sets to satisfy different requirements on accuracy, efficiency in different simulation scenario. Extending ECP to PBC code enriches the pseudo potential database. PySCF PBC program supports both the PBC conventional pseudo potential and ECP and the mix of the two kinds of potentials in the same calculation.

#!/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.

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'})