Merge branch 'master' into 2519

.github/workflows/ci_linux/build_pyscf.sh

@@ -3,7 +3,7 @@
 set -e
 
 cd ./pyscf/lib
-curl -L "https://github.com/pyscf/pyscf-build-deps/blob/master/pyscf-2.4a-deps.tar.gz?raw=true" | tar xzf -
+curl -L "https://github.com/pyscf/pyscf-build-deps/blob/master/pyscf-2.8a-deps.tar.gz?raw=true" | tar xzf -
 mkdir build; cd build
 cmake -DBUILD_LIBXC=OFF -DBUILD_XCFUN=ON -DBUILD_LIBCINT=OFF -DXCFUN_MAX_ORDER=4 ..
 make -j4

.github/workflows/publish.yml

@@ -24,15 +24,11 @@ jobs:
         run: |
           ls ${{ github.workspace }}/linux-wheels
       - name: Publish to PyPI
-        # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags')
-        uses: pypa/gh-action-pypi-publish@release/v1
-        with:
-          user: __token__
-          #password: ${{ secrets.PYPI_TEST_API_TOKEN }}
-          #repository_url: https://test.pypi.org/legacy/
-          password: ${{ secrets.PYPI_API_TOKEN }}
-          packages-dir: ${{ github.workspace }}/linux-wheels
-          verbose: true
+        run: |
+            pip3 install twine
+            export TWINE_USERNAME=__token__
+            export TWINE_PASSWORD="${{ secrets.PYPI_API_TOKEN }}"
+            twine upload --verbose linux-wheels/*
 
   release-pypi-aarch64:
     runs-on: ubuntu-latest
@@ -67,15 +63,11 @@ jobs:
       run: |
         ls ${{ github.workspace }}/linux-wheels
     - name: Publish to PyPI
-      # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags')
-      uses: pypa/gh-action-pypi-publish@release/v1
-      with:
-        user: __token__
-        #password: ${{ secrets.PYPI_TEST_API_TOKEN }}
-        #repository_url: https://test.pypi.org/legacy/
-        password: ${{ secrets.PYPI_API_TOKEN }}
-        packages-dir: ${{ github.workspace }}/linux-wheels
-        verbose: true
+      run: |
+          pip3 install twine
+          export TWINE_USERNAME=__token__
+          export TWINE_PASSWORD="${{ secrets.PYPI_API_TOKEN }}"
+          twine upload --verbose linux-wheels/*
 
   release-pypi-macos-x86:
     name: Build wheels for macos
@@ -148,15 +140,11 @@ jobs:
         run: |
           ls ${{ github.workspace }}/dist
       - name: Publish to PyPI
-        # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags')
-        uses: pypa/gh-action-pypi-publish@release/v1
-        with:
-          user: __token__
-          #password: ${{ secrets.PYPI_TEST_API_TOKEN }}
-          #repository_url: https://test.pypi.org/legacy/
-          password: ${{ secrets.PYPI_API_TOKEN }}
-          packages-dir: ${{ github.workspace }}/dist
-          verbose: true
+        run: |
+            pip3 install twine
+            export TWINE_USERNAME=__token__
+            export TWINE_PASSWORD="${{ secrets.PYPI_API_TOKEN }}"
+            twine upload --verbose linux-wheels/*
 
   release-conda-linux:
     runs-on: ubuntu-latest

CHANGELOG

@@ -366,6 +366,7 @@ PySCF 1.7.6 (2021-03-28)
   - Auxiliary second-order Green's function perturbation theory (AGF2)
   - Smearing for molecules
   - Visscher small component correction approximation for DHF
+  - DFT+U
 * Improved
   - Threading safety in Mole temporary context
   - Basis parser to support arithmetic expressions in basis data

examples/dft/14-non_collinear_gks.py

@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+#
+# Author: Qiming Sun <osirpt.sun@gmail.com>
+#
+
+'''
+GKS with non-collinear functional
+'''
+
+from pyscf import gto
+
+mol = gto.M(atom="H 0 0 0; F 0 0 1", basis='unc-sto3g', verbose=4)
+mf = mol.GKS()
+mf.xc = 'pbe'
+# Enable non-collinear functional. GKS calls collinear functional by default
+# mcol is short for multi-collinear functional. This is one treatment of
+# non-collinear method which avoids the singularity issue in functional.
+# For more details of multi-collinear method, please see
+#   Noncollinear density functional theory, Zhichen Pu, et. al., Rev. Research 5, 013036
+mf.collinear = 'mcol'
+mf.kernel()

examples/dft/16-dft_d3.py

@@ -6,18 +6,21 @@
 '''
 D3 and D4 Dispersion.
 
-This is a simplified dispersion interface to
-d3 (https://github.com/dftd3/simple-dftd3) and
-d4 (https://github.com/dftd4/dftd4) libraries.
-This interface can automatically configure the necessary settings including
-dispersion, xc, and nlc attributes of PySCF mean-field objects. However,
-advanced features of d3 and d4 program are not available.
+This example shows the functionality of the pyscf-dispersion extension, which
+implements a simplified interface to d3 (https://github.com/dftd3/simple-dftd3)
+and d4 (https://github.com/dftd4/dftd4) libraries. This interface can
+automatically configure the necessary settings including dispersion, xc, and nlc
+attributes of PySCF mean-field objects. However, advanced features of d3 and d4
+program are not available. To execute the code in this example, you would need
+to install the pyscf-dispersion extension:
+
+pip install pyscf-dispersion
 
 If you need to access more features d3 and d4 libraries, such as overwriting the
 dispersion parameters, you can use the wrapper provided by the simple-dftd3 and
-dftd4 libraries. When using these libraries, please disable the .disp attribute
-of the underlying mean-field object, and properly set the .xc and .nlc attributes
-following this example.
+dftd4 libraries. When accessing dispersion through the APIs of these libraries,
+please disable the .disp attribute of the mean-field object, and properly set
+the .xc and .nlc attributes as shown in this example.
 '''
 
 mol = pyscf.M(
@@ -134,7 +137,7 @@
 mf.disp = 'd4'
 mf.kernel() # Crash
 
-# Please note, not every xc functional can be used for D3/D4 corrections.
+# Please note, NOT every xc functional can be used for D3/D4 corrections.
 # For the valid xc and D3/D4 combinations, please refer to
 # https://github.com/dftd3/simple-dftd3/blob/main/assets/parameters.toml
 # https://github.com/dftd4/dftd4/blob/main/assets/parameters.toml

examples/dft/24-define_xc_functional.py

@@ -12,6 +12,7 @@
   functionals provided by Libxc or XcFun library.
 '''
 
+import numpy as np
 from pyscf import gto
 from pyscf import dft
 
@@ -25,28 +26,19 @@
 # half-half exact exchange and GGA functional
 hybrid_coeff = 0.5
 
-def eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None):
+def eval_gga_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None):
     # A fictitious XC functional to demonstrate the usage
-    rho0, dx, dy, dz = rho[:4]
+    rho0, dx, dy, dz = rho
     gamma = (dx**2 + dy**2 + dz**2)
     exc = .01 * rho0**2 + .02 * (gamma+.001)**.5
     vrho = .01 * 3 * rho0**2 + .02 * (gamma+.001)**.5
     vgamma = .02 * .5 * (gamma+.001)**(-.5)
-    vlapl = None
-    vtau = None
-    vxc = (vrho, vgamma, vlapl, vtau)
+    vxc = (vrho, vgamma)
     v2rho2 = 0.01 * 6 * rho0
-    v2rhosigma = .02 * .5 * (gamma+.001)**(-.5)
+    v2rhosigma = np.zeros(gamma.shape)
     v2sigma2 = 0.02 * .5 * -.5 * (gamma+.001)**(-1.5)
-    v2lapl2 = None
-    vtau2 = None
-    v2rholapl = None
-    v2rhotau = None
-    v2lapltau = None
-    v2sigmalapl = None
-    v2sigmatau = None
     # 2nd order functional derivative
-    fxc = (v2rho2, v2rhosigma, v2sigma2, v2lapl2, vtau2, v2rholapl, v2rhotau, v2lapltau, v2sigmalapl, v2sigmatau)
+    fxc = (v2rho2, v2rhosigma, v2sigma2)
     kxc = None  # 3rd order functional derivative
 
     # Mix with existing functionals
@@ -58,15 +50,42 @@ def eval_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=Non
     return exc, vxc, fxc, kxc
 
 mf = dft.RKS(mol)
-mf = mf.define_xc_(eval_xc, 'GGA', hyb=hybrid_coeff)
+mf = mf.define_xc_(eval_gga_xc, 'GGA', hyb=hybrid_coeff)
 mf.verbose = 4
 mf.kernel()
 
+def eval_mgga_xc(xc_code, rho, spin=0, relativity=0, deriv=1, omega=None, verbose=None):
+    # A fictitious XC functional to demonstrate the usage
+    rho0, dx, dy, dz, tau = rho
+    gamma = (dx**2 + dy**2 + dz**2)
+    exc = .01 * rho0**2 + .02 * (gamma+.001)**.5 +  0.01 * tau**2
+    vrho = .01 * 3 * rho0**2 + .02 * (gamma+.001)**.5
+    vgamma = .02 * .5 * (gamma+.001)**(-.5)
+    vtau = 0.02 * tau
+    vxc = (vrho, vgamma, vtau)
+    v2rho2 = 0.01 * 6 * rho0
+    v2rhosigma = np.zeros(gamma.shape)
+    v2sigma2 = 0.02 * .5 * -.5 * (gamma+.001)**(-1.5)
+    v2tau2 = np.full(tau.shape, 0.02)
+    v2rhotau = np.zeros(tau.shape)
+    v2sigmatau = np.zeros(tau.shape)
+    # 2nd order functional derivative
+    fxc = (v2rho2, v2rhosigma, v2sigma2, v2tau2, v2rhotau, v2sigmatau)
+    kxc = None  # 3rd order functional derivative
+
+    # Mix with existing functionals
+    pbe_xc = dft.libxc.eval_xc('pbe,pbe', rho[:4], spin, relativity, deriv, verbose)
+    exc += pbe_xc[0] * 0.5
+    vrho += pbe_xc[1][0] * 0.5
+    vgamma += pbe_xc[1][1] * 0.5
+    # The output follows the libxc.eval_xc API convention
+    return exc, vxc, fxc, kxc
+
 # half exact exchange in which 40% of the exchange is computed with short
 # range part of the range-separation Coulomb operator (omega = 0.8)
 beta = 0.2
 rsh_coeff = (0.8, hybrid_coeff-beta, beta)
 mf = dft.RKS(mol)
-mf = mf.define_xc_(eval_xc, 'GGA', rsh=rsh_coeff)
+mf = mf.define_xc_(eval_mgga_xc, 'MGGA', rsh=rsh_coeff)
 mf.verbose = 4
 mf.kernel()

examples/geomopt/16-apply_soscf.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+
+'''
+Automatically apply SOSCF for unconverged SCF calculations in geometry optimization.
+'''
+
+import pyscf
+from pyscf.geomopt import geometric_solver, as_pyscf_method
+
+# Force SCF to stop early at an unconverged state
+pyscf.scf.hf.RHF.max_cycle = 5
+
+mol = pyscf.M(
+    atom='''
+    O 0 0 0
+    H 0 .7 .8
+    H 0 -.7 .8
+    ''')
+mf = mol.RHF()
+try:
+    mol_eq = geometric_solver.optimize(mf)
+except RuntimeError:
+    print('geometry optimization should stop for unconverged calculation')
+
+# Apply SOSCF when SCF is not converged
+mf_scanner = mf.as_scanner()
+def apply_soscf_as_needed(mol):
+    mf_scanner(mol)
+    if not mf_scanner.converged:
+        mf_soscf = mf_scanner.newton().run()
+        for key in ['converged', 'mo_energy', 'mo_coeff', 'mo_occ', 'e_tot', 'scf_summary']:
+            setattr(mf_scanner, key, getattr(mf_soscf, key))
+    grad = mf_scanner.Gradients().kernel()
+    return mf_scanner.e_tot, grad
+
+mol_eq = geometric_solver.optimize(as_pyscf_method(mol, apply_soscf_as_needed))

examples/pbc/22-dft+u.py

@@ -0,0 +1,40 @@
+#!/usr/bin/env python
+
+'''
+DFT+U with kpoint sampling
+'''
+
+from pyscf.pbc import gto, dft
+cell = gto.Cell()
+cell.unit = 'A'
+cell.atom = 'C 0.,  0.,  0.; C 0.8917,  0.8917,  0.8917'
+cell.a = '''0.      1.7834  1.7834
+            1.7834  0.      1.7834
+            1.7834  1.7834  0.    '''
+
+cell.basis = 'gth-dzvp'
+cell.pseudo = 'gth-pade'
+cell.verbose = 4
+cell.build()
+
+kmesh = [2, 2, 2]
+kpts = cell.make_kpts(kmesh)
+# Add U term to the 2p orbital of the second Carbon atom
+U_idx = ['1 C 2p']
+U_val = [5.0]
+mf = dft.KRKSpU(cell, kpts, U_idx=U_idx, U_val=U_val, minao_ref='gth-szv')
+print(mf.U_idx)
+print(mf.U_val)
+mf.run()
+
+# When space group symmetry in k-point samples is enabled, the symmetry adapted
+# DFT+U method will be invoked automatically.
+kpts = cell.make_kpts(
+    kmesh, wrap_around=True, space_group_symmetry=True, time_reversal_symmetry=True)
+# Add U term to 2s and 2p orbitals
+U_idx = ['2p', '2s']
+U_val = [5.0, 2.0]
+mf = dft.KUKSpU(cell, kpts, U_idx=U_idx, U_val=U_val, minao_ref='gth-szv')
+print(mf.U_idx)
+print(mf.U_val)
+mf.run()

examples/scf/14-restart.py

@@ -58,12 +58,12 @@
 mf.kernel(dm)
 
 # 4. Last method can be simplified.  Calling .__dict__.update function to
-# upgrade the SCF object which assign mf.mo_coeff and mf.mo_occ etc. default
-# values. Then calling make_rdm1 function to generate a density matrix.
+# upgrade the SCF objec. Attributes such as mf.mo_coeff and mf.mo_occ etc. are
+# assigned by this function. These attributes will be used to construct an
+# initial guess automatically.
 mf = scf.RHF(mol)
 mf.__dict__.update(scf.chkfile.load('h2o.chk', 'scf'))
-dm = mf.make_rdm1()
-mf.kernel(dm)
+mf.kernel()
 
 
 #

pyscf/ao2mo/_ao2mo.py

@@ -137,6 +137,9 @@ def nr_e1(eri, mo_coeff, orbs_slice, aosym='s1', mosym='s1', out=None):
     if out.size == 0:
         return out
 
+    if eri.dtype != numpy.double:
+        raise TypeError('_ao2mo.nr_e1 is for double precision only')
+
     fdrv = getattr(libao2mo, 'AO2MOnr_e2_drv')
     pao_loc = ctypes.POINTER(ctypes.c_void_p)()
     c_nbas = ctypes.c_int(0)
@@ -184,6 +187,9 @@ def nr_e2(eri, mo_coeff, orbs_slice, aosym='s1', mosym='s1', out=None,
     if out.size == 0:
         return out
 
+    if eri.dtype != numpy.double:
+        raise TypeError('_ao2mo.nr_e2 is for double precision only')
+
     if ao_loc is None:
         pao_loc = ctypes.POINTER(ctypes.c_void_p)()
         c_nbas = ctypes.c_int(0)
@@ -283,6 +289,9 @@ def r_e2(eri, mo_coeff, orbs_slice, tao, ao_loc, aosym='s1', out=None):
     if out.size == 0:
         return out
 
+    if eri.dtype != numpy.complex128:
+        raise TypeError('_ao2mo.r_e2 is for complex double precision only')
+
     tao = numpy.asarray(tao, dtype=numpy.int32)
     if ao_loc is None:
         c_ao_loc = ctypes.POINTER(ctypes.c_void_p)()

pyscf/ao2mo/incore.py

@@ -202,6 +202,9 @@ def half_e1(eri_ao, mo_coeffs, compact=True):
     if nij_pair == 0:
         return eri1
 
+    if eri_ao.dtype != numpy.double:
+        raise TypeError('ao2mo.incore.half_e1 is for double precision only')
+
     if eri_ao.size == nao_pair**2: # 4-fold symmetry
         # half_e1 first transforms the indices which are contiguous in memory
         # transpose the 4-fold integrals to make ij the contiguous indices

pyscf/cc/ccsd.py

@@ -942,7 +942,7 @@ class CCSDBase(lib.StreamObject):
         'async_io', 'incore_complete', 'cc2', 'callback',
         'mol', 'verbose', 'stdout', 'frozen', 'level_shift',
         'mo_coeff', 'mo_occ', 'cycles', 'converged_lambda', 'emp2', 'e_hf',
-        'e_corr', 't1', 't2', 'l1', 'l2', 'chkfile',
+        'converged', 'e_corr', 't1', 't2', 'l1', 'l2', 'chkfile',
     }
 
     def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None):