Merge remote-tracking branch 'pyscf/master'

.github/workflows/ci_linux/build_pyscf.sh

@@ -5,7 +5,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 -
 mkdir build; cd build
-cmake -DBUILD_LIBXC=OFF -DBUILD_XCFUN=OFF -DBUILD_LIBCINT=OFF -DXCFUN_MAX_ORDER=4 ..
+cmake -DBUILD_LIBXC=OFF -DBUILD_XCFUN=ON -DBUILD_LIBCINT=OFF -DXCFUN_MAX_ORDER=4 ..
 make -j4
 cd ..
 rm -Rf build

FEATURES

@@ -20,6 +20,7 @@
   - Scalar relativistic DFT
   - Density fitting DFT
   - General XC functional evaluator (for Libxc or XcFun)
+  - Multi-collinear functional
   - General AO evaluator
   - VV10 NLC functional for finite size systems
   - range-separated hybrid features for RKS and UKS, including
@@ -79,6 +80,7 @@
   - EOM-IP/EA/EE-RCCSD and EOM-IP/EA/EE-UCCSD
   - RCC2
   - Density-fitting RCCSD
+  - BCCD
 
 
 * CCSD(T)
@@ -95,6 +97,7 @@
   - Selected-CI
   - Selected-CI 1, 2-particle density matrices
   - RCISD, UCISD and GCISD 1-particle transition density matrices
+  - Density-fitting RCISD, UCISD
 
 
 * Full CI
@@ -240,6 +243,8 @@
   - ddCOSMO
   - ddPCM
   - ddCOSMO analytical nuclear gradients
+  - SMD, and its nuclear gradients and hessian
+  - PCM, and its nuclear gradients and hessian
 
 
 * Tools

examples/scf/02-ghf.py

@@ -30,6 +30,6 @@
 #
 mf = scf.GHF(mol)
 dm = mf.get_init_guess() + 0j
-dm[0,:] += .1j
-dm[:,0] -= .1j
+dm[0,:] += .05j
+dm[:,0] -= .05j
 mf.kernel(dm0=dm)

pyscf/df/df.py

@@ -357,7 +357,12 @@ def get_jk(self, dm, hermi=1, with_j=True, with_k=True,
         with self.range_coulomb(omega) as rsh_df:
             return df_jk.r_get_jk(rsh_df, dm, hermi, with_j, with_k)
 
+    def get_eri(self):
+        raise NotImplementedError
+    get_ao_eri = get_eri
+
     def ao2mo(self, mo_coeffs):
         raise NotImplementedError
+    get_mo_eri = ao2mo
 
 GDF4C = DF4C

pyscf/dft/libxc.py

@@ -767,14 +767,7 @@ def _xc_key_without_underscore(xc_keys):
     'TPSS0'         : '.25*HF + .75*TPSS, TPSS',
 })  # noqa: E501
 
-# Issue 1480
-if not hasattr(__config__, 'B3LYP_WITH_VWN5'):
-    warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, '
-                  'corresponding to the original definition by Stephens et al. (issue 1480) '
-                  'and the same as the B3LYP functional in Gaussian. '
-                  'To restore the VWN5 definition, you can put the setting '
-                  '"B3LYP_WITH_VWN5 = True" in pyscf_conf.py')
-elif getattr(__config__, 'B3LYP_WITH_VWN5', False):
+if getattr(__config__, 'B3LYP_WITH_VWN5', False):
     XC_CODES['B3P86' ] = 'B3P86V5'
     XC_CODES['B3LYP' ] = 'B3LYP5'
     XC_CODES['X3LYP' ] = 'X3LYP5'
@@ -1102,6 +1095,16 @@ def parse_xc(description):
     elif not isinstance(description, str): #isinstance(description, (tuple,list)):
         return parse_xc('%s,%s' % tuple(description))
 
+    if (description.upper() in ('B3P86', 'B3LYP', 'X3LYP') and
+        not getattr(parse_xc, 'b3lyp5_warned', False) and
+        not hasattr(__config__, 'B3LYP_WITH_VWN5')):
+        parse_xc.b3lyp5_warned = True
+        warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, '
+                      'corresponding to the original definition by Stephens et al. (issue 1480) '
+                      'and the same as the B3LYP functional in Gaussian. '
+                      'To restore the VWN5 definition, you can put the setting '
+                      '"B3LYP_WITH_VWN5 = True" in pyscf_conf.py')
+
     def assign_omega(omega, hyb_or_sr, lr=0):
         if hyb[2] == omega or omega == 0:
             hyb[0] += hyb_or_sr

pyscf/dft/rks.py

@@ -461,12 +461,11 @@ def initialize_grids(self, mol=None, dm=None):
         '''Initialize self.grids the first time call get_veff'''
         if mol is None: mol = self.mol
 
+        ground_state = getattr(dm, 'ndim', 0) == 2
         if self.grids.coords is None:
             t0 = (logger.process_clock(), logger.perf_counter())
             self.grids.build(with_non0tab=True)
-            if (self.small_rho_cutoff > 1e-20 and
-                # dm.ndim == 2 indicates ground state
-                isinstance(dm, numpy.ndarray) and dm.ndim == 2):
+            if self.small_rho_cutoff > 1e-20 and ground_state:
                 # Filter grids the first time setup grids
                 self.grids = prune_small_rho_grids_(self, self.mol, dm,
                                                     self.grids)
@@ -476,9 +475,7 @@ def initialize_grids(self, mol=None, dm=None):
         if is_nlc and self.nlcgrids.coords is None:
             t0 = (logger.process_clock(), logger.perf_counter())
             self.nlcgrids.build(with_non0tab=True)
-            if (self.small_rho_cutoff > 1e-20 and
-                # dm.ndim == 2 indicates ground state
-                isinstance(dm, numpy.ndarray) and dm.ndim == 2):
+            if self.small_rho_cutoff > 1e-20 and ground_state:
                 # Filter grids the first time setup grids
                 self.nlcgrids = prune_small_rho_grids_(self, self.mol, dm,
                                                        self.nlcgrids)

pyscf/dft/test/test_he.py

@@ -14,6 +14,7 @@
 # limitations under the License.
 
 import unittest
+import tempfile
 import numpy
 from pyscf import gto
 from pyscf import lib
@@ -191,6 +192,21 @@ def test_convert(self):
         self.assertTrue(isinstance(udks.to_dhf(), scf.dhf.DHF))
         self.assertTrue(isinstance(udks.to_dks('pbe'), dft.dks.DKS))
 
+    # issue 1986
+    def test_init_guess_chkfile(self):
+        with tempfile.NamedTemporaryFile() as tmpf:
+            mol = gto.M(atom='He 0 0 0', basis='631g', charge=1, spin=1)
+            mf = dft.RKS(mol)
+            mf.chkfile = tmpf.name
+            e1 = mf.kernel()
+            mf = dft.RKS(mol)
+            mf.init_guess = 'chkfile'
+            mf.chkfile = tmpf.name
+            mf.max_cycle = 1
+            e2 = mf.kernel()
+            self.assertAlmostEqual(e1, e2, 9)
+
+
 if __name__ == "__main__":
     print("Full Tests for He")
     unittest.main()

pyscf/dft/test/test_xc_deriv.py

@@ -415,7 +415,7 @@ def test_libxc_mgga_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), 19.977512805950784, 7)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(dft.libxc.max_deriv_order('pbe,') <= 3)
+    @unittest.skipIf(dft.libxc.max_deriv_order('pbe,') <= 3, 'libxc order')
     def test_libxc_gga_deriv4(self):
         rho1 = rho[:,:4].copy()
         xc1 = dft.libxc.eval_xc_eff('PBE', rho1, deriv=4)
@@ -425,7 +425,7 @@ def test_libxc_gga_deriv4(self):
         xc1 = dft.libxc.eval_xc_eff('PBE', rho1, deriv=4)
         self.assertAlmostEqual(xc1.sum(), -869.6617638095072, 4)
 
-    @unittest.skip(not hasattr(dft, 'xcfun'))
+    @unittest.skipIf(not hasattr(dft, 'xcfun'), 'xcfun order')
     def test_xcfun_lda_deriv3(self):
         rho1 = rho[:,0].copy()
         ref = eval_xc_eff('LDA,', rho1, 3, dft.xcfun)
@@ -457,7 +457,7 @@ def test_xcfun_lda_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), -2.7477984980958627, 9)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(not hasattr(dft, 'xcfun'))
+    @unittest.skipIf(not hasattr(dft, 'xcfun'), 'xcfun order')
     def test_xcfun_gga_deriv3(self):
         rho1 = rho[:,:4].copy()
         ref = eval_xc_eff('PBE', rho1, 3, dft.xcfun)
@@ -489,7 +489,7 @@ def test_xcfun_gga_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), -3.0715856471099032, 9)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(not hasattr(dft, 'xcfun'))
+    @unittest.skipIf(not hasattr(dft, 'xcfun'), 'xcfun order')
     def test_xcfun_mgga_deriv3(self):
         rho1 = rho
         ref = eval_xc_eff('M06', rho1, 3, dft.xcfun)
@@ -521,7 +521,7 @@ def test_xcfun_mgga_deriv3(self):
         self.assertAlmostEqual(xc1.sum(), 19.977512805950784, 9)
         self.assertAlmostEqual(abs(ref[1] - xc1).max(), 0, 9)
 
-    @unittest.skip(not (hasattr(dft, 'xcfun') and dft.xcfun.MAX_DERIV_ORDER > 3))
+    @unittest.skipIf(not (hasattr(dft, 'xcfun') and dft.xcfun.MAX_DERIV_ORDER > 3), 'xcfun order')
     def test_xcfun_gga_deriv4(self):
         rho1 = rho[:,:4].copy()
         xc1 = dft.xcfun.eval_xc_eff('PBE', rho1, deriv=4)

pyscf/dft/xcfun.py

@@ -252,7 +252,12 @@
           'B97XC'   , 'B97_1XC' , 'B97_2XC' , 'M05XC'   , 'TPSSH'   ,
           'HFLYP'}
 RSH_XC = {'CAMB3LYP'}
-MAX_DERIV_ORDER = ctypes.c_int.in_dll(_itrf, 'XCFUN_max_deriv_order').value
+
+# The compatibility with the old libxcfun_itrf.so library
+try:
+    MAX_DERIV_ORDER = ctypes.c_int.in_dll(_itrf, 'XCFUN_max_deriv_order').value
+except ValueError:
+    MAX_DERIV_ORDER = 3
 
 VV10_XC = {
     'B97M_V'    : (6.0, 0.01),

pyscf/fci/cistring.py

@@ -343,12 +343,13 @@ def addr2str(norb, nelec, addr):
     '''Convert CI determinant address to string'''
     if norb >= 64:
         raise NotImplementedError('norb >= 64')
-    assert num_strings(norb, nelec) > addr
+    max_addr = num_strings(norb, nelec)
+    assert max_addr > addr
 
-    if addr < 2**31:
+    if max_addr < 2**31:
         return addrs2str(norb, nelec, [addr])[0]
-
-    return _addr2str(norb, nelec, addr)
+    else:
+        return _addr2str(norb, nelec, addr)
 
 def _addr2str(norb, nelec, addr):
     if addr == 0 or nelec == norb or nelec == 0:

pyscf/fci/test/test_cistring.py

@@ -156,6 +156,12 @@ def test_addr2str(self):
         self.assertEqual(bin(cistring.addr2str(6, 3, 8)), '0b11010')
         self.assertEqual(bin(cistring.addr2str(7, 4, 9)), '0b110011')
 
+        # Test large addresss
+        string = 0b101101111101101111001110111111110111111100
+        address = cistring.str2addr(norb=63, nelec=32, string=string)
+        string2 = cistring.addr2str(norb=63, nelec=32, addr=address)
+        self.assertEqual(string, string2)
+
     def test_str2addr(self):
         self.assertEqual(str2addr(6, 3, int('0b11001' ,2)), cistring.str2addr(6, 3, int('0b11001' ,2)))
         self.assertEqual(str2addr(6, 3, int('0b11010' ,2)), cistring.str2addr(6, 3, int('0b11010' ,2)))

pyscf/gto/basis/__init__.py

@@ -629,7 +629,13 @@ def load(filename_or_basisname, symb, optimize=OPTIMIZE_CONTRACTION):
 
         # Last, a trial to access Basis Set Exchange database
         from pyscf.gto.basis import bse
-        if bse.basis_set_exchange is not None:
+        if bse.basis_set_exchange is None:
+            import warnings
+            warnings.warn(
+                'Basis may be available in basis-set-exchange. '
+                'It is recommended to install basis-set-exchange: '
+                'pip install basis-set-exchange')
+        else:
             try:
                 bse_obj = bse.basis_set_exchange.api.get_basis(
                     filename_or_basisname, elements=symb)

pyscf/gto/mole.py

@@ -1209,6 +1209,7 @@ def copy(mol, deep=True):
     newmol._ecp    = copy.deepcopy(mol._ecp)
     newmol.pseudo  = copy.deepcopy(mol.pseudo)
     newmol._pseudo = copy.deepcopy(mol._pseudo)
+    newmol.magmom  = list(mol.magmom)
     return newmol
 
 def pack(mol):
@@ -1229,6 +1230,7 @@ def pack(mol):
             'ecp'     : mol.ecp,
             'pseudo'  : mol.pseudo,
             '_nelectron': mol._nelectron,
+            'magmom'  : mol.magmom,
             'verbose' : mol.verbose}
     return mdic
 def unpack(moldic):
@@ -2581,11 +2583,16 @@ def build(self, dump_input=True, parse_arg=ARGPARSE,
             # number of electrons are consistent.
             self.nelec
 
-        if self.magmom is None or len(self.magmom) != self.natm:
+        if self.magmom is None:
+            self.magmom = [0,] * self.natm
+        elif len(self.magmom) != self.natm:
+            logger.warn(self, 'len(magmom) != natm. Set magmom to zero')
             self.magmom = [0,] * self.natm
         if self.spin == 0 and abs(numpy.sum(self.magmom) - self.spin) > 1e-6:
             #don't check for unrestricted calcs.
             raise ValueError("mol.magmom is set incorrectly.")
+        if isinstance(self.magmom, np.ndarray):
+            self.magmom = self.magmom.tolist()
 
         if self.symmetry:
             self._build_symmetry()

pyscf/lib/CMakeLists.txt

@@ -29,6 +29,7 @@ if (CMAKE_COMPILER_IS_GNUCC) # Does it skip the link flag on old OsX?
   endif()
 endif()
 set(CMAKE_C_FLAGS "-Wall ${CMAKE_C_FLAGS}")
+set(CMAKE_C_STANDARD 99)
 
 option(BUILD_MARCH_NATIVE "gcc flag -march=native" off)
 if (BUILD_MARCH_NATIVE)
@@ -159,16 +160,18 @@ option(USE_QCINT "Using the qcint library (optimized for x86-64) for gaussian in
 
 if(BUILD_LIBCINT)
   set(LIBCINT_GIT https://github.com/sunqm/libcint.git) # libcint is a portable, cross-platform implementation
+  set(LIBCINT_VERSION v6.1.1)
   if (USE_QCINT)
     set(LIBCINT_GIT https://github.com/sunqm/qcint.git) # qcint is an optimized implementation for x86-64 architecture
+    set(LIBCINT_VERSION v6.1.2)
     if(NOT BUILD_MARCH_NATIVE)
       message(WARNING "The BUILD_MARCH_NATIVE option is not specified! qcint may not compile unless you explicitly pass compiler flags that turn on vectorization!")    
     endif()
   endif()
 
   ExternalProject_Add(libcint
     GIT_REPOSITORY ${LIBCINT_GIT}
-    GIT_TAG v6.0.0
+    GIT_TAG ${LIBCINT_VERSION}
     PREFIX ${PROJECT_BINARY_DIR}/deps
     INSTALL_DIR ${PROJECT_SOURCE_DIR}/deps
     CMAKE_ARGS -DCMAKE_BUILD_TYPE=RELEASE

pyscf/lib/ao2mo/nr_ao2mo.c

@@ -17,6 +17,7 @@
  */
 
 #include <stdlib.h>
+#include <stdio.h>
 #include <math.h>
 #include <assert.h>
 //#define NDEBUG
@@ -1223,7 +1224,11 @@ void AO2MOnr_e1_drv(int (*intor)(), void (*fill)(), void (*ftrans)(), int (*fmmm
 {
         int nao = ao_loc[nbas];
         double *eri_ao = malloc(sizeof(double) * nao*nao*nkl*ncomp);
-        assert(eri_ao);
+        if (eri_ao == NULL) {
+                fprintf(stderr, "malloc(%zu) falied in AO2MOnr_e1_drv\n",
+                        sizeof(double) * nao*nao*nkl*ncomp);
+                exit(1);
+        }
         AO2MOnr_e1fill_drv(intor, fill, eri_ao, klsh_start, klsh_count,
                            nkl, ncomp, ao_loc, cintopt, vhfopt,
                            atm, natm, bas, nbas, env);

pyscf/lib/ao2mo/nrr_ao2mo.c

@@ -18,6 +18,7 @@
  */
 
 #include <stdlib.h>
+#include <stdio.h>
 #include <complex.h>
 #include <math.h>
 #include <assert.h>
@@ -237,7 +238,11 @@ void AO2MOnrr_e1_drv(int (*intor)(), void (*fill)(),
 
 
         double *eri_ao = malloc(sizeof(double)* nao*nao*nkl*ncomp);
-        assert(eri_ao);
+        if (eri_ao == NULL) {
+                fprintf(stderr, "malloc(%zu) falied in AO2MOnrr_e1_drv\n",
+                        sizeof(double) * nao*nao*nkl*ncomp);
+                exit(1);
+        }
         int ish, kl;
         int (*fprescreen)();
         if (vhfopt) {

pyscf/lib/ao2mo/r_ao2mo.c

@@ -17,6 +17,7 @@
  */
 
 #include <stdlib.h>
+#include <stdio.h>
 #include <complex.h>
 #include <math.h>
 #include <assert.h>
@@ -833,7 +834,11 @@ void AO2MOr_e1_drv(int (*intor)(), void (*fill)(),
 
         double complex *eri_ao = malloc(sizeof(double complex)
                                         * nao*nao*nkl*ncomp);
-        assert(eri_ao);
+        if (eri_ao == NULL) {
+                fprintf(stderr, "malloc(%zu) falied in AO2MOr_e1_drv\n",
+                        sizeof(double complex) * nao*nao*nkl*ncomp);
+                exit(1);
+        }
         int ish, kl;
         int (*fprescreen)();
         if (vhfopt != NULL) {