diff --git a/pyscf/cc/ccsd.py b/pyscf/cc/ccsd.py index 2e2d8482ae..15df586dc4 100644 --- a/pyscf/cc/ccsd.py +++ b/pyscf/cc/ccsd.py @@ -821,7 +821,7 @@ def __init__(self, cc): self._scf = cc._scf.as_scanner() def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) @@ -934,12 +934,6 @@ class CCSDBase(lib.StreamObject): )) def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None): - if isinstance(mf, gto.Mole): - raise RuntimeError(''' -You see this error message because of the API updates in pyscf v0.10. -In the new API, the first argument of CC class is HF objects. Please see -http://sunqm.net/pyscf/code-rule.html#api-rules for the details of API conventions''') - from pyscf.scf import hf if isinstance(mf, hf.KohnShamDFT): raise RuntimeError( diff --git a/pyscf/ci/cisd.py b/pyscf/ci/cisd.py index 05e9924fc0..7c9fd08e27 100644 --- a/pyscf/ci/cisd.py +++ b/pyscf/ci/cisd.py @@ -23,7 +23,7 @@ from functools import reduce import numpy -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf.lib import logger from pyscf.cc import ccsd @@ -786,7 +786,7 @@ def __init__(self, ci): self._scf = ci._scf.as_scanner() def __call__(self, mol_or_geom, ci0=None, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.Mole): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/df/grad/casscf.py b/pyscf/df/grad/casscf.py index b7bd0fb87a..9ed69841eb 100644 --- a/pyscf/df/grad/casscf.py +++ b/pyscf/df/grad/casscf.py @@ -34,7 +34,7 @@ from itertools import product import numpy from scipy import linalg -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf import ao2mo from pyscf.lib import logger @@ -160,7 +160,8 @@ def __init__(self, g): lib.GradScanner.__init__(self, g) def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/df/grad/sacasscf.py b/pyscf/df/grad/sacasscf.py index cef6373d2c..828b13febe 100644 --- a/pyscf/df/grad/sacasscf.py +++ b/pyscf/df/grad/sacasscf.py @@ -21,7 +21,7 @@ import gc from functools import reduce from scipy import linalg -from pyscf.gto import Mole +from pyscf import gto from pyscf import mcscf, lib, ao2mo from pyscf.grad import lagrange from pyscf.grad import rhf as rhf_grad @@ -331,7 +331,8 @@ def __init__(self, g, state): self.state = state def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/df/incore.py b/pyscf/df/incore.py index f89215b8c3..cc26b1ed55 100644 --- a/pyscf/df/incore.py +++ b/pyscf/df/incore.py @@ -46,7 +46,7 @@ def aux_e2(mol, auxmol_or_auxbasis, intor='int3c2e', aosym='s1', comp=None, out= cintopt = gto.moleintor.make_cintopt(mol._atm, mol._bas, mol._env, 'int3c2e') ''' - if isinstance(auxmol_or_auxbasis, gto.Mole): + if isinstance(auxmol_or_auxbasis, gto.MoleBase): auxmol = auxmol_or_auxbasis else: auxbasis = auxmol_or_auxbasis @@ -76,7 +76,7 @@ def aux_e1(mol, auxmol_or_auxbasis, intor='int3c2e', aosym='s1', comp=None, out= The same arguments as function aux_e2 can be input to aux_e1. ''' - if isinstance(auxmol_or_auxbasis, gto.Mole): + if isinstance(auxmol_or_auxbasis, gto.MoleBase): auxmol = auxmol_or_auxbasis else: auxbasis = auxmol_or_auxbasis @@ -96,7 +96,7 @@ def aux_e1(mol, auxmol_or_auxbasis, intor='int3c2e', aosym='s1', comp=None, out= def fill_2c2e(mol, auxmol_or_auxbasis, intor='int2c2e', comp=None, hermi=1, out=None): '''2-center 2-electron AO integrals for auxiliary basis (auxmol) ''' - if isinstance(auxmol_or_auxbasis, gto.Mole): + if isinstance(auxmol_or_auxbasis, gto.MoleBase): auxmol = auxmol_or_auxbasis else: auxbasis = auxmol_or_auxbasis diff --git a/pyscf/grad/casci.py b/pyscf/grad/casci.py index 144f202436..715694e3c3 100644 --- a/pyscf/grad/casci.py +++ b/pyscf/grad/casci.py @@ -27,7 +27,7 @@ from functools import reduce import numpy -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf import ao2mo from pyscf.lib import logger @@ -226,7 +226,8 @@ def __init__(self, g, state): self.state = state def __call__(self, mol_or_geom, state=None, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/casscf.py b/pyscf/grad/casscf.py index 844ea13227..0ce78a82b3 100644 --- a/pyscf/grad/casscf.py +++ b/pyscf/grad/casscf.py @@ -30,7 +30,7 @@ from functools import reduce import numpy -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf import ao2mo from pyscf.lib import logger @@ -171,7 +171,8 @@ def __init__(self, g): lib.GradScanner.__init__(self, g) def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/ccsd.py b/pyscf/grad/ccsd.py index 2d8841c0dc..6ab4a5373b 100644 --- a/pyscf/grad/ccsd.py +++ b/pyscf/grad/ccsd.py @@ -25,7 +25,7 @@ import numpy from functools import reduce from pyscf import lib -from pyscf.gto import Mole +from pyscf import gto from pyscf.lib import logger from pyscf.cc import ccsd from pyscf.cc import _ccsd @@ -231,7 +231,8 @@ def __init__(self, g): lib.GradScanner.__init__(self, g) def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/cisd.py b/pyscf/grad/cisd.py index 692db23ba3..1a57c86e8f 100644 --- a/pyscf/grad/cisd.py +++ b/pyscf/grad/cisd.py @@ -21,7 +21,7 @@ ''' import numpy -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf.lib import logger from pyscf.ci import cisd @@ -89,7 +89,8 @@ def __init__(self, g, state): self.state = state def __call__(self, mol_or_geom, state=None, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/mp2.py b/pyscf/grad/mp2.py index 2e72892926..58dcc92cf5 100644 --- a/pyscf/grad/mp2.py +++ b/pyscf/grad/mp2.py @@ -24,7 +24,7 @@ import numpy from pyscf import lib from functools import reduce -from pyscf.gto import Mole +from pyscf import gto from pyscf.lib import logger from pyscf.grad import rhf as rhf_grad from pyscf.scf import cphf @@ -224,7 +224,8 @@ def __init__(self, g): lib.GradScanner.__init__(self, g) def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/rhf.py b/pyscf/grad/rhf.py index 2ec791aeee..708a23ae21 100644 --- a/pyscf/grad/rhf.py +++ b/pyscf/grad/rhf.py @@ -254,7 +254,8 @@ def __init__(self, g): lib.GradScanner.__init__(self, g) def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/sacasscf.py b/pyscf/grad/sacasscf.py index 13f086713b..190fddbb83 100644 --- a/pyscf/grad/sacasscf.py +++ b/pyscf/grad/sacasscf.py @@ -19,7 +19,7 @@ from functools import reduce from itertools import product from scipy import linalg -from pyscf.gto import Mole +from pyscf import gto from pyscf.grad import lagrange from pyscf.mcscf.addons import StateAverageMCSCFSolver, StateAverageFCISolver from pyscf.mcscf.addons import StateAverageMixFCISolver, state_average_mix_ @@ -374,7 +374,8 @@ def __init__(self, g, state): self._converged = False def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/grad/tdrhf.py b/pyscf/grad/tdrhf.py index e4c73735d3..9a3694b86e 100644 --- a/pyscf/grad/tdrhf.py +++ b/pyscf/grad/tdrhf.py @@ -22,7 +22,7 @@ from functools import reduce import numpy -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf.lib import logger from pyscf.grad import rhf as rhf_grad @@ -207,7 +207,8 @@ def __init__(self, g, state): self.state = state def __call__(self, mol_or_geom, state=None, **kwargs): - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): + assert mol_or_geom.__class__ == gto.Mole mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/gto/basis/__init__.py b/pyscf/gto/basis/__init__.py index d8e22ec5ee..9d0706d359 100644 --- a/pyscf/gto/basis/__init__.py +++ b/pyscf/gto/basis/__init__.py @@ -608,27 +608,27 @@ def load(filename_or_basisname, symb, optimize=OPTIMIZE_CONTRACTION): try: return parse_nwchem.parse(filename_or_basisname, symb, optimize=optimize) - except IndexError: - raise BasisNotFoundError(filename_or_basisname) except BasisNotFoundError as basis_err: pass + except Exception: + raise BasisNotFoundError(filename_or_basisname) try: return parse_nwchem.parse(filename_or_basisname, optimize=optimize) - except IndexError: - raise BasisNotFoundError(f'Invalid basis {filename_or_basisname}') except BasisNotFoundError: pass + except Exception: + raise BasisNotFoundError(f'Invalid basis {filename_or_basisname}') try: return parse_cp2k.parse(filename_or_basisname, optimize=optimize) - except IndexError: - raise BasisNotFoundError(f'Invalid basis {filename_or_basisname}') except BasisNotFoundError: pass + except Exception: + raise BasisNotFoundError(f'Invalid basis {filename_or_basisname}') # Last, a trial to access Basis Set Exchange database - from pyscf.basis import bse + from pyscf.gto.basis import bse if bse.basis_set_exchange is not None: try: bse_obj = bse.basis_set_exchange.api.get_basis( diff --git a/pyscf/gto/basis/bse.py b/pyscf/gto/basis/bse.py index b9b6e9561e..061d838a88 100644 --- a/pyscf/gto/basis/bse.py +++ b/pyscf/gto/basis/bse.py @@ -5,6 +5,7 @@ ''' try: + import basis_set_exchange from basis_set_exchange import lut, manip, sort except ImportError: basis_set_exchange = None diff --git a/pyscf/gto/basis/parse_cp2k.py b/pyscf/gto/basis/parse_cp2k.py index 4701fb0bad..bc2b534fe3 100644 --- a/pyscf/gto/basis/parse_cp2k.py +++ b/pyscf/gto/basis/parse_cp2k.py @@ -61,33 +61,49 @@ def load(basisfile, symb, optimize=False): return _parse(search_seg(basisfile, symb), optimize) def _parse(blines, optimize=False): - header_ln = blines.pop(0) # noqa: F841 - nsets = int(blines.pop(0)) + blines_iter = iter(blines) + try: + header_ln = next(blines_iter) # noqa: F841 + nsets = int(next(blines_iter)) + except Exception: + raise BasisNotFoundError('Not basis data') + basis = [] - for n in range(nsets): - comp = [int(p) for p in blines.pop(0).split()] - lmin, lmax, nexps, ncontractions = comp[1], comp[2], comp[3], comp[4:] - basis_n = [[l] for l in range(lmin,lmax+1)] - for nexp in range(nexps): - line = blines.pop(0) - dat = line.split() - try: - bfun = [float(x) for x in dat] - except ValueError: - if DISABLE_EVAL: - raise ValueError('Failed to parse basis %s' % line) - else: - bfun = eval(','.join(dat)) - exp = bfun.pop(0) - for i,l in enumerate(range(lmin,lmax+1)): - cl = [exp] - for c in range(ncontractions[i]): - cl.append(bfun.pop(0)) - basis_n[i].append(tuple(cl)) - basis.extend(basis_n) + try: + for n in range(nsets): + comp = [int(p) for p in next(blines_iter).split()] + lmin, lmax, nexps, ncontractions = comp[1], comp[2], comp[3], comp[4:] + basis_n = [[l] for l in range(lmin,lmax+1)] + for nexp in range(nexps): + line = next(blines_iter) + dat = line.split() + try: + bfun = [float(x) for x in dat] + except ValueError: + if DISABLE_EVAL: + raise ValueError('Failed to parse %s' % line) + else: + bfun = eval(','.join(dat)) + + if len(bfun) != sum(ncontractions) + 1: + raise ValueError('Basis data incomplete') + + bfun_iter = iter(bfun) + exp = next(bfun_iter) + for i,l in enumerate(range(lmin,lmax+1)): + cl = [exp] + for c in range(ncontractions[i]): + cl.append(next(bfun_iter)) + basis_n[i].append(cl) + basis.extend(basis_n) + except StopIteration: + raise ValueError('Basis data incomplete') + basis_sorted = [] for l in range(MAXL): basis_sorted.extend([b for b in basis if b[0] == l]) + if not basis_sorted: + raise BasisNotFoundError('Basis data not found') if optimize: basis_sorted = parse_nwchem.optimize_contraction(basis_sorted) diff --git a/pyscf/gto/basis/parse_molpro.py b/pyscf/gto/basis/parse_molpro.py index 575598056b..13e1fa2995 100644 --- a/pyscf/gto/basis/parse_molpro.py +++ b/pyscf/gto/basis/parse_molpro.py @@ -24,6 +24,7 @@ import re import numpy +import numpy as np try: from pyscf.gto.basis.parse_nwchem import optimize_contraction diff --git a/pyscf/gto/basis/parse_nwchem.py b/pyscf/gto/basis/parse_nwchem.py index 0bf5f59bea..243b8b1c78 100644 --- a/pyscf/gto/basis/parse_nwchem.py +++ b/pyscf/gto/basis/parse_nwchem.py @@ -26,6 +26,7 @@ import re import numpy +import numpy as np import scipy.linalg from pyscf.data.elements import _std_symbol from pyscf.lib.exceptions import BasisNotFoundError @@ -105,28 +106,34 @@ def _parse(raw_basis, optimize=True): if not dat or dat.startswith('#'): continue elif dat[0].isalpha(): - key = dat.split()[1].upper() + keys = dat.split() + if len(keys) == 1: + key = keys[0].upper() + else: + key = keys[1].upper() if key == 'SP': basis_parsed[0].append([0]) basis_parsed[1].append([1]) - else: + elif key in MAPSPDF: l = MAPSPDF[key] current_basis = [l] basis_parsed[l].append(current_basis) + else: + raise BasisNotFoundError('Not basis data') else: dat = dat.replace('D','e').split() try: dat = [float(x) for x in dat] except ValueError: if DISABLE_EVAL: - raise ValueError('Failed to parse basis %s' % line) + raise ValueError('Failed to parse %s' % line) else: dat = list(eval(','.join(dat))) except Exception as e: raise BasisNotFoundError('\n' + str(e) + '\nor the required basis file not existed.') if key is None: - raise RuntimeError('Failed to parse basis') + raise BasisNotFoundError('Not basis data') elif key == 'SP': basis_parsed[0][-1].append([dat[0], dat[1]]) basis_parsed[1][-1].append([dat[0], dat[2]]) diff --git a/pyscf/gto/basis/parse_nwchem_ecp.py b/pyscf/gto/basis/parse_nwchem_ecp.py index f768abb524..f0b93854c6 100644 --- a/pyscf/gto/basis/parse_nwchem_ecp.py +++ b/pyscf/gto/basis/parse_nwchem_ecp.py @@ -22,6 +22,8 @@ __all__ = ['parse', 'load', 'convert_ecp_to_nwchem'] +import numpy +import numpy as np import re from pyscf.data.elements import _std_symbol from pyscf.lib.exceptions import BasisNotFoundError @@ -78,14 +80,20 @@ def _parse_ecp(raw_ecp): if not dat or dat.startswith('#'): # comment line continue elif dat[0].isalpha(): - key = dat.split()[1].upper() + keys = dat.split() + if len(keys) == 1: + key = keys[0].upper() + else: + key = keys[1].upper() if key == 'NELEC': nelec = int(dat.split()[2]) continue elif key == 'UL': ecp_add.append([-1]) - else: + elif key in MAPSPDF: ecp_add.append([MAPSPDF[key]]) + else: + raise BasisNotFoundError('Not basis data') # up to r^6 by_ang = [[] for i in range(7)] ecp_add[-1].append(by_ang) @@ -99,6 +107,8 @@ def _parse_ecp(raw_ecp): raise ValueError('Failed to parse ecp %s' % line) else: coef = list(eval(','.join(line[1:]))) + except Exception: + raise BasisNotFoundError('Not basis data') by_ang[l].append(coef) if nelec is None: diff --git a/pyscf/gto/mole.py b/pyscf/gto/mole.py index fb3f9632fc..05a58a2c88 100644 --- a/pyscf/gto/mole.py +++ b/pyscf/gto/mole.py @@ -2157,7 +2157,7 @@ def is_au(unit): return unit.upper().startswith(('B', 'AU')) # -# GTOIntegralMixin handed three layers of basis data: input, internal format, libcint arguments. +# MoleBase handles three layers of basis data: input, internal format, libcint arguments. # The relationship of the three layers are # .atom (input) <=> ._atom (for python) <=> ._atm (for libcint) # .basis (input) <=> ._basis (for python) <=> ._bas (for libcint) @@ -2166,7 +2166,132 @@ def is_au(unit): # on the internal format. Exceptions are make_env, make_atm_env, make_bas_env, # set_common_orig_, set_rinv_orig_ which are used to manipulate the libcint arguments. # -class GTOIntegralMixin: +class MoleBase(lib.StreamObject): + '''Basic class to hold molecular structure, integrals and global options + + Attributes: + verbose : int + Print level + output : str or None + Output file, default is None which dumps msg to sys.stdout + max_memory : int, float + Allowed memory in MB + charge : int + Charge of molecule. It affects the electron numbers + spin : int or None + 2S, num. alpha electrons - num. beta electrons to control + multiplicity. If spin = None is set, multiplicity will be guessed + based on the neutral molecule. + symmetry : bool or str + Whether to use symmetry. When this variable is set to True, the + molecule will be rotated and the highest rotation axis will be + placed z-axis. + If a string is given as the name of point group, the given point + group symmetry will be used. Note that the input molecular + coordinates will not be changed in this case. + symmetry_subgroup : str + subgroup + + atom : list or str + To define molecluar structure. The internal format is + + | atom = [[atom1, (x, y, z)], + | [atom2, (x, y, z)], + | ... + | [atomN, (x, y, z)]] + + unit : str + Angstrom or Bohr + basis : dict or str + To define basis set. + nucmod : dict or str or [function(nuc_charge, nucprop) => zeta] + Nuclear model. 0 or None means point nuclear model. Other + values will enable Gaussian nuclear model. If a function is + assigned to this attribute, the function will be called to + generate the nuclear charge distribution value "zeta" and the + relevant nuclear model will be set to Gaussian model. + Default is point nuclear model. + nucprop : dict + Nuclear properties (like g-factor 'g', quadrupole moments 'Q'). + It is needed by pyscf.prop module and submodules. + cart : boolean + Using Cartesian GTO basis and integrals (6d,10f,15g) + magmom : list + Collinear spin of each atom. Default is [0,]*natm + + ** Following attributes are generated by :func:`Mole.build` ** + + stdout : file object + Default is sys.stdout if :attr:`Mole.output` is not set + topgroup : str + Point group of the system. + groupname : str + The supported subgroup of the point group. It can be one of SO3, + Dooh, Coov, D2h, C2h, C2v, D2, Cs, Ci, C2, C1 + nelectron : int + sum of nuclear charges - :attr:`Mole.charge` + symm_orb : a list of numpy.ndarray + Symmetry adapted basis. Each element is a set of symm-adapted orbitals + for one irreducible representation. The list index does **not** correspond + to the id of irreducible representation. + irrep_id : a list of int + Each element is one irreducible representation id associated with the basis + stored in symm_orb. One irrep id stands for one irreducible representation + symbol. The irrep symbol and the relevant id are defined in + :attr:`symm.param.IRREP_ID_TABLE` + irrep_name : a list of str + Each element is one irreducible representation symbol associated with the basis + stored in symm_orb. The irrep symbols are defined in + :attr:`symm.param.IRREP_ID_TABLE` + _built : bool + To label whether :func:`Mole.build` has been called. It is to + ensure certain functions being initialized only once. + _basis : dict + like :attr:`Mole.basis`, the internal format which is returned from the + parser :func:`format_basis` + + ** Following attributes are arguments used by ``libcint`` library ** + + _atm : + :code:`[[charge, ptr-of-coord, nuc-model, ptr-zeta, 0, 0], [...]]` + each element represents one atom + natm : + number of atoms + _bas : + :code:`[[atom-id, angular-momentum, num-primitive-GTO, num-contracted-GTO, 0, ptr-of-exps, ptr-of-contract-coeff, 0], [...]]` + each element represents one shell + nbas : + number of shells + _env : + list of floats to store the coordinates, GTO exponents, contract-coefficients + + Examples: + + >>> mol = Mole(atom='H^2 0 0 0; H 0 0 1.1', basis='sto3g').build() + >>> print(mol.atom_symbol(0)) + H^2 + >>> print(mol.atom_pure_symbol(0)) + H + >>> print(mol.nao_nr()) + 2 + >>> print(mol.intor('int1e_ovlp_sph')) + [[ 0.99999999 0.43958641] + [ 0.43958641 0.99999999]] + >>> mol.charge = 1 + >>> mol.build() + has no attributes Charge + ''' # noqa: E501 + + output = None + max_memory = param.MAX_MEMORY + + verbose = getattr(__config__, 'VERBOSE', logger.NOTE) + + # the unit (angstrom/bohr) of the coordinates defined by the input self.atom + unit = getattr(__config__, 'UNIT', 'angstrom') + + # Whether to hold everything in memory + incore_anyway = getattr(__config__, 'INCORE_ANYWAY', False) # Using cartesian GTO (6d,10f,15g) cart = getattr(__config__, 'gto_mole_Mole_cart', False) @@ -2177,6 +2302,7 @@ class GTOIntegralMixin: # Store the keys appeared in the module. It is used to check misinput attributes _keys = { + 'verbose', 'unit', 'incore_anyway', 'output', 'max_memory', 'cart', 'charge', 'spin', 'symmetry', 'symmetry_subgroup', 'atom', 'basis', 'nucmod', 'ecp', 'nucprop', 'magmom', 'pseudo', 'groupname', 'topgroup', 'symm_orb', 'irrep_id', 'irrep_name', @@ -2203,6 +2329,7 @@ def __init__(self): self._env = numpy.zeros(PTR_ENV_START) self._ecpbas = numpy.zeros((0,8), dtype=numpy.int32) + self.stdout = sys.stdout self.groupname = 'C1' self.topgroup = 'C1' self.symm_orb = None @@ -2218,6 +2345,11 @@ def __init__(self): self._ecp = {} self._pseudo = {} + self._built = False + # Some methods modify ._env. These method are executed in the context + # _TemporaryMoleContext which is protected by the _ctx_lock. + self._ctx_lock = None + @property def natm(self): return len(self._atm) @@ -2276,227 +2408,592 @@ def ms(self, x): else: self.spin = int(round(2*x, 4)) - @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) + copy = copy - @lib.with_doc(format_pseudo.__doc__) - def format_pseudo(self, pseudo_tab): - return format_pseudo(pseudo_tab) + pack = pack - @lib.with_doc(format_ecp.__doc__) - def format_ecp(self, ecp_tab): - return format_ecp(ecp_tab) + @classmethod + @lib.with_doc(unpack.__doc__) + def unpack(cls, moldic): + return unpack(moldic) - @lib.with_doc(expand_etb.__doc__) - def expand_etb(self, l, n, alpha, beta): - return expand_etb(l, n, alpha, beta) + @lib.with_doc(unpack.__doc__) + def unpack_(self, moldic): + self.__dict__.update(moldic) + return self - @lib.with_doc(expand_etbs.__doc__) - def expand_etbs(self, etbs): - return expand_etbs(etbs) - etbs = expand_etbs + dumps = dumps - @lib.with_doc(make_env.__doc__) - def make_env(self, atoms, basis, pre_env=[], nucmod={}, nucprop=None): - if nucprop is None: - nucprop = self.nucprop - return make_env(atoms, basis, pre_env, nucmod, nucprop) + @classmethod + @lib.with_doc(loads.__doc__) + def loads(cls, molstr): + return loads(molstr) - @lib.with_doc(make_atm_env.__doc__) - def make_atm_env(self, atom, ptr=0, nucmod=NUC_POINT, nucprop=None): - if nucprop is None: - nucprop = self.nucprop.get(atom[0], {}) - return make_atm_env(atom, ptr, nucmod, nucprop) + @lib.with_doc(loads.__doc__) + def loads_(self, molstr): + self.__dict__.update(loads(molstr).__dict__) + return self - @lib.with_doc(make_bas_env.__doc__) - def make_bas_env(self, basis_add, atom_id=0, ptr=0): - return make_bas_env(basis_add, atom_id, ptr) + def build(self, dump_input=True, parse_arg=ARGPARSE, + verbose=None, output=None, max_memory=None, + atom=None, basis=None, unit=None, nucmod=None, ecp=None, pseudo=None, + charge=None, spin=0, symmetry=None, symmetry_subgroup=None, + cart=None, magmom=None): + '''Setup moleclue and initialize some control parameters. Whenever you + change the value of the attributes of :class:`Mole`, you need call + this function to refresh the internal data of Mole. - @lib.with_doc(make_ecp_env.__doc__) - def make_ecp_env(self, _atm, _ecp, pre_env=[]): - if _ecp: - _atm, _ecpbas, _env = make_ecp_env(self, _atm, _ecp, pre_env) - else: - _atm, _ecpbas, _env = _atm, numpy.zeros((0,BAS_SLOTS)), pre_env - return _atm, _ecpbas, _env + Kwargs: + dump_input : bool + whether to dump the contents of input file in the output file + parse_arg : bool + whether to read the sys.argv and overwrite the relevant parameters + verbose : int + Print level. If given, overwrite :attr:`Mole.verbose` + output : str or None + Output file. If given, overwrite :attr:`Mole.output` + max_memory : int, float + Allowd memory in MB. If given, overwrite :attr:`Mole.max_memory` + atom : list or str + To define molecluar structure. + basis : dict or str + To define basis set. + nucmod : dict or str + Nuclear model. If given, overwrite :attr:`Mole.nucmod` + charge : int + Charge of molecule. It affects the electron numbers + If given, overwrite :attr:`Mole.charge` + spin : int + 2S, num. alpha electrons - num. beta electrons to control + multiplicity. If setting spin = None , multiplicity will be + guessed based on the neutral molecule. + If given, overwrite :attr:`Mole.spin` + symmetry : bool or str + Whether to use symmetry. If given a string of point group + name, the given point group symmetry will be used. + magmom : list + Collinear spin of each atom. Default is [0.0,]*natm - tot_electrons = tot_electrons + ''' + if not DISABLE_GC: + gc.collect() # To release circular referred objects - @lib.with_doc(gto_norm.__doc__) - def gto_norm(self, l, expnt): - return gto_norm(l, expnt) + if isinstance(dump_input, str): + sys.stderr.write('Assigning the first argument %s to mol.atom\n' % + dump_input) + dump_input, atom = True, dump_input - def set_common_origin(self, coord): - '''Update common origin for integrals of dipole, rxp etc. - **Note** the unit of the coordinates needs to be Bohr + if verbose is not None: self.verbose = verbose + if output is not None: self.output = output + if max_memory is not None: self.max_memory = max_memory + if atom is not None: self.atom = atom + if basis is not None: self.basis = basis + if unit is not None: self.unit = unit + if nucmod is not None: self.nucmod = nucmod + if ecp is not None: self.ecp = ecp + if pseudo is not None: self.pseudo = pseudo + if charge is not None: self.charge = charge + if spin != 0: self.spin = spin + if symmetry is not None: self.symmetry = symmetry + if symmetry_subgroup is not None: self.symmetry_subgroup = symmetry_subgroup + if cart is not None: self.cart = cart + if magmom is not None: self.magmom = magmom - Examples: + if parse_arg: + _update_from_cmdargs_(self) - >>> mol.set_common_origin(0) - >>> mol.set_common_origin((1,0,0)) - ''' - self._env[PTR_COMMON_ORIG:PTR_COMMON_ORIG+3] = coord - return self - set_common_orig = set_common_origin - set_common_orig_ = set_common_orig # for backward compatibility - set_common_origin_ = set_common_orig # for backward compatibility - - def with_common_origin(self, coord): - '''Return a temporary mol context which has the rquired common origin. - The required common origin has no effects out of the temporary context. - See also :func:`mol.set_common_origin` - - Examples: + # avoid opening output file twice + if (self.output is not None + # StringIO() does not have attribute 'name' + and getattr(self.stdout, 'name', None) != self.output): - >>> with mol.with_common_origin((1,0,0)): - ... mol.intor('int1e_r', comp=3) - ''' - coord0 = self._env[PTR_COMMON_ORIG:PTR_COMMON_ORIG+3].copy() - return self._TemporaryMoleContext(self.set_common_origin, (coord,), (coord0,)) - with_common_orig = with_common_origin + if self.verbose > logger.QUIET: + if os.path.isfile(self.output): + print('overwrite output file: %s' % self.output) + else: + print('output file: %s' % self.output) - def set_rinv_origin(self, coord): - r'''Update origin for operator :math:`\frac{1}{|r-R_O|}`. - **Note** the unit is Bohr + if self.output == '/dev/null': + self.stdout = open(os.devnull, 'w') + else: + self.stdout = open(self.output, 'w') - Examples: + if self.verbose >= logger.WARN: + self.check_sanity() - >>> mol.set_rinv_origin(0) - >>> mol.set_rinv_origin((0,1,0)) - ''' - self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3] = coord[:3] - return self - set_rinv_orig = set_rinv_origin - set_rinv_orig_ = set_rinv_orig # for backward compatibility - set_rinv_origin_ = set_rinv_orig # for backward compatibility + self._atom = self.format_atom(self.atom, unit=self.unit) + uniq_atoms = set([a[0] for a in self._atom]) - def with_rinv_origin(self, coord): - '''Return a temporary mol context which has the rquired origin of 1/r - operator. The required origin has no effects out of the temporary - context. See also :func:`mol.set_rinv_origin` + _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, + self.nucprop) - Examples: + if self.pseudo: + self.ecp, self.pseudo = classify_ecp_pseudo(self, self.ecp, self.pseudo) - >>> with mol.with_rinv_origin((1,0,0)): - ... mol.intor('int1e_rinv') - ''' - coord0 = self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3].copy() - return self._TemporaryMoleContext(self.set_rinv_origin, (coord,), (coord0,)) - with_rinv_orig = with_rinv_origin + if self.ecp: + # Unless explicitly input, ECP should not be assigned to ghost atoms + atoms_wo_ghost = [a for a in uniq_atoms if not is_ghost_atom(a)] + _ecp = _parse_default_basis(self.ecp, atoms_wo_ghost) + self._ecp = self.format_ecp(_ecp) + if self._ecp: + self._atm, self._ecpbas, self._env = \ + self.make_ecp_env(self._atm, self._ecp, self._env) - def set_range_coulomb(self, omega): - '''Switch on range-separated Coulomb operator for **all** 2e integrals + if self.pseudo: + # Unless explicitly input, PP should not be assigned to ghost atoms + atoms_wo_ghost = [a for a in uniq_atoms if not is_ghost_atom(a)] + _pseudo = _parse_default_basis(self.pseudo, atoms_wo_ghost) + self._pseudo = _pseudo = self.format_pseudo(_pseudo) + if _pseudo: + conflicts = set(_pseudo).intersection(self._ecp) + if conflicts: + raise RuntimeError('Pseudo potential for atoms %s are defined ' + 'in both .ecp and .pseudo.' % list(conflicts)) - Args: - omega : double + for ia, atom in enumerate(self._atom): + symb = atom[0] + if (symb in _pseudo and + # skip ghost atoms + self._atm[ia,0] != 0): + self._atm[ia,0] = sum(_pseudo[symb][0]) - | = 0 : Regular electron repulsion integral - | > 0 : Long-range operator erf(omega r12) / r12 - | < 0 : Short-range operator erfc(omega r12) /r12 - ''' - if omega is None: - self._env[PTR_RANGE_OMEGA] = 0 + if self.spin is None: + self.spin = self.nelectron % 2 else: - self._env[PTR_RANGE_OMEGA] = omega - set_range_coulomb_ = set_range_coulomb # for backward compatibility + # Access self.nelec in which the code checks whether the spin and + # number of electrons are consistent. + self.nelec - @property - def omega(self): - return self._env[PTR_RANGE_OMEGA] - omega = omega.setter(set_range_coulomb) + if self.magmom is None or len(self.magmom) != self.natm: + self.magmom = [0,] * self.natm + if self.spin == 0 and abs(numpy.sum(self.magmom) - self.spin) > 1e-6: + #don't check for unrestricted calcs. + raise ValueError("mol.magmom is set incorrectly.") - def with_range_coulomb(self, omega): - '''Return a temporary mol context which sets the required parameter - omega for range-separated Coulomb operator. - If omega = None, return the context for regular Coulomb integrals. - See also :func:`mol.set_range_coulomb` + if self.symmetry: + self._build_symmetry() - Examples: + if dump_input and not self._built and self.verbose > logger.NOTE: + self.dump_input() - >>> with mol.with_range_coulomb(omega=1.5): - ... mol.intor('int2e') - ''' - omega0 = self._env[PTR_RANGE_OMEGA].copy() - return self._TemporaryMoleContext(self.set_range_coulomb, (omega,), (omega0,)) + if self.verbose >= logger.DEBUG3: + logger.debug3(self, 'arg.atm = %s', self._atm) + logger.debug3(self, 'arg.bas = %s', self._bas) + logger.debug3(self, 'arg.env = %s', self._env) + logger.debug3(self, 'ecpbas = %s', self._ecpbas) - def with_long_range_coulomb(self, omega): - '''Return a temporary mol context for long-range part of - range-separated Coulomb operator. - ''' - return self.with_range_coulomb(abs(omega)) + self._built = True + return self + kernel = build - def with_short_range_coulomb(self, omega): - '''Return a temporary mol context for short-range part of - range-separated Coulomb operator. + def _build_symmetry(self, *args, **kwargs): ''' - return self.with_range_coulomb(-abs(omega)) - - def set_f12_zeta(self, zeta): - '''Set zeta for YP exp(-zeta r12)/r12 or STG exp(-zeta r12) type integrals + Update symmetry related attributes: topgroup, groupname, _symm_orig, + _symm_axes, irrep_id, irrep_name, symm_orb ''' - self._env[PTR_F12_ZETA] = zeta - - def set_nuc_mod(self, atm_id, zeta): - '''Change the nuclear charge distribution of the given atom ID. The charge - distribution is defined as: rho(r) = nuc_charge * Norm * exp(-zeta * r^2). - This function can **only** be called after .build() method is executed. + from pyscf import symm - Examples: + # TODO: Consider ECP info in point group symmetry initialization + self.topgroup, orig, axes = symm.detect_symm(self._atom, self._basis) - >>> for ia in range(mol.natm): - ... zeta = gto.filatov_nuc_mod(mol.atom_charge(ia)) - ... mol.set_nuc_mod(ia, zeta) - ''' - ptr = self._atm[atm_id,PTR_ZETA] - self._env[ptr] = zeta - if zeta == 0: - self._atm[atm_id,NUC_MOD_OF] = NUC_POINT + if isinstance(self.symmetry, str): + self.symmetry = str(symm.std_symb(self.symmetry)) + groupname = None + if abs(axes - np.eye(3)).max() < symm.TOLERANCE: + if symm.check_symm(self.symmetry, self._atom, self._basis): + # Try to use original axes (issue #1209) + groupname = self.symmetry + axes = np.eye(3) + else: + logger.warn(self, 'Unable to to identify input symmetry using original axes.\n' + 'Different symmetry axes will be used.') + if groupname is None: + try: + groupname, axes = symm.as_subgroup(self.topgroup, axes, + self.symmetry) + except PointGroupSymmetryError as e: + raise PointGroupSymmetryError( + 'Unable to identify input symmetry %s. Try symmetry="%s"' % + (self.symmetry, self.topgroup)) from e else: - self._atm[atm_id,NUC_MOD_OF] = NUC_GAUSS - return self - set_nuc_mod_ = set_nuc_mod # for backward compatibility + groupname, axes = symm.as_subgroup(self.topgroup, axes, + self.symmetry_subgroup) + self._symm_orig = orig + self._symm_axes = axes - def set_rinv_zeta(self, zeta): - '''Assume the charge distribution on the "rinv_origin". zeta is the parameter - to control the charge distribution: rho(r) = Norm * exp(-zeta * r^2). - **Be careful** when call this function. It affects the behavior of - int1e_rinv_* functions. Make sure to set it back to 0 after using it! - ''' - self._env[PTR_RINV_ZETA] = zeta + if self.cart and groupname in ('Dooh', 'Coov', 'SO3'): + if groupname == 'Coov': + groupname, lgroup = 'C2v', groupname + else: + groupname, lgroup = 'D2h', groupname + logger.warn(self, 'This version does not support symmetry %s ' + 'for cartesian GTO basis. Its subgroup %s is used', + lgroup, groupname) + self.groupname = groupname + + self.symm_orb, self.irrep_id = \ + symm.symm_adapted_basis(self, groupname, orig, axes) + self.irrep_name = [symm.irrep_id2name(groupname, ir) + for ir in self.irrep_id] return self - set_rinv_zeta_ = set_rinv_zeta # for backward compatibility - def with_rinv_zeta(self, zeta): - '''Return a temporary mol context which has the rquired Gaussian charge - distribution placed at "rinv_origin": rho(r) = Norm * exp(-zeta * r^2). - See also :func:`mol.set_rinv_zeta` + @lib.with_doc(format_atom.__doc__) + def format_atom(self, atom, origin=0, axes=None, unit='Ang'): + return format_atom(atom, origin, axes, unit) - Examples: + @lib.with_doc(format_basis.__doc__) + def format_basis(self, basis_tab): + return format_basis(basis_tab) - >>> with mol.with_rinv_zeta(zeta=1.5), mol.with_rinv_origin((1.,0,0)): - ... mol.intor('int1e_rinv') - ''' - zeta0 = self._env[PTR_RINV_ZETA].copy() - return self._TemporaryMoleContext(self.set_rinv_zeta, (zeta,), (zeta0,)) + @lib.with_doc(format_pseudo.__doc__) + def format_pseudo(self, pseudo_tab): + return format_pseudo(pseudo_tab) - def with_rinv_at_nucleus(self, atm_id): - '''Return a temporary mol context in which the rinv operator (1/r) is - treated like the Coulomb potential of a Gaussian charge distribution - rho(r) = Norm * exp(-zeta * r^2) at the place of the input atm_id. + @lib.with_doc(format_ecp.__doc__) + def format_ecp(self, ecp_tab): + return format_ecp(ecp_tab) - Examples: + @lib.with_doc(expand_etb.__doc__) + def expand_etb(self, l, n, alpha, beta): + return expand_etb(l, n, alpha, beta) - >>> with mol.with_rinv_at_nucleus(3): - ... mol.intor('int1e_rinv') - ''' - zeta = self._env[self._atm[atm_id,PTR_ZETA]] - rinv = self.atom_coord(atm_id) - if zeta == 0: + @lib.with_doc(expand_etbs.__doc__) + def expand_etbs(self, etbs): + return expand_etbs(etbs) + etbs = expand_etbs + + @lib.with_doc(make_env.__doc__) + def make_env(self, atoms, basis, pre_env=[], nucmod={}, nucprop=None): + if nucprop is None: + nucprop = self.nucprop + return make_env(atoms, basis, pre_env, nucmod, nucprop) + + @lib.with_doc(make_atm_env.__doc__) + def make_atm_env(self, atom, ptr=0, nucmod=NUC_POINT, nucprop=None): + if nucprop is None: + nucprop = self.nucprop.get(atom[0], {}) + return make_atm_env(atom, ptr, nucmod, nucprop) + + @lib.with_doc(make_bas_env.__doc__) + def make_bas_env(self, basis_add, atom_id=0, ptr=0): + return make_bas_env(basis_add, atom_id, ptr) + + @lib.with_doc(make_ecp_env.__doc__) + def make_ecp_env(self, _atm, _ecp, pre_env=[]): + if _ecp: + _atm, _ecpbas, _env = make_ecp_env(self, _atm, _ecp, pre_env) + else: + _atm, _ecpbas, _env = _atm, numpy.zeros((0,BAS_SLOTS)), pre_env + return _atm, _ecpbas, _env + + tot_electrons = tot_electrons + + @lib.with_doc(gto_norm.__doc__) + def gto_norm(self, l, expnt): + return gto_norm(l, expnt) + + def dump_input(self): + import __main__ + import pyscf + if hasattr(__main__, '__file__'): + try: + filename = os.path.abspath(__main__.__file__) + finput = open(filename, 'r') + self.stdout.write('#INFO: **** input file is %s ****\n' % filename) + self.stdout.write(finput.read()) + self.stdout.write('#INFO: ******************** input file end ********************\n') + self.stdout.write('\n') + self.stdout.write('\n') + finput.close() + except IOError: + logger.warn(self, 'input file does not exist') + + self.stdout.write('System: %s Threads %s\n' % + (str(platform.uname()), lib.num_threads())) + self.stdout.write('Python %s\n' % sys.version) + self.stdout.write('numpy %s scipy %s\n' % + (numpy.__version__, scipy.__version__)) + self.stdout.write('Date: %s\n' % time.ctime()) + self.stdout.write('PySCF version %s\n' % pyscf.__version__) + info = lib.repo_info(os.path.join(__file__, '..', '..')) + self.stdout.write('PySCF path %s\n' % info['path']) + if 'git' in info: + self.stdout.write(info['git'] + '\n') + + self.stdout.write('\n') + for key in os.environ: + if 'PYSCF' in key: + self.stdout.write('[ENV] %s %s\n' % (key, os.environ[key])) + if self.verbose >= logger.DEBUG2: + for key in dir(__config__): + if key[:2] != '__': + self.stdout.write('[CONFIG] %s = %s\n' % + (key, getattr(__config__, key))) + else: + conf_file = getattr(__config__, 'conf_file', None) + self.stdout.write('[CONFIG] conf_file %s\n' % conf_file) + + self.stdout.write('[INPUT] verbose = %d\n' % self.verbose) + if self.verbose >= logger.DEBUG: + self.stdout.write('[INPUT] max_memory = %s \n' % self.max_memory) + self.stdout.write('[INPUT] num. atoms = %d\n' % self.natm) + self.stdout.write('[INPUT] num. electrons = %d\n' % self.nelectron) + self.stdout.write('[INPUT] charge = %d\n' % self.charge) + self.stdout.write('[INPUT] spin (= nelec alpha-beta = 2S) = %d\n' % self.spin) + self.stdout.write('[INPUT] symmetry %s subgroup %s\n' % + (self.symmetry, self.symmetry_subgroup)) + self.stdout.write('[INPUT] Mole.unit = %s\n' % self.unit) + if self.cart: + self.stdout.write('[INPUT] Cartesian GTO integrals (6d 10f)\n') + + self.stdout.write('[INPUT] Symbol X Y Z unit' + ' X Y Z unit Magmom\n') + for ia,atom in enumerate(self._atom): + coorda = tuple([x * param.BOHR for x in atom[1]]) + coordb = tuple([x for x in atom[1]]) + magmom = self.magmom[ia] + self.stdout.write('[INPUT]%3d %-4s %16.12f %16.12f %16.12f AA ' + '%16.12f %16.12f %16.12f Bohr %4.1f\n' + % ((ia+1, _symbol(atom[0])) + coorda + coordb + (magmom,))) + if self.nucmod: + if isinstance(self.nucmod, (int, str, types.FunctionType)): + nucatms = [_symbol(atom[0]) for atom in self._atom] + else: + nucatms = self.nucmod.keys() + self.stdout.write('[INPUT] Gaussian nuclear model for atoms %s\n' % + nucatms) + + if self.nucprop: + self.stdout.write('[INPUT] nucprop %s\n' % self.nucprop) + + if self.verbose >= logger.DEBUG: + self.stdout.write('[INPUT] ---------------- BASIS SET ---------------- \n') + self.stdout.write('[INPUT] l, kappa, [nprim/nctr], ' + 'expnt, c_1 c_2 ...\n') + for atom, basis_set in self._basis.items(): + self.stdout.write('[INPUT] %s\n' % atom) + for b in basis_set: + if isinstance(b[1], int): + kappa = b[1] + b_coeff = b[2:] + else: + kappa = 0 + b_coeff = b[1:] + nprim = len(b_coeff) + nctr = len(b_coeff[0])-1 + if nprim < nctr: + logger.warn(self, 'num. primitives smaller than num. contracted basis') + self.stdout.write('[INPUT] %d %2d [%-5d/%-4d] ' + % (b[0], kappa, nprim, nctr)) + for k, x in enumerate(b_coeff): + if k == 0: + self.stdout.write('%-15.12g ' % x[0]) + else: + self.stdout.write(' '*32+'%-15.12g ' % x[0]) + for c in x[1:]: + self.stdout.write(' %4.12g' % c) + self.stdout.write('\n') + + if self.verbose >= logger.INFO: + self.stdout.write('\n') + logger.info(self, 'nuclear repulsion = %.15g', self.energy_nuc()) + if self.symmetry: + if self.topgroup == self.groupname: + logger.info(self, 'point group symmetry = %s', self.topgroup) + else: + logger.info(self, 'point group symmetry = %s, use subgroup %s', + self.topgroup, self.groupname) + logger.info(self, "symmetry origin: %s", self._symm_orig) + logger.info(self, "symmetry axis x: %s", self._symm_axes[0]) + logger.info(self, "symmetry axis y: %s", self._symm_axes[1]) + logger.info(self, "symmetry axis z: %s", self._symm_axes[2]) + for ir in range(self.symm_orb.__len__()): + logger.info(self, 'num. orbitals of irrep %s = %d', + self.irrep_name[ir], self.symm_orb[ir].shape[1]) + logger.info(self, 'number of shells = %d', self.nbas) + logger.info(self, 'number of NR pGTOs = %d', self.npgto_nr()) + logger.info(self, 'number of NR cGTOs = %d', self.nao_nr()) + logger.info(self, 'basis = %s', self.basis) + logger.info(self, 'ecp = %s', self.ecp) + if self.verbose >= logger.DEBUG2: + for i in range(len(self._bas)): + exps = self.bas_exp(i) + logger.debug1(self, 'bas %d, expnt(s) = %s', i, str(exps)) + + logger.info(self, 'CPU time: %12.2f', logger.process_clock()) + return self + + def set_common_origin(self, coord): + '''Update common origin for integrals of dipole, rxp etc. + **Note** the unit of the coordinates needs to be Bohr + + Examples: + + >>> mol.set_common_origin(0) + >>> mol.set_common_origin((1,0,0)) + ''' + self._env[PTR_COMMON_ORIG:PTR_COMMON_ORIG+3] = coord + return self + set_common_orig = set_common_origin + set_common_orig_ = set_common_orig # for backward compatibility + set_common_origin_ = set_common_orig # for backward compatibility + + def with_common_origin(self, coord): + '''Return a temporary mol context which has the rquired common origin. + The required common origin has no effects out of the temporary context. + See also :func:`mol.set_common_origin` + + Examples: + + >>> with mol.with_common_origin((1,0,0)): + ... mol.intor('int1e_r', comp=3) + ''' + coord0 = self._env[PTR_COMMON_ORIG:PTR_COMMON_ORIG+3].copy() + return self._TemporaryMoleContext(self.set_common_origin, (coord,), (coord0,)) + with_common_orig = with_common_origin + + def set_rinv_origin(self, coord): + r'''Update origin for operator :math:`\frac{1}{|r-R_O|}`. + **Note** the unit is Bohr + + Examples: + + >>> mol.set_rinv_origin(0) + >>> mol.set_rinv_origin((0,1,0)) + ''' + self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3] = coord[:3] + return self + set_rinv_orig = set_rinv_origin + set_rinv_orig_ = set_rinv_orig # for backward compatibility + set_rinv_origin_ = set_rinv_orig # for backward compatibility + + def with_rinv_origin(self, coord): + '''Return a temporary mol context which has the rquired origin of 1/r + operator. The required origin has no effects out of the temporary + context. See also :func:`mol.set_rinv_origin` + + Examples: + + >>> with mol.with_rinv_origin((1,0,0)): + ... mol.intor('int1e_rinv') + ''' + coord0 = self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3].copy() + return self._TemporaryMoleContext(self.set_rinv_origin, (coord,), (coord0,)) + with_rinv_orig = with_rinv_origin + + def set_range_coulomb(self, omega): + '''Switch on range-separated Coulomb operator for **all** 2e integrals + + Args: + omega : double + + | = 0 : Regular electron repulsion integral + | > 0 : Long-range operator erf(omega r12) / r12 + | < 0 : Short-range operator erfc(omega r12) /r12 + ''' + if omega is None: + self._env[PTR_RANGE_OMEGA] = 0 + else: + self._env[PTR_RANGE_OMEGA] = omega + set_range_coulomb_ = set_range_coulomb # for backward compatibility + + @property + def omega(self): + return self._env[PTR_RANGE_OMEGA] + omega = omega.setter(set_range_coulomb) + + def with_range_coulomb(self, omega): + '''Return a temporary mol context which sets the required parameter + omega for range-separated Coulomb operator. + If omega = None, return the context for regular Coulomb integrals. + See also :func:`mol.set_range_coulomb` + + Examples: + + >>> with mol.with_range_coulomb(omega=1.5): + ... mol.intor('int2e') + ''' + omega0 = self._env[PTR_RANGE_OMEGA].copy() + return self._TemporaryMoleContext(self.set_range_coulomb, (omega,), (omega0,)) + + def with_long_range_coulomb(self, omega): + '''Return a temporary mol context for long-range part of + range-separated Coulomb operator. + ''' + return self.with_range_coulomb(abs(omega)) + + def with_short_range_coulomb(self, omega): + '''Return a temporary mol context for short-range part of + range-separated Coulomb operator. + ''' + return self.with_range_coulomb(-abs(omega)) + + def set_f12_zeta(self, zeta): + '''Set zeta for YP exp(-zeta r12)/r12 or STG exp(-zeta r12) type integrals + ''' + self._env[PTR_F12_ZETA] = zeta + + def set_nuc_mod(self, atm_id, zeta): + '''Change the nuclear charge distribution of the given atom ID. The charge + distribution is defined as: rho(r) = nuc_charge * Norm * exp(-zeta * r^2). + This function can **only** be called after .build() method is executed. + + Examples: + + >>> for ia in range(mol.natm): + ... zeta = gto.filatov_nuc_mod(mol.atom_charge(ia)) + ... mol.set_nuc_mod(ia, zeta) + ''' + ptr = self._atm[atm_id,PTR_ZETA] + self._env[ptr] = zeta + if zeta == 0: + self._atm[atm_id,NUC_MOD_OF] = NUC_POINT + else: + self._atm[atm_id,NUC_MOD_OF] = NUC_GAUSS + return self + set_nuc_mod_ = set_nuc_mod # for backward compatibility + + def set_rinv_zeta(self, zeta): + '''Assume the charge distribution on the "rinv_origin". zeta is the parameter + to control the charge distribution: rho(r) = Norm * exp(-zeta * r^2). + **Be careful** when call this function. It affects the behavior of + int1e_rinv_* functions. Make sure to set it back to 0 after using it! + ''' + self._env[PTR_RINV_ZETA] = zeta + return self + set_rinv_zeta_ = set_rinv_zeta # for backward compatibility + + def with_rinv_zeta(self, zeta): + '''Return a temporary mol context which has the rquired Gaussian charge + distribution placed at "rinv_origin": rho(r) = Norm * exp(-zeta * r^2). + See also :func:`mol.set_rinv_zeta` + + Examples: + + >>> with mol.with_rinv_zeta(zeta=1.5), mol.with_rinv_origin((1.,0,0)): + ... mol.intor('int1e_rinv') + ''' + zeta0 = self._env[PTR_RINV_ZETA].copy() + return self._TemporaryMoleContext(self.set_rinv_zeta, (zeta,), (zeta0,)) + + def with_rinv_at_nucleus(self, atm_id): + '''Return a temporary mol context in which the rinv operator (1/r) is + treated like the Coulomb potential of a Gaussian charge distribution + rho(r) = Norm * exp(-zeta * r^2) at the place of the input atm_id. + + Examples: + + >>> with mol.with_rinv_at_nucleus(3): + ... mol.intor('int1e_rinv') + ''' + zeta = self._env[self._atm[atm_id,PTR_ZETA]] + rinv = self.atom_coord(atm_id) + if zeta == 0: self._env[AS_RINV_ORIG_ATOM] = atm_id # required by ecp gradients return self.with_rinv_origin(rinv) else: @@ -2577,6 +3074,14 @@ def set_geom_(self, atoms_or_coords, unit=None, symmetry=None, ia+1, mol.atom_symbol(ia), *coords) return mol + def update(self, chkfile): + return self.update_from_chk(chkfile) + def update_from_chk(self, chkfile): + with h5py.File(chkfile, 'r') as fh5: + mol = loads(fh5['mol'][()]) + self.__dict__.update(mol.__dict__) + return self + def has_ecp(self): '''Whether pseudo potential is used in the system.''' return len(self._ecpbas) > 0 or self._pseudo @@ -2659,824 +3164,479 @@ def atom_charges(self): def atom_nelec_core(self, atm_id): '''Number of core electrons for pseudo potential. ''' - return charge(self.atom_symbol(atm_id)) - self.atom_charge(atm_id) - - def atom_coord(self, atm_id, unit='Bohr'): - r'''Coordinates (ndarray) of the given atom id - - Args: - atm_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') - >>> mol.atom_coord(1) - [ 0. 0. 2.07869874] - ''' - ptr = self._atm[atm_id,PTR_COORD] - if not is_au(unit): - return self._env[ptr:ptr+3] * param.BOHR - else: - return self._env[ptr:ptr+3].copy() - - def atom_coords(self, unit='Bohr'): - '''np.asarray([mol.atom_coords(i) for i in range(mol.natm)])''' - ptr = self._atm[:,PTR_COORD] - c = self._env[numpy.vstack((ptr,ptr+1,ptr+2)).T].copy() - if not is_au(unit): - c *= param.BOHR - return c - - atom_mass_list = atom_mass_list - - def atom_nshells(self, atm_id): - r'''Number of basis/shells of the given atom - - Args: - atm_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') - >>> mol.atom_nshells(1) - 5 - ''' - return (self._bas[:,ATOM_OF] == atm_id).sum() - - def atom_shell_ids(self, atm_id): - r'''A list of the shell-ids of the given atom - - Args: - atm_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.atom_shell_ids(1) - [3, 4, 5, 6, 7] - ''' - return numpy.where(self._bas[:,ATOM_OF] == atm_id)[0] - - def bas_coord(self, bas_id): - r'''Coordinates (ndarray) associated with the given basis id - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') - >>> mol.bas_coord(1) - [ 0. 0. 2.07869874] - ''' - atm_id = self.bas_atom(bas_id) - ptr = self._atm[atm_id,PTR_COORD] - return self._env[ptr:ptr+3].copy() - - def bas_atom(self, bas_id): - r'''The atom (0-based id) that the given basis sits on - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_atom(7) - 1 - ''' - return self._bas[bas_id,ATOM_OF].copy() - - def bas_angular(self, bas_id): - r'''The angular momentum associated with the given basis - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_atom(7) - 2 - ''' - return self._bas[bas_id,ANG_OF].copy() - - def bas_nctr(self, bas_id): - r'''The number of contracted GTOs for the given shell - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_atom(3) - 3 - ''' - return self._bas[bas_id,NCTR_OF].copy() - - def bas_nprim(self, bas_id): - r'''The number of primitive GTOs for the given shell - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_atom(3) - 11 - ''' - return self._bas[bas_id,NPRIM_OF].copy() - - def bas_kappa(self, bas_id): - r'''Kappa (if l < j, -l-1, else l) of the given shell - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_kappa(3) - 0 - ''' - return self._bas[bas_id,KAPPA_OF].copy() - - def bas_exp(self, bas_id): - r'''exponents (ndarray) of the given shell - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_kappa(0) - [ 13.01 1.962 0.4446] - ''' - nprim = self.bas_nprim(bas_id) - ptr = self._bas[bas_id,PTR_EXP] - return self._env[ptr:ptr+nprim].copy() - - def bas_exps(self): - '''exponents of all basis - return [mol.bas_exp(i) for i in range(self.nbas)] - ''' - nprims = self._bas[:,NPRIM_OF] - pexps = self._bas[:,PTR_EXP] - exps = [self._env[i0:i1] for i0, i1 in zip(pexps, pexps + nprims)] - return exps - - def _libcint_ctr_coeff(self, bas_id): - nprim = self.bas_nprim(bas_id) - nctr = self.bas_nctr(bas_id) - ptr = self._bas[bas_id,PTR_COEFF] - return self._env[ptr:ptr+nprim*nctr].reshape(nctr,nprim).T - - def bas_ctr_coeff(self, bas_id): - r'''Contract coefficients (ndarray) of the given shell - - Args: - bas_id : int - 0-based - - Examples: - - >>> mol.M(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') - >>> mol.bas_ctr_coeff(0) - [[ 10.03400444] - [ 4.1188704 ] - [ 1.53971186]] - ''' - l = self.bas_angular(bas_id) - es = self.bas_exp(bas_id) - cs = self._libcint_ctr_coeff(bas_id) - cs = numpy.einsum('pi,p->pi', cs, 1/gto_norm(l, es)) - return cs - - def bas_len_spinor(self, bas_id): - '''The number of spinor associated with given basis - If kappa is 0, return 4l+2 - ''' - l = self.bas_angular(bas_id) - k = self.bas_kappa(bas_id) - return len_spinor(l, k) - - def bas_len_cart(self, bas_id): - '''The number of Cartesian function associated with given basis - ''' - return len_cart(self._bas[bas_id,ANG_OF]) - - npgto_nr = npgto_nr + return charge(self.atom_symbol(atm_id)) - self.atom_charge(atm_id) - nao_nr = nao_nr - nao_2c = nao_2c - nao_cart = nao_cart + def atom_coord(self, atm_id, unit='Bohr'): + r'''Coordinates (ndarray) of the given atom id - nao_nr_range = nao_nr_range - nao_2c_range = nao_2c_range + Args: + atm_id : int + 0-based - ao_loc_nr = ao_loc_nr - ao_loc_2c = ao_loc_2c + Examples: - @property - def nao(self): - if self._nao is None: - return self.nao_nr() + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') + >>> mol.atom_coord(1) + [ 0. 0. 2.07869874] + ''' + ptr = self._atm[atm_id,PTR_COORD] + if not is_au(unit): + return self._env[ptr:ptr+3] * param.BOHR else: - return self._nao - @nao.setter - def nao(self, x): - self._nao = x + return self._env[ptr:ptr+3].copy() - ao_loc = property(ao_loc_nr) + def atom_coords(self, unit='Bohr'): + '''np.asarray([mol.atom_coords(i) for i in range(mol.natm)])''' + ptr = self._atm[:,PTR_COORD] + c = self._env[numpy.vstack((ptr,ptr+1,ptr+2)).T].copy() + if not is_au(unit): + c *= param.BOHR + return c - tmap = time_reversal_map = time_reversal_map + atom_mass_list = atom_mass_list - def intor(self, intor, comp=None, hermi=0, aosym='s1', out=None, - shls_slice=None, grids=None): - '''Integral generator. + def atom_nshells(self, atm_id): + r'''Number of basis/shells of the given atom Args: - intor : str - Name of the 1e or 2e AO integrals. Ref to :func:`getints` for the - complete list of available 1-electron integral names + atm_id : int + 0-based - Kwargs: - comp : int - Components of the integrals, e.g. int1e_ipovlp_sph has 3 components. - hermi : int - Symmetry of the integrals + Examples: - | 0 : no symmetry assumed (default) - | 1 : hermitian - | 2 : anti-hermitian + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') + >>> mol.atom_nshells(1) + 5 + ''' + return (self._bas[:,ATOM_OF] == atm_id).sum() - grids : ndarray - Coordinates of grids for the int1e_grids integrals + def atom_shell_ids(self, atm_id): + r'''A list of the shell-ids of the given atom - Returns: - ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp + Args: + atm_id : int + 0-based Examples: - >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') - >>> mol.intor('int1e_ipnuc_sph', comp=3) # - [[[ 0. 0. ] - [ 0. 0. ]] - [[ 0. 0. ] - [ 0. 0. ]] - [[ 0.10289944 0.48176097] - [-0.48176097 -0.10289944]]] - >>> mol.intor('int1e_nuc_spinor') - [[-1.69771092+0.j 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j] - [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j -0.67146312+0.j] - [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] - [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.atom_shell_ids(1) + [3, 4, 5, 6, 7] ''' - if not self._built: - logger.warn(self, 'Warning: intor envs of %s not initialized.', self) - # FIXME: Whether to check _built and call build? ._bas and .basis - # may not be consistent. calling .build() may leads to wrong intor env. - #self.build(False, False) - intor = self._add_suffix(intor) - bas = self._bas - env = self._env - if 'ECP' in intor: - assert (self._ecp is not None) - bas = numpy.vstack((self._bas, self._ecpbas)) - env[AS_ECPBAS_OFFSET] = len(self._bas) - env[AS_NECPBAS] = len(self._ecpbas) - if shls_slice is None: - shls_slice = (0, self.nbas, 0, self.nbas) - elif '_grids' in intor: - assert grids is not None - env = numpy.append(env, grids.ravel()) - env[NGRIDS] = grids.shape[0] - env[PTR_GRIDS] = env.size - grids.size - return moleintor.getints(intor, self._atm, bas, env, - shls_slice, comp, hermi, aosym, out=out) - - def _add_suffix(self, intor, cart=None): - if not (intor[:4] == 'cint' or - intor.endswith(('_sph', '_cart', '_spinor', '_ssc'))): - if cart is None: - cart = self.cart - if cart: - intor = intor + '_cart' - else: - intor = intor + '_sph' - return intor + return numpy.where(self._bas[:,ATOM_OF] == atm_id)[0] - def intor_symmetric(self, intor, comp=None, grids=None): - '''One-electron integral generator. The integrals are assumed to be hermitian + def bas_coord(self, bas_id): + r'''Coordinates (ndarray) associated with the given basis id Args: - intor : str - Name of the 1-electron integral. Ref to :func:`getints` for the - complete list of available 1-electron integral names - - Kwargs: - comp : int - Components of the integrals, e.g. int1e_ipovlp_sph has 3 components. - grids : ndarray - Coordinates of grids for the int1e_grids integrals - - Returns: - ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp + bas_id : int + 0-based Examples: - >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') - >>> mol.intor_symmetric('int1e_nuc_spinor') - [[-1.69771092+0.j 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j] - [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j -0.67146312+0.j] - [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] - [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1') + >>> mol.bas_coord(1) + [ 0. 0. 2.07869874] ''' - return self.intor(intor, comp, 1, aosym='s4', grids=grids) + atm_id = self.bas_atom(bas_id) + ptr = self._atm[atm_id,PTR_COORD] + return self._env[ptr:ptr+3].copy() - def intor_asymmetric(self, intor, comp=None, grids=None): - '''One-electron integral generator. The integrals are assumed to be anti-hermitian + def bas_atom(self, bas_id): + r'''The atom (0-based id) that the given basis sits on Args: - intor : str - Name of the 1-electron integral. Ref to :func:`getints` for the - complete list of available 1-electron integral names + bas_id : int + 0-based - Kwargs: - comp : int - Components of the integrals, e.g. int1e_ipovlp has 3 components. - grids : ndarray - Coordinates of grids for the int1e_grids integrals + Examples: - Returns: - ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_atom(7) + 1 + ''' + return self._bas[bas_id,ATOM_OF].copy() + + def bas_angular(self, bas_id): + r'''The angular momentum associated with the given basis + + Args: + bas_id : int + 0-based Examples: - >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') - >>> mol.intor_asymmetric('int1e_nuc_spinor') - [[-1.69771092+0.j 0.00000000+0.j 0.67146312+0.j 0.00000000+0.j] - [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j 0.67146312+0.j] - [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] - [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_atom(7) + 2 ''' - return self.intor(intor, comp, 2, aosym='a4', grids=grids) + return self._bas[bas_id,ANG_OF].copy() - @lib.with_doc(moleintor.getints_by_shell.__doc__) - def intor_by_shell(self, intor, shells, comp=None, grids=None): - intor = self._add_suffix(intor) - if 'ECP' in intor: - assert (self._ecp is not None) - bas = numpy.vstack((self._bas, self._ecpbas)) - self._env[AS_ECPBAS_OFFSET] = len(self._bas) - self._env[AS_NECPBAS] = len(self._ecpbas) - else: - bas = self._bas - return moleintor.getints_by_shell(intor, shells, self._atm, bas, - self._env, comp) + def bas_nctr(self, bas_id): + r'''The number of contracted GTOs for the given shell - eval_ao = eval_gto = eval_gto + Args: + bas_id : int + 0-based - def get_ao_indices(self, bas_list, ao_loc=None): - ''' - Generate (dis-continued) AO indices for basis specified in bas_list - ''' - if ao_loc is None: - ao_loc = self.ao_loc - return lib.locs_to_indices(ao_loc, bas_list) + Examples: - sph_labels = spheric_labels = sph_labels - cart_labels = cart_labels - ao_labels = ao_labels + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_atom(3) + 3 + ''' + return self._bas[bas_id,NCTR_OF].copy() - spinor_labels = spinor_labels + def bas_nprim(self, bas_id): + r'''The number of primitive GTOs for the given shell - search_ao_label = search_ao_label + Args: + bas_id : int + 0-based - def search_shell_id(self, atm_id, l): - return search_shell_id(self, atm_id, l) + Examples: - search_ao_nr = search_ao_nr - search_ao_r = search_ao_r + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_atom(3) + 11 + ''' + return self._bas[bas_id,NPRIM_OF].copy() - aoslice_by_atom = aoslice_nr_by_atom = offset_ao_by_atom = offset_nr_by_atom = aoslice_by_atom - aoslice_2c_by_atom = offset_2c_by_atom = offset_2c_by_atom + def bas_kappa(self, bas_id): + r'''Kappa (if l < j, -l-1, else l) of the given shell - condense_to_shell = condense_to_shell - get_overlap_cond = get_overlap_cond + Args: + bas_id : int + 0-based - to_uncontracted_cartesian_basis = to_uncontracted_cartesian_basis - decontract_basis = decontract_basis + Examples: - ao_rotation_matrix = ao_rotation_matrix + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_kappa(3) + 0 + ''' + return self._bas[bas_id,KAPPA_OF].copy() - def cart2sph_coeff(self, normalized='sp'): - '''Transformation matrix that transforms Cartesian GTOs to spherical - GTOs for all basis functions + def bas_exp(self, bas_id): + r'''exponents (ndarray) of the given shell - Kwargs: - normalized : string or boolean - How the Cartesian GTOs are normalized. Except s and p functions, - Cartesian GTOs do not have the universal normalization coefficients - for the different components of the same shell. The value of this - argument can be one of 'sp', 'all', None. 'sp' means the Cartesian s - and p basis are normalized. 'all' means all Cartesian functions are - normalized. None means none of the Cartesian functions are normalized. - The default value 'sp' is the convention used by libcint library. + Args: + bas_id : int + 0-based Examples: - >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') - >>> c = mol.cart2sph_coeff() - >>> s0 = mol.intor('int1e_ovlp_sph') - >>> s1 = c.T.dot(mol.intor('int1e_ovlp_cart')).dot(c) - >>> print(abs(s1-s0).sum()) - >>> 4.58676826646e-15 + >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_kappa(0) + [ 13.01 1.962 0.4446] ''' - c2s_l = [cart2sph(l, normalized=normalized) for l in range(12)] - c2s = [] - for ib in range(self.nbas): - l = self.bas_angular(ib) - for n in range(self.bas_nctr(ib)): - c2s.append(c2s_l[l]) - return scipy.linalg.block_diag(*c2s) + nprim = self.bas_nprim(bas_id) + ptr = self._bas[bas_id,PTR_EXP] + return self._env[ptr:ptr+nprim].copy() - def sph2spinor_coeff(self): - '''Transformation matrix that transforms real-spherical GTOs to spinor - GTOs for all basis functions + def bas_exps(self): + '''exponents of all basis + return [mol.bas_exp(i) for i in range(self.nbas)] + ''' + nprims = self._bas[:,NPRIM_OF] + pexps = self._bas[:,PTR_EXP] + exps = [self._env[i0:i1] for i0, i1 in zip(pexps, pexps + nprims)] + return exps - Examples: + def _libcint_ctr_coeff(self, bas_id): + nprim = self.bas_nprim(bas_id) + nctr = self.bas_nctr(bas_id) + ptr = self._bas[bas_id,PTR_COEFF] + return self._env[ptr:ptr+nprim*nctr].reshape(nctr,nprim).T - >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') - >>> ca, cb = mol.sph2spinor_coeff() - >>> 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 - ''' - from pyscf.symm import sph - return sph.sph2spinor_coeff(self) + def bas_ctr_coeff(self, bas_id): + r'''Contract coefficients (ndarray) of the given shell + Args: + bas_id : int + 0-based -class GlobalOptionMixin: + Examples: - verbose = getattr(__config__, 'VERBOSE', logger.NOTE) + >>> mol.M(atom='H 0 0 0; Cl 0 0 1.1', basis='cc-pvdz') + >>> mol.bas_ctr_coeff(0) + [[ 10.03400444] + [ 4.1188704 ] + [ 1.53971186]] + ''' + l = self.bas_angular(bas_id) + es = self.bas_exp(bas_id) + cs = self._libcint_ctr_coeff(bas_id) + cs = numpy.einsum('pi,p->pi', cs, 1/gto_norm(l, es)) + return cs - # the unit (angstrom/bohr) of the coordinates defined by the input self.atom - unit = getattr(__config__, 'UNIT', 'angstrom') + def bas_len_spinor(self, bas_id): + '''The number of spinor associated with given basis + If kappa is 0, return 4l+2 + ''' + l = self.bas_angular(bas_id) + k = self.bas_kappa(bas_id) + return len_spinor(l, k) - # Whether to hold everything in memory - incore_anyway = getattr(__config__, 'INCORE_ANYWAY', False) + def bas_len_cart(self, bas_id): + '''The number of Cartesian function associated with given basis + ''' + return len_cart(self._bas[bas_id,ANG_OF]) - output = None - max_memory = param.MAX_MEMORY + npgto_nr = npgto_nr - _keys = { - 'verbose', 'unit', 'incore_anyway', 'output', 'max_memory', - } + nao_nr = nao_nr + nao_2c = nao_2c + nao_cart = nao_cart - def __init__(self): - self.stdout = sys.stdout - self._built = False + nao_nr_range = nao_nr_range + nao_2c_range = nao_2c_range - # Some methods modify ._env. These method are executed in the context - # _TemporaryMoleContext which is protected by the _ctx_lock. - self._ctx_lock = None + ao_loc_nr = ao_loc_nr + ao_loc_2c = ao_loc_2c - copy = copy + @property + def nao(self): + if self._nao is None: + return self.nao_nr() + else: + return self._nao + @nao.setter + def nao(self, x): + self._nao = x - pack = pack + ao_loc = property(ao_loc_nr) - @classmethod - @lib.with_doc(unpack.__doc__) - def unpack(cls, moldic): - return unpack(moldic) + tmap = time_reversal_map = time_reversal_map - @lib.with_doc(unpack.__doc__) - def unpack_(self, moldic): - self.__dict__.update(moldic) - return self + def intor(self, intor, comp=None, hermi=0, aosym='s1', out=None, + shls_slice=None, grids=None): + '''Integral generator. - dumps = dumps + Args: + intor : str + Name of the 1e or 2e AO integrals. Ref to :func:`getints` for the + complete list of available 1-electron integral names - @classmethod - @lib.with_doc(loads.__doc__) - def loads(cls, molstr): - return loads(molstr) + Kwargs: + comp : int + Components of the integrals, e.g. int1e_ipovlp_sph has 3 components. + hermi : int + Symmetry of the integrals - @lib.with_doc(loads.__doc__) - def loads_(self, molstr): - self.__dict__.update(loads(molstr).__dict__) - return self + | 0 : no symmetry assumed (default) + | 1 : hermitian + | 2 : anti-hermitian - def update(self, chkfile): - return self.update_from_chk(chkfile) - def update_from_chk(self, chkfile): - with h5py.File(chkfile, 'r') as fh5: - mol = loads(fh5['mol'][()]) - self.__dict__.update(mol.__dict__) - return self + grids : ndarray + Coordinates of grids for the int1e_grids integrals - def build(self, dump_input=True, parse_arg=ARGPARSE, - verbose=None, output=None, max_memory=None, - atom=None, basis=None, unit=None, nucmod=None, ecp=None, pseudo=None, - charge=None, spin=0, symmetry=None, symmetry_subgroup=None, - cart=None, magmom=None): - '''Setup moleclue and initialize some control parameters. Whenever you - change the value of the attributes of :class:`Mole`, you need call - this function to refresh the internal data of Mole. + Returns: + ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp - Kwargs: - dump_input : bool - whether to dump the contents of input file in the output file - parse_arg : bool - whether to read the sys.argv and overwrite the relevant parameters - verbose : int - Print level. If given, overwrite :attr:`Mole.verbose` - output : str or None - Output file. If given, overwrite :attr:`Mole.output` - max_memory : int, float - Allowd memory in MB. If given, overwrite :attr:`Mole.max_memory` - atom : list or str - To define molecluar structure. - basis : dict or str - To define basis set. - nucmod : dict or str - Nuclear model. If given, overwrite :attr:`Mole.nucmod` - charge : int - Charge of molecule. It affects the electron numbers - If given, overwrite :attr:`Mole.charge` - spin : int - 2S, num. alpha electrons - num. beta electrons to control - multiplicity. If setting spin = None , multiplicity will be - guessed based on the neutral molecule. - If given, overwrite :attr:`Mole.spin` - symmetry : bool or str - Whether to use symmetry. If given a string of point group - name, the given point group symmetry will be used. - magmom : list - Collinear spin of each atom. Default is [0.0,]*natm + Examples: + >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') + >>> mol.intor('int1e_ipnuc_sph', comp=3) # + [[[ 0. 0. ] + [ 0. 0. ]] + [[ 0. 0. ] + [ 0. 0. ]] + [[ 0.10289944 0.48176097] + [-0.48176097 -0.10289944]]] + >>> mol.intor('int1e_nuc_spinor') + [[-1.69771092+0.j 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j] + [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j -0.67146312+0.j] + [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] + [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] ''' - if not DISABLE_GC: - gc.collect() # To release circular referred objects + if not self._built: + logger.warn(self, 'Warning: intor envs of %s not initialized.', self) + # FIXME: Whether to check _built and call build? ._bas and .basis + # may not be consistent. calling .build() may leads to wrong intor env. + #self.build(False, False) + intor = self._add_suffix(intor) + bas = self._bas + env = self._env + if 'ECP' in intor: + assert (self._ecp is not None) + bas = numpy.vstack((self._bas, self._ecpbas)) + env[AS_ECPBAS_OFFSET] = len(self._bas) + env[AS_NECPBAS] = len(self._ecpbas) + if shls_slice is None: + shls_slice = (0, self.nbas, 0, self.nbas) + elif '_grids' in intor: + assert grids is not None + env = numpy.append(env, grids.ravel()) + env[NGRIDS] = grids.shape[0] + env[PTR_GRIDS] = env.size - grids.size + return moleintor.getints(intor, self._atm, bas, env, + shls_slice, comp, hermi, aosym, out=out) - if isinstance(dump_input, str): - sys.stderr.write('Assigning the first argument %s to mol.atom\n' % - dump_input) - dump_input, atom = True, dump_input + def _add_suffix(self, intor, cart=None): + if not (intor[:4] == 'cint' or + intor.endswith(('_sph', '_cart', '_spinor', '_ssc'))): + if cart is None: + cart = self.cart + if cart: + intor = intor + '_cart' + else: + intor = intor + '_sph' + return intor - if verbose is not None: self.verbose = verbose - if output is not None: self.output = output - if max_memory is not None: self.max_memory = max_memory - if atom is not None: self.atom = atom - if basis is not None: self.basis = basis - if unit is not None: self.unit = unit - if nucmod is not None: self.nucmod = nucmod - if ecp is not None: self.ecp = ecp - if pseudo is not None: self.pseudo = pseudo - if charge is not None: self.charge = charge - if spin != 0: self.spin = spin - if symmetry is not None: self.symmetry = symmetry - if symmetry_subgroup is not None: self.symmetry_subgroup = symmetry_subgroup - if cart is not None: self.cart = cart - if magmom is not None: self.magmom = magmom + def intor_symmetric(self, intor, comp=None, grids=None): + '''One-electron integral generator. The integrals are assumed to be hermitian - if parse_arg: - _update_from_cmdargs_(self) + Args: + intor : str + Name of the 1-electron integral. Ref to :func:`getints` for the + complete list of available 1-electron integral names - # avoid opening output file twice - if (self.output is not None - # StringIO() does not have attribute 'name' - and getattr(self.stdout, 'name', None) != self.output): + Kwargs: + comp : int + Components of the integrals, e.g. int1e_ipovlp_sph has 3 components. + grids : ndarray + Coordinates of grids for the int1e_grids integrals - if self.verbose > logger.QUIET: - if os.path.isfile(self.output): - print('overwrite output file: %s' % self.output) - else: - print('output file: %s' % self.output) + Returns: + ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp - if self.output == '/dev/null': - self.stdout = open(os.devnull, 'w') - else: - self.stdout = open(self.output, 'w') + Examples: - if self.verbose >= logger.WARN: - self.check_sanity() + >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') + >>> mol.intor_symmetric('int1e_nuc_spinor') + [[-1.69771092+0.j 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j] + [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j -0.67146312+0.j] + [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] + [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] + ''' + return self.intor(intor, comp, 1, aosym='s4', grids=grids) - self._atom = self.format_atom(self.atom, unit=self.unit) - uniq_atoms = set([a[0] for a in self._atom]) + def intor_asymmetric(self, intor, comp=None, grids=None): + '''One-electron integral generator. The integrals are assumed to be anti-hermitian - _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, - self.nucprop) + Args: + intor : str + Name of the 1-electron integral. Ref to :func:`getints` for the + complete list of available 1-electron integral names - if self.pseudo: - self.ecp, self.pseudo = classify_ecp_pseudo(self, self.ecp, self.pseudo) + Kwargs: + comp : int + Components of the integrals, e.g. int1e_ipovlp has 3 components. + grids : ndarray + Coordinates of grids for the int1e_grids integrals - if self.ecp: - # Unless explicitly input, ECP should not be assigned to ghost atoms - atoms_wo_ghost = [a for a in uniq_atoms if not is_ghost_atom(a)] - _ecp = _parse_default_basis(self.ecp, atoms_wo_ghost) - self._ecp = self.format_ecp(_ecp) - if self._ecp: - self._atm, self._ecpbas, self._env = \ - self.make_ecp_env(self._atm, self._ecp, self._env) + Returns: + ndarray of 1-electron integrals, can be either 2-dim or 3-dim, depending on comp - if self.pseudo: - # Unless explicitly input, PP should not be assigned to ghost atoms - atoms_wo_ghost = [a for a in uniq_atoms if not is_ghost_atom(a)] - _pseudo = _parse_default_basis(self.pseudo, atoms_wo_ghost) - self._pseudo = _pseudo = self.format_pseudo(_pseudo) - if _pseudo: - conflicts = set(_pseudo).intersection(self._ecp) - if conflicts: - raise RuntimeError('Pseudo potential for atoms %s are defined ' - 'in both .ecp and .pseudo.' % list(conflicts)) + Examples: - for ia, atom in enumerate(self._atom): - symb = atom[0] - if (symb in _pseudo and - # skip ghost atoms - self._atm[ia,0] != 0): - self._atm[ia,0] = sum(_pseudo[symb][0]) + >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g') + >>> mol.intor_asymmetric('int1e_nuc_spinor') + [[-1.69771092+0.j 0.00000000+0.j 0.67146312+0.j 0.00000000+0.j] + [ 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j 0.67146312+0.j] + [-0.67146312+0.j 0.00000000+0.j -1.69771092+0.j 0.00000000+0.j] + [ 0.00000000+0.j -0.67146312+0.j 0.00000000+0.j -1.69771092+0.j]] + ''' + return self.intor(intor, comp, 2, aosym='a4', grids=grids) - if self.spin is None: - self.spin = self.nelectron % 2 + @lib.with_doc(moleintor.getints_by_shell.__doc__) + def intor_by_shell(self, intor, shells, comp=None, grids=None): + intor = self._add_suffix(intor) + if 'ECP' in intor: + assert (self._ecp is not None) + bas = numpy.vstack((self._bas, self._ecpbas)) + self._env[AS_ECPBAS_OFFSET] = len(self._bas) + self._env[AS_NECPBAS] = len(self._ecpbas) else: - # Access self.nelec in which the code checks whether the spin and - # number of electrons are consistent. - self.nelec + bas = self._bas + return moleintor.getints_by_shell(intor, shells, self._atm, bas, + self._env, comp) - if not self.magmom: - self.magmom = [0.,]*self.natm - if self.spin == 0 and abs(numpy.sum(numpy.asarray(self.magmom)) - self.spin) > 1e-6: - #don't check for unrestricted calcs. - raise ValueError("mol.magmom is set incorrectly.") + eval_ao = eval_gto = eval_gto - if self.symmetry: - self._build_symmetry() + energy_nuc = get_enuc = energy_nuc - if dump_input and not self._built and self.verbose > logger.NOTE: - self.dump_input() + def get_ao_indices(self, bas_list, ao_loc=None): + ''' + Generate (dis-continued) AO indices for basis specified in bas_list + ''' + if ao_loc is None: + ao_loc = self.ao_loc + return lib.locs_to_indices(ao_loc, bas_list) - if self.verbose >= logger.DEBUG3: - logger.debug3(self, 'arg.atm = %s', self._atm) - logger.debug3(self, 'arg.bas = %s', self._bas) - logger.debug3(self, 'arg.env = %s', self._env) - logger.debug3(self, 'ecpbas = %s', self._ecpbas) + sph_labels = spheric_labels = sph_labels + cart_labels = cart_labels + ao_labels = ao_labels - self._built = True - return self - kernel = build + spinor_labels = spinor_labels - def dump_input(self): - import __main__ - import pyscf - if hasattr(__main__, '__file__'): - try: - filename = os.path.abspath(__main__.__file__) - finput = open(filename, 'r') - self.stdout.write('#INFO: **** input file is %s ****\n' % filename) - self.stdout.write(finput.read()) - self.stdout.write('#INFO: ******************** input file end ********************\n') - self.stdout.write('\n') - self.stdout.write('\n') - finput.close() - except IOError: - logger.warn(self, 'input file does not exist') + search_ao_label = search_ao_label - self.stdout.write('System: %s Threads %s\n' % - (str(platform.uname()), lib.num_threads())) - self.stdout.write('Python %s\n' % sys.version) - self.stdout.write('numpy %s scipy %s\n' % - (numpy.__version__, scipy.__version__)) - self.stdout.write('Date: %s\n' % time.ctime()) - self.stdout.write('PySCF version %s\n' % pyscf.__version__) - info = lib.repo_info(os.path.join(__file__, '..', '..')) - self.stdout.write('PySCF path %s\n' % info['path']) - if 'git' in info: - self.stdout.write(info['git'] + '\n') + def search_shell_id(self, atm_id, l): + return search_shell_id(self, atm_id, l) - self.stdout.write('\n') - for key in os.environ: - if 'PYSCF' in key: - self.stdout.write('[ENV] %s %s\n' % (key, os.environ[key])) - if self.verbose >= logger.DEBUG2: - for key in dir(__config__): - if key[:2] != '__': - self.stdout.write('[CONFIG] %s = %s\n' % - (key, getattr(__config__, key))) - else: - conf_file = getattr(__config__, 'conf_file', None) - self.stdout.write('[CONFIG] conf_file %s\n' % conf_file) + search_ao_nr = search_ao_nr + search_ao_r = search_ao_r - self.stdout.write('[INPUT] verbose = %d\n' % self.verbose) - if self.verbose >= logger.DEBUG: - self.stdout.write('[INPUT] max_memory = %s \n' % self.max_memory) - self.stdout.write('[INPUT] num. atoms = %d\n' % self.natm) - self.stdout.write('[INPUT] num. electrons = %d\n' % self.nelectron) - self.stdout.write('[INPUT] charge = %d\n' % self.charge) - self.stdout.write('[INPUT] spin (= nelec alpha-beta = 2S) = %d\n' % self.spin) - self.stdout.write('[INPUT] symmetry %s subgroup %s\n' % - (self.symmetry, self.symmetry_subgroup)) - self.stdout.write('[INPUT] Mole.unit = %s\n' % self.unit) - if self.cart: - self.stdout.write('[INPUT] Cartesian GTO integrals (6d 10f)\n') + aoslice_by_atom = aoslice_nr_by_atom = offset_ao_by_atom = offset_nr_by_atom = aoslice_by_atom + aoslice_2c_by_atom = offset_2c_by_atom = offset_2c_by_atom - self.stdout.write('[INPUT] Symbol X Y Z unit' - ' X Y Z unit Magmom\n') - for ia,atom in enumerate(self._atom): - coorda = tuple([x * param.BOHR for x in atom[1]]) - coordb = tuple([x for x in atom[1]]) - magmom = self.magmom[ia] - self.stdout.write('[INPUT]%3d %-4s %16.12f %16.12f %16.12f AA ' - '%16.12f %16.12f %16.12f Bohr %4.1f\n' - % ((ia+1, _symbol(atom[0])) + coorda + coordb + (magmom,))) - if self.nucmod: - if isinstance(self.nucmod, (int, str, types.FunctionType)): - nucatms = [_symbol(atom[0]) for atom in self._atom] - else: - nucatms = self.nucmod.keys() - self.stdout.write('[INPUT] Gaussian nuclear model for atoms %s\n' % - nucatms) + condense_to_shell = condense_to_shell + get_overlap_cond = get_overlap_cond + + to_uncontracted_cartesian_basis = to_uncontracted_cartesian_basis + decontract_basis = decontract_basis + + ao_rotation_matrix = ao_rotation_matrix + + def cart2sph_coeff(self, normalized='sp'): + '''Transformation matrix that transforms Cartesian GTOs to spherical + GTOs for all basis functions + + Kwargs: + normalized : string or boolean + How the Cartesian GTOs are normalized. Except s and p functions, + Cartesian GTOs do not have the universal normalization coefficients + for the different components of the same shell. The value of this + argument can be one of 'sp', 'all', None. 'sp' means the Cartesian s + and p basis are normalized. 'all' means all Cartesian functions are + normalized. None means none of the Cartesian functions are normalized. + The default value 'sp' is the convention used by libcint library. - if self.nucprop: - self.stdout.write('[INPUT] nucprop %s\n' % self.nucprop) + Examples: - if self.verbose >= logger.DEBUG: - self.stdout.write('[INPUT] ---------------- BASIS SET ---------------- \n') - self.stdout.write('[INPUT] l, kappa, [nprim/nctr], ' - 'expnt, c_1 c_2 ...\n') - for atom, basis_set in self._basis.items(): - self.stdout.write('[INPUT] %s\n' % atom) - for b in basis_set: - if isinstance(b[1], int): - kappa = b[1] - b_coeff = b[2:] - else: - kappa = 0 - b_coeff = b[1:] - nprim = len(b_coeff) - nctr = len(b_coeff[0])-1 - if nprim < nctr: - logger.warn(self, 'num. primitives smaller than num. contracted basis') - self.stdout.write('[INPUT] %d %2d [%-5d/%-4d] ' - % (b[0], kappa, nprim, nctr)) - for k, x in enumerate(b_coeff): - if k == 0: - self.stdout.write('%-15.12g ' % x[0]) - else: - self.stdout.write(' '*32+'%-15.12g ' % x[0]) - for c in x[1:]: - self.stdout.write(' %4.12g' % c) - self.stdout.write('\n') + >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') + >>> c = mol.cart2sph_coeff() + >>> s0 = mol.intor('int1e_ovlp_sph') + >>> s1 = c.T.dot(mol.intor('int1e_ovlp_cart')).dot(c) + >>> print(abs(s1-s0).sum()) + >>> 4.58676826646e-15 + ''' + c2s_l = [cart2sph(l, normalized=normalized) for l in range(12)] + c2s = [] + for ib in range(self.nbas): + l = self.bas_angular(ib) + for n in range(self.bas_nctr(ib)): + c2s.append(c2s_l[l]) + return scipy.linalg.block_diag(*c2s) - if self.verbose >= logger.INFO: - self.stdout.write('\n') - logger.info(self, 'nuclear repulsion = %.15g', self.energy_nuc()) - if self.symmetry: - if self.topgroup == self.groupname: - logger.info(self, 'point group symmetry = %s', self.topgroup) - else: - logger.info(self, 'point group symmetry = %s, use subgroup %s', - self.topgroup, self.groupname) - logger.info(self, "symmetry origin: %s", self._symm_orig) - logger.info(self, "symmetry axis x: %s", self._symm_axes[0]) - logger.info(self, "symmetry axis y: %s", self._symm_axes[1]) - logger.info(self, "symmetry axis z: %s", self._symm_axes[2]) - for ir in range(self.symm_orb.__len__()): - logger.info(self, 'num. orbitals of irrep %s = %d', - self.irrep_name[ir], self.symm_orb[ir].shape[1]) - logger.info(self, 'number of shells = %d', self.nbas) - logger.info(self, 'number of NR pGTOs = %d', self.npgto_nr()) - logger.info(self, 'number of NR cGTOs = %d', self.nao_nr()) - logger.info(self, 'basis = %s', self.basis) - logger.info(self, 'ecp = %s', self.ecp) - if self.verbose >= logger.DEBUG2: - for i in range(len(self._bas)): - exps = self.bas_exp(i) - logger.debug1(self, 'bas %d, expnt(s) = %s', i, str(exps)) + def sph2spinor_coeff(self): + '''Transformation matrix that transforms real-spherical GTOs to spinor + GTOs for all basis functions - logger.info(self, 'CPU time: %12.2f', logger.process_clock()) - return self + Examples: + + >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='ccpvtz') + >>> ca, cb = mol.sph2spinor_coeff() + >>> 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 + ''' + from pyscf.symm import sph + return sph.sph2spinor_coeff(self) def apply(self, fn, *args, **kwargs): if callable(fn): @@ -3507,133 +3667,19 @@ def _TemporaryMoleContext(self, method, args, args_bak): self._ctx_lock = None -class Mole(GTOIntegralMixin, GlobalOptionMixin, lib.StreamObject): - '''Molecule geometry data, GTO basis paramaters and global options - - Attributes: - verbose : int - Print level - output : str or None - Output file, default is None which dumps msg to sys.stdout - max_memory : int, float - Allowed memory in MB - charge : int - Charge of molecule. It affects the electron numbers - spin : int or None - 2S, num. alpha electrons - num. beta electrons to control - multiplicity. If spin = None is set, multiplicity will be guessed - based on the neutral molecule. - symmetry : bool or str - Whether to use symmetry. When this variable is set to True, the - molecule will be rotated and the highest rotation axis will be - placed z-axis. - If a string is given as the name of point group, the given point - group symmetry will be used. Note that the input molecular - coordinates will not be changed in this case. - symmetry_subgroup : str - subgroup - - atom : list or str - To define molecluar structure. The internal format is - - | atom = [[atom1, (x, y, z)], - | [atom2, (x, y, z)], - | ... - | [atomN, (x, y, z)]] - - unit : str - Angstrom or Bohr - basis : dict or str - To define basis set. - nucmod : dict or str or [function(nuc_charge, nucprop) => zeta] - Nuclear model. 0 or None means point nuclear model. Other - values will enable Gaussian nuclear model. If a function is - assigned to this attribute, the function will be called to - generate the nuclear charge distribution value "zeta" and the - relevant nuclear model will be set to Gaussian model. - Default is point nuclear model. - nucprop : dict - Nuclear properties (like g-factor 'g', quadrupole moments 'Q'). - It is needed by pyscf.prop module and submodules. - cart : boolean - Using Cartesian GTO basis and integrals (6d,10f,15g) - magmom : list - Collinear spin of each atom. Default is [0.0,]*natm - - ** Following attributes are generated by :func:`Mole.build` ** - - stdout : file object - Default is sys.stdout if :attr:`Mole.output` is not set - topgroup : str - Point group of the system. - groupname : str - The supported subgroup of the point group. It can be one of SO3, - Dooh, Coov, D2h, C2h, C2v, D2, Cs, Ci, C2, C1 - nelectron : int - sum of nuclear charges - :attr:`Mole.charge` - symm_orb : a list of numpy.ndarray - Symmetry adapted basis. Each element is a set of symm-adapted orbitals - for one irreducible representation. The list index does **not** correspond - to the id of irreducible representation. - irrep_id : a list of int - Each element is one irreducible representation id associated with the basis - stored in symm_orb. One irrep id stands for one irreducible representation - symbol. The irrep symbol and the relevant id are defined in - :attr:`symm.param.IRREP_ID_TABLE` - irrep_name : a list of str - Each element is one irreducible representation symbol associated with the basis - stored in symm_orb. The irrep symbols are defined in - :attr:`symm.param.IRREP_ID_TABLE` - _built : bool - To label whether :func:`Mole.build` has been called. It is to - ensure certain functions being initialized only once. - _basis : dict - like :attr:`Mole.basis`, the internal format which is returned from the - parser :func:`format_basis` - - ** Following attributes are arguments used by ``libcint`` library ** - - _atm : - :code:`[[charge, ptr-of-coord, nuc-model, ptr-zeta, 0, 0], [...]]` - each element reperesents one atom - natm : - number of atoms - _bas : - :code:`[[atom-id, angular-momentum, num-primitive-GTO, num-contracted-GTO, 0, ptr-of-exps, ptr-of-contract-coeff, 0], [...]]` - each element reperesents one shell - nbas : - number of shells - _env : - list of floats to store the coordinates, GTO exponents, contract-coefficients - - Examples: - - >>> mol = Mole(atom='H^2 0 0 0; H 0 0 1.1', basis='sto3g').build() - >>> print(mol.atom_symbol(0)) - H^2 - >>> print(mol.atom_pure_symbol(0)) - H - >>> print(mol.nao_nr()) - 2 - >>> print(mol.intor('int1e_ovlp_sph')) - [[ 0.99999999 0.43958641] - [ 0.43958641 0.99999999]] - >>> mol.charge = 1 - >>> mol.build() - has no attributes Charge - ''' # noqa: E501 +class Mole(MoleBase): + '''A Mole object to hold the basic information of a molecule. + ''' __add__ = conc_mol - energy_nuc = get_enuc = energy_nuc inertia_moment = inertia_moment tostring = tostring tofile = tofile def __init__(self, **kwargs): - GTOIntegralMixin.__init__(self) - GlobalOptionMixin.__init__(self) - if kwargs: - self.__dict__.update(kwargs) + MoleBase.__init__(self) + for key, val in kwargs.items(): + setattr(self, key, val) def fromstring(self, string, format='xyz'): '''Update the Mole object based on the input geometry string''' @@ -3669,6 +3715,11 @@ def __getattr__(self, key): from pyscf import __all__ # noqa from pyscf import scf, dft + for mod in (scf, dft): + method = getattr(mod, key, None) + if callable(method): + return method(self) + if 'TD' in key[:3]: if key in ('TDHF', 'TDA'): mf = scf.HF(self) @@ -3688,57 +3739,6 @@ def __getattr__(self, key): mf.run() return method - def _build_symmetry(self, *args, **kwargs): - ''' - Update symmetry related attributes: topgroup, groupname, _symm_orig, - _symm_axes, irrep_id, irrep_name, symm_orb - ''' - from pyscf import symm - - # TODO: Consider ECP info in point group symmetry initialization - self.topgroup, orig, axes = symm.detect_symm(self._atom, self._basis) - - if isinstance(self.symmetry, str): - self.symmetry = str(symm.std_symb(self.symmetry)) - groupname = None - if abs(axes - np.eye(3)).max() < symm.TOLERANCE: - if symm.check_symm(self.symmetry, self._atom, self._basis): - # Try to use original axes (issue #1209) - groupname = self.symmetry - axes = np.eye(3) - else: - logger.warn(self, 'Unable to to identify input symmetry using original axes.\n' - 'Different symmetry axes will be used.') - if groupname is None: - try: - groupname, axes = symm.as_subgroup(self.topgroup, axes, - self.symmetry) - except PointGroupSymmetryError as e: - raise PointGroupSymmetryError( - 'Unable to identify input symmetry %s. Try symmetry="%s"' % - (self.symmetry, self.topgroup)) from e - else: - groupname, axes = symm.as_subgroup(self.topgroup, axes, - self.symmetry_subgroup) - self._symm_orig = orig - self._symm_axes = axes - - if self.cart and groupname in ('Dooh', 'Coov', 'SO3'): - if groupname == 'Coov': - groupname, lgroup = 'C2v', groupname - else: - groupname, lgroup = 'D2h', groupname - logger.warn(self, 'This version does not support symmetry %s ' - 'for cartesian GTO basis. Its subgroup %s is used', - lgroup, groupname) - self.groupname = groupname - - self.symm_orb, self.irrep_id = \ - symm.symm_adapted_basis(self, groupname, orig, axes) - self.irrep_name = [symm.irrep_id2name(groupname, ir) - for ir in self.irrep_id] - return self - def ao2mo(self, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e', **kwargs): '''Integral transformation for arbitrary orbitals and arbitrary diff --git a/pyscf/lo/orth.py b/pyscf/lo/orth.py index c463038fe1..bf6fb7ac5f 100644 --- a/pyscf/lo/orth.py +++ b/pyscf/lo/orth.py @@ -287,7 +287,8 @@ def orth_ao(mf_or_mol, method=ORTH_METHOD, pre_orth_ao=REF_BASIS, s=None): ''' from pyscf import scf from pyscf.lo import nao - if isinstance(mf_or_mol, gto.Mole): + if isinstance(mf_or_mol, gto.MoleBase): + assert mf_or_mol.__class__ == gto.Mole mol = mf_or_mol mf = None else: diff --git a/pyscf/mcscf/__init__.py b/pyscf/mcscf/__init__.py index e404c234d2..8fecd54896 100644 --- a/pyscf/mcscf/__init__.py +++ b/pyscf/mcscf/__init__.py @@ -193,8 +193,8 @@ def CASSCF(mf_or_mol, ncas, nelecas, ncore=None, frozen=None): from pyscf import gto from pyscf import scf - if isinstance(mf_or_mol, gto.Mole): - mf = scf.RHF(mf_or_mol) + if isinstance(mf_or_mol, gto.MoleBase): + mf = mf_or_mol.RHF() else: mf = mf_or_mol @@ -215,8 +215,8 @@ def CASSCF(mf_or_mol, ncas, nelecas, ncore=None, frozen=None): def CASCI(mf_or_mol, ncas, nelecas, ncore=None): from pyscf import gto from pyscf import scf - if isinstance(mf_or_mol, gto.Mole): - mf = scf.RHF(mf_or_mol) + if isinstance(mf_or_mol, gto.MoleBase): + mf = mf_or_mol.RHF() else: mf = mf_or_mol @@ -238,8 +238,8 @@ def CASCI(mf_or_mol, ncas, nelecas, ncore=None): def UCASCI(mf_or_mol, ncas, nelecas, ncore=None): from pyscf import gto from pyscf import scf - if isinstance(mf_or_mol, gto.Mole): - mf = scf.UHF(mf_or_mol) + if isinstance(mf_or_mol, gto.MoleBase): + mf = mf_or_mol.UHF() else: mf = mf_or_mol @@ -252,8 +252,8 @@ def UCASCI(mf_or_mol, ncas, nelecas, ncore=None): def UCASSCF(mf_or_mol, ncas, nelecas, ncore=None, frozen=None): from pyscf import gto from pyscf import scf - if isinstance(mf_or_mol, gto.Mole): - mf = scf.UHF(mf_or_mol) + if isinstance(mf_or_mol, gto.MoleBase): + mf = mf_or_mol.UHF() else: mf = mf_or_mol @@ -271,8 +271,8 @@ def DFCASSCF(mf_or_mol, ncas, nelecas, auxbasis=None, ncore=None, frozen=None): from pyscf import gto from pyscf import scf - if isinstance(mf_or_mol, gto.Mole): - mf = scf.RHF(mf_or_mol).density_fit() + if isinstance(mf_or_mol, gto.MoleBase): + mf = mf_or_mol.RHF().density_fit() else: mf = mf_or_mol @@ -288,8 +288,8 @@ def DFCASSCF(mf_or_mol, ncas, nelecas, auxbasis=None, ncore=None, def DFCASCI(mf_or_mol, ncas, nelecas, auxbasis=None, ncore=None): from pyscf import gto from pyscf import scf - if isinstance(mf_or_mol, gto.Mole): - mf = scf.RHF(mf_or_mol).density_fit() + if isinstance(mf_or_mol, gto.MoleBase): + mf = mf_or_mol.RHF().density_fit() else: mf = mf_or_mol diff --git a/pyscf/mcscf/casci.py b/pyscf/mcscf/casci.py index 0bbea9f61d..88869aa0b2 100644 --- a/pyscf/mcscf/casci.py +++ b/pyscf/mcscf/casci.py @@ -19,6 +19,7 @@ import sys from functools import reduce +import warnings import numpy from pyscf import lib from pyscf.lib import logger @@ -648,21 +649,12 @@ def __init__(self, mc): self._scf = mc._scf.as_scanner() def __call__(self, mol_or_geom, mo_coeff=None, ci0=None): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) - # These properties can be updated when calling mf_scanner(mol) if - # they are shared with mc._scf. In certain scenario the properties - # may be created for mc separately, e.g. when mcscf.approx_hessian is - # called. For safety, the code below explicitly resets these - # properties. self.reset (mol) - for key in ('with_df', 'with_x2c', 'with_solvent', 'with_dftd3'): - sub_mod = getattr(self, key, None) - if sub_mod: - sub_mod.reset(mol) if mo_coeff is None: mf_scanner = self._scf @@ -915,7 +907,8 @@ def _eig(self, h, *args): return scf.hf.eig(h, None) def get_h2cas(self, mo_coeff=None): - raise DeprecationWarning + '''An alias of get_h2eff method''' + return self.get_h2eff(mo_coeff) def get_h2eff(self, mo_coeff=None): '''Compute the active space two-particle Hamiltonian. @@ -928,9 +921,10 @@ def ao2mo(self, mo_coeff=None): raise NotImplementedError def get_h1cas(self, mo_coeff=None, ncas=None, ncore=None): - raise DeprecationWarning + '''An alias of get_h1eff method''' + return self.get_h1eff(mo_coeff, ncas, ncore) - get_h1eff = h1e_for_cas + get_h1eff = h1e_for_cas = h1e_for_cas def casci(self, mo_coeff=None, ci0=None, verbose=None): raise NotImplementedError diff --git a/pyscf/mcscf/mc1step.py b/pyscf/mcscf/mc1step.py index 9ea4e58f70..9ce98e9317 100644 --- a/pyscf/mcscf/mc1step.py +++ b/pyscf/mcscf/mc1step.py @@ -21,7 +21,7 @@ from functools import reduce import numpy import scipy.linalg -from pyscf.gto import Mole +from pyscf import gto from pyscf import lib from pyscf.lib import logger from pyscf.mcscf import casci @@ -538,7 +538,7 @@ def __init__(self, mc): def __call__(self, mol_or_geom, **kwargs): from pyscf.mcscf.addons import project_init_guess - if isinstance(mol_or_geom, Mole): + if isinstance(mol_or_geom, gto.MoleBase): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/mp/mp2.py b/pyscf/mp/mp2.py index b6adc6ebc4..a90548f1d2 100644 --- a/pyscf/mp/mp2.py +++ b/pyscf/mp/mp2.py @@ -419,7 +419,7 @@ def __init__(self, mp): self._scf = mp._scf.as_scanner() def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) diff --git a/pyscf/pbc/df/incore.py b/pyscf/pbc/df/incore.py index ac4ddf1bd0..016afca89d 100644 --- a/pyscf/pbc/df/incore.py +++ b/pyscf/pbc/df/incore.py @@ -66,7 +66,7 @@ def aux_e2(cell, auxcell_or_auxbasis, intor='int3c2e', aosym='s1', comp=None, Returns: (nao_pair, naux) array ''' - if isinstance(auxcell_or_auxbasis, (gto.Mole, pbcgto.Cell)): + if isinstance(auxcell_or_auxbasis, gto.MoleBase): auxcell = auxcell_or_auxbasis else: assert isinstance(auxcell_or_auxbasis, str) @@ -132,7 +132,7 @@ def wrap_int3c(cell, auxcell, intor='int3c2e', aosym='s1', comp=1, def fill_2c2e(cell, auxcell_or_auxbasis, intor='int2c2e', hermi=0, kpt=np.zeros(3)): '''2-center 2-electron AO integrals (L|ij), where L is the auxiliary basis. ''' - if isinstance(auxcell_or_auxbasis, gto.Mole): + if isinstance(auxcell_or_auxbasis, gto.MoleBase): auxcell = auxcell_or_auxbasis else: auxcell = make_auxcell(cell, auxcell_or_auxbasis) diff --git a/pyscf/pbc/df/outcore.py b/pyscf/pbc/df/outcore.py index bddeac2132..b4d5ad67ca 100644 --- a/pyscf/pbc/df/outcore.py +++ b/pyscf/pbc/df/outcore.py @@ -38,7 +38,7 @@ def aux_e1(cell, auxcell_or_auxbasis, erifile, intor='int3c2e', aosym='s2ij', co kptij_lst : (*,2,3) array A list of (kpti, kptj) ''' - if isinstance(auxcell_or_auxbasis, gto.Mole): + if isinstance(auxcell_or_auxbasis, gto.MoleBase): auxcell = auxcell_or_auxbasis else: auxcell = make_auxcell(cell, auxcell_or_auxbasis) @@ -159,7 +159,7 @@ def _aux_e2(cell, auxcell_or_auxbasis, erifile, intor='int3c2e', aosym='s2ij', c kptij_lst : (*,2,3) array A list of (kpti, kptj) ''' - if isinstance(auxcell_or_auxbasis, gto.Mole): + if isinstance(auxcell_or_auxbasis, gto.MoleBase): auxcell = auxcell_or_auxbasis else: auxcell = make_auxcell(cell, auxcell_or_auxbasis) diff --git a/pyscf/pbc/df/rsdf_helper.py b/pyscf/pbc/df/rsdf_helper.py index bc93221ecc..dac891ef3b 100644 --- a/pyscf/pbc/df/rsdf_helper.py +++ b/pyscf/pbc/df/rsdf_helper.py @@ -673,14 +673,14 @@ def _get_3c2e_Rcuts(bas_lst_or_mol, auxbas_lst_or_auxmol, dijs_lst, omega, where i and j shls are separated by d specified by "dijs_lst". """ - if isinstance(bas_lst_or_mol, mol_gto.mole.Mole): + if isinstance(bas_lst_or_mol, mol_gto.mole.MoleBase): mol = bas_lst_or_mol else: bas_lst = bas_lst_or_mol mol = MoleNoBasSort() mol.build(dump_input=False, parse_arg=False, atom="H 0 0 0", basis=bas_lst, spin=None) - if isinstance(auxbas_lst_or_auxmol, mol_gto.mole.Mole): + if isinstance(auxbas_lst_or_auxmol, mol_gto.mole.MoleBase): auxmol = auxbas_lst_or_auxmol else: auxbas_lst = auxbas_lst_or_auxmol @@ -1002,7 +1002,7 @@ def _aux_e2_nospltbas(cell, auxcell_or_auxbasis, omega, erifile, ''' log = logger.Logger(cell.stdout, cell.verbose) - if isinstance(auxcell_or_auxbasis, mol_gto.Mole): + if isinstance(auxcell_or_auxbasis, mol_gto.MoleBase): auxcell = auxcell_or_auxbasis else: auxcell = make_auxcell(cell, auxcell_or_auxbasis) diff --git a/pyscf/pbc/dft/__init__.py b/pyscf/pbc/dft/__init__.py index 0877830c44..9731c164d6 100644 --- a/pyscf/pbc/dft/__init__.py +++ b/pyscf/pbc/dft/__init__.py @@ -34,9 +34,12 @@ from pyscf.pbc.dft.rks import KohnShamDFT from pyscf.pbc.lib import kpts as libkpts -gto.Cell.GKS = GKS = gks.GKS -gto.Cell.UKS = UKS = uks.UKS -gto.Cell.ROKS = ROKS = roks.ROKS +GKS = gks.GKS +UKS = uks.UKS +ROKS = roks.ROKS +gto.Cell.GKS = property(GKS) +gto.Cell.UKS = property(UKS) +gto.Cell.ROKS = property(ROKS) def KRKS(cell, *args, **kwargs): for arg in args: @@ -46,7 +49,7 @@ def KRKS(cell, *args, **kwargs): if isinstance(kwargs['kpts'], libkpts.KPoints): return krks_ksymm.KRKS(cell, *args, **kwargs) return krks.KRKS(cell, *args, **kwargs) -gto.Cell.KRKS = KRKS +gto.Cell.KRKS = property(KRKS) def KUKS(cell, *args, **kwargs): for arg in args: @@ -56,10 +59,12 @@ def KUKS(cell, *args, **kwargs): if isinstance(kwargs['kpts'], libkpts.KPoints): return kuks_ksymm.KUKS(cell, *args, **kwargs) return kuks.KUKS(cell, *args, **kwargs) -gto.Cell.KUKS = KUKS +gto.Cell.KUKS = property(KUKS) -gto.Cell.KROKS = KROKS = kroks.KROKS -gto.Cell.KGKS = KGKS = kgks.KGKS +KROKS = kroks.KROKS +KGKS = kgks.KGKS +gto.Cell.KROKS = property(KROKS) +gto.Cell.KGKS = property(KGKS) def KRKSpU(cell, *args, **kwargs): for arg in args: @@ -69,7 +74,7 @@ def KRKSpU(cell, *args, **kwargs): if isinstance(kwargs['kpts'], libkpts.KPoints): return krkspu_ksymm.KRKSpU(cell, *args, **kwargs) return krkspu.KRKSpU(cell, *args, **kwargs) -gto.Cell.KRKSpU = KRKSpU +gto.Cell.KRKSpU = property(KRKSpU) def KUKSpU(cell, *args, **kwargs): for arg in args: @@ -79,7 +84,7 @@ def KUKSpU(cell, *args, **kwargs): if isinstance(kwargs['kpts'], libkpts.KPoints): return kukspu_ksymm.KUKSpU(cell, *args, **kwargs) return kukspu.KUKSpU(cell, *args, **kwargs) -gto.Cell.KUKSpU = KUKSpU +gto.Cell.KUKSpU = property(KUKSpU) def RKS(cell, *args, **kwargs): if cell.spin == 0: @@ -87,7 +92,7 @@ def RKS(cell, *args, **kwargs): else: return roks.ROKS(cell, *args, **kwargs) RKS.__doc__ = rks.RKS.__doc__ -gto.Cell.RKS = RKS +gto.Cell.RKS = property(RKS) def KS(cell, *args, **kwargs): if cell.spin == 0: @@ -97,7 +102,7 @@ def KS(cell, *args, **kwargs): KS.__doc__ = ''' A wrap function to create DFT object (RKS or UKS) for PBC systems.\n ''' + rks.RKS.__doc__ -gto.Cell.KS = KS +gto.Cell.KS = property(KS) def KKS(cell, *args, **kwargs): if cell.spin == 0: @@ -107,4 +112,4 @@ def KKS(cell, *args, **kwargs): KKS.__doc__ = ''' A wrap function to create DFT object with k-point sampling (KRKS or KUKS).\n ''' + krks.KRKS.__doc__ -gto.Cell.KKS = KKS +gto.Cell.KKS = property(KKS) diff --git a/pyscf/pbc/gto/cell.py b/pyscf/pbc/gto/cell.py index b9354fa8ff..ddab942d05 100644 --- a/pyscf/pbc/gto/cell.py +++ b/pyscf/pbc/gto/cell.py @@ -152,8 +152,8 @@ def conc_cell(cell1, cell2): cell3 = Cell() cell3.__dict__.update(mol3.__dict__) - # Whether to update the lattice_vectors? - cell3.a = cell1.a + # lattice_vectors needs to be consistent with cell3.unit (Bohr) + cell3.a = cell1.lattice_vectors() cell3.mesh = np.max((cell1.mesh, cell2.mesh), axis=0) ke_cutoff1 = cell1.ke_cutoff @@ -833,7 +833,7 @@ def _mesh_inf_vaccum(cell): return int(meshz*.5 + .999) * 2 -class Cell(mole.GTOIntegralMixin, mole.GlobalOptionMixin, lib.StreamObject): +class Cell(mole.MoleBase): '''A Cell object holds the basic information of a crystal. Attributes: @@ -892,8 +892,7 @@ class Cell(mole.GTOIntegralMixin, mole.GlobalOptionMixin, lib.StreamObject): } def __init__(self, **kwargs): - mole.GTOIntegralMixin.__init__(self) - mole.GlobalOptionMixin.__init__(self) + mole.MoleBase.__init__(self) self.a = None # lattice vectors, (a1,a2,a3) # if set, defines a spherical cutoff # of fourier components, with .5 * G**2 < ke_cutoff @@ -980,6 +979,11 @@ def __getattr__(self, key): from pyscf.pbc import scf, dft from pyscf.dft import XC + for mod in (scf, dft): + method = getattr(mod, key, None) + if callable(method): + return method(self) + if key[0] == 'K': # with k-point sampling if 'TD' in key[:4]: if key in ('KTDHF', 'KTDA'): @@ -1128,7 +1132,7 @@ def build(self, dump_input=True, parse_arg=mole.ARGPARSE, self.symmorphic = symmorphic dump_input = dump_input and not self._built and self.verbose > logger.NOTE - mole.GlobalOptionMixin.build(self, dump_input, parse_arg, *args, **kwargs) + mole.MoleBase.build(self, dump_input, parse_arg, *args, **kwargs) exp_min = np.array([self.bas_exp(ib).min() for ib in range(self.nbas)]) if self.exp_to_discard is None: diff --git a/pyscf/pbc/gto/pseudo/__init__.py b/pyscf/pbc/gto/pseudo/__init__.py index 55b9095095..64166afeb3 100644 --- a/pyscf/pbc/gto/pseudo/__init__.py +++ b/pyscf/pbc/gto/pseudo/__init__.py @@ -16,7 +16,7 @@ # Author: Qiming Sun # Timothy Berkelbach -from pyscf.gto.basis import parse +from pyscf.gto.basis.parse_cp2k_pp import parse from pyscf.gto.basis import load_pseudo as load from pyscf.gto.basis import PP_ALIAS as ALIAS from pyscf.pbc.gto.pseudo.pp import * diff --git a/pyscf/pbc/gto/test/test_cell.py b/pyscf/pbc/gto/test/test_cell.py index 61dc9f6e72..7d9f23528b 100644 --- a/pyscf/pbc/gto/test/test_cell.py +++ b/pyscf/pbc/gto/test/test_cell.py @@ -466,6 +466,7 @@ def test_conc_cell(self): self.assertTrue(len(cl3._ecpbas), 20) self.assertTrue(len(cl3._bas), 12) self.assertTrue(len(cl3._atm), 8) + self.assertAlmostEqual(abs(cl3.lattice_vectors() - cl1.lattice_vectors()).max(), 0, 12) def test_eval_gto(self): cell = pgto.M(a=np.eye(3)*4, atom='He 1 1 1', basis=[[2,(1,.5),(.5,.5)]], precision=1e-10) diff --git a/pyscf/scf/hf.py b/pyscf/scf/hf.py index 335becbeef..9009ceb8d4 100644 --- a/pyscf/scf/hf.py +++ b/pyscf/scf/hf.py @@ -1349,7 +1349,7 @@ def __init__(self, mf_obj): self._last_mol_fp = mf_obj.mol.ao_loc def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False) @@ -1638,10 +1638,6 @@ def init_guess_by_1e(self, mol=None): @lib.with_doc(init_guess_by_chkfile.__doc__) def init_guess_by_chkfile(self, chkfile=None, project=None): - if isinstance(chkfile, gto.Mole): - raise TypeError(''' - You see this error message because of the API updates. - The first argument needs to be the name of a chkfile.''') if chkfile is None: chkfile = self.chkfile return init_guess_by_chkfile(self.mol, chkfile, project=project) def from_chk(self, chkfile=None, project=None): diff --git a/pyscf/tdscf/rhf.py b/pyscf/tdscf/rhf.py index 6f5d69bec2..d3d62d812f 100644 --- a/pyscf/tdscf/rhf.py +++ b/pyscf/tdscf/rhf.py @@ -626,7 +626,7 @@ def __init__(self, td): self._scf = td._scf.as_scanner() def __call__(self, mol_or_geom, **kwargs): - if isinstance(mol_or_geom, gto.Mole): + if isinstance(mol_or_geom, gto.MoleBase): mol = mol_or_geom else: mol = self.mol.set_geom_(mol_or_geom, inplace=False)