xref: /dports/science/py-pyscf/pyscf-2.0.1/pyscf/gto/mole.py

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

'''
Mole class and helper functions to handle paramters and attributes for GTO
integrals. This module serves the interface to the integral library libcint.
'''

import os, sys
import types
import re
import platform
import gc
import time

import json
import ctypes
import numpy
import numpy as np
import h5py
import scipy.special
import scipy.linalg
import contextlib
import threading
from pyscf import lib
from pyscf.lib import param
from pyscf.data import elements
from pyscf.lib import logger
from pyscf.gto import cmd_args
from pyscf.gto import basis
from pyscf.gto import moleintor
from pyscf.gto.eval_gto import eval_gto
from pyscf.gto.ecp import core_configuration
from pyscf import __config__

from pyscf.data.elements import ELEMENTS, ELEMENTS_PROTON, \
        _rm_digit, charge, _symbol, _std_symbol, _atom_symbol, is_ghost_atom, \
        _std_symbol_without_ghost

from pyscf.lib.exceptions import BasisNotFoundError, PointGroupSymmetryError
import warnings

# For code compatibility in python-2 and python-3
if sys.version_info >= (3,):
    unicode = str


# for _atm, _bas, _env
CHARGE_OF  = 0
PTR_COORD  = 1
NUC_MOD_OF = 2
PTR_ZETA   = PTR_FRAC_CHARGE = 3
ATM_SLOTS  = 6
ATOM_OF    = 0
ANG_OF     = 1
NPRIM_OF   = 2
NCTR_OF    = 3
RADI_POWER = 3 # for ECP
KAPPA_OF   = 4
SO_TYPE_OF = 4 # for ECP
PTR_EXP    = 5
PTR_COEFF  = 6
BAS_SLOTS  = 8
# pointer to env
PTR_EXPCUTOFF   = 0
PTR_COMMON_ORIG = 1
PTR_RINV_ORIG   = 4
PTR_RINV_ZETA   = 7
PTR_RANGE_OMEGA = 8
PTR_F12_ZETA    = 9
PTR_GTG_ZETA    = 10
NGRIDS          = 11
PTR_GRIDS       = 12
AS_RINV_ORIG_ATOM = 17
AS_ECPBAS_OFFSET = 18
AS_NECPBAS      = 19
PTR_ENV_START   = 20
# parameters from libcint
NUC_POINT = 1
NUC_GAUSS = 2
# nucleus with fractional charges. It can be used to mimic MM particles
NUC_FRAC_CHARGE = 3
NUC_ECP = 4  # atoms with pseudo potential

BASE = getattr(__config__, 'BASE', 0)
NORMALIZE_GTO = getattr(__config__, 'NORMALIZE_GTO', True)
DISABLE_EVAL = getattr(__config__, 'DISABLE_EVAL', False)
# Whether to disable the explicit call to gc.collect(). gc.collect() may cause
# non-negligible overhead (https://github.com/pyscf/pyscf/issues/1038).
DISABLE_GC = getattr(__config__, 'DISABLE_GC', False)

def M(**kwargs):
    r'''This is a shortcut to build up Mole object.

    Args: Same to :func:`Mole.build`

    Examples:

    >>> from pyscf import gto
    >>> mol = gto.M(atom='H 0 0 0; F 0 0 1', basis='6-31g')
    '''
    mol = Mole()
    mol.build(**kwargs)
    return mol

def gaussian_int(n, alpha):
    r'''int_0^inf x^n exp(-alpha x^2) dx'''
    n1 = (n + 1) * .5
    return scipy.special.gamma(n1) / (2. * alpha**n1)

def gto_norm(l, expnt):
    r'''Normalized factor for GTO radial part   :math:`g=r^l e^{-\alpha r^2}`

    .. math::

        \frac{1}{\sqrt{\int g^2 r^2 dr}}
        = \sqrt{\frac{2^{2l+3} (l+1)! (2a)^{l+1.5}}{(2l+2)!\sqrt{\pi}}}

    Ref: H. B. Schlegel and M. J. Frisch, Int. J. Quant.  Chem., 54(1995), 83-87.

    Args:
        l (int):
            angular momentum
        expnt :
            exponent :math:`\alpha`

    Returns:
        normalization factor

    Examples:

    >>> print(gto_norm(0, 1))
    2.5264751109842591
    '''
    if l >= 0:
        #f = 2**(2*l+3) * math.factorial(l+1) * (2*expnt)**(l+1.5) \
        #        / (math.factorial(2*l+2) * math.sqrt(math.pi))
        #return math.sqrt(f)
        return 1/numpy.sqrt(gaussian_int(l*2+2, 2*expnt))
    else:
        raise ValueError('l should be > 0')

def cart2sph(l, c_tensor=None, normalized=None):
    '''
    Cartesian to real spherical transformation matrix

    Kwargs:
        normalized :
            How the Cartesian GTOs are normalized.  'sp' means the s and p
            functions are normalized (this is the convention used by libcint
            library).
    '''
    nf = (l+1)*(l+2)//2
    if c_tensor is None:
        c_tensor = numpy.eye(nf)
    else:
        c_tensor = numpy.asarray(c_tensor, order='F').reshape(-1,nf)
    if l == 0 or l == 1:
        if normalized == 'sp':
            return c_tensor
        elif l == 0:
            return c_tensor * 0.282094791773878143
        else:
            return c_tensor * 0.488602511902919921
    else:
        assert(l <= 12)
        nd = l * 2 + 1
        ngrid = c_tensor.shape[0]
        c2sph = numpy.zeros((ngrid,nd), order='F')
        fn = moleintor.libcgto.CINTc2s_ket_sph
        fn(c2sph.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(ngrid),
           c_tensor.ctypes.data_as(ctypes.c_void_p), ctypes.c_int(l))
        return c2sph

def cart2spinor_kappa(kappa, l=None, normalized=None):
    '''Cartesian to spinor transformation matrix for kappa

    Kwargs:
        normalized :
            How the Cartesian GTOs are normalized.  'sp' means the s and p
            functions are normalized (this is the convention used by libcint
            library).
    '''
    if kappa < 0:
        l = -kappa - 1
        nd = l * 2 + 2
    elif kappa > 0:
        l = kappa
        nd = l * 2
    else:
        assert(l is not None)
        assert(l <= 12)
        nd = l * 4 + 2
    nf = (l+1)*(l+2)//2
    c2smat = numpy.zeros((nf*2,nd), order='F', dtype=numpy.complex128)
    cmat = numpy.eye(nf)
    fn = moleintor.libcgto.CINTc2s_ket_spinor_sf1
    fn(c2smat.ctypes.data_as(ctypes.c_void_p),
       c2smat[nf:].ctypes.data_as(ctypes.c_void_p),
       cmat.ctypes.data_as(ctypes.c_void_p),
       ctypes.c_int(nf*2), ctypes.c_int(nf),
       ctypes.c_int(1), ctypes.c_int(kappa), ctypes.c_int(l))
    if normalized != 'sp':
        if l == 0:
            c2smat *= 0.282094791773878143
        elif l == 1:
            c2smat *= 0.488602511902919921
    # c2smat[0] is the transformation for spin up
    # c2smat[1] is the transformation for spin down
    c2smat = c2smat.reshape(2,nf,nd)
    return c2smat
cart2j_kappa = cart2spinor_kappa

def cart2spinor_l(l, normalized=None):
    '''Cartesian to spinor transformation matrix for angular moment l

    Kwargs:
        normalized :
            How the Cartesian GTOs are normalized.  'sp' means the s and p
            functions are normalized (this is the convention used by libcint
            library).
    '''
    return cart2spinor_kappa(0, l, normalized)
cart2j_l = cart2spinor_l

def sph2spinor_kappa(kappa, l=None):
    '''Real spherical to spinor transformation matrix for kappa'''
    from pyscf.symm.sph import sph2spinor
    ua, ub = sph2spinor(l)
    if kappa < 0:
        l = -kappa - 1
        ua = ua[:,l*2:]
        ub = ub[:,l*2:]
    elif kappa > 0:
        l = kappa
        ua = ua[:,:l*2]
        ub = ub[:,:l*2]
    else:
        assert(l is not None)
        assert(l <= 12)
    return ua, ub

def sph2spinor_l(l):
    '''Real spherical to spinor transformation matrix for angular moment l'''
    return sph2spinor_kappa(0, l)

def atom_types(atoms, basis=None):
    '''symmetry inequivalent atoms'''
    atmgroup = {}
    for ia, a in enumerate(atoms):
        if 'GHOST' in a[0].upper():
            a = ['X'+a[0][5:]] + list(a[1:])
        if a[0] in atmgroup:
            atmgroup[a[0]].append(ia)
        elif basis is None:
            atmgroup[a[0]] = [ia]
        else:
            stdsymb = _std_symbol(a[0])
            if a[0] in basis:
                if stdsymb in basis and basis[a[0]] == basis[stdsymb]:
                    if stdsymb in atmgroup:
                        atmgroup[stdsymb].append(ia)
                    else:
                        atmgroup[stdsymb] = [ia]
                else:
                    atmgroup[a[0]] = [ia]
            elif stdsymb in atmgroup:
                atmgroup[stdsymb].append(ia)
            else:
                atmgroup[stdsymb] = [ia]
    return atmgroup


def format_atom(atoms, origin=0, axes=None,
                unit=getattr(__config__, 'UNIT', 'Ang')):
    '''Convert the input :attr:`Mole.atom` to the internal data format.
    Including, changing the nuclear charge to atom symbol, converting the
    coordinates to AU, rotate and shift the molecule.
    If the :attr:`~Mole.atom` is a string, it takes ";" and "\\n"
    for the mark to separate atoms;  "," and arbitrary length of blank space
    to spearate the individual terms for an atom.  Blank line will be ignored.

    Args:
        atoms : list or str
            the same to :attr:`Mole.atom`

    Kwargs:
        origin : ndarray
            new axis origin.
        axes : ndarray
            (new_x, new_y, new_z), new coordinates
        unit : str or number
            If unit is one of strings (B, b, Bohr, bohr, AU, au), the coordinates
            of the input atoms are the atomic unit;  If unit is one of strings
            (A, a, Angstrom, angstrom, Ang, ang), the coordinates are in the
            unit of angstrom;  If a number is given, the number are considered
            as the Bohr value (in angstrom), which should be around 0.53.
            Set unit=1 if wishing to preserve the unit of the coordinates.

    Returns:
        "atoms" in the internal format. The internal format is
            | atom = [[atom1, (x, y, z)],
            |         [atom2, (x, y, z)],
            |         ...
            |         [atomN, (x, y, z)]]

    Examples:

    >>> gto.format_atom('9,0,0,0; h@1 0 0 1', origin=(1,1,1))
    [['F', [-1.0, -1.0, -1.0]], ['H@1', [-1.0, -1.0, 0.0]]]
    >>> gto.format_atom(['9,0,0,0', (1, (0, 0, 1))], origin=(1,1,1))
    [['F', [-1.0, -1.0, -1.0]], ['H', [-1, -1, 0]]]
    '''
    def str2atm(line):
        dat = line.split()
        try:
            coords = [float(x) for x in dat[1:4]]
        except ValueError:
            if DISABLE_EVAL:
                raise ValueError('Failed to parse geometry %s' % line)
            else:
                coords = list(eval(','.join(dat[1:4])))
        if len(coords) != 3:
            raise ValueError('Coordinates error in %s' % line)
        return [_atom_symbol(dat[0]), coords]

    if isinstance(atoms, (str, unicode)):
        # The input atoms points to a geometry file
        if os.path.isfile(atoms):
            try:
                atoms = fromfile(atoms)
            except ValueError:
                sys.stderr.write('\nFailed to parse geometry file  %s\n\n' % atoms)
                raise

        atoms = str(atoms.replace(';','\n').replace(',',' ').replace('\t',' '))
        fmt_atoms = []
        for dat in atoms.split('\n'):
            dat = dat.strip()
            if dat and dat[0] != '#':
                fmt_atoms.append(dat)

        if len(fmt_atoms[0].split()) < 4:
            fmt_atoms = from_zmatrix('\n'.join(fmt_atoms))
        else:
            fmt_atoms = [str2atm(line) for line in fmt_atoms]
    else:
        fmt_atoms = []
        for atom in atoms:
            if isinstance(atom, (str, unicode)):
                if atom.lstrip()[0] != '#':
                    fmt_atoms.append(str2atm(atom.replace(',',' ')))
            else:
                if isinstance(atom[1], (int, float)):
                    fmt_atoms.append([_atom_symbol(atom[0]), atom[1:4]])
                else:
                    fmt_atoms.append([_atom_symbol(atom[0]), atom[1]])

    if len(fmt_atoms) == 0:
        return []

    if axes is None:
        axes = numpy.eye(3)

    if isinstance(unit, (str, unicode)):
        if unit.upper().startswith(('B', 'AU')):
            unit = 1.
        else: #unit[:3].upper() == 'ANG':
            unit = 1./param.BOHR
    else:
        unit = 1./unit

    c = numpy.array([a[1] for a in fmt_atoms], dtype=numpy.double)
    c = numpy.einsum('ix,kx->ki', axes * unit, c - origin)
    z = [a[0] for a in fmt_atoms]
    return list(zip(z, c.tolist()))

#TODO: sort exponents
def format_basis(basis_tab):
    '''Convert the input :attr:`Mole.basis` to the internal data format.

    ``{ atom: [(l, ((-exp, c_1, c_2, ..),
                    (-exp, c_1, c_2, ..))),
               (l, ((-exp, c_1, c_2, ..),
                    (-exp, c_1, c_2, ..)))], ... }``

    Args:
        basis_tab : dict
            Similar to :attr:`Mole.basis`, it **cannot** be a str

    Returns:
        Formated :attr:`~Mole.basis`

    Examples:

    >>> gto.format_basis({'H':'sto-3g', 'H^2': '3-21g'})
    {'H': [[0,
        [3.4252509099999999, 0.15432897000000001],
        [0.62391373000000006, 0.53532813999999995],
        [0.16885539999999999, 0.44463454000000002]]],
     'H^2': [[0,
        [5.4471780000000001, 0.15628500000000001],
        [0.82454700000000003, 0.90469100000000002]],
        [0, [0.18319199999999999, 1.0]]]}
    '''
    basis_converter = _generate_basis_converter()
    fmt_basis = {}
    for atom, atom_basis in basis_tab.items():
        symb = _atom_symbol(atom)
        fmt_basis[symb] = basis_converter(symb, atom_basis)

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

def _generate_basis_converter():
    def nparray_to_list(item):
        val = []
        for x in item:
            if isinstance(x, (tuple, list)):
                val.append(nparray_to_list(x))
            elif isinstance(x, numpy.ndarray):
                val.append(x.tolist())
            else:
                val.append(x)
        return val

    def load(basis_name, symb):
        if basis_name.lower().startswith('unc'):
            return uncontract(basis.load(basis_name[3:], symb))
        else:
            return basis.load(basis_name, symb)

    def converter(symb, raw_basis):
        if isinstance(raw_basis, (str, unicode)):
            bset = load(str(raw_basis), _std_symbol_without_ghost(symb))
        elif (any(isinstance(x, (str, unicode)) for x in raw_basis)
              # The first element is the basis of internal format
              or not isinstance(raw_basis[0][0], int)):
            stdsymb = _std_symbol_without_ghost(symb)
            bset = []
            for rawb in raw_basis:
                if isinstance(rawb, (str, unicode)):
                    bset += load(str(rawb), stdsymb)
                else:
                    bset += nparray_to_list(rawb)
        else:
            bset = nparray_to_list(raw_basis)
        return bset
    return converter

def uncontracted_basis(_basis):
    '''Uncontract internal format _basis

    Examples:

    >>> gto.uncontract(gto.load('sto3g', 'He'))
    [[0, [6.36242139, 1]], [0, [1.158923, 1]], [0, [0.31364979, 1]]]
    '''
    MAXL = 10
    ubasis_raw = [[] for l in range(MAXL)]
    ubasis_exp = [[] for l in range(MAXL)]
    for b in _basis:
        angl = b[0]
        kappa = b[1]
        if isinstance(kappa, int):
            coeffs = b[2:]
        else:
            coeffs = b[1:]

        if isinstance(kappa, int) and kappa != 0:
            warnings.warn('For basis with kappa != 0, the uncontract basis might be wrong. '
                          'Please double check the resultant attribute mol._basis')
            for p in coeffs:
                ubasis_raw[angl].append([angl, kappa, [p[0], 1]])
                ubasis_exp[angl].append(p[0])
        else:
            for p in coeffs:
                ubasis_raw[angl].append([angl, [p[0], 1]])
                ubasis_exp[angl].append(p[0])

    # Check linear dependency
    ubasis = []
    for l in range(MAXL):
        basis_l = ubasis_raw[l]
        if basis_l:
            es = numpy.array(ubasis_exp[l])
            # Remove duplicated primitive basis functions
            es, e_idx = numpy.unique(es.round(9), True)
            # from large exponent to small exponent
            for i in reversed(e_idx):
                ubasis.append(basis_l[i])
    return ubasis
uncontract = uncontracted_basis
contract = contracted_basis = basis.to_general_contraction

def to_uncontracted_cartesian_basis(mol):
    '''Decontract the basis of a Mole or a Cell.  Returns a Mole (Cell) object
    with the uncontracted basis environment and a list of coefficients that
    transform the uncontracted cartesian basis to the original basis.  Each
    element in the list corresponds to one shell of the original Mole (Cell).

    Examples:

    >>> mol = gto.M(atom='Ne', basis='ccpvdz')
    >>> pmol, ctr_coeff = mol.to_uncontracted_cartesian_basis()
    >>> c = scipy.linalg.block_diag(*ctr_coeff)
    >>> s = reduce(numpy.dot, (c.T, pmol.intor('int1e_ovlp'), c))
    >>> abs(s-mol.intor('int1e_ovlp')).max()
    0.0
    '''
    import copy
    lmax = mol._bas[:,ANG_OF].max()
    if mol.cart:
        c2s = [numpy.eye((l+1)*(l+2)//2) for l in range(lmax+1)]
    else:
        c2s = [cart2sph(l, normalized='sp') for l in range(lmax+1)]

    pmol = copy.copy(mol)
    pmol.cart = True
    _bas = []
    _env = mol._env.copy()
    contr_coeff = []
    for ib in range(mol.nbas):
        l = mol._bas[ib,ANG_OF]
        ncart = (l+1)*(l+2)//2
        es = mol.bas_exp(ib)
        cs = mol._libcint_ctr_coeff(ib)
        nprim = cs.shape[0]
        norm = gto_norm(l, es)
        c = numpy.einsum('pi,p,xm->pxim', cs, 1./norm, c2s[l])
        contr_coeff.append(c.reshape(nprim*ncart,-1))

        pexp = mol._bas[ib,PTR_EXP]
        pc = mol._bas[ib,PTR_COEFF]
        bs = numpy.empty((nprim,8), dtype=numpy.int32)
        bs[:] = mol._bas[ib]
        bs[:,NCTR_OF] = bs[:,NPRIM_OF] = 1
        bs[:,PTR_EXP] = numpy.arange(pexp, pexp+nprim)
        bs[:,PTR_COEFF] = numpy.arange(pc, pc+nprim)
        _env[pc:pc+nprim] = norm
        _bas.append(bs)

    pmol._bas = numpy.asarray(numpy.vstack(_bas), dtype=numpy.int32)
    pmol._env = _env
    return pmol, contr_coeff

def format_ecp(ecp_tab):
    r'''Convert the input :attr:`ecp` (dict) to the internal data format::

      { atom: (nelec,  # core electrons
               ((l,  # l=-1 for UL, l>=0 for Ul to indicate |l><l|
                 (((exp_1, c_1),  # for r^0
                   (exp_2, c_2),
                   ...),
                  ((exp_1, c_1),  # for r^1
                   (exp_2, c_2),
                   ...),
                  ((exp_1, c_1),  # for r^2
                   ...))))),
       ...}
    '''
    fmt_ecp = {}
    for atom, atom_ecp in ecp_tab.items():
        symb = _atom_symbol(atom)

        if isinstance(atom_ecp, (str, unicode)):
            stdsymb = _std_symbol_without_ghost(symb)
            ecp_dat = basis.load_ecp(str(atom_ecp), stdsymb)
            if ecp_dat is None or len(ecp_dat) == 0:
                #raise BasisNotFoundError('ECP not found for  %s' % symb)
                sys.stderr.write('ECP %s not found for  %s\n' %
                                 (atom_ecp, symb))
            else:
                fmt_ecp[symb] = ecp_dat
        else:
            fmt_ecp[symb] = atom_ecp
    return fmt_ecp

# transform etb to basis format
def expand_etb(l, n, alpha, beta):
    r'''Generate the exponents of even tempered basis for :attr:`Mole.basis`.
    .. math::

        e = e^{-\alpha * \beta^{i-1}} for i = 1 .. n

    Args:
        l : int
            Angular momentum
        n : int
            Number of GTOs

    Returns:
        Formated :attr:`~Mole.basis`

    Examples:

    >>> gto.expand_etb(1, 3, 1.5, 2)
    [[1, [6.0, 1]], [1, [3.0, 1]], [1, [1.5, 1]]]
    '''
    return [[l, [alpha*beta**i, 1]] for i in reversed(range(n))]
def expand_etbs(etbs):
    r'''Generate even tempered basis.  See also :func:`expand_etb`

    Args:
        etbs = [(l, n, alpha, beta), (l, n, alpha, beta),...]

    Returns:
        Formated :attr:`~Mole.basis`

    Examples:

    >>> gto.expand_etbs([(0, 2, 1.5, 2.), (1, 2, 1, 2.)])
    [[0, [6.0, 1]], [0, [3.0, 1]], [1, [1., 1]], [1, [2., 1]]]
    '''
    return lib.flatten([expand_etb(*etb) for etb in etbs])
etbs = expand_etbs

# concatenate two mol
def conc_env(atm1, bas1, env1, atm2, bas2, env2):
    r'''Concatenate two Mole's integral parameters.  This function can be used
    to construct the environment for cross integrals like

    .. math::

        \langle \mu | \nu \rangle, \mu \in mol1, \nu \in mol2

    Returns:
        Concatenated atm, bas, env

    Examples:
        Compute the overlap between H2 molecule and O atom

    >>> mol1 = gto.M(atom='H 0 1 0; H 0 0 1', basis='sto3g')
    >>> mol2 = gto.M(atom='O 0 0 0', basis='sto3g')
    >>> atm3, bas3, env3 = gto.conc_env(mol1._atm, mol1._bas, mol1._env,
    ...                                 mol2._atm, mol2._bas, mol2._env)
    >>> gto.moleintor.getints('int1e_ovlp_sph', atm3, bas3, env3, range(2), range(2,5))
    [[ 0.04875181  0.44714688  0.          0.37820346  0.        ]
     [ 0.04875181  0.44714688  0.          0.          0.37820346]]
    '''
    off = len(env1)
    natm_off = len(atm1)
    atm2 = numpy.copy(atm2)
    bas2 = numpy.copy(bas2)
    atm2[:,PTR_COORD] += off
    atm2[:,PTR_ZETA ] += off
    bas2[:,ATOM_OF  ] += natm_off
    bas2[:,PTR_EXP  ] += off
    bas2[:,PTR_COEFF] += off
    return (numpy.asarray(numpy.vstack((atm1,atm2)), dtype=numpy.int32),
            numpy.asarray(numpy.vstack((bas1,bas2)), dtype=numpy.int32),
            numpy.hstack((env1,env2)))

def conc_mol(mol1, mol2):
    '''Concatenate two Mole objects.
    '''
    if not mol1._built:
        logger.warn(mol1, 'Warning: object %s not initialized. Initializing %s',
                    mol1, mol1)
        mol1.build()
    if not mol2._built:
        logger.warn(mol2, 'Warning: object %s not initialized. Initializing %s',
                    mol2, mol2)
        mol2.build()

    mol3 = Mole()
    mol3._built = True

    mol3._atm, mol3._bas, mol3._env = \
            conc_env(mol1._atm, mol1._bas, mol1._env,
                     mol2._atm, mol2._bas, mol2._env)
    off = len(mol1._env)
    natm_off = len(mol1._atm)
    if len(mol2._ecpbas) == 0:
        mol3._ecpbas = mol1._ecpbas
    else:
        ecpbas2 = numpy.copy(mol2._ecpbas)
        ecpbas2[:,ATOM_OF  ] += natm_off
        ecpbas2[:,PTR_EXP  ] += off
        ecpbas2[:,PTR_COEFF] += off
        if len(mol1._ecpbas) == 0:
            mol3._ecpbas = ecpbas2
        else:
            mol3._ecpbas = numpy.hstack((mol1._ecpbas, ecpbas2))

    mol3.verbose = mol1.verbose
    mol3.output = mol1.output
    mol3.max_memory = mol1.max_memory
    mol3.charge = mol1.charge + mol2.charge
    mol3.spin = abs(mol1.spin - mol2.spin)
    mol3.symmetry = False
    mol3.symmetry_subgroup = None
    mol3.cart = mol1.cart and mol2.cart

    mol3._atom = mol1._atom + mol2._atom
    mol3.atom = mol3._atom
    mol3.unit = 'Bohr'

    mol3._basis.update(mol2._basis)
    mol3._basis.update(mol1._basis)
    mol3._pseudo.update(mol2._pseudo)
    mol3._pseudo.update(mol1._pseudo)
    mol3._ecp.update(mol2._ecp)
    mol3._ecp.update(mol1._ecp)
    mol3.basis = mol3._basis
    mol3.ecp = mol3._ecp

    mol3.nucprop.update(mol1.nucprop)
    mol3.nucprop.update(mol2.nucprop)
    return mol3

# <bas-of-mol1|intor|bas-of-mol2>
def intor_cross(intor, mol1, mol2, comp=None, grids=None):
    r'''1-electron integrals from two molecules like

    .. math::

        \langle \mu | intor | \nu \rangle, \mu \in mol1, \nu \in mol2

    Args:
        intor : str
            Name of the 1-electron integral, such as int1e_ovlp_sph (spherical overlap),
            int1e_nuc_cart (cartesian nuclear attraction), int1e_ipovlp_spinor
            (spinor overlap gradients), etc.  Ref to :func:`getints` for the
            full list of available 1-electron integral names
        mol1, mol2:
            :class:`Mole` objects

    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

    Examples:
        Compute the overlap between H2 molecule and O atom

    >>> mol1 = gto.M(atom='H 0 1 0; H 0 0 1', basis='sto3g')
    >>> mol2 = gto.M(atom='O 0 0 0', basis='sto3g')
    >>> gto.intor_cross('int1e_ovlp_sph', mol1, mol2)
    [[ 0.04875181  0.44714688  0.          0.37820346  0.        ]
     [ 0.04875181  0.44714688  0.          0.          0.37820346]]
    '''
    nbas1 = len(mol1._bas)
    nbas2 = len(mol2._bas)
    atmc, basc, envc = conc_env(mol1._atm, mol1._bas, mol1._env,
                                mol2._atm, mol2._bas, mol2._env)
    if '_grids' in intor:
        assert grids is not None
        envc = numpy.append(envc, grids.ravel())
        envc[NGRIDS] = grids.shape[0]
        envc[PTR_GRIDS] = envc.size - grids.size

    shls_slice = (0, nbas1, nbas1, nbas1+nbas2)

    if (intor.endswith('_sph') or intor.startswith('cint') or
        intor.endswith('_spinor') or intor.endswith('_cart')):
        return moleintor.getints(intor, atmc, basc, envc, shls_slice, comp, 0)
    elif mol1.cart == mol2.cart:
        intor = mol1._add_suffix(intor)
        return moleintor.getints(intor, atmc, basc, envc, shls_slice, comp, 0)
    elif mol1.cart:
        mat = moleintor.getints(intor+'_cart', atmc, basc, envc, shls_slice, comp, 0)
        return numpy.dot(mat, mol2.cart2sph_coeff())
    else:
        mat = moleintor.getints(intor+'_cart', atmc, basc, envc, shls_slice, comp, 0)
        return numpy.dot(mol1.cart2sph_coeff().T, mat)

# append (charge, pointer to coordinates, nuc_mod) to _atm
def make_atm_env(atom, ptr=0, nuclear_model=NUC_POINT, nucprop={}):
    '''Convert the internal format :attr:`Mole._atom` to the format required
    by ``libcint`` integrals
    '''
    nuc_charge = charge(atom[0])
    if nuclear_model == NUC_POINT:
        zeta = 0
    elif nuclear_model == NUC_GAUSS:
        zeta = dyall_nuc_mod(nuc_charge, nucprop)
    else:  # callable(nuclear_model)
        zeta = nuclear_model(nuc_charge, nucprop)
        nuclear_model = NUC_GAUSS

    _env = numpy.hstack((atom[1], zeta))
    _atm = numpy.zeros(6, dtype=numpy.int32)
    _atm[CHARGE_OF] = nuc_charge
    _atm[PTR_COORD] = ptr
    _atm[NUC_MOD_OF] = nuclear_model
    _atm[PTR_ZETA ] = ptr + 3
    return _atm, _env

# append (atom, l, nprim, nctr, kappa, ptr_exp, ptr_coeff, 0) to bas
# absorb normalization into GTO contraction coefficients
def make_bas_env(basis_add, atom_id=0, ptr=0):
    '''Convert :attr:`Mole.basis` to the argument ``bas`` for ``libcint`` integrals
    '''
    _bas = []
    _env = []
    for b in basis_add:
        if not b:  # == []
            continue
        angl = b[0]
        if angl > 14:
            sys.stderr.write('Warning: integral library does not support basis '
                             'with angular momentum > 14\n')

        if isinstance(b[1], int):
            kappa = b[1]
            b_coeff = numpy.array(sorted(list(b[2:]), reverse=True))
        else:
            kappa = 0
            b_coeff = numpy.array(sorted(list(b[1:]), reverse=True))
        es = b_coeff[:,0]
        cs = b_coeff[:,1:]
        nprim, nctr = cs.shape
        cs = numpy.einsum('pi,p->pi', cs, gto_norm(angl, es))
        if NORMALIZE_GTO:
            cs = _nomalize_contracted_ao(angl, es, cs)

        _env.append(es)
        _env.append(cs.T.reshape(-1))
        ptr_exp = ptr
        ptr_coeff = ptr_exp + nprim
        ptr = ptr_coeff + nprim * nctr
        _bas.append([atom_id, angl, nprim, nctr, kappa, ptr_exp, ptr_coeff, 0])
    _env = lib.flatten(_env) # flatten nested lists
    return (numpy.array(_bas, numpy.int32).reshape(-1,BAS_SLOTS),
            numpy.array(_env, numpy.double))

def _nomalize_contracted_ao(l, es, cs):
    #ee = numpy.empty((nprim,nprim))
    #for i in range(nprim):
    #    for j in range(i+1):
    #        ee[i,j] = ee[j,i] = gaussian_int(angl*2+2, es[i]+es[j])
    #s1 = 1/numpy.sqrt(numpy.einsum('pi,pq,qi->i', cs, ee, cs))
    ee = es.reshape(-1,1) + es.reshape(1,-1)
    ee = gaussian_int(l*2+2, ee)
    s1 = 1. / numpy.sqrt(numpy.einsum('pi,pq,qi->i', cs, ee, cs))
    return numpy.einsum('pi,i->pi', cs, s1)

def make_env(atoms, basis, pre_env=[], nucmod={}, nucprop={}):
    '''Generate the input arguments for ``libcint`` library based on internal
    format :attr:`Mole._atom` and :attr:`Mole._basis`
    '''
    _atm = []
    _bas = []
    _env = []
    ptr_env = len(pre_env)

    for ia, atom in enumerate(atoms):
        symb = atom[0]
        stdsymb = _rm_digit(symb)

        if ia+1 in nucprop:
            prop = nucprop[ia+1]
        elif symb in nucprop:
            prop = nucprop[symb]
        else:
            prop = nucprop.get(stdsymb, {})

        nuclear_model = NUC_POINT
        if nucmod:
            if nucmod is None:
                nuclear_model = NUC_POINT
            elif isinstance(nucmod, (int, str, unicode, types.FunctionType)):
                nuclear_model = _parse_nuc_mod(nucmod)
            elif ia+1 in nucmod:
                nuclear_model = _parse_nuc_mod(nucmod[ia+1])
            elif symb in nucmod:
                nuclear_model = _parse_nuc_mod(nucmod[symb])
            elif stdsymb in nucmod:
                nuclear_model = _parse_nuc_mod(nucmod[stdsymb])
        atm0, env0 = make_atm_env(atom, ptr_env, nuclear_model, prop)
        ptr_env = ptr_env + len(env0)
        _atm.append(atm0)
        _env.append(env0)

    _basdic = {}
    for symb, basis_add in basis.items():
        bas0, env0 = make_bas_env(basis_add, 0, ptr_env)
        ptr_env = ptr_env + len(env0)
        _basdic[symb] = bas0
        _env.append(env0)

    for ia, atom in enumerate(atoms):
        symb = atom[0]
        puresymb = _rm_digit(symb)
        if symb in _basdic:
            b = _basdic[symb].copy()
        elif puresymb in _basdic:
            b = _basdic[puresymb].copy()
        else:
            if symb[:2].upper() == 'X-':
                symb = symb[2:]
            elif symb[:6].upper() == 'GHOST-':
                symb = symb[6:]
            puresymb = _rm_digit(symb)
            if symb in _basdic:
                b = _basdic[symb].copy()
            elif puresymb in _basdic:
                b = _basdic[puresymb].copy()
            else:
                sys.stderr.write('Warning: Basis not found for atom %d %s\n' % (ia, symb))
                continue
        b[:,ATOM_OF] = ia
        _bas.append(b)

    if _atm:
        _atm = numpy.asarray(numpy.vstack(_atm), numpy.int32).reshape(-1, ATM_SLOTS)
    else:
        _atm = numpy.zeros((0,ATM_SLOTS), numpy.int32)
    if _bas:
        _bas = numpy.asarray(numpy.vstack(_bas), numpy.int32).reshape(-1, BAS_SLOTS)
    else:
        _bas = numpy.zeros((0,BAS_SLOTS), numpy.int32)
    if _env:
        _env = numpy.hstack((pre_env,numpy.hstack(_env)))
    else:
        _env = numpy.array(pre_env, copy=False)
    return _atm, _bas, _env

def make_ecp_env(mol, _atm, ecp, pre_env=[]):
    '''Generate the input arguments _ecpbas for ECP integrals
    '''
    _env = []
    ptr_env = len(pre_env)

    _ecpdic = {}
    for symb, ecp_add in ecp.items():
        ecp0 = []
        nelec = ecp_add[0]
        for lb in ecp_add[1]:
            for rorder, bi in enumerate(lb[1]):
                if len(bi) > 0:
                    ec = numpy.array(sorted(bi, reverse=True))
                    nexp, ncol = ec.shape
                    _env.append(ec[:,0])
                    _env.append(ec[:,1])
                    ptr_exp, ptr_coeff = ptr_env, ptr_env + nexp
                    ecp0.append([0, lb[0], nexp, rorder, 0,
                                 ptr_exp, ptr_coeff, 0])
                    ptr_env += nexp * 2

                    if ncol == 3:  # Has SO-ECP
                        _env.append(ec[:,2])
                        ptr_coeff, ptr_env = ptr_env, ptr_env + nexp
                        ecp0.append([0, lb[0], nexp, rorder, 1,
                                     ptr_exp, ptr_coeff, 0])

        _ecpdic[symb] = (nelec, numpy.asarray(ecp0, dtype=numpy.int32))

    _ecpbas = []
    if _ecpdic:
        _atm = _atm.copy()
        for ia, atom in enumerate(mol._atom):
            symb = atom[0]
            if symb in _ecpdic:
                ecp0 = _ecpdic[symb]
            elif _rm_digit(symb) in _ecpdic:
                ecp0 = _ecpdic[_rm_digit(symb)]
            else:
                ecp0 = None
            if ecp0 is not None:
                _atm[ia,CHARGE_OF ] = charge(symb) - ecp0[0]
                _atm[ia,NUC_MOD_OF] = NUC_ECP
                b = ecp0[1].copy().reshape(-1,BAS_SLOTS)
                b[:,ATOM_OF] = ia
                _ecpbas.append(b)
    if _ecpbas:
        _ecpbas = numpy.asarray(numpy.vstack(_ecpbas), numpy.int32)
        _env = numpy.hstack([pre_env] + _env)
    else:
        _ecpbas = numpy.zeros((0,BAS_SLOTS), numpy.int32)
        _env = pre_env
    return _atm, _ecpbas, _env

def tot_electrons(mol):
    '''Total number of electrons for the given molecule

    Returns:
        electron number in integer

    Examples:

    >>> mol = gto.M(atom='H 0 1 0; C 0 0 1', charge=1)
    >>> mol.tot_electrons()
    6
    '''
    if mol._atm.size != 0:
        nelectron = mol.atom_charges().sum()
    elif mol._atom:
        nelectron = sum(charge(a[0]) for a in mol._atom)
    else:
        nelectron = sum(charge(a[0]) for a in format_atom(mol.atom))
    nelectron -= mol.charge
    nelectron_int = int(round(nelectron))

    if abs(nelectron - nelectron_int) > 1e-4:
        logger.warn(mol, 'Found fractional number of electrons %f. Round it to %d',
                    nelectron, nelectron_int)
    return nelectron_int

def copy(mol):
    '''Deepcopy of the given :class:`Mole` object
    '''
    import copy
    newmol = copy.copy(mol)
    newmol._atm    = numpy.copy(mol._atm)
    newmol._bas    = numpy.copy(mol._bas)
    newmol._env    = numpy.copy(mol._env)
    newmol._ecpbas = numpy.copy(mol._ecpbas)

    newmol.atom    = copy.deepcopy(mol.atom)
    newmol._atom   = copy.deepcopy(mol._atom)
    newmol.basis   = copy.deepcopy(mol.basis)
    newmol._basis  = copy.deepcopy(mol._basis)
    newmol.ecp     = copy.deepcopy(mol.ecp)
    newmol._ecp    = copy.deepcopy(mol._ecp)
    return newmol

def pack(mol):
    '''Pack the input args of :class:`Mole` to a dict.

    Note this function only pack the input arguments (not the entire Mole
    class).  Modifications to mol._atm, mol._bas, mol._env are not tracked.
    Use :func:`dumps` to serialize the entire Mole object.
    '''
    mdic = {'atom'    : mol.atom,
            'unit'    : mol.unit,
            'basis'   : mol.basis,
            'charge'  : mol.charge,
            'spin'    : mol.spin,
            'symmetry': mol.symmetry,
            'nucmod'  : mol.nucmod,
            'nucprop' : mol.nucprop,
            'ecp'     : mol.ecp,
            '_nelectron': mol._nelectron,
            'verbose' : mol.verbose}
    return mdic
def unpack(moldic):
    '''Unpack a dict which is packed by :func:`pack`, to generate the input
    arguments for :class:`Mole` object.
    '''
    mol = Mole()
    mol.__dict__.update(moldic)
    return mol


def dumps(mol):
    '''Serialize Mole object to a JSON formatted str.
    '''
    exclude_keys = set(('output', 'stdout', '_keys',
                        # Constructing in function loads
                        'symm_orb', 'irrep_id', 'irrep_name'))
    nparray_keys = set(('_atm', '_bas', '_env', '_ecpbas',
                        '_symm_orig', '_symm_axes'))

    moldic = dict(mol.__dict__)
    for k in exclude_keys:
        if k in moldic:
            del(moldic[k])
    for k in nparray_keys:
        if isinstance(moldic[k], numpy.ndarray):
            moldic[k] = moldic[k].tolist()
    moldic['atom'] = repr(mol.atom)
    moldic['basis']= repr(mol.basis)
    moldic['ecp' ] = repr(mol.ecp)

    try:
        return json.dumps(moldic)
    except TypeError:
        def skip_value(dic):
            dic1 = {}
            for k,v in dic.items():
                if (v is None or
                    isinstance(v, (str, unicode, bool, int, float))):
                    dic1[k] = v
                elif isinstance(v, (list, tuple)):
                    dic1[k] = v   # Should I recursively skip_vaule?
                elif isinstance(v, set):
                    dic1[k] = list(v)
                elif isinstance(v, dict):
                    dic1[k] = skip_value(v)
                else:
                    msg =('Function mol.dumps drops attribute %s because '
                          'it is not JSON-serializable' % k)
                    warnings.warn(msg)
            return dic1
        return json.dumps(skip_value(moldic), skipkeys=True)

def loads(molstr):
    '''Deserialize a str containing a JSON document to a Mole object.
    '''
    # the numpy function array is used by eval function
    from numpy import array  # noqa
    moldic = json.loads(molstr)
    if sys.version_info < (3,):
        # Convert to utf8 because JSON loads fucntion returns unicode.
        def byteify(inp):
            if isinstance(inp, dict):
                return dict([(byteify(k), byteify(v)) for k, v in inp.iteritems()])
            elif isinstance(inp, (tuple, list)):
                return [byteify(x) for x in inp]
            elif isinstance(inp, unicode):
                return inp.encode('utf-8')
            else:
                return inp
        moldic = byteify(moldic)
    mol = Mole()
    mol.__dict__.update(moldic)
    mol.atom = eval(mol.atom)
    mol.basis= eval(mol.basis)
    mol.ecp  = eval(mol.ecp)
    mol._atm = numpy.array(mol._atm, dtype=numpy.int32)
    mol._bas = numpy.array(mol._bas, dtype=numpy.int32)
    mol._env = numpy.array(mol._env, dtype=numpy.double)
    mol._ecpbas = numpy.array(mol._ecpbas, dtype=numpy.int32)

    if mol.symmetry and mol._symm_orig is not None:
        from pyscf import symm
        mol._symm_orig = numpy.array(mol._symm_orig)
        mol._symm_axes = numpy.array(mol._symm_axes)
        mol.symm_orb, mol.irrep_id = \
                symm.symm_adapted_basis(mol, mol.groupname,
                                        mol._symm_orig, mol._symm_axes)
        mol.irrep_name = [symm.irrep_id2name(mol.groupname, ir)
                           for ir in mol.irrep_id]

    elif mol.symmetry and mol.symm_orb is not None:
        # Backward compatibility. To load symm_orb from chkfile of pyscf-1.6
        # and earlier.
        symm_orb = []

        # decompress symm_orb
        for val, x, y, shape in mol.symm_orb:
            if isinstance(val[0], list):
                # backward compatibility for chkfile of pyscf-1.4 in which val
                # is an array of real floats. In pyscf-1.5, val can be a list
                # of list, to include the imaginary part
                val_real, val_imag = val
            else:
                val_real, val_imag = val, None
            if val_imag is None:
                c = numpy.zeros(shape)
                c[numpy.array(x),numpy.array(y)] = numpy.array(val_real)
            else:
                c = numpy.zeros(shape, dtype=numpy.complex128)
                val = numpy.array(val_real) + numpy.array(val_imag) * 1j
                c[numpy.array(x),numpy.array(y)] = val
            symm_orb.append(c)
        mol.symm_orb = symm_orb

    return mol


def len_spinor(l, kappa):
    '''The number of spinor associated with given angular momentum and kappa.  If kappa is 0,
    return 4l+2
    '''
    if kappa == 0:
        n = (l * 4 + 2)
    elif kappa < 0:
        n = (l * 2 + 2)
    else:
        n = (l * 2)
    return n

def len_cart(l):
    '''The number of Cartesian function associated with given angular momentum.
    '''
    return (l + 1) * (l + 2) // 2

def npgto_nr(mol, cart=None):
    '''Total number of primitive spherical GTOs for the given :class:`Mole` object'''
    if cart is None:
        cart = mol.cart
    l = mol._bas[:,ANG_OF]
    if cart:
        return ((l+1)*(l+2)//2 * mol._bas[:,NPRIM_OF]).sum()
    else:
        return ((l*2+1) * mol._bas[:,NPRIM_OF]).sum()
def nao_nr(mol, cart=None):
    '''Total number of contracted GTOs for the given :class:`Mole` object'''
    if cart is None:
        cart = mol.cart
    if cart:
        return nao_cart(mol)
    else:
        return ((mol._bas[:,ANG_OF]*2+1) * mol._bas[:,NCTR_OF]).sum()
def nao_cart(mol):
    '''Total number of contracted cartesian GTOs for the given :class:`Mole` object'''
    l = mol._bas[:,ANG_OF]
    return ((l+1)*(l+2)//2 * mol._bas[:,NCTR_OF]).sum()

# nao_id0:nao_id1 corresponding to bas_id0:bas_id1
def nao_nr_range(mol, bas_id0, bas_id1):
    '''Lower and upper boundary of contracted spherical basis functions associated
    with the given shell range

    Args:
        mol :
            :class:`Mole` object
        bas_id0 : int
            start shell id
        bas_id1 : int
            stop shell id

    Returns:
        tupel of start basis function id and the stop function id

    Examples:

    >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g')
    >>> gto.nao_nr_range(mol, 2, 4)
    (2, 6)
    '''
    ao_loc = moleintor.make_loc(mol._bas[:bas_id1], 'sph')
    nao_id0 = ao_loc[bas_id0]
    nao_id1 = ao_loc[-1]
    return nao_id0, nao_id1

def nao_2c(mol):
    '''Total number of contracted spinor GTOs for the given :class:`Mole` object'''
    l = mol._bas[:,ANG_OF]
    kappa = mol._bas[:,KAPPA_OF]
    dims = (l*4+2) * mol._bas[:,NCTR_OF]
    dims[kappa<0] = l[kappa<0] * 2 + 2
    dims[kappa>0] = l[kappa>0] * 2
    return dims.sum()

# nao_id0:nao_id1 corresponding to bas_id0:bas_id1
def nao_2c_range(mol, bas_id0, bas_id1):
    '''Lower and upper boundary of contracted spinor basis functions associated
    with the given shell range

    Args:
        mol :
            :class:`Mole` object
        bas_id0 : int
            start shell id, 0-based
        bas_id1 : int
            stop shell id, 0-based

    Returns:
        tupel of start basis function id and the stop function id

    Examples:

    >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g')
    >>> gto.nao_2c_range(mol, 2, 4)
    (4, 12)
    '''
    ao_loc = moleintor.make_loc(mol._bas[:bas_id1], '')
    nao_id0 = ao_loc[bas_id0]
    nao_id1 = ao_loc[-1]
    return nao_id0, nao_id1

def ao_loc_nr(mol, cart=None):
    '''Offset of every shell in the spherical basis function spectrum

    Returns:
        list, each entry is the corresponding start basis function id

    Examples:

    >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g')
    >>> gto.ao_loc_nr(mol)
    [0, 1, 2, 3, 6, 9, 10, 11, 12, 15, 18]
    '''
    if cart is None:
        cart = mol.cart
    if cart:
        return moleintor.make_loc(mol._bas, 'cart')
    else:
        return moleintor.make_loc(mol._bas, 'sph')

def ao_loc_2c(mol):
    '''Offset of every shell in the spinor basis function spectrum

    Returns:
        list, each entry is the corresponding start id of spinor function

    Examples:

    >>> mol = gto.M(atom='O 0 0 0; C 0 0 1', basis='6-31g')
    >>> gto.ao_loc_2c(mol)
    [0, 2, 4, 6, 12, 18, 20, 22, 24, 30, 36]
    '''
    return moleintor.make_loc(mol._bas, 'spinor')

def time_reversal_map(mol):
    r'''The index to map the spinor functions and its time reversal counterpart.
    The returned indices have postive or negative values.  For the i-th basis function,
    if the returned j = idx[i] < 0, it means :math:`T|i\rangle = -|j\rangle`,
    otherwise :math:`T|i\rangle = |j\rangle`
    '''
    tao = []
    i = 0
    for b in mol._bas:
        l = b[ANG_OF]
        if b[KAPPA_OF] == 0:
            djs = (l * 2, l * 2 + 2)
        elif b[KAPPA_OF] > 0:
            djs = (l * 2,)
        else:
            djs = (l * 2 + 2,)
        if l % 2 == 0:
            for n in range(b[NCTR_OF]):
                for dj in djs:
                    for m in range(0, dj, 2):
                        tao.append(-(i + dj - m))
                        tao.append(  i + dj - m - 1)
                    i += dj
        else:
            for n in range(b[NCTR_OF]):
                for dj in djs:
                    for m in range(0, dj, 2):
                        tao.append(  i + dj - m)
                        tao.append(-(i + dj - m - 1))
                    i += dj
    return numpy.asarray(tao, dtype=numpy.int32)


CHECK_GEOM = getattr(__config__, 'gto_mole_check_geom', True)

def energy_nuc(mol, charges=None, coords=None):
    '''Compute nuclear repulsion energy (AU) or static Coulomb energy

    Returns
        float
    '''
    if charges is None: charges = mol.atom_charges()
    if len(charges) <= 1:
        return 0
    rr = inter_distance(mol, coords)
    rr[numpy.diag_indices_from(rr)] = 1e200
    if CHECK_GEOM and numpy.any(rr < 1e-5):
        for atm_idx in numpy.argwhere(rr<1e-5):
            logger.warn(mol, 'Atoms %s have the same coordinates', atm_idx)
        raise RuntimeError('Ill geometry')
    e = numpy.einsum('i,ij,j->', charges, 1./rr, charges) * .5
    return e

def inter_distance(mol, coords=None):
    '''
    Inter-particle distance array
    '''
    if coords is None: coords = mol.atom_coords()
    rr = numpy.linalg.norm(coords.reshape(-1,1,3) - coords, axis=2)
    rr[numpy.diag_indices_from(rr)] = 0
    return rr

def sph_labels(mol, fmt=True, base=BASE):
    '''Labels for spherical GTO functions

    Kwargs:
        fmt : str or bool
        if fmt is boolean, it controls whether to format the labels and the
        default format is "%d%3s %s%-4s".  if fmt is string, the string will
        be used as the print format.

    Returns:
        List of [(atom-id, symbol-str, nl-str, str-of-real-spherical-notation]
        or formatted strings based on the argument "fmt"

    Examples:

    >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='sto-3g')
    >>> gto.sph_labels(mol)
    [(0, 'H', '1s', ''), (1, 'Cl', '1s', ''), (1, 'Cl', '2s', ''), (1, 'Cl', '3s', ''),
     (1, 'Cl', '2p', 'x'), (1, 'Cl', '2p', 'y'), (1, 'Cl', '2p', 'z'), (1, 'Cl', '3p', 'x'),
     (1, 'Cl', '3p', 'y'), (1, 'Cl', '3p', 'z')]
    '''
    count = numpy.zeros((mol.natm, 9), dtype=int)
    label = []
    for ib in range(mol.nbas):
        ia = mol.bas_atom(ib)
        l = mol.bas_angular(ib)
        strl = param.ANGULAR[l]
        nc = mol.bas_nctr(ib)
        symb = mol.atom_symbol(ia)
        nelec_ecp = mol.atom_nelec_core(ia)
        if nelec_ecp == 0 or l > 3:
            shl_start = count[ia,l]+l+1
        else:
            coreshl = core_configuration(nelec_ecp)
            shl_start = coreshl[l]+count[ia,l]+l+1
        count[ia,l] += nc
        for n in range(shl_start, shl_start+nc):
            for m in range(-l, l+1):
                label.append((ia+base, symb, '%d%s' % (n, strl),
                              str(param.REAL_SPHERIC[l][l+m])))

    if isinstance(fmt, (str, unicode)):
        return [(fmt % x) for x in label]
    elif fmt:
        return ['%d %s %s%-4s' % x for x in label]
    else:
        return label
spheric_labels = spherical_labels = sph_labels

def cart_labels(mol, fmt=True, base=BASE):
    '''Labels of Cartesian GTO functions

    Kwargs:
        fmt : str or bool
        if fmt is boolean, it controls whether to format the labels and the
        default format is "%d%3s %s%-4s".  if fmt is string, the string will
        be used as the print format.

    Returns:
        List of [(atom-id, symbol-str, nl-str, str-of-xyz-notation)]
        or formatted strings based on the argument "fmt"
    '''
    cartxyz = []
    for l in range(max(mol._bas[:,ANG_OF])+1):
        xyz = []
        for x in range(l, -1, -1):
            for y in range(l-x, -1, -1):
                z = l-x-y
                xyz.append('x'*x + 'y'*y + 'z'*z)
        cartxyz.append(xyz)

    count = numpy.zeros((mol.natm, 9), dtype=int)
    label = []
    for ib in range(len(mol._bas)):
        ia = mol.bas_atom(ib)
        l = mol.bas_angular(ib)
        strl = param.ANGULAR[l]
        nc = mol.bas_nctr(ib)
        symb = mol.atom_symbol(ia)
        nelec_ecp = mol.atom_nelec_core(ia)
        if nelec_ecp == 0 or l > 3:
            shl_start = count[ia,l]+l+1
        else:
            coreshl = core_configuration(nelec_ecp)
            shl_start = coreshl[l]+count[ia,l]+l+1
        count[ia,l] += nc
        ncart = (l + 1) * (l + 2) // 2
        for n in range(shl_start, shl_start+nc):
            for m in range(ncart):
                label.append((ia+base, symb, '%d%s' % (n, strl), cartxyz[l][m]))

    if isinstance(fmt, (str, unicode)):
        return [(fmt % x) for x in label]
    elif fmt:
        return ['%d%3s %s%-4s' % x for x in label]
    else:
        return label

def ao_labels(mol, fmt=True, base=BASE):
    '''Labels of AO basis functions

    Kwargs:
        fmt : str or bool
            if fmt is boolean, it controls whether to format the labels and the
            default format is "%d%3s %s%-4s".  if fmt is string, the string will
            be used as the print format.

    Returns:
        List of [(atom-id, symbol-str, nl-str, str-of-AO-notation)]
        or formatted strings based on the argument "fmt"
    '''
    if mol.cart:
        return mol.cart_labels(fmt, base)
    else:
        return mol.sph_labels(fmt, base)

def spinor_labels(mol, fmt=True, base=BASE):
    '''
    Labels of spinor GTO functions
    '''
    count = numpy.zeros((mol.natm, 9), dtype=int)
    label = []
    for ib in range(mol.nbas):
        ia = mol.bas_atom(ib)
        l = mol.bas_angular(ib)
        kappa = mol.bas_kappa(ib)
        strl = param.ANGULAR[l]
        nc = mol.bas_nctr(ib)
        symb = mol.atom_symbol(ia)
        nelec_ecp = mol.atom_nelec_core(ia)
        if nelec_ecp == 0 or l > 3:
            shl_start = count[ia,l]+l+1
        else:
            coreshl = core_configuration(nelec_ecp)
            shl_start = coreshl[l]+count[ia,l]+l+1
        count[ia,l] += nc
        for n in range(shl_start, shl_start+nc):
            if kappa >= 0:
                for m in range(-l*2+1, l*2, 2):
                    label.append((ia+base, symb, '%d%s%d/2' % (n, strl, l*2-1),
                                  '%d/2'%m))
            if kappa <= 0:
                for m in range(-l*2-1, l*2+2, 2):
                    label.append((ia+base, symb, '%d%s%d/2' % (n, strl, l*2+1),
                                  '%d/2'%m))

    if isinstance(fmt, (str, unicode)):
        return [(fmt % x) for x in label]
    elif fmt:
        return ['%d %s %s,%-5s' % x for x in label]
    else:
        return label

def search_ao_label(mol, label):
    '''Find the index of the AO basis function based on the given ao_label

    Args:
        ao_label : string or a list of strings
            The regular expression pattern to match the orbital labels
            returned by mol.ao_labels()

    Returns:
        A list of index for the AOs that matches the given ao_label RE pattern

    Examples:

    >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='ccpvtz')
    >>> mol.search_ao_label('Cl.*p')
    [19 20 21 22 23 24 25 26 27 28 29 30]
    >>> mol.search_ao_label('Cl 2p')
    [19 20 21]
    >>> mol.search_ao_label(['Cl.*d', 'Cl 4p'])
    [25 26 27 31 32 33 34 35 36 37 38 39 40]
    '''
    return _aolabels2baslst(mol, label)

def _aolabels2baslst(mol, aolabels_or_baslst, base=BASE):
    if callable(aolabels_or_baslst):
        baslst = [i for i,x in enumerate(mol.ao_labels(base=base))
                  if aolabels_or_baslst(x)]
    elif isinstance(aolabels_or_baslst, str):
        aolabels = re.sub(' +', ' ', aolabels_or_baslst.strip(), count=1)
        aolabels = re.compile(aolabels)
        baslst = [i for i,s in enumerate(mol.ao_labels(base=base))
                  if re.search(aolabels, s)]
    elif len(aolabels_or_baslst) > 0 and isinstance(aolabels_or_baslst[0], str):
        aolabels = [re.compile(re.sub(' +', ' ', x.strip(), count=1))
                    for x in aolabels_or_baslst]
        baslst = [i for i,t in enumerate(mol.ao_labels(base=base))
                  if any(re.search(x, t) for x in aolabels)]
    else:
        baslst = [i-base for i in aolabels_or_baslst]
    return numpy.asarray(baslst, dtype=int)

def search_shell_id(mol, atm_id, l):
    '''Search the first basis/shell id (**not** the basis function id) which
    matches the given atom-id and angular momentum

    Args:
        atm_id : int
            atom id, 0-based
        l : int
            angular momentum

    Returns:
        basis id, 0-based.  If not found, return None

    Examples:

    >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='sto-3g')
    >>> mol.search_shell_id(1, 1) # Cl p shell
    4
    >>> mol.search_shell_id(1, 2) # Cl d shell
    None
    '''
    for ib in range(len(mol._bas)):
        ia = mol.bas_atom(ib)
        l1 = mol.bas_angular(ib)
        if ia == atm_id and l1 == l:
            return ib

def search_ao_nr(mol, atm_id, l, m, atmshell):
    '''Search the first basis function id (**not** the shell id) which matches
    the given atom-id, angular momentum magnetic angular momentum, principal shell.

    Args:
        atm_id : int
            atom id, 0-based
        l : int
            angular momentum
        m : int
            magnetic angular momentum
        atmshell : int
            principal quantum number

    Returns:
        basis function id, 0-based.  If not found, return None

    Examples:

    >>> mol = gto.M(atom='H 0 0 0; Cl 0 0 1', basis='sto-3g')
    >>> mol.search_ao_nr(1, 1, -1, 3) # Cl 3px
    7
    '''
    ibf = 0
    for ib in range(len(mol._bas)):
        ia = mol.bas_atom(ib)
        l1 = mol.bas_angular(ib)
        if mol.cart:
            degen = (l1 + 1) * (l1 + 2) // 2
        else:
            degen = l1 * 2 + 1
        nc = mol.bas_nctr(ib)
        if ia == atm_id and l1 == l:
            if atmshell > nc+l1:
                atmshell = atmshell - nc
            else:
                return ibf + (atmshell-l1-1)*degen + (l1+m)
        ibf += degen * nc
    raise RuntimeError('Required AO not found')

def search_ao_r(mol, atm_id, l, j, m, atmshell):
    raise RuntimeError('TODO')
#TODO:    ibf = 0
#TODO:    for ib in range(len(mol._bas)):
#TODO:        ia = mol.bas_atom(ib)
#TODO:        l1 = mol.bas_angular(ib)
#TODO:        nc = mol.bas_nctr(ib)
#TODO:        k = mol.bas_kappa(bas_id)
#TODO:        degen = len_spinor(l1, k)
#TODO:        if ia == atm_id and l1 == l and k == kappa:
#TODO:            if atmshell > nc+l1:
#TODO:                atmshell = atmshell - nc
#TODO:            else:
#TODO:                return ibf + (atmshell-l1-1)*degen + (degen+m)
#TODO:        ibf += degen

def offset_2c_by_atom(mol):
    '''2-component AO offset for each atom.  Return a list, each item
    of the list gives (start-shell-id, stop-shell-id, start-AO-id, stop-AO-id)
    '''
    return aoslice_by_atom(mol, mol.ao_loc_2c())

def aoslice_by_atom(mol, ao_loc=None):
    '''AO offsets for each atom.  Return a list, each item of the list gives
    (start-shell-id, stop-shell-id, start-AO-id, stop-AO-id)
    '''
    if ao_loc is None:
        ao_loc = mol.ao_loc_nr()

    aorange = numpy.empty((mol.natm,4), dtype=int)

    if mol.natm == 0:
        return aorange

    bas_atom = mol._bas[:,ATOM_OF]
    delimiter = numpy.where(bas_atom[0:-1] != bas_atom[1:])[0] + 1

    if mol.natm == len(delimiter) + 1:
        aorange[:,0] = shell_start = numpy.append(0, delimiter)
        aorange[:,1] = shell_end = numpy.append(delimiter, mol.nbas)

    else:  # Some atoms miss basis
        shell_start = numpy.empty(mol.natm, dtype=int)
        shell_start[:] = -1
        shell_start[0] = 0
        shell_start[bas_atom[0]] = 0
        shell_start[bas_atom[delimiter]] = delimiter

        shell_end = numpy.empty(mol.natm, dtype=int)
        shell_end[0] = 0
        shell_end[bas_atom[delimiter-1]] = delimiter
        shell_end[bas_atom[-1]] = mol.nbas

        for i in range(1, mol.natm):
            if shell_start[i] == -1:
                shell_start[i] = shell_end[i] = shell_end[i-1]

    aorange[:,0] = shell_start
    aorange[:,1] = shell_end
    aorange[:,2] = ao_loc[shell_start]
    aorange[:,3] = ao_loc[shell_end]
    return aorange
offset_nr_by_atom = aoslice_by_atom

def same_basis_set(mol1, mol2):
    '''Check whether two molecules use the same basis sets.
    The two molecules can have different geometry.
    '''
    atomtypes1 = atom_types(mol1._atom, mol1._basis)
    atomtypes2 = atom_types(mol2._atom, mol2._basis)
    if set(atomtypes1.keys()) != set(atomtypes2.keys()):
        return False
    for k in atomtypes1:
        if len(atomtypes1[k]) != len(atomtypes2[k]):
            return False
        elif mol1._basis.get(k, None) != mol2._basis.get(k, None):
            return False
    return True

def same_mol(mol1, mol2, tol=1e-5, cmp_basis=True, ignore_chiral=False):
    '''Compare the two molecules whether they have the same structure.

    Kwargs:
        tol : float
            In Bohr
        cmp_basis : bool
            Whether to compare basis functions for the two molecules
    '''
    from pyscf import symm

    if mol1._atom.__len__() != mol2._atom.__len__():
        return False

    chg1 = mol1._atm[:,CHARGE_OF]
    chg2 = mol2._atm[:,CHARGE_OF]
    if not numpy.all(numpy.sort(chg1) == numpy.sort(chg2)):
        return False

    if cmp_basis and not same_basis_set(mol1, mol2):
        return False

    def finger(mol, chgs, coord):
        center = numpy.einsum('z,zr->r', chgs, coord) / chgs.sum()
        im = inertia_moment(mol, chgs, coord)
        # Divid im by chgs.sum(), to normalize im. Otherwise the input tol may
        # not reflect the actual deviation.
        im /= chgs.sum()

        w, v = scipy.linalg.eigh(im)
        axes = v.T
        if numpy.linalg.det(axes) < 0:
            axes *= -1
        r = numpy.dot(coord-center, axes.T)
        return w, r

    coord1 = mol1.atom_coords()
    coord2 = mol2.atom_coords()
    w1, r1 = finger(mol1, chg1, coord1)
    w2, r2 = finger(mol2, chg2, coord2)
    if not (numpy.allclose(w1, w2, atol=tol)):
        return False

    rotate_xy  = numpy.array([[-1., 0., 0.],
                              [ 0.,-1., 0.],
                              [ 0., 0., 1.]])
    rotate_yz  = numpy.array([[ 1., 0., 0.],
                              [ 0.,-1., 0.],
                              [ 0., 0.,-1.]])
    rotate_zx  = numpy.array([[-1., 0., 0.],
                              [ 0., 1., 0.],
                              [ 0., 0.,-1.]])

    def inspect(z1, r1, z2, r2):
        place = int(-numpy.log10(tol)) - 1
        idx = symm.argsort_coords(r2, place)
        z2 = z2[idx]
        r2 = r2[idx]
        for rot in (1, rotate_xy, rotate_yz, rotate_zx):
            r1new = numpy.dot(r1, rot)
            idx = symm.argsort_coords(r1new, place)
            if (numpy.all(z1[idx] == z2) and
                numpy.allclose(r1new[idx], r2, atol=tol)):
                return True
        return False

    return (inspect(chg1, r1, chg2, r2) or
            (ignore_chiral and inspect(chg1, r1, chg2, -r2)))
is_same_mol = same_mol

def chiral_mol(mol1, mol2=None):
    '''Detect whether the given molelcule is chiral molecule or two molecules
    are chiral isomers.
    '''
    if mol2 is None:
        mol2 = mol1.copy()
        ptr_coord = mol2._atm[:,PTR_COORD]
        mol2._env[ptr_coord  ] *= -1
        mol2._env[ptr_coord+1] *= -1
        mol2._env[ptr_coord+2] *= -1
    return (not same_mol(mol1, mol2, ignore_chiral=False) and
            same_mol(mol1, mol2, ignore_chiral=True))

def inertia_moment(mol, mass=None, coords=None):
    if mass is None:
        mass = mol.atom_mass_list()
    if coords is None:
        coords = mol.atom_coords()
    mass_center = numpy.einsum('i,ij->j', mass, coords)/mass.sum()
    coords = coords - mass_center
    im = numpy.einsum('i,ij,ik->jk', mass, coords, coords)
    im = numpy.eye(3) * im.trace() - im
    return im

def atom_mass_list(mol, isotope_avg=False):
    '''A list of mass for all atoms in the molecule

    Kwargs:
        isotope_avg : boolean
            Whether to use the isotope average mass as the atomic mass
    '''
    if isotope_avg:
        mass_table = elements.MASSES
    else:
        mass_table = elements.ISOTOPE_MAIN

    nucprop = mol.nucprop
    if nucprop:
        mass = []
        for ia in range(mol.natm):
            z = mol.atom_charge(ia)
            symb = mol.atom_symbol(ia)
            stdsymb = _std_symbol(symb)
            if ia+1 in nucprop:
                prop = nucprop[ia+1]
            elif symb in nucprop:
                prop = nucprop[symb]
            else:
                prop = nucprop.get(stdsymb, {})

            mass.append(prop.get('mass', mass_table[z]))
    else:
        #mass = [mass_table[z] for z in mol.atom_charges()]
        mass = []
        for ia in range(mol.natm):
            z = charge(mol.atom_symbol(ia))
            mass.append(mass_table[z])

    return numpy.array(mass)

def condense_to_shell(mol, mat, compressor='max'):
    '''The given matrix is first partitioned to blocks, based on AO shell as
    delimiter. Then call compressor function to abstract each block.
    '''
    ao_loc = mol.ao_loc_nr()
    if callable(compressor):
        abstract = numpy.empty((mol.nbas, mol.nbas))
        for i, i0 in enumerate(ao_loc[:mol.nbas]):
            for j, j0 in enumerate(ao_loc[:mol.nbas]):
                abstract[i,j] = compressor(mat[i0:ao_loc[i+1],j0:ao_loc[j+1]])
    else:
        abstract = lib.condense(compressor, mat, ao_loc)
    return abstract


def tostring(mol, format='raw'):
    '''Convert molecular geometry to a string of the required format.

    Supported output formats:
        | raw: Each line is  <symobl> <x> <y> <z>
        | xyz: XYZ cartesian coordinates format
        | zmat: Z-matrix format
    '''
    format = format.lower()
    if format == 'xyz' or format == 'raw':
        coords = mol.atom_coords() * param.BOHR
        output = []
        if format == 'xyz':
            output.append('%d' % mol.natm)
            output.append('XYZ from PySCF')

        for i in range(mol.natm):
            symb = mol.atom_pure_symbol(i)
            x, y, z = coords[i]
            output.append('%-4s %14.5f %14.5f %14.5f' %
                          (symb, x, y, z))
        return '\n'.join(output)
    elif format == 'zmat':
        coords = mol.atom_coords() * param.BOHR
        zmat = cart2zmat(coords).splitlines()
        output = []
        for i, line in enumerate(zmat):
            symb = mol.atom_pure_symbol(i)
            output.append('%-4s   %s' % (symb, line))
        return '\n'.join(output)
    else:
        raise NotImplementedError

def tofile(mol, filename, format=None):
    '''Write molecular geometry to a file of the required format.

    Supported output formats:
        | raw: Each line is  <symobl> <x> <y> <z>
        | xyz: XYZ cartesian coordinates format
        | zmat: Z-matrix format
    '''
    if format is None:  # Guess format based on filename
        format = os.path.splitext(filename)[1][1:]
    string = tostring(mol, format)
    with open(filename, 'w') as f:
        f.write(string)
        f.write('\n')
    return string

def fromfile(filename, format=None):
    '''Read molecular geometry from a file
    (in testing)

    Supported formats:
        | raw: Each line is  <symobl> <x> <y> <z>
        | xyz: XYZ cartesian coordinates format
        | zmat: Z-matrix format
    '''
    if format is None:  # Guess format based on filename
        format = os.path.splitext(filename)[1][1:].lower()
        if format not in ('xyz', 'zmat', 'sdf'):
            format = 'raw'
    with open(filename, 'r') as f:
        return fromstring(f.read(), format)


def fromstring(string, format='xyz'):
    '''Convert the string of the specified format to internal format
    (in testing)

    Supported formats:
        | raw: Each line is  <symobl> <x> <y> <z>
        | xyz: XYZ cartesian coordinates format
        | zmat: Z-matrix format
    '''
    format = format.lower()
    if format == 'zmat':
        return from_zmatrix(string)
    elif format == 'xyz':
        dat = string.splitlines()
        natm = int(dat[0])
        return '\n'.join(dat[2:natm+2])
    elif format == 'sdf':
        raw = string.splitlines()
        natoms, nbonds = raw[3].split()[:2]
        atoms = []
        for line in raw[4:4+int(natoms)]:
            d = line.split()
            atoms.append('%s %s %s %s' % (d[3], d[0], d[1], d[2]))
        return '\n'.join(atoms)
    elif format == 'raw':
        return string
    else:
        raise NotImplementedError


#
# Mole class handles three layers: input, internal format, libcint arguments.
# The relationship of the three layers are, eg
#    .atom (input) <=>  ._atom (for python) <=> ._atm (for libcint)
#   .basis (input) <=> ._basis (for python) <=> ._bas (for libcint)
# input layer does not talk to libcint directly.  Data are held in python
# internal fomrat layer.  Most of methods defined in this class only operates
# 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 Mole(lib.StreamObject):
    '''Basic class to hold molecular structure 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)

        ** 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`
        _keys : a set of str
            Store the keys appeared in the module.  It is used to check misinput attributes

        ** 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()
    <class 'pyscf.gto.mole.Mole'> has no attributes Charge
    '''  # noqa: E501

    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)

    def __init__(self, **kwargs):
        self.output = None
        self.max_memory = param.MAX_MEMORY

        self.charge = 0
        self.spin = 0 # 2j == nelec_alpha - nelec_beta
        self.symmetry = False
        self.symmetry_subgroup = None
        self.cart = False

# Save inputs
# self.atom = [(symb/nuc_charge, (coord(Angstrom):0.,0.,0.)), ...]
        self.atom = []
# self.basis = {atom_type/nuc_charge: [l, kappa, (expnt, c_1, c_2,..),..]}
        self.basis = 'sto-3g'
# self.nucmod = {atom_symbol: nuclear_model, atom_id: nuc_mod}, atom_id is 1-based
        self.nucmod = {}
# self.ecp = {atom_symbol: [[l, (r_order, expnt, c),...]]}
        self.ecp = {}
# Nuclear property. self.nucprop = {atom_symbol: {key: value}}
        self.nucprop = {}
##################################################
# don't modify the following private variables, they are not input options
        self._atm = numpy.zeros((0,6), dtype=numpy.int32)
        self._bas = numpy.zeros((0,8), dtype=numpy.int32)
        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
        self.irrep_id = None
        self.irrep_name = None
        self._symm_orig = None
        self._symm_axes = None
        self._nelectron = None
        self._nao = None
        self._enuc = None
        self._atom = []
        self._basis = {}
        self._ecp = {}
        self._built = False

        # _pseudo is created to make the mol object consistenet with the mol
        # object converted from Cell.to_mol(). It is initialized in the
        # Cell.build() method only. Assigning _pseudo to mol object basically
        # has no effects. Mole.build() method does not have code to access the
        # contents of _pseudo.
        self._pseudo = {}

        keys = set(('verbose', 'unit', 'cart', 'incore_anyway'))
        self._keys = set(self.__dict__.keys()).union(keys)
        self.__dict__.update(kwargs)

    @property
    def natm(self):
        return len(self._atm)
    @property
    def nbas(self):
        return len(self._bas)

    @property
    def nelec(self):
        ne = self.nelectron
        nalpha = (ne + self.spin) // 2
        nbeta = nalpha - self.spin
        assert(nalpha >= 0 and nbeta >= 0)
        if nalpha + nbeta != ne:
            raise RuntimeError('Electron number %d and spin %d are not consistent\n'
                               'Note mol.spin = 2S = Nalpha - Nbeta, not 2S+1' %
                               (ne, self.spin))
        return nalpha, nbeta
    @nelec.setter
    def nelec(self, neleca_nelecb):
        neleca, nelecb = neleca_nelecb
        self._nelectron = neleca + nelecb
        self.spin = neleca - nelecb

    @property
    def nelectron(self):
        if self._nelectron is None:
            return self.tot_electrons()
        else:
            return self._nelectron
    @nelectron.setter
    def nelectron(self, n):
        self._nelectron = n

    @property
    def multiplicity(self):
        return self.spin + 1
    @multiplicity.setter
    def multiplicity(self, x):
        if x is None:
            self.spin = None
        else:
            self.spin = x - 1

    @property
    def ms(self):
        '''Spin quantum number. multiplicity = ms*2+1'''
        if self.spin % 2 == 0:
            return self.spin // 2
        else:
            return self.spin * .5
    @ms.setter
    def ms(self, x):
        if x is None:
            self.spin = None
        else:
            self.spin = int(round(2*x, 4))

    def __getattr__(self, key):
        '''To support accessing methods (mol.HF, mol.KS, mol.CCSD, mol.CASSCF, ...)
        from Mole object.
        '''
        if key[:2] == '__':  # Skip Python builtins
            raise AttributeError('Mole object has no attribute %s' % key)
        elif key in ('_ipython_canary_method_should_not_exist_',
                     '_repr_mimebundle_'):
            # https://github.com/mewwts/addict/issues/26
            # https://github.com/jupyter/notebook/issues/2014
            raise AttributeError

        # Import all available modules. Some methods are registered to other
        # classes/modules when importing modules in __all__.
        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)
            else:
                mf = dft.KS(self)
                xc = key.split('TD', 1)[1]
                if xc in dft.XC:
                    mf.xc = xc
                    key = 'TDDFT'
        else:
            mf = scf.HF(self)

        method = getattr(mf, key, None)
        if method is None:
            raise AttributeError('Mole object has no attribute %s' % key)

        # Initialize SCF object for post-SCF methods if applicable
        if self.nelectron != 0:
            mf.run()
        return method

# need "deepcopy" here because in shallow copy, _env may get new elements but
# with ptr_env unchanged
# def __copy__(self):
#        cls = self.__class__
#        newmol = cls.__new__(cls)
#        newmol = ...
# do not use __copy__ to aovid iteratively call copy.copy
    copy = copy

    pack = pack

    @classmethod
    @lib.with_doc(unpack.__doc__)
    def unpack(cls, moldic):
        return unpack(moldic)

    @lib.with_doc(unpack.__doc__)
    def unpack_(self, moldic):
        self.__dict__.update(moldic)
        return self

    dumps = dumps

    @classmethod
    @lib.with_doc(loads.__doc__)
    def loads(cls, molstr):
        return loads(molstr)

    @lib.with_doc(loads.__doc__)
    def loads_(self, molstr):
        self.__dict__.update(loads(molstr).__dict__)
        return self

    def build(self, dump_input=True, parse_arg=True,
              verbose=None, output=None, max_memory=None,
              atom=None, basis=None, unit=None, nucmod=None, ecp=None,
              charge=None, spin=0, symmetry=None, symmetry_subgroup=None,
              cart=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.

        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.

        '''
        if not DISABLE_GC:
            gc.collect()  # To release circular referred objects

        if isinstance(dump_input, (str, unicode)):
            sys.stderr.write('Assigning the first argument %s to mol.atom\n' %
                             dump_input)
            dump_input, atom = True, dump_input

        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 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 parse_arg:
            _update_from_cmdargs_(self)

        # 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):

            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)

            if self.output == '/dev/null':
                self.stdout = open(os.devnull, 'w')
            else:
                self.stdout = open(self.output, 'w')

        if self.verbose >= logger.WARN:
            self.check_sanity()

        self._atom = self.format_atom(self.atom, unit=self.unit)
        uniq_atoms = set([a[0] for a in self._atom])

        if isinstance(self.basis, (str, unicode, tuple, list)):
            # specify global basis for whole molecule
            _basis = dict(((a, self.basis) for a in uniq_atoms))
        elif 'default' in self.basis:
            default_basis = self.basis['default']
            _basis = dict(((a, default_basis) for a in uniq_atoms))
            _basis.update(self.basis)
            del(_basis['default'])
        else:
            _basis = self.basis
        self._basis = self.format_basis(_basis)

# TODO: Consider ECP info in point group symmetry initialization
        if self.ecp:
            # Unless explicitly input, ECP should not be assigned to ghost atoms
            if isinstance(self.ecp, (str, unicode)):
                _ecp = dict([(a, str(self.ecp))
                             for a in uniq_atoms if not is_ghost_atom(a)])
            elif 'default' in self.ecp:
                default_ecp = self.ecp['default']
                _ecp = dict(((a, default_ecp)
                             for a in uniq_atoms if not is_ghost_atom(a)))
                _ecp.update(self.ecp)
                del(_ecp['default'])
            else:
                _ecp = self.ecp
            self._ecp = self.format_ecp(_ecp)

        env = self._env[:PTR_ENV_START]
        self._atm, self._bas, self._env = \
                self.make_env(self._atom, self._basis, env, self.nucmod,
                              self.nucprop)
        self._atm, self._ecpbas, self._env = \
                self.make_ecp_env(self._atm, self._ecp, self._env)

        if self.spin is None:
            self.spin = self.nelectron % 2
        else:
            # Access self.nelec in which the code checks whether the spin and
            # number of electrons are consistent.
            self.nelec

        if self.symmetry:
            self._build_symmetry()

        if dump_input and not self._built and self.verbose > logger.NOTE:
            self.dump_input()

        logger.debug3(self, 'arg.atm = %s', str(self._atm))
        logger.debug3(self, 'arg.bas = %s', str(self._bas))
        logger.debug3(self, 'arg.env = %s', str(self._env))
        logger.debug3(self, 'ecpbas  = %s', str(self._ecpbas))

        self._built = True
        return self
    kernel = build

    def _build_symmetry(self):
        '''
        Update symmetry related attributes: topgroup, groupname, _symm_orig,
        _symm_axes, irrep_id, irrep_name, symm_orb
        '''
        from pyscf import symm
        self.topgroup, orig, axes = symm.detect_symm(self._atom, self._basis)

        if isinstance(self.symmetry, (str, unicode)):
            self.symmetry = str(symm.std_symb(self.symmetry))
            try:
                self.groupname, axes = symm.as_subgroup(self.topgroup, axes,
                                                        self.symmetry)
            except PointGroupSymmetryError:
                raise PointGroupSymmetryError(
                    'Unable to identify input symmetry %s. Try symmetry="%s"' %
                    (self.symmetry, self.topgroup))
        else:
            self.groupname, axes = symm.as_subgroup(self.topgroup, axes,
                                                    self.symmetry_subgroup)
        self._symm_orig = orig
        self._symm_axes = axes

        if self.cart and self.groupname in ('Dooh', 'Coov', 'SO3'):
            if self.groupname == 'Coov':
                self.groupname, lgroup = 'C2v', self.groupname
            else:
                self.groupname, lgroup = 'D2h', self.groupname
            logger.warn(self, 'This version does not support symmetry %s '
                        'for cartesian GTO basis. Its subgroup %s is used',
                        lgroup, self.groupname)

        self.symm_orb, self.irrep_id = \
                symm.symm_adapted_basis(self, self.groupname, orig, axes)
        self.irrep_name = [symm.irrep_id2name(self.groupname, ir)
                           for ir in self.irrep_id]
        return self

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

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

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

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

    @lib.with_doc(expand_etbs.__doc__)
    def expand_etbs(self, etbs):
        return expand_etbs(etbs)
    etbs = expand_etbs

    @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__
        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())
        import pyscf
        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')

        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]])
            self.stdout.write('[INPUT]%3d %-4s %16.12f %16.12f %16.12f AA  '
                              '%16.12f %16.12f %16.12f Bohr\n'
                              % ((ia+1, _symbol(atom[0])) + coorda + coordb))
        if self.nucmod:
            if isinstance(self.nucmod, (int, str, unicode,
                                        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)
                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):
        '''Retuen 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):
        '''Retuen 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):
        '''Retuen 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):
        '''Retuen 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):
        '''Retuen 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):
        '''Retuen 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):
        '''Retuen 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:
            self._env[AS_RINV_ORIG_ATOM] = atm_id  # required by ecp gradients
            rinv0 = self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3].copy()
            zeta0 = self._env[PTR_RINV_ZETA].copy()

            def set_rinv(z, r):
                self._env[PTR_RINV_ZETA] = z
                self._env[PTR_RINV_ORIG:PTR_RINV_ORIG+3] = r
            return self._TemporaryMoleContext(set_rinv, (zeta,rinv), (zeta0,rinv0))
    with_rinv_as_nucleus = with_rinv_at_nucleus  # For backward compatibility

    def with_integral_screen(self, threshold):
        '''Retuen a temporary mol context which has the rquired integral
        screen threshold
        '''
        expcutoff0 = self._env[PTR_EXPCUTOFF]
        expcutoff = abs(numpy.log(threshold))
        def set_cutoff(cut):
            self._env[PTR_EXPCUTOFF] = cut
        return self._TemporaryMoleContext(set_cutoff, (expcutoff,), (expcutoff0,))

    def set_geom_(self, atoms_or_coords, unit=None, symmetry=None,
                  inplace=True):
        '''Update geometry
        '''
        import copy
        if inplace:
            mol = self
        else:
            mol = copy.copy(self)
            mol._env = mol._env.copy()
        if unit is None:
            unit = mol.unit
        else:
            mol.unit = unit
        if symmetry is None:
            symmetry = mol.symmetry

        if isinstance(atoms_or_coords, numpy.ndarray):
            mol.atom = list(zip([x[0] for x in mol._atom],
                                atoms_or_coords.tolist()))
        else:
            mol.atom = atoms_or_coords

        if isinstance(atoms_or_coords, numpy.ndarray) and not symmetry:
            mol._atom = mol.atom

            if isinstance(unit, (str, unicode)):
                if unit.upper().startswith(('B', 'AU')):
                    unit = 1.
                else: #unit[:3].upper() == 'ANG':
                    unit = 1./param.BOHR
            else:
                unit = 1./unit
            ptr = mol._atm[:,PTR_COORD]
            mol._env[ptr+0] = unit * atoms_or_coords[:,0]
            mol._env[ptr+1] = unit * atoms_or_coords[:,1]
            mol._env[ptr+2] = unit * atoms_or_coords[:,2]
        else:
            mol.symmetry = symmetry
            mol.build(False, False)

        if mol.verbose >= logger.INFO:
            logger.info(mol, 'New geometry')
            for ia, atom in enumerate(mol._atom):
                coorda = tuple([x * param.BOHR for x in atom[1]])
                coordb = tuple([x for x in atom[1]])
                coords = coorda + coordb
                logger.info(mol, ' %3d %-4s %16.12f %16.12f %16.12f AA  '
                            '%16.12f %16.12f %16.12f Bohr\n',
                            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

    def has_ecp_soc(self):
        '''Whether spin-orbit coupling is enabled in ECP.'''
        return (len(self._ecpbas) > 0 and
                numpy.any(self._ecpbas[:,SO_TYPE_OF] == 1))


#######################################################
#NOTE: atm_id or bas_id start from 0
    def atom_symbol(self, atm_id):
        r'''For the given atom id, return the input symbol (without striping special characters)

        Args:
            atm_id : int
                0-based

        Examples:

        >>> mol.build(atom='H^2 0 0 0; H 0 0 1.1')
        >>> mol.atom_symbol(0)
        H^2
        '''
        return _symbol(self._atom[atm_id][0])

    def atom_pure_symbol(self, atm_id):
        r'''For the given atom id, return the standard symbol (striping special characters)

        Args:
            atm_id : int
                0-based

        Examples:

        >>> mol.build(atom='H^2 0 0 0; H 0 0 1.1')
        >>> mol.atom_symbol(0)
        H
        '''
        return _std_symbol(self._atom[atm_id][0])

    @property
    def elements(self):
        '''A list of elements in the molecule'''
        return [self.atom_pure_symbol(i) for i in range(self.natm)]

    def atom_charge(self, atm_id):
        r'''Nuclear effective charge of the given atom id
        Note "atom_charge /= charge(atom_symbol)" when ECP is enabled.
        Number of electrons screened by ECP can be obtained by charge(atom_symbol)-atom_charge

        Args:
            atm_id : int
                0-based

        Examples:

        >>> mol.build(atom='H 0 0 0; Cl 0 0 1.1')
        >>> mol.atom_charge(1)
        17
        '''
        if self._atm[atm_id,NUC_MOD_OF] != NUC_FRAC_CHARGE:
            # regular QM atoms
            return self._atm[atm_id,CHARGE_OF]
        else:
            # MM atoms with fractional charges
            return self._env[self._atm[atm_id,PTR_FRAC_CHARGE]]

    def atom_charges(self):
        '''np.asarray([mol.atom_charge(i) for i in range(mol.natm)])'''
        z = self._atm[:,CHARGE_OF]
        if numpy.any(self._atm[:,NUC_MOD_OF] == NUC_FRAC_CHARGE):
            # Create the integer nuclear charges first then replace the MM
            # particles with the MM charges that saved in _env[PTR_FRAC_CHARGE]
            z = numpy.array(z, dtype=numpy.double)
            idx = self._atm[:,NUC_MOD_OF] == NUC_FRAC_CHARGE
            # MM fractional charges can be positive or negative
            z[idx] = self._env[self._atm[idx,PTR_FRAC_CHARGE]]
        return z

    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 unit[:3].upper() == 'ANG':
            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 unit[:3].upper() == 'ANG':
            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 _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

    nao_nr = nao_nr
    nao_2c = nao_2c
    nao_cart = nao_cart

    nao_nr_range = nao_nr_range
    nao_2c_range = nao_2c_range

    ao_loc_nr = ao_loc_nr
    ao_loc_2c = ao_loc_2c

    @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

    ao_loc = property(ao_loc_nr)

    tmap = time_reversal_map = time_reversal_map

    inertia_moment = inertia_moment

    tostring = tostring
    tofile = tofile

    def fromstring(self, string, format='xyz'):
        '''Update the Mole object based on the input geometry string'''
        atom = self.format_atom(fromstring(string, format), unit=1)
        self.set_geom_(atom, unit='Angstrom', inplace=True)
        if format == 'sdf' and 'M  CHG' in string:
            raise NotImplementedError
            #FIXME self.charge = 0
        return self

    def fromfile(self, filename, format=None):
        '''Update the Mole object based on the input geometry file'''
        atom = self.format_atom(fromfile(filename, format), unit=1)
        self.set_geom_(atom, unit='Angstrom', inplace=True)
        if format == 'sdf':
            raise NotImplementedError
        return self

    def intor(self, intor, comp=None, hermi=0, aosym='s1', out=None,
              shls_slice=None, grids=None):
        '''Integral generator.

        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

        Kwargs:
            comp : int
                Components of the integrals, e.g. int1e_ipovlp_sph has 3 components.
            hermi : int
                Symmetry of the integrals

                | 0 : no symmetry assumed (default)
                | 1 : hermitian
                | 2 : anti-hermitian

            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

        Examples:

        >>> mol.build(atom='H 0 0 0; H 0 0 1.1', basis='sto-3g')
        >>> mol.intor('int1e_ipnuc_sph', comp=3) # <nabla i | V_nuc | j>
        [[[ 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 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

    def intor_symmetric(self, intor, comp=None, grids=None):
        '''One-electron integral generator. The integrals are assumed to be hermitian

        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

        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]]
        '''
        return self.intor(intor, comp, 1, aosym='s4', grids=grids)

    def intor_asymmetric(self, intor, comp=None, grids=None):
        '''One-electron integral generator. The integrals are assumed to be anti-hermitian

        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 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

        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]]
        '''
        return self.intor(intor, comp, 2, aosym='a4', grids=grids)

    @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)

    eval_ao = eval_gto = eval_gto

    energy_nuc = energy_nuc
    def get_enuc(self):
        return self.energy_nuc()

    sph_labels = spheric_labels = sph_labels
    cart_labels = cart_labels
    ao_labels = ao_labels

    spinor_labels = spinor_labels

    search_ao_label = search_ao_label

    def search_shell_id(self, atm_id, l):
        return search_shell_id(self, atm_id, l)

    search_ao_nr = search_ao_nr
    search_ao_r = search_ao_r

    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

    condense_to_shell = condense_to_shell

    to_uncontracted_cartesian_basis = to_uncontracted_cartesian_basis

    __add__ = conc_mol

    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.

        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
        '''
        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)

    def sph2spinor_coeff(self):
        '''Transformation matrix that transforms real-spherical GTOs to spinor
        GTOs for all basis functions

        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):
            return lib.StreamObject.apply(self, fn, *args, **kwargs)
        elif isinstance(fn, (str, unicode)):
            method = getattr(self, fn.upper())
            return method(*args, **kwargs)
        else:
            raise TypeError('First argument of .apply method must be a '
                            'function/class or a name (string) of a method.')

    def ao2mo(self, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e',
              **kwargs):
        '''Integral transformation for arbitrary orbitals and arbitrary
        integrals.  See more detalied documentation in func:`ao2mo.kernel`.

        Args:
            mo_coeffs (an np array or a list of arrays) : A matrix of orbital
                coefficients if it is a numpy ndarray, or four sets of orbital
                coefficients, corresponding to the four indices of (ij|kl).

        Kwargs:
            erifile (str or h5py File or h5py Group object) : The file/object
                to store the transformed integrals.  If not given, the return
                value is an array (in memory) of the transformed integrals.
            dataname : str
                *Note* this argument is effective if erifile is given.
                The dataset name in the erifile (ref the hierarchy of HDF5 format
                http://www.hdfgroup.org/HDF5/doc1.6/UG/09_Groups.html).  By assigning
                different dataname, the existed integral file can be reused.  If
                the erifile contains the specified dataname, the old integrals
                will be replaced by the new one under the key dataname.
            intor (str) : integral name Name of the 2-electron integral.  Ref
                to :func:`getints_by_shell`
                for the complete list of available 2-electron integral names

        Returns:
            An array of transformed integrals if erifile is not given.
            Otherwise, return the file/fileobject if erifile is assigned.


        Examples:

        >>> import pyscf
        >>> mol = pyscf.M(atom='O 0 0 0; H 0 1 0; H 0 0 1', basis='sto3g')
        >>> mo1 = numpy.random.random((mol.nao_nr(), 10))
        >>> mo2 = numpy.random.random((mol.nao_nr(), 8))

        >>> eri1 = mol.ao2mo(mo1)
        >>> print(eri1.shape)
        (55, 55)

        >>> eri1 = mol.ao2mo(mo1, compact=False)
        >>> print(eri1.shape)
        (100, 100)

        >>> eri1 = mol.ao2mo(eri, (mo1,mo2,mo2,mo2))
        >>> print(eri1.shape)
        (80, 36)

        >>> eri1 = mol.ao2mo(eri, (mo1,mo2,mo2,mo2), erifile='water.h5')
        '''
        from pyscf import ao2mo
        return ao2mo.kernel(self, mo_coeffs, erifile, dataname, intor, **kwargs)

    @contextlib.contextmanager
    def _TemporaryMoleContext(self, method, args, args_bak):
        '''Almost every method depends on the Mole environment. Ensure the
        modification in temporary environment being thread safe
        '''
        haslock = hasattr(self, '_lock')
        if not haslock:
            self._lock = threading.RLock()

        with self._lock:
            method(*args)
            try:
                yield
            finally:
                method(*args_bak)
                if not haslock:
                    del self._lock

def _parse_nuc_mod(str_or_int_or_fn):
    nucmod = NUC_POINT
    if callable(str_or_int_or_fn):
        nucmod = str_or_int_or_fn
    elif (isinstance(str_or_int_or_fn, (str, unicode)) and
          str_or_int_or_fn[0].upper() == 'G'): # 'gauss_nuc'
        nucmod = NUC_GAUSS
    elif str_or_int_or_fn != 0:
        nucmod = NUC_GAUSS
    return nucmod

def _update_from_cmdargs_(mol):
    try:
        # Detect whether in Ipython shell
        __IPYTHON__  # noqa:
        return
    except Exception:
        pass

    if not mol._built: # parse cmdline args only once
        opts = cmd_args.cmd_args()

        if opts.verbose:
            mol.verbose = opts.verbose
        if opts.max_memory:
            mol.max_memory = opts.max_memory

        if opts.output:
            mol.output = opts.output


def from_zmatrix(atomstr):
    '''>>> a = """H
    H 1 2.67247631453057
    H 1 4.22555607338457 2 50.7684795164077
    H 1 2.90305235726773 2 79.3904651036893 3 6.20854462618583"""
    >>> for x in zmat2cart(a): print(x)
    ['H', array([ 0.,  0.,  0.])]
    ['H', array([ 2.67247631,  0.        ,  0.        ])]
    ['H', array([ 2.67247631,  0.        ,  3.27310166])]
    ['H', array([ 0.53449526,  0.30859098,  2.83668811])]
    '''
    from pyscf.symm import rotation_mat
    atomstr = atomstr.replace(';','\n').replace(',',' ')
    symbols = []
    coord = []
    min_items_per_line = 1
    for line_id, line in enumerate(atomstr.splitlines()):
        line = line.strip()
        if line and line[0] != '#':
            rawd = line.split()
            if len(rawd) < min_items_per_line:
                raise ValueError('Zmatrix format error at L%d %s' % (line_id, line))

            symbols.append(rawd[0])
            if len(rawd) < 3:
                coord.append(numpy.zeros(3))
                min_items_per_line = 3
            elif len(rawd) == 3:
                if DISABLE_EVAL:
                    coord.append(numpy.array((float(rawd[2]), 0, 0)))
                else:
                    coord.append(numpy.array((eval(rawd[2]), 0, 0)))
                min_items_per_line = 5
            elif len(rawd) == 5:
                if DISABLE_EVAL:
                    vals = rawd[1:]
                else:
                    vals = eval(','.join(rawd[1:]))
                bonda = int(vals[0]) - 1
                bond  = float(vals[1])
                anga  = int(vals[2]) - 1
                ang   = float(vals[3])/180*numpy.pi
                assert(ang >= 0)
                v1 = coord[anga] - coord[bonda]
                if not numpy.allclose(v1[:2], 0):
                    vecn = numpy.cross(v1, numpy.array((0.,0.,1.)))
                else: # on z
                    vecn = numpy.array((0.,0.,1.))
                rmat = rotation_mat(vecn, ang)
                c = numpy.dot(rmat, v1) * (bond/numpy.linalg.norm(v1))
                coord.append(coord[bonda]+c)
                min_items_per_line = 7
            else:
                if DISABLE_EVAL:
                    vals = rawd[1:]
                else:
                    vals = eval(','.join(rawd[1:]))
                bonda = int(vals[0]) - 1
                bond  = float(vals[1])
                anga  = int(vals[2]) - 1
                ang   = float(vals[3])/180*numpy.pi
                assert(ang >= 0 and ang <= numpy.pi)
                v1 = coord[anga] - coord[bonda]
                v1 /= numpy.linalg.norm(v1)
                if ang < 1e-7:
                    c = v1 * bond
                elif numpy.pi-ang < 1e-7:
                    c = -v1 * bond
                else:
                    diha  = int(vals[4]) - 1
                    dih   = float(vals[5])/180*numpy.pi
                    v2 = coord[diha] - coord[anga]
                    vecn = numpy.cross(v2, -v1)
                    vecn_norm = numpy.linalg.norm(vecn)
                    if vecn_norm < 1e-7:
                        if not numpy.allclose(v1[:2], 0):
                            vecn = numpy.cross(v1, numpy.array((0.,0.,1.)))
                        else: # on z
                            vecn = numpy.array((0.,0.,1.))
                        rmat = rotation_mat(vecn, ang)
                        c = numpy.dot(rmat, v1) * bond
                    else:
                        rmat = rotation_mat(v1, -dih)
                        vecn = numpy.dot(rmat, vecn) / vecn_norm
                        rmat = rotation_mat(vecn, ang)
                        c = numpy.dot(rmat, v1) * bond
                coord.append(coord[bonda]+c)
    atoms = list(zip([_atom_symbol(x) for x in symbols], coord))
    return atoms
zmat2cart = zmat = from_zmatrix

def cart2zmat(coord):
    '''>>> c = numpy.array((
    (0.000000000000,  1.889726124565,  0.000000000000),
    (0.000000000000,  0.000000000000, -1.889726124565),
    (1.889726124565, -1.889726124565,  0.000000000000),
    (1.889726124565,  0.000000000000,  1.133835674739)))
    >>> print(cart2zmat(c))
    1
    1 2.67247631453057
    1 4.22555607338457 2 50.7684795164077
    1 2.90305235726773 2 79.3904651036893 3 6.20854462618583
    '''
    zstr = []
    zstr.append('1')
    if len(coord) > 1:
        r1 = coord[1] - coord[0]
        nr1 = numpy.linalg.norm(r1)
        zstr.append('1 %.15g' % nr1)
    if len(coord) > 2:
        r2 = coord[2] - coord[0]
        nr2 = numpy.linalg.norm(r2)
        a = numpy.arccos(numpy.dot(r1,r2)/(nr1*nr2))
        zstr.append('1 %.15g 2 %.15g' % (nr2, a*180/numpy.pi))
    if len(coord) > 3:
        o0, o1, o2 = coord[:3]
        p0, p1, p2 = 1, 2, 3
        for k, c in enumerate(coord[3:]):
            r0 = c - o0
            nr0 = numpy.linalg.norm(r0)
            r1 = o1 - o0
            nr1 = numpy.linalg.norm(r1)
            a1 = numpy.arccos(numpy.dot(r0,r1)/(nr0*nr1))
            b0 = numpy.cross(r0, r1)
            nb0 = numpy.linalg.norm(b0)

            if abs(nb0) < 1e-7: # o0, o1, c in line
                a2 = 0
                zstr.append('%d %.15g %d %.15g %d %.15g' %
                            (p0, nr0, p1, a1*180/numpy.pi, p2, a2))
            else:
                b1 = numpy.cross(o2-o0, r1)
                nb1 = numpy.linalg.norm(b1)

                if abs(nb1) < 1e-7:  # o0 o1 o2 in line
                    a2 = 0
                    zstr.append('%d %.15g %d %.15g %d %.15g' %
                                (p0, nr0, p1, a1*180/numpy.pi, p2, a2))
                    o2 = c
                    p2 = 4 + k
                else:
                    if numpy.dot(numpy.cross(b1, b0), r1) < 0:
                        a2 = numpy.arccos(numpy.dot(b1, b0) / (nb0*nb1))
                    else:
                        a2 =-numpy.arccos(numpy.dot(b1, b0) / (nb0*nb1))
                    zstr.append('%d %.15g %d %.15g %d %.15g' %
                                (p0, nr0, p1, a1*180/numpy.pi, p2, a2*180/numpy.pi))

    return '\n'.join(zstr)

def dyall_nuc_mod(nuc_charge, nucprop={}):
    ''' Generate the nuclear charge distribution parameter zeta
    rho(r) = nuc_charge * Norm * exp(-zeta * r^2)

    Ref. L. Visscher and K. Dyall, At. Data Nucl. Data Tables, 67, 207 (1997)
    '''
    mass = nucprop.get('mass', elements.ISOTOPE_MAIN[nuc_charge])
    r = (0.836 * mass**(1./3) + 0.570) / 52917.7249
    zeta = 1.5 / (r**2)
    return zeta

def filatov_nuc_mod(nuc_charge, nucprop={}):
    ''' Generate the nuclear charge distribution parameter zeta
    rho(r) = nuc_charge * Norm * exp(-zeta * r^2)

    Ref. M. Filatov and D. Cremer, Theor. Chem. Acc. 108, 168 (2002)
         M. Filatov and D. Cremer, Chem. Phys. Lett. 351, 259 (2002)
    '''
    c = param.LIGHT_SPEED
    nuc_charge = charge(nuc_charge)
    r = (-0.263188*nuc_charge + 106.016974 + 138.985999/nuc_charge) / c**2
    zeta = 1 / (r**2)
    return zeta

def fakemol_for_charges(coords, expnt=1e16):
    '''Construct a fake Mole object that holds the charges on the given
    coordinates (coords).  The shape of the charge can be a normal
    distribution with the Gaussian exponent (expnt).
    '''
    nbas = coords.shape[0]
    fakeatm = numpy.zeros((nbas,ATM_SLOTS), dtype=numpy.int32)
    fakebas = numpy.zeros((nbas,BAS_SLOTS), dtype=numpy.int32)
    fakeenv = [0] * PTR_ENV_START
    ptr = PTR_ENV_START
    fakeatm[:,PTR_COORD] = numpy.arange(ptr, ptr+nbas*3, 3)
    fakeenv.append(coords.ravel())
    ptr += nbas*3
    fakebas[:,ATOM_OF] = numpy.arange(nbas)
    fakebas[:,NPRIM_OF] = 1
    fakebas[:,NCTR_OF] = 1
# approximate point charge with gaussian distribution exp(-1e16*r^2)
    fakebas[:,PTR_EXP] = ptr
    fakebas[:,PTR_COEFF] = ptr+1
    fakeenv.append([expnt, 1/(2*numpy.sqrt(numpy.pi)*gaussian_int(2,expnt))])
    ptr += 2
    fakemol = Mole()
    fakemol._atm = fakeatm
    fakemol._bas = fakebas
    fakemol._env = numpy.hstack(fakeenv)
    fakemol._built = True
    return fakemol

del(BASE)