Source code for pyscf.lo.iao

#!/usr/bin/env python
# Copyright 2014-2018 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
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.
# Author: Qiming Sun <>
#         Paul J. Robinson <>

Intrinsic Atomic Orbitals
ref. JCTC, 9, 4834

from functools import reduce
import numpy
import scipy.linalg
from pyscf.lib import logger
from pyscf import gto
from pyscf import scf
from pyscf import __config__

# Alternately, use ANO for minao
# orthogonalize iao by orth.lowdin(c.T*mol.intor(ovlp)*c)
MINAO = getattr(__config__, 'lo_iao_minao', 'minao')

[docs]def iao(mol, orbocc, minao=MINAO): '''Intrinsic Atomic Orbitals. [Ref. JCTC, 9, 4834] Args: mol : the molecule or cell object orbocc : 2D array occupied orbitals Returns: non-orthogonal IAO orbitals. Orthogonalize them as C (C^T S C)^{-1/2}, eg using :func:`orth.lowdin` >>> orbocc = mf.mo_coeff[:,mf.mo_occ>0] >>> c = iao(mol, orbocc) >>>, orth.lowdin(reduce(, (c.T,s,c)))) ''' if mol.has_ecp(): logger.warn(mol, 'ECP/PP is used. MINAO is not a good reference AO basis in IAO.') pmol = reference_mol(mol, minao) # For PBC, we must use the pbc code for evaluating the integrals lest the # pbc conditions be ignored. # DO NOT import pbcgto early and check whether mol is a cell object. # "from pyscf.pbc import gto as pbcgto and isinstance(mol, pbcgto.Cell)" # The code should work even pbc module is not availabe. if hasattr(mol, 'pbc_intor'): # cell object has pbc_intor method from pyscf.pbc import gto as pbcgto s1 = mol.pbc_intor('int1e_ovlp', hermi=1) s2 = pmol.pbc_intor('int1e_ovlp', hermi=1) s12 = pbcgto.cell.intor_cross('int1e_ovlp', mol, pmol) else: #s1 is the one electron overlap integrals (coulomb integrals) s1 = mol.intor_symmetric('int1e_ovlp') #s2 is the same as s1 except in minao s2 = pmol.intor_symmetric('int1e_ovlp') #overlap integrals of the two molecules s12 = gto.mole.intor_cross('int1e_ovlp', mol, pmol) #transpose of overlap s21 = s12.T s1cd = scipy.linalg.cho_factor(s1) s2cd = scipy.linalg.cho_factor(s2) p12 = scipy.linalg.cho_solve(s1cd, s12) ctild = scipy.linalg.cho_solve(s2cd,, orbocc)) ctild = scipy.linalg.cho_solve(s1cd,, ctild)) ccs1 = reduce(, (orbocc, orbocc.conj().T, s1)) ccs2 = reduce(, (ctild, ctild.conj().T, s1)) #a is the set of IAOs in the original basis a = (p12 + reduce(, (ccs1, ccs2, p12)) * 2 -, p12) -, p12)) return a
[docs]def reference_mol(mol, minao=MINAO): '''Create a molecule which uses reference minimal basis''' pmol = mol.copy() if hasattr(pmol, 'rcut'): pmol.rcut = None, False, basis=minao) return pmol
[docs]def fast_iao_mullikan_pop(mol, dm, iaos, verbose=logger.DEBUG): ''' Args: mol : the molecule or cell object iaos : 2D array (orthogonal or non-orthogonal) IAO orbitals Returns: mullikan population analysis in the basis IAO ''' pmol = reference_mol(mol) if hasattr(mol, 'pbc_intor'): # whether mol object is a cell ovlpS = mol.pbc_intor('int1e_ovlp') else: ovlpS = mol.intor_symmetric('int1e_ovlp') # Transform DM in big basis to IAO basis # |IAO> = |big> C # DM_{IAO} = C^{-1} DM (C^{-1})^T = S_{IAO}^{-1} C^T S DM S C S_{IAO}^{-1} cs =, ovlpS) s_iao =, iaos) iao_inv = numpy.linalg.solve(s_iao, cs) if isinstance(dm, numpy.ndarray) and dm.ndim == 2: dm = reduce(, (iao_inv, dm, iao_inv.conj().T)) return scf.hf.mulliken_pop(pmol, dm, s_iao, verbose) else: dm = [reduce(, (iao_inv, dm[0], iao_inv.conj().T)), reduce(, (iao_inv, dm[1], iao_inv.conj().T))] return scf.uhf.mulliken_pop(pmol, dm, s_iao, verbose)