Merge branch 'master' into dsrg-mrpt

Author: shuhangli98

.github/workflows/ci_conda.yml

@@ -0,0 +1,25 @@
+name: CI_conda
+
+on:
+  workflow_dispatch:
+
+jobs:
+  build-conda-linux:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+    steps:
+      - uses: actions/checkout@v4
+      - name: Setup conda
+        uses: s-weigand/setup-conda@v1
+        with:
+          update-conda: true
+          conda-channels: conda-forge
+      - run: conda --version
+      - run: which python
+      - name: Build conda package
+        run: |
+          export CMAKE_BUILD_PARALLEL_LEVEL=4
+          conda install -y conda-build
+          conda config --set anaconda_upload False
+          conda build --output-folder . conda

.github/workflows/publish.yml

@@ -39,7 +39,7 @@ jobs:
     strategy:
       fail-fast: false
     env:
-      img: quay.io/pypa/manylinux2014_aarch64:2023-03-12-25fd859
+      img: quay.io/pypa/manylinux2014_aarch64:latest
     steps:
     - name: Checkout
       uses: actions/checkout@v4

README.md

@@ -14,6 +14,7 @@ Python-based Simulations of Chemistry Framework
 * [Documentation](http://www.pyscf.org)
 * [Installation](#installation)
 * [Features](../master/FEATURES)
+* [News](https://pyscf.org/news.html): **2nd PySCF Developers Meeting!**
 
 
 # Installation

examples/gto/04-input_basis.py

@@ -65,7 +65,7 @@
 mol = gto.M(
     atom = '''O 0 0 0; H1 0 1 0; H2 0 0 1''',
     basis = {'O': gto.parse('''
-# Parse NWChem format basis string (see https://bse.pnl.gov/bse/portal).
+# Parse NWChem format basis string (see https://www.basissetexchange.org/).
 # Comment lines are ignored
 #BASIS SET: (6s,3p) -> [2s,1p]
 O    S

examples/gto/31-basis_set_exchange.py

@@ -14,6 +14,22 @@
 a local copy of the Basis Set Exchange.
 '''
 
+# PySCF has a native interface to the Basis Set Exchange, which is
+# used as a fall-through if the basis set was not found within PySCF
+# itself
+mol = gto.M(
+    atom = '''
+N          0.6683566134        0.2004327755        0.0000000000
+H          0.9668193796       -0.3441960976        0.8071193402
+H          0.9668193796       -0.3441960976       -0.8071193402
+F         -0.7347916126       -0.0467759204        0.0000000000
+''',
+    basis = 'HGBSP1-7',
+    verbose = 4
+)
+
+
+# One can also use the Basis Set Exchange in a more direct fashion
 mol = gto.M(
     atom = '''
 N          0.6683566134        0.2004327755        0.0000000000

examples/pbc/01-input_output_geometry.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python
+
+'''
+This example demonstrates some ways to input the Cell geometry, including
+lattice vectors, and to write the Cell geometry to file.
+
+Example solid is wurtzite BN.
+'''
+
+from pyscf.pbc import gto
+
+# Input Cartesian coordinates for the lattice vectors a and the atomic
+# positions. Coordinates are in Angstrom by default
+
+cell = gto.Cell()
+cell.a = [[  2.5539395809,  0.0000000000,  0.0000000000],
+          [ -1.2769697905,  2.2117765568,  0.0000000000],
+          [  0.0000000000,  0.0000000000,  4.2268548012]]
+cell.atom = '''
+            B     1.276969829         0.737258874         4.225688066
+            N     1.276969829         0.737258874         2.642950986
+            B     0.000000000         1.474517748         2.112260792
+            N     0.000000000         1.474517748         0.529523459
+            '''
+cell.pseudo = 'gth-pade'
+cell.basis = 'gth-szv'
+cell.build()
+
+# Write cell geometry to file
+# - Format is guessed from filename
+# - These can be read by VESTA, Avogadro, etc.
+# - XYZ is Extended XYZ file format, which includes lattice vectors in the
+#   comment line
+cell.tofile('bn.vasp')
+cell.tofile('POSCAR')
+cell.tofile('bn.xyz')
+
+# Read a and atom from file
+from pyscf.pbc.gto.cell import fromfile
+a, atom = fromfile('bn.vasp')
+a, atom = fromfile('bn.xyz')
+
+# Read a and atom from file directly into Cell
+cell = gto.Cell()
+cell.fromfile('bn.vasp')
+cell.pseudo = 'gth-pade'
+cell.basis = 'gth-szv'
+cell.build()

examples/pbc/09-band_ase.py

@@ -4,7 +4,7 @@
 
 import matplotlib.pyplot as plt
 
-from ase.lattice import bulk
+from ase.build import bulk
 from ase.dft.kpoints import sc_special_points as special_points, get_bandpath
 
 c = bulk('C', 'diamond', a=3.5668)

examples/solvent/05-pcm.py

@@ -63,4 +63,13 @@
 td.kernel()
 
 mf = mol.RHF().PCM().run()
-td = mf.TDA().run()
+td = mf.TDA().run()
+
+# infinite epsilon with any PCM computations
+cm = pcm.PCM(mol)
+cm.eps = float('inf') # ideal conductor
+mf = dft.RKS(mol, xc='b3lyp')
+mf = mf.PCM(cm)
+mf.kernel()
+
+

examples/solvent/07-cosmors.py

@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+'''
+An example of using COSMO-RS functionality.
+'''
+
+#%% Imports
+
+import io
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pyscf import gto, dft
+from pyscf.solvent import pcm
+from pyscf.solvent.cosmors import get_sas_volume
+from pyscf.solvent.cosmors import write_cosmo_file
+from pyscf.solvent.cosmors import get_cosmors_parameters
+
+
+#%% Set parameters
+
+# get molecule
+coords = '''
+C        0.000000    0.000000             -0.542500
+O        0.000000    0.000000              0.677500
+H        0.000000    0.9353074360871938   -1.082500
+H        0.000000   -0.9353074360871938   -1.082500
+'''
+mol = gto.M(atom=coords, basis='6-31G*', verbose=1)
+
+# configure PCM
+cm = pcm.PCM(mol)
+cm.eps = float('inf') # f_epsilon = 1 is required for COSMO-RS
+cm.method = 'C-PCM' # or COSMO, IEF-PCM, SS(V)PE, see https://manual.q-chem.com/5.4/topic_pcm-em.html
+cm.lebedev_order = 29 # lebedev grids on the cavity surface, lebedev_order=29  <--> # of grids = 302
+
+
+#%% COSMO-files
+
+# run DFT SCF (any level of theory is OK, though DFT is optimal)
+mf = dft.RKS(mol, xc='b3lyp')
+mf = mf.PCM(cm)
+mf.kernel()
+
+# generate COSMO-file
+with io.StringIO() as outp: # with open('formaldehyde.cosmo', 'w') as outf:
+    write_cosmo_file(outp, mf)
+    print(outp.getvalue())
+
+
+# if PCM DFT were computed with fepsilon < 1 <=> eps != inf the ValueError will be raised
+# use ignore_low_feps=True to overrule it, but please be sure you know what you're doing
+
+# run DFT SCF
+cm = pcm.PCM(mol)
+cm.eps = 32.613 # methanol
+mf = dft.RKS(mol, xc='b3lyp')
+mf = mf.PCM(cm)
+mf.kernel()
+
+# try to get COSMO-file
+with io.StringIO() as outp:
+    try:
+        write_cosmo_file(outp, mf)
+        print(outp.getvalue())
+    except ValueError as e:
+        print(e)
+
+# overruling
+with io.StringIO() as outp:
+    write_cosmo_file(outp, mf, ignore_low_feps=True)
+    print(outp.getvalue())
+
+
+# The molecular volume is computed for the solvent-accessible surface generated
+# within pcm.PCM object. Please note that r_sol is assumed to be equal to zero,
+# so that SAS is a union of vdW spheres.
+# Increasing integration step increases accuracy of integration, but for most COSMO-RS
+# modelling and ML step=0.2 should be enough:
+
+# compute volumes with different step values
+steps = [1.0, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01]
+Vs = [get_sas_volume(mf.with_solvent.surface, step) for step in steps]
+# plot
+ax = plt.gca()
+ax.plot(Vs)
+ax.set_xticks(range(len(steps)))
+_ = ax.set_xticklabels(steps)
+plt.show()
+
+
+#%% Sigma-profiles
+
+# compute SCF PCM
+cm = pcm.PCM(mol)
+cm.eps = float('inf')
+mf = dft.RKS(mol, xc='b3lyp')
+mf = mf.PCM(cm)
+mf.kernel()
+
+# compute sigma-profile and related parameters
+params = get_cosmors_parameters(mf)
+print(params)
+plt.plot(params['Screening charge, e/A**2'],
+         params['Screening charge density, A**2'])
+plt.show()
+
+# custom sigma grid
+sigmas_grid = np.linspace(-0.025, 0.025, 101)
+params = get_cosmors_parameters(mf, sigmas_grid)
+plt.plot(params['Screening charge, e/A**2'],
+         params['Screening charge density, A**2'])
+plt.show()
+
+

pyproject.toml

@@ -57,8 +57,9 @@ cppe = ["cppe"]
 pyqmc = ["pyqmc"]
 bse = ["basis-set-exchange"]
 dispersion = ["pyscf-dispersion"]
+ccpy = ["coupled-cluster-py"]
 
-all = ["pyscf[forge,geomopt,doci,properties,semiempirical,cppe,pyqmc,bse,dispersion]"]
+all = ["pyscf[forge,geomopt,doci,properties,semiempirical,cppe,pyqmc,bse,dispersion,ccpy]"]
 
 # extras which should not be installed by "all" components
 cornell_shci = ["pyscf-cornell-shci"]

pyscf/ao2mo/__init__.py

@@ -31,6 +31,7 @@
 import tempfile
 import numpy
 import h5py
+from pyscf import gto
 from pyscf.ao2mo import incore
 from pyscf.ao2mo import outcore
 from pyscf.ao2mo import r_outcore
@@ -143,7 +144,7 @@ def full(eri_or_mol, mo_coeff, erifile=None, dataname='eri_mo', intor='int2e',
     '''
     if isinstance(eri_or_mol, numpy.ndarray):
         return incore.full(eri_or_mol, mo_coeff, *args, **kwargs)
-    else:
+    elif isinstance(eri_or_mol, gto.MoleBase):
         if '_spinor' in intor:
             mod = r_outcore
         else:
@@ -157,6 +158,11 @@ def full(eri_or_mol, mo_coeff, erifile=None, dataname='eri_mo', intor='int2e',
                             *args, **kwargs)
         else:
             return mod.full_iofree(eri_or_mol, mo_coeff, intor, *args, **kwargs)
+    else:
+        raise RuntimeError('ERI is not available. If this is generated by mf._eri, '
+                           'the integral tensor is too big to store in memory. '
+                           'You should either increase mol.max_memory, or set '
+                           'mol.incore_anyway. See issue #2473.')
 
 def general(eri_or_mol, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e',
             *args, **kwargs):
@@ -293,7 +299,7 @@ def general(eri_or_mol, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e
     '''
     if isinstance(eri_or_mol, numpy.ndarray):
         return incore.general(eri_or_mol, mo_coeffs, *args, **kwargs)
-    else:
+    elif isinstance(eri_or_mol, gto.MoleBase):
         if '_spinor' in intor:
             mod = r_outcore
         else:
@@ -307,6 +313,11 @@ def general(eri_or_mol, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e
                                *args, **kwargs)
         else:
             return mod.general_iofree(eri_or_mol, mo_coeffs, intor, *args, **kwargs)
+    else:
+        raise RuntimeError('ERI is not available. If this is generated by mf._eri, '
+                           'the integral tensor is too big to store in memory. '
+                           'You should either increase mol.max_memory, or set '
+                           'mol.incore_anyway. See issue #2473.')
 
 def kernel(eri_or_mol, mo_coeffs, erifile=None, dataname='eri_mo', intor='int2e',
            *args, **kwargs):

pyscf/ao2mo/_ao2mo.py

@@ -73,7 +73,7 @@ def nr_e1fill(intor, sh_range, atm, bas, env,
     natm = ctypes.c_int(c_atm.shape[0])
     nbas = ctypes.c_int(c_bas.shape[0])
     ao_loc = make_loc(bas, intor)
-    nao = ao_loc[-1]
+    nao = int(ao_loc[-1])
 
     klsh0, klsh1, nkl = sh_range
 

pyscf/ao2mo/nrr_outcore.py

@@ -380,9 +380,9 @@ def _count_naopair(mol, nao):
     ao_loc = mol.ao_loc_2c()
     nao_pair = 0
     for i in range(mol.nbas):
-        di = ao_loc[i+1] - ao_loc[i]
+        di = int(ao_loc[i+1] - ao_loc[i])
         for j in range(i+1):
-            dj = ao_loc[j+1] - ao_loc[j]
+            dj = int(ao_loc[j+1] - ao_loc[j])
             nao_pair += di * dj
     return nao_pair
 

pyscf/ao2mo/outcore.py

@@ -710,8 +710,8 @@ def guess_shell_ranges(mol, aosym, max_iobuf, max_aobuf=None, ao_loc=None,
                        compress_diag=True):
     if ao_loc is None: ao_loc = mol.ao_loc_nr()
     max_iobuf = max(1, max_iobuf)
-
-    dims = ao_loc[1:] - ao_loc[:-1]
+    ao_loc_long = ao_loc.astype(numpy.int64)
+    dims = ao_loc_long[1:] - ao_loc_long[:-1]
     dijs = (dims.reshape(-1,1) * dims)
     nbas = dijs.shape[0]
 
@@ -773,7 +773,7 @@ def balance_partition(ao_loc, blksize, start_id=0, stop_id=None):
     displs = [i+start_id for i in displs]
     tasks = []
     for i0, i1 in zip(displs[:-1],displs[1:]):
-        tasks.append((i0, i1, ao_loc[i1]-ao_loc[i0]))
+        tasks.append((i0, i1, int(ao_loc[i1]-ao_loc[i0])))
     return tasks
 
 del (MAX_MEMORY)

pyscf/ao2mo/r_outcore.py

@@ -302,9 +302,9 @@ def _count_naopair(mol, nao):
     ao_loc = mol.ao_loc_2c()
     nao_pair = 0
     for i in range(mol.nbas):
-        di = ao_loc[i+1] - ao_loc[i]
+        di = int(ao_loc[i+1] - ao_loc[i])
         for j in range(i+1):
-            dj = ao_loc[j+1] - ao_loc[j]
+            dj = int(ao_loc[j+1] - ao_loc[j])
             nao_pair += di * dj
     return nao_pair