pyscf · sunqm · Apr 15, 2020 · Mar 23, 2020
diff --git a/pyscf/dft/numint.py b/pyscf/dft/numint.py
@@ -22,6 +22,7 @@
 import scipy.linalg
 from pyscf import lib
 from pyscf.lib import logger
+from pyscf.dft.sap import sap_effective_charge
 try:
     from pyscf.dft import libxc
 except (ImportError, OSError):
@@ -717,6 +718,68 @@ def nr_vxc(mol, grids, xc_code, dms, spin=0, relativity=0, hermi=0,
     return ni.nr_vxc(mol, grids, xc_code, dms, spin, relativity,
                      hermi, max_memory, verbose)
 
+def nr_sap_vxc(ni, mol, grids, max_memory=2000, verbose=None):
+    '''Calculate superposition of atomic potentials matrix on given meshgrids.
+
+    Args:
+        ni : an instance of :class:`NumInt`
+
+        mol : an instance of :class:`Mole`
+
+        grids : an instance of :class:`Grids`
+            grids.coords and grids.weights are needed for coordinates and weights of meshgrids.
+
+    Kwargs:
+        max_memory : int or float
+            The maximum size of cache to use (in MB).
+
+    Returns:
+        vmat is the XC potential matrix in 2D array of shape (nao,nao)
+        where nao is the number of AO functions.
+
+    Examples:
+    >>> import numpy
+    >>> from pyscf import gto, dft
+    >>> mol = gto.M(atom='H 0 0 0; H 0 0 1.1')
+    >>> grids = dft.gen_grid.Grids(mol)
+    >>> ni = dft.numint.NumInt()
+    >>> vsap = ni.nr_sap(mol, grids)
+    '''
+    shls_slice = (0, mol.nbas)
+    ao_loc = mol.ao_loc_nr()
+    nao = mol.nao
+    vmat = numpy.zeros((nao,nao))
+    aow = None
+    vxcw = None
+    ao_deriv = 0
+
+    atom_coords = mol.atom_coords()
+    atom_charges = mol.atom_charges()
+    eps = numpy.finfo(float).eps
+
+    for ao, mask, weight, coords in ni.block_loop(mol, grids, nao, ao_deriv, max_memory):
+        aow = numpy.ndarray(ao.shape, order='F', buffer=aow)
+        vxc = numpy.ndarray(coords.shape[0], buffer=vxcw)
+        vxc.fill(0.0)
+        # Form potential
+        for igrid in range(vxc.size):
+            for iatom in range(len(mol._atm)):
+                # Distance from nucleus
+                rnuc = numpy.linalg.norm(atom_coords[iatom,:] - coords[igrid,:])
+                if rnuc > eps:
+                    # Zeff(R)
+                    Zeff = sap_effective_charge(atom_charges[iatom], rnuc)
+                    vxc[igrid] -= Zeff/rnuc
+
+        # *.5 because vmat + vmat.T
+        #:aow = numpy.einsum('pi,p->pi', ao, .5*weight*vrho, out=aow)
+        aow = _scale_ao(ao, .5*weight*vxc, out=aow)
+        vmat += _dot_ao_ao(mol, ao, aow, mask, shls_slice, ao_loc)
+        vxc = None
+
+    vmat = vmat + vmat.conj().T
+    return vmat
+
 def nr_rks(ni, mol, grids, xc_code, dms, relativity=0, hermi=0,
            max_memory=2000, verbose=None):
     '''Calculate RKS XC functional and potential matrix on given meshgrids
@@ -1887,8 +1950,13 @@ def nr_fxc(self, mol, grids, xc_code, dm0, dms, spin=0, relativity=0, hermi=0,
             return self.nr_uks_fxc(mol, grids, xc_code, dm0, dms, relativity,
                                    hermi, rho0, vxc, fxc, max_memory, verbose)
 
+        #@lib.with_doc(nr_sap.__doc__)
+    def nr_sap(self, mol, grids, max_memory=2000, verbose=None):
+        return self.nr_sap_vxc(mol, grids, max_memory, verbose)
+
     nr_rks = nr_rks
     nr_uks = nr_uks
+    nr_sap_vxc = nr_sap_vxc
     nr_rks_fxc = nr_rks_fxc
     nr_uks_fxc = nr_uks_fxc
     cache_xc_kernel  = cache_xc_kernel

diff --git a/pyscf/dft/rks.py b/pyscf/dft/rks.py
@@ -144,6 +144,60 @@ def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
     vxc = lib.tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)
     return vxc
 
+def get_vsap(ks, mol=None):
+    '''Superposition of atomic potentials
+
+    S. Lehtola, Assessment of initial guesses for self-consistent
+    field calculations. Superposition of Atomic Potentials: simple yet
+    efficient, J. Chem. Theory Comput. 15, 1593 (2019). DOI:
+    10.1021/acs.jctc.8b01089. arXiv:1810.11659.
+
+    This function evaluates the effective charge of a neutral atom,
+    given by exchange-only LDA on top of spherically symmetric
+    unrestricted Hartree-Fock calculations as described in
+
+    S. Lehtola, L. Visscher, E. Engel, Efficient implementation of the
+    superposition of atomic potentials initial guess for electronic
+    structure calculations in Gaussian basis sets, J. Chem. Phys., in
+    press (2020).
+
+    The potentials have been calculated for the ground-states of
+    spherically symmetric atoms at the non-relativistic level of theory
+    as described in
+
+    S. Lehtola, "Fully numerical calculations on atoms with fractional
+    occupations and range-separated exchange functionals", Phys. Rev. A
+    101, 012516 (2020). DOI: 10.1103/PhysRevA.101.012516
+
+    using accurate finite-element calculations as described in
+
+    S. Lehtola, "Fully numerical Hartree-Fock and density functional
+    calculations. I. Atoms", Int. J. Quantum Chem. e25945 (2019).
+    DOI: 10.1002/qua.25945
+
+    .. note::
+        This function will modify the input ks object.
+
+    Args:
+        ks : an instance of :class:`RKS`
+            XC functional are controlled by ks.xc attribute.  Attribute
+            ks.grids might be initialized.
+
+    Returns:
+        matrix Vsap = Vnuc + J + Vxc.
+    '''
+    if mol is None: mol = ks.mol
+    t0 = (time.clock(), time.time())
+
+    if ks.grids.coords is None:
+        ks.grids.build(with_non0tab=True)
+        t0 = logger.timer(ks, 'setting up grids', *t0)
+
+    ni = ks._numint
+    max_memory = ks.max_memory - lib.current_memory()[0]
+    vsap = ni.nr_sap(mol, ks.grids, max_memory=max_memory)
+    return vsap
+
 # The vhfopt of standard Coulomb operator can be used here as an approximate
 # opt since long-range part Coulomb is always smaller than standard Coulomb.
 # It's safe to prescreen LR integrals with the integral estimation from
@@ -416,6 +470,7 @@ def dump_flags(self, verbose=None):
         return KohnShamDFT.dump_flags(self, verbose)
 
     get_veff = get_veff
+    get_vsap = get_vsap
     energy_elec = energy_elec
 
     def nuc_grad_method(self):

diff --git a/pyscf/dft/sap.py b/pyscf/dft/sap.py
@@ -0,0 +1,111 @@
+# Copyright (c) 2020, Susi Lehtola
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# * Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+# * Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+# * Neither the name of the <organization> nor the
+# names of its contributors may be used to endorse or promote products
+# derived from this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+# COPYRIGHT HOLDER BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+# USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+# OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+# SUCH DAMAGE.
+#
+# Routines for the implementation of the superposition of atomic
+# potentials guess for electronic structure calculations, see
+#
+# S. Lehtola, "Assessment of Initial Guesses for Self-Consistent Field
+# Calculations. Superposition of Atomic Potentials: Simple yet
+# Efficient", J. Chem. Theory Comput. 15, 1593 (2019).
+# DOI: 10.1021/acs.jctc.8b01089
+#
+# This function evaluates the effective charge of a neutral atom,
+# given by exchange-only LDA on top of spherically symmetric
+# unrestricted Hartree-Fock calculations as described in
+#
+# S. Lehtola, L. Visscher, E. Engel, Efficient implementation of the
+# superposition of atomic potentials initial guess for electronic
+# structure calculations in Gaussian basis sets, J. Chem. Phys., in
+# press (2020).
+#
+# The potentials have been calculated for the ground-states of
+# spherically symmetric atoms at the non-relativistic level of theory
+# as described in
+#
+# S. Lehtola, "Fully numerical calculations on atoms with fractional
+# occupations and range-separated exchange functionals", Phys. Rev. A
+# 101, 012516 (2020). DOI: 10.1103/PhysRevA.101.012516
+#
+# using accurate finite-element calculations as described in
+#
+# S. Lehtola, "Fully numerical Hartree-Fock and density functional
+# calculations. I. Atoms", Int. J. Quantum Chem. e25945 (2019).
+# DOI: 10.1002/qua.25945
+
+import numpy
+from pyscf.dft.sap_data import sap_Zeff
+
+def sap_effective_charge(Z, r):
+    '''
+    Calculates the effective charge for the superposition of atomic potentials.
+
+    S. Lehtola, "Assessment of Initial Guesses for Self-Consistent Field
+    Calculations. Superposition of Atomic Potentials: Simple yet
+    Efficient", J. Chem. Theory Comput. 15, 1593 (2019).
+    DOI: 10.1021/acs.jctc.8b01089
+
+    This function evaluates the effective charge of a neutral atom,
+    given by exchange-only LDA on top of spherically symmetric
+    unrestricted Hartree-Fock calculations as described in
+
+    S. Lehtola, L. Visscher, E. Engel, Efficient implementation of the
+    superposition of atomic potentials initial guess for electronic
+    structure calculations in Gaussian basis sets, J. Chem. Phys., in
+    press (2020).
+
+    The potentials have been calculated for the ground-states of
+    spherically symmetric atoms at the non-relativistic level of theory
+    as described in
+
+    S. Lehtola, "Fully numerical calculations on atoms with fractional
+    occupations and range-separated exchange functionals", Phys. Rev. A
+    101, 012516 (2020). DOI: 10.1103/PhysRevA.101.012516
+
+    using accurate finite-element calculations as described in
+
+    S. Lehtola, "Fully numerical Hartree-Fock and density functional
+    calculations. I. Atoms", Int. J. Quantum Chem. e25945 (2019).
+    DOI: 10.1002/qua.25945
+
+    Input:
+       Z: atomic charge
+       r: distance from nucleus
+    Output:
+       Z(r): screened charge
+    '''
+
+    if Z < 1:
+        return 0.0
+    if Z >= sap_Zeff.shape[1]:
+        raise ValueError('Atoms beyond Og are not supported')
+    if r < 0.0:
+        raise ValueError('Distance cannot be negative')
+    if r >= sap_Zeff.shape[0]:
+        return 0.0
+
+    # Linear interpolation
+    return numpy.interp(r, sap_Zeff[0,:], sap_Zeff[Z,:])
diff --git a/pyscf/dft/sap_data.py b/pyscf/dft/sap_data.py
diff --git a/pyscf/dft/uks.py b/pyscf/dft/uks.py
@@ -27,7 +27,6 @@
 from pyscf.scf import uhf
 from pyscf.dft import rks
 
-
 def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
     '''Coulomb + XC functional for UKS.  See pyscf/dft/rks.py
     :func:`get_veff` fore more details.
@@ -115,6 +114,50 @@ def get_veff(ks, mol=None, dm=None, dm_last=0, vhf_last=0, hermi=1):
     vxc = lib.tag_array(vxc, ecoul=ecoul, exc=exc, vj=vj, vk=vk)
     return vxc
 
+def get_vsap(ks, mol=None):
+    '''Superposition of atomic potentials
+
+    S. Lehtola, Assessment of initial guesses for self-consistent
+    field calculations. Superposition of Atomic Potentials: simple yet
+    efficient, J. Chem. Theory Comput. 15, 1593 (2019). DOI:
+    10.1021/acs.jctc.8b01089. arXiv:1810.11659.
+
+    This function evaluates the effective charge of a neutral atom,
+    given by exchange-only LDA on top of spherically symmetric
+    unrestricted Hartree-Fock calculations as described in
+
+    S. Lehtola, L. Visscher, E. Engel, Efficient implementation of the
+    superposition of atomic potentials initial guess for electronic
+    structure calculations in Gaussian basis sets, J. Chem. Phys., in
+    press (2020).
+
+    The potentials have been calculated for the ground-states of
+    spherically symmetric atoms at the non-relativistic level of theory
+    as described in
+
+    S. Lehtola, "Fully numerical calculations on atoms with fractional
+    occupations and range-separated exchange functionals", Phys. Rev. A
+    101, 012516 (2020). DOI: 10.1103/PhysRevA.101.012516
+
+    using accurate finite-element calculations as described in
+
+    S. Lehtola, "Fully numerical Hartree-Fock and density functional
+    calculations. I. Atoms", Int. J. Quantum Chem. e25945 (2019).
+    DOI: 10.1002/qua.25945
+
+    .. note::
+        This function will modify the input ks object.
+
+    Args:
+        ks : an instance of :class:`RKS`
+            XC functional are controlled by ks.xc attribute.  Attribute
+            ks.grids might be initialized.
+
+    Returns:
+        matrix Vsap = Vnuc + J + Vxc.
+    '''
+    Vsap = rks.get_vsap(ks, mol)
+    return [Vsap, Vsap]
 
 def energy_elec(ks, dm=None, h1e=None, vhf=None):
     if dm is None: dm = ks.make_rdm1()
@@ -139,6 +182,7 @@ def dump_flags(self, verbose=None):
         return self
 
     get_veff = get_veff
+    get_vsap = get_vsap
     energy_elec = energy_elec
 
     def nuc_grad_method(self):