Source code for pyscf.symm.addons

#!/usr/bin/env python
# Copyright 2014-2020 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Qiming Sun <osirpt.sun@gmail.com>
#

from functools import reduce
import numpy
import scipy.linalg
from pyscf import lib
from pyscf.lib import logger
from pyscf.lib.exceptions import PointGroupSymmetryError
from pyscf.symm import basis
from pyscf.symm import param
from pyscf import __config__


MULTI_IRREPS = -1

[docs] def label_orb_symm(mol, irrep_name, symm_orb, mo, s=None, check=getattr(__config__, 'symm_addons_label_orb_symm_check', True), tol=getattr(__config__, 'symm_addons_label_orb_symm_tol', 1e-9)): '''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 :py:meth:`symm.param.IRREP_ID_TABLE` Args: mol : an instance of :class:`Mole` irrep_name : list 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_orb : list of 2d array the symmetry adapted basis mo : 2d 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] ''' nmo = mo.shape[1] if s is None: s = mol.intor_symmetric('int1e_ovlp') s_mo = numpy.dot(s, mo) norm = numpy.zeros((len(irrep_name), nmo)) for i, csym in enumerate(symm_orb): moso = numpy.dot(csym.T.conj(), s_mo) ovlpso = reduce(numpy.dot, (csym.T.conj(), s, csym)) try: s_moso = lib.cho_solve(ovlpso, moso) except numpy.linalg.LinAlgError: ovlpso[numpy.diag_indices(csym.shape[1])] += 1e-12 s_moso = lib.cho_solve(ovlpso, moso) norm[i] = numpy.einsum('ki,ki->i', moso.conj(), s_moso).real norm /= numpy.sum(norm, axis=0) # for orbitals which are not normalized iridx = numpy.argmax(norm, axis=0) orbsym = numpy.asarray([irrep_name[i] for i in iridx]) logger.debug(mol, 'irreps of each MO %s', orbsym) if check: largest_norm = norm[iridx,numpy.arange(nmo)] orbidx = numpy.where(largest_norm < 1-tol)[0] if orbidx.size > 0: idx = numpy.where(largest_norm < 1-tol*1e2)[0] if idx.size > 0: raise ValueError('orbitals %s not symmetrized, norm = %s' % (idx, largest_norm[idx])) else: logger.warn(mol, 'orbitals %s not strictly symmetrized.', numpy.unique(orbidx)) logger.warn(mol, 'They can be symmetrized with ' 'pyscf.symm.symmetrize_space function.') logger.debug(mol, 'norm = %s', largest_norm[orbidx]) return orbsym
[docs] def symmetrize_orb(mol, mo, orbsym=None, s=None, check=getattr(__config__, 'symm_addons_symmetrize_orb_check', False)): '''Symmetrize the given orbitals. This function is different to the :func:`symmetrize_space`: In this function, each orbital is symmetrized by removing non-symmetric components. :func:`symmetrize_space` symmetrizes the entire space by mixing different orbitals. Note this function might return non-orthorgonal orbitals. Call :func:`symmetrize_space` to find the symmetrized orbitals that are close to the given orbitals. Args: mo : 2D float array The orbital space to symmetrize Kwargs: orbsym : integer list Irrep id for each orbital. If not given, the irreps are guessed by calling :func:`label_orb_symm`. s : 2D 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'] ''' if s is None: s = mol.intor_symmetric('int1e_ovlp') if orbsym is None: orbsym = label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo, s=s, check=check) orbsym = numpy.asarray(orbsym) s_mo = numpy.dot(s, mo) mo1 = numpy.empty_like(mo) if orbsym[0] in mol.irrep_name: irrep_id = mol.irrep_name else: irrep_id = mol.irrep_id for i, ir in enumerate(irrep_id): idx = orbsym == ir csym = mol.symm_orb[i] ovlpso = reduce(numpy.dot, (csym.T.conj(), s, csym)) sc = lib.cho_solve(ovlpso, numpy.dot(csym.T.conj(), s_mo[:,idx])) mo1[:,idx] = numpy.dot(csym, sc) return mo1
[docs] def find_symmetric_mo(moso, ovlpso, thr=1e-8): '''Find the list of MO of that thransform like particular irrep Args: moso : 2D float array Overlap matrix of symmetry-adapted AO and MO, it can be obtained by reduce(numpy.dot, (csym.T.conj(), s, mo)), where csym taken from mol.symm_orb, and s is the AO overlap matrix ovlpso : 2D float array Overlap matrix between symmetry-adapted AO, it can be obtained by reduce(numpy.dot, (csym.T.conj(), s, csym)) Kwargs: thr : float Threshold to consider MO symmetry-adapted Returns: 1D bool array to select symmetry adapted MO ''' irrep_dim = ovlpso.shape[0] try: diag = numpy.einsum('ki,ki->i', moso.conj(), lib.cho_solve(ovlpso, moso)) except numpy.linalg.LinAlgError: ovlpso[numpy.diag_indices(irrep_dim)] += 1e-12 diag = numpy.einsum('ki,ki->i', moso.conj(), lib.cho_solve(ovlpso, moso)) idx = abs(1-diag) < thr return idx
[docs] def symmetrize_space(mol, mo, s=None, check=getattr(__config__, 'symm_addons_symmetrize_space_check', True), tol=getattr(__config__, 'symm_addons_symmetrize_space_tol', 1e-9), clean=getattr(__config__, 'symm_addons_symmetrize_space_clean', False)): '''Symmetrize the given orbital space. This function is different to the :func:`symmetrize_orb`: In this function, the given orbitals are mixed to reveal the symmtery; :func:`symmetrize_orb` projects out non-symmetric components for each orbital. Args: mol : an instance of :class:`Mole` mo : 2D float array The orbital space to symmetrize Kwargs: s : 2D float array Overlap matrix. If not given, overlap is computed with the input mol. check : bool Whether to check orthogonality of input orbitals and try to fix it tol : float Orthogonality tolerance clean : bool Whether to zero out symmetry forbidden orbital coefficients 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'] ''' from pyscf.tools import mo_mapping from pyscf.lo import orth if s is None: s = mol.intor_symmetric('int1e_ovlp') nmo = mo.shape[1] s_mo = numpy.dot(s, mo) if check: moso = reduce(numpy.dot, (mo.conj().T, s, mo)) max_non_orth = abs(moso - numpy.eye(nmo)).max() logger.debug(mol, 'Non-orthogonality of input orbitals before symmetrization %8.2e', max_non_orth) if max_non_orth > tol: logger.info(mol, 'Input orbitals are not orthogonal, perform Lowdin orthogonalization') mo = numpy.dot(mo, orth.lowdin(moso)) s_mo = numpy.dot(s, mo) max_non_orth = abs(reduce(numpy.dot, (mo.conj().T, s, mo)) - numpy.eye(nmo)).max() logger.debug(mol, 'Non-orthogonality of orbitals before symmetrization %8.2e', max_non_orth) if max_non_orth > tol: raise ValueError('Input orbitals are not orthogonalized') mo1 = [] for i, csym in enumerate(mol.symm_orb): moso = numpy.dot(csym.T.conj(), s_mo) ovlpso = reduce(numpy.dot, (csym.T.conj(), s, csym)) # excluding orbitals which are already symmetrized idx = find_symmetric_mo(moso, ovlpso) orb_irrep = mo[:, idx] if clean: # rotate MO to symm-adapted AO and back to zero-out non-symmetric part SALC_mo = numpy.dot(csym.T.conj(), orb_irrep) orb_irrep = numpy.dot(csym, SALC_mo) moso1 = moso[:, ~idx] dm = numpy.dot(moso1, moso1.T.conj()) if dm.trace() > 1e-8 and sum(idx) < csym.shape[1]: logger.debug(mol, '%5d orbital(s) symmetrized in irrep %3d', csym.shape[1] - sum(idx), i) e, u = scipy.linalg.eigh(dm, ovlpso) orb_symmetrized = numpy.dot(csym, u[:, abs(1-e) < 1e-6]) orb_irrep = numpy.hstack([orb_irrep, orb_symmetrized]) moso = numpy.dot(orb_irrep.T.conj(), numpy.dot(s, orb_irrep)) if moso.shape[0] == 0: continue max_non_orth = abs(moso - numpy.eye(moso.shape[0])).max() logger.debug(mol, 'Non-orthogonality in irrep %3d after symmetrization: %8.2e', i, max_non_orth) if max_non_orth > tol: logger.info(mol, 'Symmetrized orbitals in irrep %3d not orthogonal, perform Lowdin orthogonalization', i) orb_irrep = numpy.dot(orb_irrep, orth.lowdin(moso)) max_non_orth = abs(numpy.dot(orb_irrep.T.conj(), numpy.dot(s, orb_irrep)) - numpy.eye(orb_irrep.shape[1])).max() logger.debug(mol, 'Non-orthogonality in irrep %3d after symmetrization and orthogonalizastion: %8.2e', i, max_non_orth) mo1.append(orb_irrep) mo1 = numpy.hstack(mo1) if mo1.shape[1] != nmo: raise ValueError('mo1.shape[1] != nmo: %d != %d The input orbital space is not symmetrized.\n One ' 'possible reason is that the input mol and orbitals ' 'are of different orientation.' % (mo1.shape[1], nmo)) if check: moso1 = reduce(numpy.dot, (mo1.conj().T, s, mo1)) max_non_orth = abs(moso1 - numpy.eye(nmo)).max() logger.debug(mol, 'Non-orthogonality of output orbitals after symmetrization %8.2e', max_non_orth) if max_non_orth > tol: logger.info(mol, 'Symmetrized output orbitals are not orthogonalized, perform Lowdin orthogonalization') mo1 = numpy.dot(mo1, orth.lowdin(moso1)) max_non_orth = abs(reduce(numpy.dot, (mo1.conj().T, s, mo1)) - numpy.eye(nmo)).max() logger.debug(mol, 'Non-orthogonality of output orbitals after orthogonalization %8.2e', max_non_orth) if max_non_orth > tol: raise ValueError('Output orbitals are not orthogonalized') idx = mo_mapping.mo_1to1map(reduce(numpy.dot, (mo.T.conj(), s, mo1))) return mo1[:, idx]
[docs] def symmetrize_multidim(mol, mo, s=None, check=getattr(__config__, 'symm_addons_symmetrize_space_check', True), tol=getattr(__config__, 'symm_addons_symmetrize_space_tol', 1e-10), keep_phase=getattr(__config__, 'symm_addons_symmetrize_multidim_keep_phase', True)): '''Symmetrize orbitals with respect to multidimensional irreps. Make coefficients of partner functions of multidimensional irreps to be the same. The functions uses the convention of the libmsym interface, that introduces underscores to the labels of multidimensional irreps partners. Args: mol : an instance of :class:`Mole` Symmetry-adapted basis with multidimensional irreps should be generated by libmsym mo : 2D float array The orbital space to symmetrize Kwargs: s : 2D float array Overlap matrix. If not given, overlap is computed with the input mol. check : bool Whether to check orthogonality of input orbitals and try to fix it tol : float Orthogonality tolerance keep_phase : bool Whether to keep original orbital phases, rather then make them coherent with the first partner Returns: 2D orbital coefficients ''' from pyscf.tools import mo_mapping from pyscf.lo import orth if s is None: s = mol.intor_symmetric('int1e_ovlp') nmo = mo.shape[1] s_mo = numpy.dot(s, mo) mo1 = [] if check: moso = reduce(numpy.dot, (mo.conj().T, s, mo)) max_non_orth = abs(moso - numpy.eye(nmo)).max() logger.info(mol, 'Non-orthogonality of input orbitals %8.2e', max_non_orth) if max_non_orth > tol: mo_lowdin = numpy.dot(mo, orth.lowdin(moso)) max_non_orth_lowdin = abs(reduce(numpy.dot, (mo_lowdin.conj().T, s, mo_lowdin)) - numpy.eye(nmo)).max() logger.info(mol, 'Non-orthogonality after Lowdin orthogonalization %8.2e', max_non_orth) if (max_non_orth_lowdin - max_non_orth) > tol/100: mo = mo_lowdin s_mo = numpy.dot(s, mo) logger.info(mol, 'Use Lowdin-orthogonalizied input orbitals') else: logger.info(mol, 'Use original input orbitals') irreps_mdim = [] irreps_1dim = [] for irrep in mol.irrep_name: if "_" in irrep: added = False base = irrep.split("_")[0] for partners in irreps_mdim: if any(base in elem for elem in partners): partners.append(irrep) added = True if not added: # new base irreps_mdim.append([irrep]) else: irreps_1dim.append(irrep) for irrep in irreps_1dim: i = mol.irrep_name.index(irrep) csym = mol.symm_orb[i] moso = numpy.dot(csym.T.conj(), s_mo) ovlpso = reduce(numpy.dot, (csym.T.conj(), s, csym)) # find MO that thransform like irrep idx = find_symmetric_mo(moso, ovlpso) if sum(idx) != csym.shape[1]: raise ValueError('Number of symmetry-adapted MOs in not equal to dimensionality ' 'of irrep %s: %d != %d' % (mol.irrep_name[irrep], csym.shape[1], sum(idx))) mo1.append(mo[:, idx]) for partners in irreps_mdim: SALC_ao_partners = [] SALC_ov_partners = [] SALC_mo_partners = [] for irrep in partners: csym = mol.symm_orb[mol.irrep_name.index(irrep)] SALC_ao_partners.append(csym) moso = numpy.dot(csym.T.conj(), s_mo) ovlpso = reduce(numpy.dot, (csym.T.conj(), s, csym)) SALC_ov_partners.append(ovlpso) # find MO that thransform like irrep idx = find_symmetric_mo(moso, ovlpso) mo_irrep = mo[:, idx] if sum(idx) != csym.shape[1]: raise ValueError('Number of symmetry-adapted MOs in not equal to dimensionality ' 'of irrep %s: %d != %d' % (mol.irrep_name[irrep], csym.shape[1], sum(idx))) # rotate MO to symm-adapted AO SALC_mo = numpy.dot(csym.T.conj(), mo_irrep) SALC_mo_partners.append(SALC_mo) irrep_dim = SALC_mo_partners[0].shape[1] phases = [numpy.ones(irrep_dim)] idxes = [numpy.array(range(irrep_dim))] SALC_mo_average = SALC_mo_partners[0] for i in range(1, len(partners)): if not numpy.allclose(SALC_ov_partners[0], SALC_ov_partners[i]): raise ValueError('Symmetry-adapted AO of partner functions %s and %s are not compatible' % (partners[0], partners[i])) ov_0i = reduce(numpy.dot, (SALC_mo_partners[0].T.conj(), SALC_ov_partners[0], SALC_mo_partners[i])) idx = mo_mapping.mo_1to1map(ov_0i) idxes.append(idx) phase = numpy.rint(numpy.diagonal(ov_0i[:, idx])) phases.append(phase) partners_mo_i = phase*SALC_mo_partners[i][:, idx] SALC_mo_average += partners_mo_i SALC_mo_average /= len(partners) for csym, phase, SALC_mo in zip(SALC_ao_partners, phases, SALC_mo_partners): if keep_phase: SALC_mo = phase*SALC_mo_average else: SALC_mo = SALC_mo_average orb_done = numpy.dot(csym, SALC_mo) moso1 = reduce(numpy.dot, (orb_done.T.conj(), s, orb_done)) max_non_orth = abs(moso1 - numpy.eye(moso1.shape[0])).max() logger.debug(mol, 'Non-orthogonality after symmetrization of %s: %s', partners[i], max_non_orth) if max_non_orth > tol: orb_done = numpy.dot(orb_done, orth.lowdin(moso1)) max_non_orth = abs(numpy.dot(orb_done.T.conj(), numpy.dot(s, orb_done)) - numpy.eye(orb_done.shape[1])).max() logger.debug(mol, 'After additional orthogonalization: %s', max_non_orth) mo1.append(orb_done) mo1 = numpy.hstack(mo1) if check: moso1 = reduce(numpy.dot, (mo1.conj().T, s, mo1)) max_non_orth = abs(moso1 - numpy.eye(nmo)).max() logger.debug(mol, 'Non-orthogonality of output orbitals after symmetrization %8.2e', max_non_orth) if max_non_orth > tol: logger.info(mol, 'Output orbitals are not orthogonalized, perform Lowdin orthogonalization') mo1 = numpy.dot(mo1, orth.lowdin(moso1)) max_non_orth = abs(reduce(numpy.dot, (mo1.conj().T, s, mo1)) - numpy.eye(nmo)).max() logger.info(mol, 'Non-orthogonality of output orbitals after orthogonalization %8.2e', max_non_orth) if max_non_orth > tol: raise ValueError('Output orbitals are not orthogonalized') idx = mo_mapping.mo_1to1map(reduce(numpy.dot, (mo.T.conj(), s, mo1))) return mo1[:, idx]
[docs] def std_symb(gpname): '''std_symb('d2h') returns D2h; std_symb('D2H') returns D2h''' if gpname == 'SO3': return gpname else: return str(gpname[0].upper() + gpname[1:].lower())
[docs] def irrep_name2id(gpname, symb): '''Convert the irrep symbol to internal irrep ID Args: gpname : str The point group symbol symb : str Irrep symbol Returns: Irrep ID, int ''' gpname = std_symb(gpname) symb = std_symb(symb) if gpname == 'SO3': return basis.so3_irrep_symb2id(symb) elif gpname in ('Dooh', 'Coov'): return basis.linearmole_irrep_symb2id(gpname, symb) else: return param.IRREP_ID_TABLE[gpname][symb]
[docs] def irrep_id2name(gpname, irrep_id): '''Convert the internal irrep ID to irrep symbol Args: gpname : str The point group symbol irrep_id : int See IRREP_ID_TABLE in pyscf/symm/param.py Returns: Irrep sybmol, str ''' gpname = std_symb(gpname) if gpname == 'SO3': return basis.so3_irrep_id2symb(irrep_id) elif gpname in ('Dooh', 'Coov'): return basis.linearmole_irrep_id2symb(gpname, irrep_id) else: # irrep_id may be obtained from high symmetry (Dooh, Coov) irrep_id_in_d2h = irrep_id % 10 return param.CHARACTER_TABLE[gpname][irrep_id_in_d2h][0]
[docs] def irrep_name(pgname, irrep_id): raise PointGroupSymmetryError('This function was obsoleted. Use irrep_id2name')
[docs] def route(target, nelec, orbsym): '''Pick orbitals to form a determinant which has the right symmetry. If solution is not found, return [] ''' def riter(target, nelec, orbsym): if nelec == 1: if target in orbsym: return [orbsym.index(target)] else: return [] else: for i, ir in enumerate(orbsym): off = i + 1 orb_left = orbsym[off:] res = riter(target ^ ir, nelec-1, orb_left) if res: return [i] + [off+x for x in res] return [] if isinstance(orbsym, numpy.ndarray): orbsym = orbsym.tolist() return riter(target, nelec, orbsym)
[docs] def eigh(h, orbsym): '''Solve eigenvalue problem based on the symmetry information for basis. See also pyscf/lib/linalg_helper.py :func:`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]), ...) ''' return lib.eigh_by_blocks(h, labels=orbsym)
[docs] def direct_prod(orbsym1, orbsym2, groupname='D2h'): if groupname == 'SO3': prod = orbsym1[:,None] ^ orbsym2 orbsym1_not_s = orbsym1 != 0 orbsym2_not_s = orbsym2 != 0 prod[orbsym1_not_s[:,None] & orbsym2_not_s != 0] = MULTI_IRREPS prod[orbsym1[:,None] == orbsym2] = 0 elif groupname == 'Dooh': orbsym1_octa = (orbsym1 // 10) * 8 + orbsym1 % 10 orbsym2_octa = (orbsym2 // 10) * 8 + orbsym2 % 10 prod = orbsym1_octa[:,None] ^ orbsym2_octa prod = (prod % 8) + (prod // 8) * 10 orbsym1_irrepE = (orbsym1 >= 2) & (orbsym1 != 4) & (orbsym1 != 5) orbsym2_irrepE = (orbsym2 >= 2) & (orbsym2 != 4) & (orbsym2 != 5) prod[orbsym1_irrepE[:,None] & orbsym2_irrepE] = MULTI_IRREPS prod[orbsym1[:,None] == orbsym2] = 0 elif groupname == 'Coov': prod = orbsym1[:,None] ^ orbsym2 orbsym1_irrepE = orbsym1 >= 2 orbsym2_irrepE = orbsym2 >= 2 prod[orbsym1_irrepE[:,None] & orbsym2_irrepE] = MULTI_IRREPS prod[orbsym1[:,None] == orbsym2] = 0 else: # D2h and subgroup prod = orbsym1[:,None] ^ orbsym2 return prod
if __name__ == "__main__": from pyscf import gto from pyscf import scf mol = gto.Mole() mol.build( atom = [['H', (0,0,0)], ['H', (0,0,1)]], basis = {'H': 'cc-pvdz'}, symmetry = 1 ) mf = scf.RHF(mol) mf.scf() nao, nmo = mf.mo_coeff.shape print(label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mf.mo_coeff)) numpy.random.seed(1) u = numpy.random.random((nmo,nmo))*1e-2 u = scipy.linalg.expm(u - u.T) mo = symmetrize_orb(mol, numpy.dot(mf.mo_coeff, u)) print(label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mo)) orbsym = [0, 3, 0, 2, 5, 6] res = route(7, 3, orbsym) print(res, reduce(lambda x, y: x ^ y, [orbsym[i] for i in res]))