Source code for pyscf.symm.basis

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

'''
Generate symmetry adapted basis
'''

from functools import reduce
import numpy
from pyscf.data.elements import _symbol, _rm_digit
from pyscf import gto
from pyscf import lib
from pyscf.lib.exceptions import PointGroupSymmetryError
from pyscf.symm import geom
from pyscf.symm import param
from pyscf.symm.Dmatrix import Dmatrix, get_euler_angles

__all__ = ['tot_parity_odd',
           'symm_adapted_basis',
           'dump_symm_adapted_basis',
           'symmetrize_matrix',
           'linearmole_symm_descent',
           'linearmole_irrep_symb2id',
           'linearmole_irrep_id2symb',
           'linearmole_symm_adapted_basis',
           'so3_irrep_symb2id',
           'so3_irrep_id2symb']

OP_PARITY_ODD = {
    'E'  : (0, 0, 0),
    'C2x': (0, 1, 1),
    'C2y': (1, 0, 1),
    'C2z': (1, 1, 0),
    'i'  : (1, 1, 1),
    'sx' : (1, 0, 0),
    'sy' : (0, 1, 0),
    'sz' : (0, 0, 1),
}

[docs] def tot_parity_odd(op, l, m): if op == 'E': return 0 else: ox,oy,oz = OP_PARITY_ODD[op] gx,gy,gz = param.SPHERIC_GTO_PARITY_ODD[l][l+m] return (ox and gx) ^ (oy and gy) ^ (oz and gz)
[docs] def symm_adapted_basis(mol, gpname, orig=0, coordinates=None): if gpname == 'SO3': return so3_symm_adapted_basis(mol, gpname, orig, coordinates) elif gpname in ('Dooh', 'Coov'): return linearmole_symm_adapted_basis(mol, gpname, orig, coordinates) # prop_atoms are the atoms relocated wrt the charge center with proper # orientation if coordinates is None: coordinates = numpy.eye(3) prop_atoms = mol.format_atom(mol._atom, orig, coordinates, 'Bohr') eql_atom_ids = geom.symm_identical_atoms(gpname, prop_atoms) ops = numpy.asarray([param.D2H_OPS[op] for op in param.OPERATOR_TABLE[gpname]]) chartab = numpy.array([x[1:] for x in param.CHARACTER_TABLE[gpname]]) nirrep = chartab.__len__() aoslice = mol.aoslice_by_atom() nao = mol.nao_nr() atom_coords = numpy.array([a[1] for a in prop_atoms]) sodic = [[] for i in range(8)] for atom_ids in eql_atom_ids: r0 = atom_coords[atom_ids[0]] op_coords = numpy.einsum('x,nxy->ny', r0, ops) # Using ops to generate other atoms from atom_ids[0] coords0 = atom_coords[atom_ids] dc = abs(op_coords.reshape(-1,1,3) - coords0).max(axis=2) op_relate_idx = numpy.argwhere(dc < geom.TOLERANCE)[:,1] ao_loc = numpy.array([aoslice[atom_ids[i],2] for i in op_relate_idx]) b0, b1 = aoslice[atom_ids[0],:2] ip = 0 for ib in range(b0, b1): l = mol.bas_angular(ib) if mol.cart: degen = (l + 1) * (l + 2) // 2 cbase = numpy.zeros((degen,nirrep,nao)) for op_id, op in enumerate(ops): n = 0 for x in range(l, -1, -1): for y in range(l-x, -1, -1): z = l-x-y idx = ao_loc[op_id] + n sign = op[0,0]**x * op[1,1]**y * op[2,2]**z cbase[n,:,idx] += sign * chartab[:,op_id] n += 1 else: degen = l * 2 + 1 cbase = numpy.zeros((degen,nirrep,nao)) for op_id, op in enumerate(param.OPERATOR_TABLE[gpname]): for n, m in enumerate(range(-l, l+1)): idx = ao_loc[op_id] + n if tot_parity_odd(op, l, m): cbase[n,:,idx] -= chartab[:,op_id] else: cbase[n,:,idx] += chartab[:,op_id] norms = numpy.sqrt(numpy.einsum('mij,mij->mi', cbase, cbase)) for i in range(mol.bas_nctr(ib)): for n, ir in numpy.argwhere(norms > 1e-12): c = numpy.zeros(nao) c[ip:] = cbase[n,ir,:nao-ip] / norms[n,ir] sodic[ir].append(c) ip += degen l_idx = ao_l_dict(mol) Ds = _momentum_rotation_matrices(mol, coordinates) so = [] irrep_ids = [] for ir, c in enumerate(sodic): if len(c) > 0: irrep_ids.append(ir) c_ir = numpy.vstack(c).T nso = c_ir.shape[1] for l, idx in l_idx.items(): c = c_ir[idx].reshape(-1,Ds[l].shape[1],nso) c_ir[idx] = numpy.einsum('nm,smp->snp', Ds[l], c).reshape(-1,nso) so.append(c_ir) return so, irrep_ids
def ao_l_dict(mol): ao_loc = mol.ao_loc_nr() l_idx = {} ANG_OF = 1 for l in range(mol._bas[:,ANG_OF].max()+1): idx = [numpy.arange(ao_loc[ib], ao_loc[ib+1]) for ib in numpy.where(mol._bas[:,ANG_OF] == l)[0]] if idx: l_idx[l] = numpy.hstack(idx) return l_idx def _momentum_rotation_matrices(mol, axes): '''Cache the rotation matrices for each angular momentum''' alpha, beta, gamma = get_euler_angles(numpy.eye(3), axes) ANG_OF = 1 l_max = mol._bas[:,ANG_OF].max() if not mol.cart: return [Dmatrix(l, alpha, beta, gamma, reorder_p=True) for l in range(l_max+1)] pp = Dmatrix(1, alpha, beta, gamma, reorder_p=True) Ds = [numpy.ones((1,1))] for l in range(1, l_max+1): # All possible x,y,z combinations cidx = numpy.sort(lib.cartesian_prod([(0, 1, 2)] * l), axis=1) addr = 0 affine = numpy.ones((1,1)) for i in range(l): nd = affine.shape[0] * 3 affine = numpy.einsum('ik,jl->ijkl', affine, pp).reshape(nd, nd) addr = addr * 3 + cidx[:,i] uniq_addr, rev_addr = numpy.unique(addr, return_inverse=True) ncart = (l + 1) * (l + 2) // 2 assert ncart == uniq_addr.size trans = numpy.zeros((ncart,ncart)) for i, k in enumerate(rev_addr): trans[k] += affine[i,uniq_addr] Ds.append(trans) return Ds
[docs] def dump_symm_adapted_basis(mol, so): raise PointGroupSymmetryError('TODO')
[docs] def symmetrize_matrix(mat, so): return [reduce(numpy.dot, (c.conj().T,mat,c)) for c in so]
def _basis_offset_for_atoms(atoms, basis_tab): basoff = [0] n = 0 for at in atoms: symb = _symbol(at[0]) if symb in basis_tab: bas0 = basis_tab[symb] else: bas0 = basis_tab[_rm_digit(symb)] for b in bas0: angl = b[0] n += _num_contract(b) * (angl*2+1) basoff.append(n) return n, basoff def _num_contract(basis): if isinstance(basis[1], int): # This branch should never be reached if basis_tab is formatted by # function mole.format_basis nctr = len(basis[2]) - 1 else: nctr = len(basis[1]) - 1 return nctr ############################### # SO3 (real spherical harmonics) # Irreps ID maps # SO3 -> Dooh # s 0 -> A1g 0 # pz 105 -> A1u 5 # py 106 -> E1uy 6 # px 107 -> E1ux 7 # dz2 200 -> A1g 0 # dyz 203 -> E1gy 3 # dxz 202 -> E1gx 2 # dxy 211 -> E2gy 11 # dx2y2 210 -> E2gx 10 # f0 305 -> A1u 5 # f-1 306 -> E1uy 6 # f+1 307 -> E1ux 7 # f-2 314 -> E2uy 14 # f+2 315 -> E2ux 15 # f-3 316 -> E3uy 16 # f+3 317 -> E3ux 17 # g0 400 -> A1g 0 # g-1 403 -> E1gy 3 # g+1 402 -> E1gx 2 # g-2 411 -> E2gy 11 # g+2 410 -> E2gx 10 # g-3 413 -> E3gy 13 # g+3 412 -> E3gx 12 # g-4 421 -> E4gy 21 # g+4 420 -> E4gx 20 _SO3_SYMB2ID = { 's+0' : 0, 'p-1': 106, 'p+0': 105, 'p+1': 107, 'd-2': 211, 'd-1': 203, 'd+0': 200, 'd+1': 202, 'd+2': 210, 'f-3': 316, 'f-2': 314, 'f-1': 306, 'f+0': 305, 'f+1': 307, 'f+2': 315, 'f+3': 317, 'g-4': 421, 'g-3': 413, 'g-2': 411, 'g-1': 403, 'g+0': 400, 'g+1': 402, 'g+2': 410, 'g+3': 412, 'g+4': 420, 'h-5': 526, 'h-4': 524, 'h-3': 516, 'h-2': 514, 'h-1': 506, 'h+0': 505, 'h+1': 507, 'h+2': 515, 'h+3': 517, 'h+4': 525, 'h+5': 527, 'i-6': 631, 'i-5': 623, 'i-4': 621, 'i-3': 613, 'i-2': 611, 'i-1': 603, 'i+0': 600, 'i+1': 602, 'i+2': 610, 'i+3': 612, 'i+4': 620, 'i+5': 622, 'i+6': 630, } _SO3_ID2SYMB = {v: k for k, v in _SO3_SYMB2ID.items()} _ANGULAR = 'spdfghiklmnortu'
[docs] def so3_irrep_symb2id(symb): symb = symb.lower() if symb in _SO3_SYMB2ID: return _SO3_SYMB2ID[symb] else: s = ''.join([i for i in symb if i.isalpha()]) n = int(''.join([i for i in symb if not i.isalpha()])) e2 = abs(n) // 2 l = _ANGULAR.index(s) if abs(n) > l: raise KeyError(symb) gu = 'u' if l % 2 else 'g' xy = 'y' if n < 0 else 'x' parity = '_odd' if n % 2 else '_even' return l * 100 + e2 * 10 + DOOH_IRREP_ID_TABLE[parity + gu + xy]
[docs] def so3_irrep_id2symb(irrep_id): if irrep_id in _SO3_ID2SYMB: return _SO3_ID2SYMB[irrep_id] else: l, e2 = divmod(irrep_id, 100) e2, n = divmod(e2, 10) if n in (0, 1, 5, 4): rn = e2 * 2 elif n in (2, 3, 6, 7): rn = e2 * 2 + 1 else: raise KeyError(f'Unknown Irrep ID {irrep_id} for SO3 group') if l <= 1 or rn > l: raise KeyError(f'Unknown Irrep ID {irrep_id} for SO3 group') if n in (6, 4, 3, 1): rn = -rn return f'{_ANGULAR[l]}{rn:+d}'
def so3_symm_adapted_basis(mol, gpname, orig=0, coordinates=None): assert gpname == 'SO3' assert mol.natm == 1 ao_loc = mol.ao_loc_nr(cart=False) nao_sph = ao_loc[-1] if mol.cart: coeff = mol.cart2sph(normalized='sp') else: coeff = numpy.eye(nao_sph) nao = coeff.shape[0] lmax = max(mol._bas[:,gto.ANG_OF]) so = [] irrep_names = [] for l in range(lmax+1): bas_idx = mol._bas[:,gto.ANG_OF] == l if sum(bas_idx) == 0: # Skip any unused angular momentum channels continue cs = [coeff[:,p0:p1] for p0, p1 in zip(ao_loc[:-1][bas_idx], ao_loc[1:][bas_idx])] c_groups = numpy.hstack(cs).reshape(nao, -1, l*2+1) if l == 1: so.extend([c_groups[:,:,1], c_groups[:,:,2], c_groups[:,:,0]]) irrep_names.extend(['p-1', 'p+0', 'p+1']) else: for m in range(-l, l+1): so.append(c_groups[:,:,l+m]) irrep_names.append('%s%+d' % (_ANGULAR[l], m)) irrep_ids = [so3_irrep_symb2id(ir) for ir in irrep_names] return so, irrep_ids ############################### # Linear molecule # Irreps ID maps # Dooh -> D2h | Coov -> C2v # A1g 0 Ag 0 | A1 0 A1 0 # A2g 1 B1g 1 | A2 1 A2 1 # A1u 5 B1u 5 | E1x 2 B1 2 # A2u 4 Au 4 | E1y 3 B2 3 # E1gx 2 B2g 2 | E2x 10 A1 0 # E1gy 3 B3g 3 | E2y 11 A2 1 # E1ux 7 B3u 7 | E3x 12 B1 2 # E1uy 6 B2u 6 | E3y 13 B2 3 # E2gx 10 Ag 0 | E4x 20 A1 0 # E2gy 11 B1g 1 | E4y 21 A2 1 # E2ux 15 B1u 5 | E5x 22 B1 2 # E2uy 14 Au 4 | E5y 23 B2 3 # E3gx 12 B2g 2 | # E3gy 13 B3g 3 | # E3ux 17 B3u 7 | # E3uy 16 B2u 6 | # E4gx 20 Ag 0 | # E4gy 21 B1g 1 | # E4ux 25 B1u 5 | # E4uy 24 Au 4 | # E5gx 22 B2g 2 | # E5gy 23 B3g 3 | # E5ux 27 B3u 7 | # E5uy 26 B2u 6 | DOOH_IRREP_ID_TABLE = { 'A1g' : 0, 'A2g' : 1, 'A1u' : 5, 'A2u' : 4, 'E1gx': 2, 'E1gy': 3, 'E1ux': 7, 'E1uy': 6, '_evengx': 0, '_evengy': 1, '_evenux': 5, '_evenuy': 4, '_oddgx': 2, '_oddgy': 3, '_oddux': 7, '_odduy': 6, } COOV_IRREP_ID_TABLE = { 'A1' : 0, 'A2' : 1, 'E1x': 2, 'E1y': 3, '_evenx': 0, '_eveny': 1, '_oddx': 2, '_oddy': 3, }
[docs] def linearmole_symm_descent(gpname, irrep_id): '''Map irreps to D2h or C2v''' if gpname in ('Dooh', 'Coov'): return irrep_id % 10 else: raise PointGroupSymmetryError('%s is not proper for linear molecule.' % gpname)
[docs] def linearmole_irrep_symb2id(gpname, symb): if gpname == 'Dooh': if symb in DOOH_IRREP_ID_TABLE: return DOOH_IRREP_ID_TABLE[symb] else: try: n = int(''.join([i for i in symb if i.isdigit()])) if n % 2: return (n//2)*10 + DOOH_IRREP_ID_TABLE['_odd'+symb[-2:]] else: return (n//2)*10 + DOOH_IRREP_ID_TABLE['_even'+symb[-2:]] except (KeyError, ValueError): raise PointGroupSymmetryError(f'Incorrect Dooh irrep {symb}') elif gpname == 'Coov': if symb in COOV_IRREP_ID_TABLE: return COOV_IRREP_ID_TABLE[symb] else: if 'g' in symb or 'u' in symb: raise PointGroupSymmetryError(f'Incorrect Coov irrep {symb}') try: n = int(''.join([i for i in symb if i.isdigit()])) if n % 2: return (n//2)*10 + COOV_IRREP_ID_TABLE['_odd'+symb[-1]] else: return (n//2)*10 + COOV_IRREP_ID_TABLE['_even'+symb[-1]] except (KeyError, ValueError): raise PointGroupSymmetryError(f'Incorrect Coov irrep {symb}') else: raise PointGroupSymmetryError(f'Incorrect cylindrical symmetry group {gpname}')
DOOH_IRREP_SYMBS = ('A1g' , 'A2g' , 'E1gx', 'E1gy' , 'A2u', 'A1u' , 'E1uy', 'E1ux') DOOH_IRREP_SYMBS_EXT = ('gx' , 'gy' , 'gx', 'gy' , 'uy', 'ux' , 'uy', 'ux') COOV_IRREP_SYMBS = ('A1' , 'A2' , 'E1x', 'E1y')
[docs] def linearmole_irrep_id2symb(gpname, irrep_id): if gpname == 'Dooh': if irrep_id < 10: return DOOH_IRREP_SYMBS[irrep_id] else: l = abs(linearmole_irrep2momentum(irrep_id)) n = irrep_id % 10 return 'E%d%s' % (l, DOOH_IRREP_SYMBS_EXT[n]) elif gpname == 'Coov': if irrep_id < 10: return COOV_IRREP_SYMBS[irrep_id] else: l = abs(linearmole_irrep2momentum(irrep_id)) n = irrep_id % 10 if n >= 4: raise PointGroupSymmetryError(f'Incorrect Coov irrep {irrep_id}') if n % 2: xy = 'y' else: xy = 'x' return 'E%d%s' % (l, xy) else: raise PointGroupSymmetryError(f'Incorrect cylindrical symmetry group {gpname}')
def linearmole_irrep2momentum(irrep_id): if irrep_id % 10 in (0, 1, 5, 4): l = irrep_id // 10 * 2 else: l = irrep_id // 10 * 2 + 1 if irrep_id % 10 in (1, 3, 4, 6): # Ey l *= -1 return l
[docs] def linearmole_symm_adapted_basis(mol, gpname, orig=0, coordinates=None): assert (gpname in ('Dooh', 'Coov')) assert (not mol.cart) if coordinates is None: coordinates = numpy.eye(3) prop_atoms = mol.format_atom(mol._atom, orig, coordinates, 'Bohr') eql_atom_ids = geom.symm_identical_atoms(gpname, prop_atoms) aoslice = mol.aoslice_by_atom() basoff = aoslice[:,2] nao = mol.nao_nr() sodic = {} shalf = numpy.sqrt(.5) def plus(i0, i1): c = numpy.zeros(nao) c[i0] = c[i1] = shalf return c def minus(i0, i1): c = numpy.zeros(nao) c[i0] = shalf c[i1] =-shalf return c def identity(i0): c = numpy.zeros(nao) c[i0] = 1 return c def add_so(irrep_name, c): if irrep_name in sodic: sodic[irrep_name].append(c) else: sodic[irrep_name] = [c] if gpname == 'Dooh': for atom_ids in eql_atom_ids: if len(atom_ids) == 2: at0 = atom_ids[0] at1 = atom_ids[1] ip = 0 b0, b1, p0, p1 = aoslice[at0] for ib in range(b0, b1): angl = mol.bas_angular(ib) nc = mol.bas_nctr(ib) degen = angl * 2 + 1 if angl == 1: for i in range(nc): aoff = ip + i*degen + angl # m = 0 idx0 = basoff[at0] + aoff + 1 idx1 = basoff[at1] + aoff + 1 add_so('A1g', minus(idx0, idx1)) add_so('A1u', plus (idx0, idx1)) # m = +/- 1 idx0 = basoff[at0] + aoff - 1 idy0 = basoff[at0] + aoff idx1 = basoff[at1] + aoff - 1 idy1 = basoff[at1] + aoff add_so('E1ux', plus (idx0, idx1)) add_so('E1uy', plus (idy0, idy1)) add_so('E1gx', minus(idx0, idx1)) add_so('E1gy', minus(idy0, idy1)) else: for i in range(nc): aoff = ip + i*degen + angl # m = 0 idx0 = basoff[at0] + aoff idx1 = basoff[at1] + aoff if angl % 2: # p-sigma, f-sigma add_so('A1g', minus(idx0, idx1)) add_so('A1u', plus (idx0, idx1)) else: # s-sigma, d-sigma add_so('A1g', plus (idx0, idx1)) add_so('A1u', minus(idx0, idx1)) # +/-m for m in range(1,angl+1): idx0 = basoff[at0] + aoff + m idy0 = basoff[at0] + aoff - m idx1 = basoff[at1] + aoff + m idy1 = basoff[at1] + aoff - m if angl % 2: # odd parity add_so('E%dux'%m, plus (idx0, idx1)) add_so('E%duy'%m, plus (idy0, idy1)) add_so('E%dgx'%m, minus(idx0, idx1)) add_so('E%dgy'%m, minus(idy0, idy1)) else: add_so('E%dgy'%m, plus (idy0, idy1)) add_so('E%dgx'%m, plus (idx0, idx1)) add_so('E%duy'%m, minus(idy0, idy1)) add_so('E%dux'%m, minus(idx0, idx1)) ip += nc * degen elif len(atom_ids) == 1: at0 = atom_ids[0] ip = 0 b0, b1, p0, p1 = aoslice[at0] for ib in range(b0, b1): angl = mol.bas_angular(ib) nc = mol.bas_nctr(ib) degen = angl * 2 + 1 if angl == 1: for i in range(nc): aoff = ip + i*degen + angl # m = 0 idx0 = basoff[at0] + aoff + 1 add_so('A1u', identity(idx0)) # m = +/- 1 idx0 = basoff[at0] + aoff - 1 idy0 = basoff[at0] + aoff add_so('E1uy', identity(idy0)) add_so('E1ux', identity(idx0)) else: for i in range(nc): aoff = ip + i*degen + angl idx0 = basoff[at0] + aoff # m = 0 if angl % 2: add_so('A1u', identity(idx0)) else: add_so('A1g', identity(idx0)) # +/-m for m in range(1,angl+1): idx0 = basoff[at0] + aoff + m idy0 = basoff[at0] + aoff - m if angl % 2: # p, f functions add_so('E%dux'%m, identity(idx0)) add_so('E%duy'%m, identity(idy0)) else: # d, g functions add_so('E%dgy'%m, identity(idy0)) add_so('E%dgx'%m, identity(idx0)) ip += nc * degen elif gpname == 'Coov': for atom_ids in eql_atom_ids: at0 = atom_ids[0] ip = 0 b0, b1, p0, p1 = aoslice[at0] for ib in range(b0, b1): angl = mol.bas_angular(ib) nc = mol.bas_nctr(ib) degen = angl * 2 + 1 if angl == 1: for i in range(nc): aoff = ip + i*degen + angl # m = 0 idx0 = basoff[at0] + aoff + 1 add_so('A1', identity(idx0)) # m = +/- 1 idx0 = basoff[at0] + aoff - 1 idy0 = basoff[at0] + aoff add_so('E1x', identity(idx0)) add_so('E1y', identity(idy0)) else: for i in range(nc): aoff = ip + i*degen + angl idx0 = basoff[at0] + aoff # m = 0 add_so('A1', identity(idx0)) # +/-m for m in range(1,angl+1): idx0 = basoff[at0] + aoff + m idy0 = basoff[at0] + aoff - m add_so('E%dx'%m, identity(idx0)) add_so('E%dy'%m, identity(idy0)) ip += nc * degen irrep_ids = [] irrep_names = list(sodic.keys()) for irname in irrep_names: irrep_ids.append(linearmole_irrep_symb2id(gpname, irname)) irrep_idx = numpy.argsort(irrep_ids) irrep_ids = [irrep_ids[i] for i in irrep_idx] ao_loc = mol.ao_loc_nr() l_idx = {} ANG_OF = 1 for l in range(mol._bas[:,ANG_OF].max()+1): idx = [numpy.arange(ao_loc[ib], ao_loc[ib+1]) for ib in numpy.where(mol._bas[:,ANG_OF] == l)[0]] if idx: l_idx[l] = numpy.hstack(idx) Ds = _momentum_rotation_matrices(mol, coordinates) so = [] for i in irrep_idx: c_ir = numpy.vstack(sodic[irrep_names[i]]).T nso = c_ir.shape[1] for l, idx in l_idx.items(): c = c_ir[idx].reshape(-1,Ds[l].shape[1],nso) c_ir[idx] = numpy.einsum('nm,smp->snp', Ds[l], c).reshape(-1,nso) so.append(c_ir) return so, irrep_ids
if __name__ == "__main__": h2o = gto.Mole() h2o.verbose = 0 h2o.output = None h2o.atom = [['O' , (1. , 0. , 0. ,)], [1 , (0. , -.757 , 0.587,)], [1 , (0. , 0.757 , 0.587,)] ] h2o.basis = {'H': 'cc-pvdz', 'O': 'cc-pvdz',} h2o.build() gpname, origin, axes = geom.detect_symm(h2o._atom) atoms = gto.format_atom(h2o._atom, origin, axes) h2o.build(False, False, atom=atoms) print(gpname) eql_atoms = geom.symm_identical_atoms(gpname, atoms) print(symm_adapted_basis(h2o, gpname, eql_atoms)[1]) mol = gto.M( atom = [['H', (0,0,0)], ['H', (0,0,-1)], ['H', (0,0,1)]], basis = 'ccpvtz', charge=1) gpname, orig, axes = geom.detect_symm(mol._atom) atoms = gto.format_atom(mol._atom, orig, axes) mol.build(False, False, atom=atoms) print(gpname) eql_atoms = geom.symm_identical_atoms(gpname, atoms) print(symm_adapted_basis(mol, gpname, eql_atoms)[1]) mol = gto.M( atom = [['H', (0,0,0)], ['H', (0,0,-1)], ['He', (0,0,1)]], basis = 'ccpvtz') gpname, orig, axes = geom.detect_symm(mol._atom) atoms = gto.format_atom(mol._atom, orig, axes) mol.build(False, False, atom=atoms) print(gpname) eql_atoms = geom.symm_identical_atoms(gpname, atoms) print(symm_adapted_basis(mol, gpname, eql_atoms)[1])