Skip to content

Commit

Permalink
Fix basis order issue in molden.load (issue pyscf#1961)
Browse files Browse the repository at this point in the history
  • Loading branch information
@sunqm
sunqm committed Dec 11, 2023
1 parent efe8f52 commit 37be001
Show file tree
Hide file tree
Showing 3 changed files with 85 additions and 38 deletions.
59 changes: 23 additions & 36 deletions pyscf/gto/mole.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -419,7 +419,7 @@ def str2atm(line):
return list(zip(z, c.tolist()))

#TODO: sort exponents
def format_basis(basis_tab):
def format_basis(basis_tab, sort_basis=True):
'''Convert the input :attr:`Mole.basis` to the internal data format.
``{ atom: [(l, ((-exp, c_1, c_2, ..),
Expand Down Expand Up @@ -457,10 +457,16 @@ def format_basis(basis_tab):
fmt_basis = {}
for atom, atom_basis in basis_tab.items():
symb = _atom_symbol(atom)
fmt_basis[symb] = basis_converter(symb, atom_basis)

if len(fmt_basis[symb]) == 0:
_basis = basis_converter(symb, atom_basis)
if len(_basis) == 0:
raise BasisNotFoundError('Basis not found for %s' % symb)

# Sort basis accroding to angular momentum. This is important for method
# decontract_basis, which assumes that basis functions with the same
# angular momentum are grouped together. Related to issue #1620 #1770
if sort_basis:
_basis = sorted([b for b in _basis if b], key=lambda b: b[0])
fmt_basis[symb] = _basis
return fmt_basis

def _generate_basis_converter():
Expand Down Expand Up @@ -503,11 +509,6 @@ def converter(symb, raw_basis):
_basis.extend(nparray_to_list(rawb))
else:
_basis = nparray_to_list(raw_basis)

# Sort basis accroding to angular momentum. This is important for method
# decontract_basis, which assumes that basis functions with the same
# angular momentum are grouped together. Related to issue #1620 #1770
_basis = sorted([b for b in _basis if b], key=lambda b: b[0])
return _basis
return converter

Expand Down Expand Up @@ -654,6 +655,8 @@ def _to_full_contraction(mol, bas_idx):
bas_idx = ib0 + numpy.where(mol._bas[ib0:ib1,ANG_OF] == l)[0]
if len(bas_idx) == 0:
continue
if bas_idx[0] + len(bas_idx) != bas_idx[-1] + 1:
raise NotImplementedError('Discontinuous bases of same angular momentum')

mol_exps, b_coeff = _to_full_contraction(mol, bas_idx)
nprim, nc = b_coeff.shape
Expand Down Expand Up @@ -2523,11 +2526,13 @@ def build(self, dump_input=True, parse_arg=ARGPARSE,
if self.verbose >= logger.WARN:
self.check_sanity()

self._atom = self.format_atom(self.atom, unit=self.unit)
if self.atom:
self._atom = self.format_atom(self.atom, unit=self.unit)
uniq_atoms = {a[0] for a in self._atom}

_basis = _parse_default_basis(self.basis, uniq_atoms)
self._basis = self.format_basis(_basis)
if self.basis:
_basis = _parse_default_basis(self.basis, uniq_atoms)
self._basis = self.format_basis(_basis)
env = self._env[:PTR_ENV_START]
self._atm, self._bas, self._env = \
self.make_env(self._atom, self._basis, env, self.nucmod,
Expand Down Expand Up @@ -2643,30 +2648,12 @@ def _build_symmetry(self, *args, **kwargs):
for ir in self.irrep_id]
return self

@lib.with_doc(format_atom.__doc__)
def format_atom(self, atom, origin=0, axes=None, unit='Ang'):
return format_atom(atom, origin, axes, unit)

@lib.with_doc(format_basis.__doc__)
def format_basis(self, basis_tab):
return format_basis(basis_tab)

@lib.with_doc(format_pseudo.__doc__)
def format_pseudo(self, pseudo_tab):
return format_pseudo(pseudo_tab)

@lib.with_doc(format_ecp.__doc__)
def format_ecp(self, ecp_tab):
return format_ecp(ecp_tab)

@lib.with_doc(expand_etb.__doc__)
def expand_etb(self, l, n, alpha, beta):
return expand_etb(l, n, alpha, beta)

@lib.with_doc(expand_etbs.__doc__)
def expand_etbs(self, etbs):
return expand_etbs(etbs)
etbs = expand_etbs
format_atom = staticmethod(format_atom)
format_basis = staticmethod(format_basis)
format_pseudo = staticmethod(format_pseudo)
format_ecp = staticmethod(format_ecp)
expand_etb = staticmethod(expand_etb)
expand_etbs = etbs = staticmethod(expand_etbs)

@lib.with_doc(make_env.__doc__)
def make_env(self, atoms, basis, pre_env=[], nucmod={}, nucprop=None):
Expand Down
17 changes: 15 additions & 2 deletions pyscf/tools/molden.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -239,8 +239,16 @@ def read_one_bas(lsym, nb, fac=1):
elif dat[0].upper() in 'SPDFGHIJ':
basis[symb].append(read_one_bas(*dat))

mol.basis = envs['basis'] = basis
mol.atom = [atoms[i] for i in atom_seq]
uniq_atoms = {a[0] for a in mol.atom}

# To avoid the mol.build() sort the basis, disable mol.basis and set the
# internal data _basis directly. It is a workaround to solve issue #1961.
# Mole.decontract_basis function should be rewritten to support
# discontinuous bases that have the same angular momentum.
mol.basis = {}
_basis = gto.mole._parse_default_basis(basis, uniq_atoms)
mol._basis = envs['basis'] = gto.format_basis(_basis, sort_basis=False)
return mol

def _parse_mo(lines, envs):
Expand Down Expand Up @@ -353,6 +361,8 @@ def load(moldenfile, verbose=0):
else:
sys.stderr.write('Unknown section %s\n' % sec_title)

mo_energy, mo_coeff, mo_occ, irrep_labels, spins = None, None, None, None, None

for sec_kind in _SEC_ORDER:
if sec_kind == 'MO' and 'MO' in sec_kinds:
if len(sec_kinds['MO']) == 1:
Expand Down Expand Up @@ -392,6 +402,8 @@ def load(moldenfile, verbose=0):
if isinstance(mo_occ, tuple):
mol.spin = int(mo_occ[0].sum() - mo_occ[1].sum())

if not mol._built:
mol.build(0, 0)
return mol, mo_energy, mo_coeff, mo_occ, irrep_labels, spins

parse = read = load
Expand Down Expand Up @@ -503,8 +515,9 @@ def remove_high_l(mol, mo_coeff=None):
'''
pmol = mol.copy()
pmol.basis = {}
pmol._basis = {}
for symb, bas in mol._basis.items():
pmol.basis[symb] = [b for b in bas if b[0] <= 4]
pmol._basis[symb] = [b for b in bas if b[0] <= 4]
pmol.build(0, 0)
if mo_coeff is None:
return pmol, None
Expand Down
47 changes: 47 additions & 0 deletions pyscf/tools/test/test_molden.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,53 @@ def test_dump_cartesian_gto_symm_orbital(self):
mo_coeff = res[2]
self.assertAlmostEqual(abs(mf.mo_coeff-mo_coeff).max(), 0, 12)

def test_basis_not_sorted(self):
with tempfile.NamedTemporaryFile('w') as ftmp:
ftmp.write('''\
[Molden Format]
made by pyscf v[2.4.0]
[Atoms] (AU)
Ne 1 10 0.00000000000000 0.00000000000000 0.00000000000000
[GTO]
1 0
s 8 1.00
17880 0.00073769817327139
2683 0.0056746782244738
611.5 0.028871187450674
173.5 0.10849560938601
56.64 0.29078802505672
20.42 0.44814064476114
7.81 0.25792047270532
1.653 0.015056839544698
s 8 1.00
17880 -0.00033214909943481
2683 -0.0026205019065875
611.5 -0.013009816761002
173.5 -0.053420003125961
56.64 -0.14716522424261
20.42 -0.33838075724805
7.81 -0.20670101921688
1.653 1.0950299234565
s 1 1.00
0.4869 1
p 3 1.00
28.39 0.066171986640049
6.27 0.34485329752845
1.695 0.73045763818875
p 1 1.00
0.4317 1
d 1 1.00
2.202 1
s 1 1.00
1 1
[5d]
[7f]
[9g]
''')
ftmp.flush()
mol = molden.load(ftmp.name)[0]
self.assertEqual(mol._bas[:,1].tolist(), [0, 0, 0, 1, 1, 2, 0])

if __name__ == "__main__":
print("Full Tests for molden")
Expand Down

0 comments on commit 37be001

Please sign in to comment.