10.30. symm – Point group symmetry and spin symmetry

This module offers the functions to detect point group symmetry, basis symmetriziation, Clebsch-Gordon coefficients. This module works as a plugin of PySCF package. Symmetry is not hard coded in each method.

PySCF supports D2h symmetry and linear molecule symmetry (Dooh and Coov). For D2h, the direct production of representations are

D2h

A1g

B1g

B2g

B3g

A1u

B1u

B2u

B3u

A1g

A1g

B1g

B1g

A1g

B2g

B2g

B3g

A1g

B3g

B3g

B2g

B1g

A1g

A1u

A1u

B1u

B2u

B3u

A1g

B1u

B1u

A1u

B3u

B2u

B1g

A1g

B2u

B2u

B3u

A1u

B1u

B2g

B3g

A1g

B3u

B3u

B2u

B1u

A1u

B3g

B2g

B1g

A1g

The multiplication table for XOR operator is

XOR

000

001

010

011

100

101

110

111

000

000

001

001

000

010

010

011

000

011

011

010

001

000

100

100

101

110

111

000

101

101

100

111

110

001

000

110

110

111

100

101

010

011

000

111

111

110

101

100

011

010

001

000

Comparing the two table, we notice that the two tables can be changed to each other with the mapping

D2h

XOR

ID

A1g

000

0

B1g

001

1

B2g

010

2

B3g

011

3

A1u

100

4

B1u

101

5

B2u

110

6

B3u

111

7

The XOR operator and the D2h subgroups have the similar relationships. We therefore use the XOR operator ID to assign the irreps (see pyscf/symm/param.py).

C2h

XOR

ID

C2v

XOR

ID

D2

XOR

ID

Ag

00

0

A1

00

0

A1

00

0

Bg

01

1

A2

01

1

B1

01

1

Au

10

2

B1

10

2

B2

10

2

Bu

11

3

B2

11

3

B3

11

3

Cs

XOR

ID

Ci

XOR

ID

C2

XOR

ID

A’

0

0

Ag

0

0

A

0

0

B”

1

1

Au

1

1

B

1

1

To easily get the relationship between the linear molecule symmetry and D2h/C2v, the ID for irreducible representations of linear molecule symmetry are chosen as (see pyscf/symm/basis.py)

\(D_{\infty h}\)

ID

\(D_{2h}\)

ID

A1g

0

Ag

0

A2g

1

B1g

1

A1u

5

B1u

5

A2u

4

Au

4

E1gx

2

B2g

2

E1gy

3

B3g

3

E1uy

6

B2u

6

E1ux

7

B3u

7

E2gx

10

Ag

0

E2gy

11

B1g

1

E2ux

15

B1u

5

E2uy

14

Au

4

E3gx

12

B2g

2

E3gy

13

B3g

3

E3uy

16

B2u

6

E3ux

17

B3u

7

E4gx

20

Ag

0

E4gy

21

B1g

1

E4ux

25

B1u

5

E4uy

24

Au

4

E5gx

22

B2g

2

E5gy

23

B3g

3

E5uy

26

B2u

6

E5ux

27

B3u

7

and

\(C_{\infty v}\)

ID

\(C_{2v}\)

ID

A1

0

A1

0

A2

1

A2

1

E1x

2

B1

2

E1y

3

B2

3

E2x

10

A1

0

E2y

11

A2

1

E3x

12

B1

2

E3y

13

B2

3

E4x

20

A1

0

E4y

21

A2

1

E5x

22

B1

2

E5y

23

B2

3

So that, the subduction from linear molecule symmetry to D2h/C2v can be achieved by the modular operation %10.

In many output messages, the irreducible representations are labeld with the IDs instead of the irreps’ symbols. We can use symm.addons.irrep_id2name() to convert the ID to irreps’ symbol, e.g.:

>>> from pyscf import symm
>>> [symm.irrep_id2name('Dooh', x) for x in [7, 6, 0, 10, 11, 0, 5, 3, 2, 5, 15, 14]]
['E1ux', 'E1uy', 'A1g', 'E2gx', 'E2gy', 'A1g', 'A1u', 'E1gy', 'E1gx', 'A1u', 'E2ux', 'E2uy']

10.30.1. Enabling symmetry in other module

  • SCF

    To control the HF determinant symmetry, one can assign occupancy for particular irreps, e.g.

#!/usr/bin/env python
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

from pyscf import gto
from pyscf import scf

'''
Specify irrep_nelec to control the wave function symmetry
'''


mol = gto.Mole()
mol.build(
    verbose = 0,
    atom = '''
        C     0.   0.   0.625
        C     0.   0.  -0.625 ''',
    basis = 'cc-pVDZ',
    spin = 0,
    symmetry = True,
)

mf = scf.RHF(mol)

# Frozen occupancy
# 'A1g': 4 electrons
# 'E1gx': 2 electrons
# 'E1gy': 2 electrons
# Rest 4 electrons are put in irreps A1u, E1ux, E1uy ... based on Aufbau principle
# The irrep names can be found in pyscf/symm/param.py
mf.irrep_nelec = {'A1g': 4, 'E1gx': 2, 'E1gy': 2}
e = mf.kernel()
print('E = %.15g  ref = -74.1112374269129' % e)


mol.symmetry = 'D2h'
mol.charge = 1
mol.spin = 1
mol.build(dump_input=False, parse_arg=False)
mf = scf.RHF(mol)

# Frozen occupancy
# 'Ag': 2 alpha, 1 beta electrons
# 'B1u': 4 electrons
# 'B2u': 2 electrons
# 'B3u': 2 electrons
mf.irrep_nelec = {'Ag': (2,1), 'B1u': 4, 'B2u': 2, 'B3u': 2,}
e = mf.kernel()
print('E = %.15g  ref = -74.4026583773135' % e)

# Frozen occupancy
# 'Ag': 4 electrons
# 'B1u': 2 alpha, 1 beta electrons
# 'B2u': 2 electrons
# 'B3u': 2 electrons
mf.irrep_nelec = {'Ag': 4, 'B1u': (2,1), 'B2u': 2, 'B3u': 2,}
e = mf.kernel()
print('E = %.15g  ref = -74.8971476600812' % e)
  • FCI

    FCI wavefunction symmetry can be controlled by initial guess. Function fci.addons.symm_initguess() can generate the FCI initial guess with the right symmetry.

  • MCSCF

    The symmetry of active space in the CASCI/CASSCF calculations can controlled

#!/usr/bin/env python
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

'''
Active space can be adjusted by specifing the number of orbitals for each irrep.
'''

from pyscf import gto, scf, mcscf

mol = gto.Mole()
mol.build(
    atom = 'N  0  0  0; N  0  0  2',
    basis = 'ccpvtz',
    symmetry = True,
)
myhf = scf.RHF(mol)
myhf.kernel()

mymc = mcscf.CASSCF(myhf, 8, 4)
# Select active orbitals which have the specified symmetry
# 2 E1gx orbitals, 2 E1gy orbitals, 2 E1ux orbitals, 2 E1uy orbitals
cas_space_symmetry = {'E1gx':2, 'E1gy':2, 'E1ux':2, 'E1uy':2}
mo = mcscf.sort_mo_by_irrep(mymc, myhf.mo_coeff, cas_space_symmetry)
mymc.kernel(mo)


mol = gto.Mole()
mol.build(
    atom = 'Ti',
    basis = 'ccpvdz',
    symmetry = True,
    spin = 2,
)
myhf = scf.RHF(mol)
myhf.kernel()
myhf.analyze()

mymc = mcscf.CASSCF(myhf, 14, 2)
# Put 3d, 4d, 4s, 4p orbtials in active space
cas_space_symmetry = {'A1g': 3,  # 4s, 3d(z^2), 4d(z^2)
                      'E2gx':2, 'E2gy':2, 'E1gx':2, 'E1gy':2,  # 3d and 4d
                      'A1u':1, 'E1ux':1, 'E1uy':1,  # 4p
                     }
mo = mcscf.sort_mo_by_irrep(mymc, myhf.mo_coeff, cas_space_symmetry)
mymc.verbose = 4
mymc.fcisolver.wfnsym = 'E3gx'
mymc.kernel(mo)
mymc.analyze()
  • MP2 and CCSD

    Point group symmetry are not supported in CCSD, MP2.

10.30.1.1. Examples

Relevant examples examples/symm/30-guess_symmetry.py examples/symm/31-symmetrize_space.py examples/symm/32-symmetrize_natural_orbital.py examples/symm/33-lz_adaption.py

10.30.1.2. Program reference

10.30.2. geom

exception pyscf.symm.geom.RotationAxisNotFound[source]
pyscf.symm.geom.alias_axes(axes, ref)[source]

Rename axes, make it as close as possible to the ref axes

pyscf.symm.geom.detect_symm(atoms, basis=None, verbose=2)[source]

Detect the point group symmetry for given molecule.

Return group name, charge center, and nex_axis (three rows for x,y,z)

pyscf.symm.geom.rotation_mat(vec, theta)[source]

rotate angle theta along vec new(x,y,z) = R * old(x,y,z)

pyscf.symm.geom.symm_identical_atoms(gpname, atoms)[source]

Requires

10.30.3. Symmtry adapted basis

Generate symmetry adapted basis

pyscf.symm.basis.linearmole_symm_descent(gpname, irrepid)[source]

Map irreps to D2h or C2v

10.30.4. addons

pyscf.symm.addons.eigh(h, orbsym)[source]

Solve eigenvalue problem based on the symmetry information for basis. See also pyscf/lib/linalg_helper.py eigh_by_blocks()

Examples:

>>> from pyscf import gto, symm
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1', basis='ccpvdz', symmetry=True)
>>> c = numpy.hstack(mol.symm_orb)
>>> vnuc_so = reduce(numpy.dot, (c.T, mol.intor('int1e_nuc_sph'), c))
>>> orbsym = symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, c)
>>> symm.eigh(vnuc_so, orbsym)
(array([-4.50766885, -1.80666351, -1.7808565 , -1.7808565 , -1.74189134,
        -0.98998583, -0.98998583, -0.40322226, -0.30242374, -0.07608981]),
 ...)
pyscf.symm.addons.irrep_id2name(gpname, irrep_id)[source]

Convert the internal irrep ID to irrep symbol

Args:
gpnamestr

The point group symbol

irrep_idint

See IRREP_ID_TABLE in pyscf/symm/param.py

Returns:

Irrep sybmol, str

pyscf.symm.addons.irrep_name2id(gpname, symb)[source]

Convert the irrep symbol to internal irrep ID

Args:
gpnamestr

The point group symbol

symbstr

Irrep symbol

Returns:

Irrep ID, int

pyscf.symm.addons.label_orb_symm(mol, irrep_name, symm_orb, mo, s=None, check=True, tol=1e-09)[source]

Label the symmetry of given orbitals

irrep_name can be either the symbol or the ID of the irreducible representation. If the ID is provided, it returns the numeric code associated with XOR operator, see symm.param.IRREP_ID_TABLE()

Args:

mol : an instance of Mole

irrep_namelist of str or int

A list of irrep ID or name, it can be either mol.irrep_id or mol.irrep_name. It can affect the return “label”.

symm_orblist of 2d array

the symmetry adapted basis

mo2d array

the orbitals to label

Returns:

list of symbols or integers to represent the irreps for the given orbitals

Examples:

>>> from pyscf import gto, scf, symm
>>> mol = gto.M(atom='H 0 0 0; H 0 0 1', basis='ccpvdz',verbose=0, symmetry=1)
>>> mf = scf.RHF(mol)
>>> mf.kernel()
>>> symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mf.mo_coeff)
['Ag', 'B1u', 'Ag', 'B1u', 'B2u', 'B3u', 'Ag', 'B2g', 'B3g', 'B1u']
>>> symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mf.mo_coeff)
[0, 5, 0, 5, 6, 7, 0, 2, 3, 5]
pyscf.symm.addons.route(target, nelec, orbsym)[source]

Pick orbitals to form a determinant which has the right symmetry. If solution is not found, return []

pyscf.symm.addons.std_symb(gpname)[source]

std_symb(‘d2h’) returns D2h; std_symb(‘D2H’) returns D2h

pyscf.symm.addons.symmetrize_orb(mol, mo, orbsym=None, s=None, check=False)[source]

Symmetrize the given orbitals.

This function is different to the symmetrize_space(): In this function, each orbital is symmetrized by removing non-symmetric components. symmetrize_space() symmetrizes the entire space by mixing different orbitals.

Note this function might return non-orthorgonal orbitals. Call symmetrize_space() to find the symmetrized orbitals that are close to the given orbitals.

Args:
mo2D float array

The orbital space to symmetrize

Kwargs:
orbsyminteger list

Irrep id for each orbital. If not given, the irreps are guessed by calling label_orb_symm().

s2D float array

Overlap matrix. If given, use this overlap than the the overlap of the input mol.

Returns:

2D orbital coefficients

Examples:

>>> from pyscf import gto, symm, scf
>>> mol = gto.M(atom = 'C  0  0  0; H  1  1  1; H -1 -1  1; H  1 -1 -1; H -1  1 -1',
...             basis = 'sto3g')
>>> mf = scf.RHF(mol).run()
>>> mol.build(0, 0, symmetry='D2')
>>> mo = symm.symmetrize_orb(mol, mf.mo_coeff)
>>> print(symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mo))
['A', 'A', 'B1', 'B2', 'B3', 'A', 'B1', 'B2', 'B3']
pyscf.symm.addons.symmetrize_space(mol, mo, s=None, check=True, tol=1e-07)[source]

Symmetrize the given orbital space.

This function is different to the symmetrize_orb(): In this function, the given orbitals are mixed to reveal the symmtery; symmetrize_orb() projects out non-symmetric components for each orbital.

Args:
mo2D float array

The orbital space to symmetrize

Kwargs:
s2D float array

Overlap matrix. If not given, overlap is computed with the input mol.

Returns:

2D orbital coefficients

Examples:

>>> from pyscf import gto, symm, scf
>>> mol = gto.M(atom = 'C  0  0  0; H  1  1  1; H -1 -1  1; H  1 -1 -1; H -1  1 -1',
...             basis = 'sto3g')
>>> mf = scf.RHF(mol).run()
>>> mol.build(0, 0, symmetry='D2')
>>> mo = symm.symmetrize_space(mol, mf.mo_coeff)
>>> print(symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mo))
['A', 'A', 'A', 'B1', 'B1', 'B2', 'B2', 'B3', 'B3']

10.30.5. Clebsch Gordon coefficients

pyscf.symm.cg.cg_spin(l, jdouble, mjdouble, spin)[source]

Clebsch Gordon coefficient of <l,m,1/2,spin|j,mj>

10.30.6. Parameters

10.30.7. Spherical harmonics

Spherical harmonics

pyscf.symm.sph.cart2spinor(l)[source]

Cartesian to spinor for angular moment l

pyscf.symm.sph.multipoles(r, lmax, reorder_dipole=True)[source]

Compute all multipoles upto lmax

rad = numpy.linalg.norm(r, axis=1) ylms = real_ylm(r/rad.reshape(-1,1), lmax) pol = [rad**l*y for l, y in enumerate(ylms)]

Kwargs:
reorder_pbool

sort dipole to the order (x,y,z)

pyscf.symm.sph.real2spinor_whole(mol)

Transformation matrix that transforms real-spherical GTOs to spinor GTOs for all basis functions

Examples:

>>> from pyscf import gto
>>> from pyscf.symm import sph
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz')
>>> ca, cb = sph.sph2spinor_coeff(mol)
>>> s0 = mol.intor('int1e_ovlp_spinor')
>>> s1 = ca.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(ca)
>>> s1+= cb.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(cb)
>>> print(abs(s1-s0).max())
>>> 6.66133814775e-16
pyscf.symm.sph.real_sph_vec(r, lmax, reorder_p=False)[source]

Computes (all) real spherical harmonics up to the angular momentum lmax

pyscf.symm.sph.sph2spinor_coeff(mol)[source]

Transformation matrix that transforms real-spherical GTOs to spinor GTOs for all basis functions

Examples:

>>> from pyscf import gto
>>> from pyscf.symm import sph
>>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz')
>>> ca, cb = sph.sph2spinor_coeff(mol)
>>> s0 = mol.intor('int1e_ovlp_spinor')
>>> s1 = ca.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(ca)
>>> s1+= cb.conj().T.dot(mol.intor('int1e_ovlp_sph')).dot(cb)
>>> print(abs(s1-s0).max())
>>> 6.66133814775e-16
pyscf.symm.sph.sph_pure2real(l, reorder_p=True)[source]

Transformation matrix: from the pure spherical harmonic functions Y_m to the real spherical harmonic functions O_m.

O_m = sum Y_m’ * U(m’,m)

Y(-1) = 1/sqrt(2){-iO(-1) + O(1)}; Y(1) = 1/sqrt(2){-iO(-1) - O(1)} Y(-2) = 1/sqrt(2){-iO(-2) + O(2)}; Y(2) = 1/sqrt(2){iO(-2) + O(2)} O(-1) = i/sqrt(2){Y(-1) + Y(1)}; O(1) = 1/sqrt(2){Y(-1) - Y(1)} O(-2) = i/sqrt(2){Y(-2) - Y(2)}; O(2) = 1/sqrt(2){Y(-2) + Y(2)}

Kwargs:

reorder_p (bool): Whether the p functions are in the (x,y,z) order.

Returns:

2D array U_{complex,real}

pyscf.symm.sph.sph_real2pure(l, reorder_p=True)[source]

Transformation matrix: from real spherical harmonic functions to the pure spherical harmonic functions.

Kwargs:

reorder_p (bool): Whether the real p functions are in the (x,y,z) order.

10.30.8. Wigner rotation Dmatrix

Wigner rotation D-matrix for real spherical harmonics

pyscf.symm.Dmatrix.Dmatrix(l, alpha, beta, gamma, reorder_p=False)[source]

Wigner rotation D-matrix

D_{mm’} = <lm|R(alpha,beta,gamma)|lm’> alpha, beta, gamma are Euler angles (in z-y-z convention)

Kwargs:

reorder_p (bool): Whether to put the p functions in the (x,y,z) order.

pyscf.symm.Dmatrix.dmatrix(l, beta, reorder_p=False)[source]

Wigner small-d matrix (in z-y-z convention)

pyscf.symm.Dmatrix.get_euler_angles(c1, c2)[source]

Find the three Euler angles (alpha, beta, gamma in z-y-z convention) that rotates coordinates c1 to coordinates c2.

yp = numpy.einsum(‘j,kj->k’, c1[1], geom.rotation_mat(c1[2], beta)) tmp = numpy.einsum(‘ij,kj->ik’, c1 , geom.rotation_mat(c1[2], alpha)) tmp = numpy.einsum(‘ij,kj->ik’, tmp, geom.rotation_mat(yp , beta )) c2 = numpy.einsum(‘ij,kj->ik’, tmp, geom.rotation_mat(c2[2], gamma))

(For backward compatibility) if c1 and c2 are two points in the real space, the Euler angles define the rotation transforms the old coordinates to the new coordinates (new_x, new_y, new_z) in which c1 is identical to c2.

tmp = numpy.einsum(‘j,kj->k’, c1 , geom.rotation_mat((0,0,1), gamma)) tmp = numpy.einsum(‘j,kj->k’, tmp, geom.rotation_mat((0,1,0), beta) ) c2 = numpy.einsum(‘j,kj->k’, tmp, geom.rotation_mat((0,0,1), alpha))