10.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. Amongst other details, the basis set and pseudopotentials are parsed in this module.

10.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 of 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'
cell.pseudo = 'gth-pade'
#
# 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',
    pseudo = 'gth-pade',
    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 mesh (a length-3 list or 1x3 array).

In the Cell class, many parameters are determined automatically according to the attribute precision, which likewise can be set manually or left to its default value (1e-8). The parameters determined by precision include ke_cutoff, mesh, rcut, ew_cut and ew_eta. These parameters can also be set manually by the user.

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 mesh 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 mesh 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.

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 mesh. mesh 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 mesh = [nx,ny,nz], the number of vectors in Gv is nx*ny*nz.

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).

10.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

pseudodict or str

To define pseudopotential.

precisionfloat

To control Ewald sums and lattice sums accuracy

rcutfloat

Cutoff radius (unit Bohr) in lattice summation. The default value is estimated based on the required precision.

ke_cutofffloat

If set, defines a spherical cutoff of planewaves, with .5 * G**2 < ke_cutoff The default value is estimated based on precision

dimensionint

Periodic dimensions. Default is 3

low_dim_ft_typestr

For semi-empirical periodic systems, whether to calculate integrals at the non-PBC dimension using the sampled mesh grids in infinity vacuum (inf_vacuum) or truncated Coulomb potential (analytic_2d_1). Unless explicitly specified, analytic_2d_1 is used for 2D system and inf_vacuum is assumed for 1D and 0D.

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

ew_eta, ew_cutfloat

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
ao2mo(mo_coeffs, intor='int2e', erifile=None, dataname='eri_mo', **kwargs)[source]

Integral transformation for arbitrary orbitals and arbitrary integrals. See more detalied documentation in func:ao2mo.kernel.

Args:
mo_coeffs (an np array or a list of arrays)A matrix of orbital

coefficients if it is a numpy ndarray, or four sets of orbital coefficients, corresponding to the four indices of (ij|kl).

Kwargs:
erifile (str or h5py File or h5py Group object)The file/object

to store the transformed integrals. If not given, the return value is an array (in memory) of the transformed integrals.

datanamestr

Note this argument is effective if erifile is given. The dataset name in the erifile (ref the hierarchy of HDF5 format http://www.hdfgroup.org/HDF5/doc1.6/UG/09_Groups.html). By assigning different dataname, the existed integral file can be reused. If the erifile contains the specified dataname, the old integrals will be replaced by the new one under the key dataname.

intor (str)integral name Name of the 2-electron integral. Ref

to getints_by_shell() for the complete list of available 2-electron integral names

Returns:

An array of transformed integrals if erifile is not given. Otherwise, return the file/fileobject if erifile is assigned.

Examples:

>>> import pyscf
>>> mol = pyscf.M(atom='O 0 0 0; H 0 1 0; H 0 0 1', basis='sto3g')
>>> mo1 = numpy.random.random((mol.nao_nr(), 10))
>>> mo2 = numpy.random.random((mol.nao_nr(), 8))
>>> eri1 = mol.ao2mo(mo1)
>>> print(eri1.shape)
(55, 55)
>>> eri1 = mol.ao2mo(mo1, compact=False)
>>> print(eri1.shape)
(100, 100)
>>> eri1 = mol.ao2mo(eri, (mo1,mo2,mo2,mo2))
>>> print(eri1.shape)
(80, 36)
>>> eri1 = mol.ao2mo(eri, (mo1,mo2,mo2,mo2), erifile='water.h5')
bas_rcut(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-Rcut)

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.

ke_cutofffloat

If set, defines a spherical cutoff of planewaves, with .5 * G**2 < ke_cutoff The default value is estimated based on precision

precisionfloat

To control Ewald sums and lattice sums accuracy

nimgs(3,) ndarray of ints

Number of repeated images in lattice summation to produce periodicity. This value can be estimated based on the required precision. It’s recommended NOT making changes to this value.

rcutfloat

Cutoff radius (unit Bohr) in lattice summation to produce periodicity. The value can be estimated based on the required precision. It’s recommended NOT making changes to this value.

ew_eta, ew_cutfloat

Parameters eta and cut to converge Ewald summation. See get_ewald_params()

pseudodict or str

To define pseudopotential.

ecpdict or str

To define ECP type pseudopotential.

h(3,3) ndarray

a.T. Deprecated

dimensionint

Default is 3

low_dim_ft_typestr

For semi-empirical periodic systems, whether to calculate integrals at the non-PBC dimension using the sampled mesh grids in infinity vacuum (inf_vacuum) or truncated Coulomb potential (analytic_2d_1). Unless explicitly specified, analytic_2d_1 is used for 2D system and inf_vacuum is assumed for 1D and 0D.

copy()[source]

Deepcopy of the given Mole object

dumps()

Serialize Cell object to a JSON formatted str.

energy_nuc(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.

See Also:

pyscf.pbc.gto.get_ewald_params

eval_ao(eval_name, coords, comp=None, kpts=None, kpt=None, shls_slice=None, non0tab=None, ao_loc=None, out=None)

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>

atmint32 ndarray

libcint integral function argument

basint32 ndarray

libcint integral function argument

envfloat64 ndarray

libcint integral function argument

coords2D array, shape (N,3)

The coordinates of the grids.

Kwargs:
shls_slice2-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.

non0tab2D 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()

outndarray

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

atmint32 ndarray

libcint integral function argument

basint32 ndarray

libcint integral function argument

envfloat64 ndarray

libcint integral function argument

coords2D array, shape (N,3)

The coordinates of the grids.

Kwargs:
shls_slice2-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.

non0tab2D 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()

outndarray

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(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.

See Also:

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_tabdict

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_tabdict

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(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(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(mesh=None, **kwargs)

Calculate G-vectors and weights.

Returns:
Gv(ngris, 3) ndarray of floats

The array of G-vectors.

get_SI(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(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:
rcutnumber

real space cut-off for interaction

Returns:

cut : ndarray of 3 ints defining N_x

get_ewald_params(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_cutfloat

The Ewald ‘eta’ and ‘cut’ parameters.

get_kpts(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_aroundbool

To ensure all kpts are in first Brillouin zone.

with_gamma_pointbool

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(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(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(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.

ke_cutofffloat

If set, defines a spherical cutoff of planewaves, with .5 * G**2 < ke_cutoff The default value is estimated based on precision

precisionfloat

To control Ewald sums and lattice sums accuracy

nimgs(3,) ndarray of ints

Number of repeated images in lattice summation to produce periodicity. This value can be estimated based on the required precision. It’s recommended NOT making changes to this value.

rcutfloat

Cutoff radius (unit Bohr) in lattice summation to produce periodicity. The value can be estimated based on the required precision. It’s recommended NOT making changes to this value.

ew_eta, ew_cutfloat

Parameters eta and cut to converge Ewald summation. See get_ewald_params()

pseudodict or str

To define pseudopotential.

ecpdict or str

To define ECP type pseudopotential.

h(3,3) ndarray

a.T. Deprecated

dimensionint

Default is 3

low_dim_ft_typestr

For semi-empirical periodic systems, whether to calculate integrals at the non-PBC dimension using the sampled mesh grids in infinity vacuum (inf_vacuum) or truncated Coulomb potential (analytic_2d_1). Unless explicitly specified, analytic_2d_1 is used for 2D system and inf_vacuum is assumed for 1D and 0D.

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_ecp_env(_atm, _ecp, pre_env=[])[source]

Generate the input arguments _ecpbas for ECP integrals

make_kpts(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_aroundbool

To ensure all kpts are in first Brillouin zone.

with_gamma_pointbool

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()

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

pbc_eval_ao(eval_name, coords, comp=None, kpts=None, kpt=None, shls_slice=None, non0tab=None, ao_loc=None, out=None)

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>

atmint32 ndarray

libcint integral function argument

basint32 ndarray

libcint integral function argument

envfloat64 ndarray

libcint integral function argument

coords2D array, shape (N,3)

The coordinates of the grids.

Kwargs:
shls_slice2-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.

non0tab2D 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()

outndarray

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_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>

atmint32 ndarray

libcint integral function argument

basint32 ndarray

libcint integral function argument

envfloat64 ndarray

libcint integral function argument

coords2D array, shape (N,3)

The coordinates of the grids.

Kwargs:
shls_slice2-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.

non0tab2D 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()

outndarray

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

See also Mole.intor

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(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.

10.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.

10.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.

10.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       
                                ''')},
    pseudo = 'gth-pade',
    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.

10.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.

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