From 84202739146cb307951c3d27b62a43c3a02f2193 Mon Sep 17 00:00:00 2001 From: Xiaojie Wu Date: Sun, 17 Nov 2024 17:19:02 -0800 Subject: [PATCH 1/4] Update the blacklist of dispersion correction --- pyscf/scf/dispersion.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyscf/scf/dispersion.py b/pyscf/scf/dispersion.py index 6bc02bb2a9..1c4a86f01c 100644 --- a/pyscf/scf/dispersion.py +++ b/pyscf/scf/dispersion.py @@ -50,7 +50,7 @@ # These xc functionals are not supported yet _black_list = { - 'wb97x-d', 'wb97x-d3', + 'wb97x-d', 'wb97x-d3', 'wb97x_d', 'wb97x_d3', 'wb97m-d3bj2b', 'wb97m-d3bjatm', 'b97m-d3bj2b', 'b97m-d3bjatm', } From 612d95656acabbae21d9ca9525ffda7e242c0adc Mon Sep 17 00:00:00 2001 From: Christopher Hillenbrand Date: Sun, 24 Nov 2024 22:56:13 -0500 Subject: [PATCH 2/4] some lightweight array helper functions (#2515) * some lightweight array helper functions * satisfy linter * make requested changes * Add out argument for entrywise_mul * restore lib.zeros for testing --- pyscf/lib/np_helper/np_helper.c | 94 ++++++++++++++++++++ pyscf/lib/np_helper/np_helper.h | 18 ++++ pyscf/lib/numpy_helper.py | 130 +++++++++++++++++++++++++++- pyscf/lib/test/test_numpy_helper.py | 36 ++++++++ 4 files changed, 274 insertions(+), 4 deletions(-) diff --git a/pyscf/lib/np_helper/np_helper.c b/pyscf/lib/np_helper/np_helper.c index a997a24223..d6c6b96848 100644 --- a/pyscf/lib/np_helper/np_helper.c +++ b/pyscf/lib/np_helper/np_helper.c @@ -14,6 +14,7 @@ */ #include +#include #include "np_helper/np_helper.h" void NPdset0(double *p, const size_t n) @@ -47,3 +48,96 @@ void NPzcopy(double complex *out, const double complex *in, const size_t n) out[i] = in[i]; } } + +/* + * These are mostly useful for first-touch array allocation on NUMA systems. + * Use with numpy.empty. + */ +void NPomp_dset0(const size_t n, double *out) +{ +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < n; i++) { + out[i] = 0.0; + } +} + +void NPomp_zset0(const size_t n, double complex *out) +{ +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < n; i++) { + out[i] = 0.0; + } +} + + +/* + * Copy a double precision matrix with multithreading. + */ +void NPomp_dcopy(const size_t m, + const size_t n, + const double *__restrict in, const size_t in_stride, + double *__restrict out, const size_t out_stride) +{ +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < m; i++) { +#pragma omp simd + for (size_t j = 0; j < n; j++) { + out[i * out_stride + j] = in[i * in_stride + j]; + } + } +} + +/* + * Copy a complex double precision matrix with multithreading. + */ +void NPomp_zcopy(const size_t m, + const size_t n, + const double complex *__restrict in, const size_t in_stride, + double complex *__restrict out, const size_t out_stride) +{ +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < m; i++) { +#pragma omp simd + for (size_t j = 0; j < n; j++) { + out[i * out_stride + j] = in[i * in_stride + j]; + } + } +} + +/* + * Elementwise multiplication of two double matrices. + * B <- A \circ B + */ +void NPomp_dmul(const size_t m, + const size_t n, + const double *__restrict a, const size_t a_stride, + double *__restrict b, const size_t b_stride, + double *__restrict out, const size_t out_stride) +{ +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < m; i++) { +#pragma omp simd + for (size_t j = 0; j < n; j++) { + out[i * out_stride + j] = b[i * b_stride + j] * a[i * a_stride + j]; + } + } +} + +/* + * Elementwise multiplication of two complex double matrices. + * B <- A \circ B + */ +void NPomp_zmul(const size_t m, + const size_t n, + const double complex *__restrict a, const size_t a_stride, + double complex *__restrict b, const size_t b_stride, + double complex *__restrict out, const size_t out_stride) +{ +#pragma omp parallel for schedule(static) + for (size_t i = 0; i < m; i++) { +#pragma omp simd + for (size_t j = 0; j < n; j++) { + out[i * out_stride + j] = b[i * b_stride + j] * a[i * a_stride + j]; + } + } +} diff --git a/pyscf/lib/np_helper/np_helper.h b/pyscf/lib/np_helper/np_helper.h index e92e8ab957..5c7f834399 100644 --- a/pyscf/lib/np_helper/np_helper.h +++ b/pyscf/lib/np_helper/np_helper.h @@ -70,6 +70,24 @@ void NPzset0(double complex *p, const size_t n); void NPdcopy(double *out, const double *in, const size_t n); void NPzcopy(double complex *out, const double complex *in, const size_t n); +void NPomp_dset0(const size_t n, double *out); +void NPomp_zset0(const size_t n, double complex *out); + +void NPomp_dcopy(const size_t m, const size_t n, + const double *in, const size_t in_stride, + double *out, const size_t out_stride); +void NPomp_zcopy(const size_t m, const size_t n, + const double complex *in, const size_t in_stride, + double complex *out, const size_t out_stride); +void NPomp_dmul(const size_t m, const size_t n, + const double *a, const size_t a_stride, + double *b, const size_t b_stride, + double *out, const size_t out_stride); +void NPomp_zmul(const size_t m, const size_t n, + const double complex *a, const size_t a_stride, + double complex *b, const size_t b_stride, + double complex *out, const size_t out_stride); + void NPdgemm(const char trans_a, const char trans_b, const int m, const int n, const int k, const int lda, const int ldb, const int ldc, diff --git a/pyscf/lib/numpy_helper.py b/pyscf/lib/numpy_helper.py index 020aa64a8d..3b4a5e531f 100644 --- a/pyscf/lib/numpy_helper.py +++ b/pyscf/lib/numpy_helper.py @@ -552,12 +552,10 @@ def inplace_transpose_scale(a, alpha=1.0): alpha : float, optional scaling factor, by default 1.0 """ - assert a.ndim == 2 + lda, order, _ = leading_dimension_order(a) assert a.shape[0] == a.shape[1] n = a.shape[0] - astrides = [s // a.itemsize for s in a.strides] - lda = max(astrides) - assert min(astrides) == 1 + assert order in ('C', 'F') if a.dtype == numpy.double: _np_helper.NPomp_d_itranspose_scale( ctypes.c_int(n), ctypes.c_double(alpha), @@ -974,6 +972,37 @@ def frompointer(pointer, count, dtype=float): a = numpy.ndarray(count, dtype=numpy.int8, buffer=buf) return a.view(dtype) +def leading_dimension_order(a): + """Return the leading dimension and the order of a matrix. + + Parameters + ---------- + a : ndarray + 2D array. + + Returns + ------- + lda : int + Leading dimension of the array -- the stride between rows or columns. + order : str + 'F' for col major, 'C' for row major, 'G' for neither. + a_cshape : tuple + If a is row major, a.shape; if a is col major, a.T.shape; otherwise None. + """ + assert a.ndim == 2 + astrides = [s//a.itemsize for s in a.strides] + lda = max(astrides) + if astrides[0] == 1: + order = 'F' + a_cshape = a.T.shape + elif astrides[1] == 1: + order = 'C' + a_cshape = a.shape + else: + order = 'G' + a_cshape = None + return lda, order, a_cshape + norm = numpy.linalg.norm def cond(x, p=None): @@ -1177,6 +1206,99 @@ def expm(a): y, buf = buf, y return y +def omatcopy(a, out=None): + """Copies a matrix. + + Parameters + ---------- + a : ndarray + Matrix to be copied. The order of the matrix is preserved. + a can be either row or column major. + out : ndarray, optional + Matrix to be overwritten. A new one is allocated if not provided. + + Returns + ------- + out : ndarray + Copy of a with the same order. + """ + lda, _, a_cshape = leading_dimension_order(a) + if out is None: + out = numpy.empty_like(a) + ld_out, _, out_cshape = leading_dimension_order(out) + assert out_cshape == a_cshape and a_cshape is not None + if a.dtype == numpy.double: + fn = _np_helper.NPomp_dcopy + elif a.dtype == numpy.complex128: + fn = _np_helper.NPomp_zcopy + else: + raise NotImplementedError + fn(ctypes.c_size_t(a_cshape[0]), + ctypes.c_size_t(a_cshape[1]), + a.ctypes.data_as(ctypes.c_void_p), + ctypes.c_size_t(lda), + out.ctypes.data_as(ctypes.c_void_p), + ctypes.c_size_t(ld_out)) + return out + +def zeros(shape, dtype=numpy.double, order='C'): + """Allocate and zero an array in parallel. Useful for multi-socket systems + due to the first touch policy. On most systems np.zeros does not count + as first touch. Arrays returned by this function will (ideally) have + pages backing them that are distributed across the sockets. + """ + dtype = numpy.dtype(dtype) + if dtype == numpy.double: + out = numpy.empty(shape, dtype=dtype, order=order) + _np_helper.NPomp_dset0(ctypes.c_size_t(out.size), + out.ctypes.data_as(ctypes.c_void_p)) + elif dtype == numpy.complex128: + out = numpy.empty(shape, dtype=dtype, order=order) + _np_helper.NPomp_zset0(ctypes.c_size_t(out.size), + out.ctypes.data_as(ctypes.c_void_p)) + else: # fallback + out = numpy.zeros(shape, dtype=dtype, order=order) + return out + +def entrywise_mul(a, b, out=None): + """Entrywise multiplication of two matrices. + + Parameters + ---------- + a : ndarray + b : ndarray + out : ndarray, optional + Output matrix. A new one is allocated if not provided. + + Returns + ------- + ndarray + a * b + """ + assert a.ndim == 2 and b.ndim == 2 + assert a.shape == b.shape and a.dtype == b.dtype + lda, _, a_cshape = leading_dimension_order(a) + ldb, _, b_cshape = leading_dimension_order(b) + if out is None: + out = numpy.empty_like(b) + ld_out, _, out_cshape = leading_dimension_order(out) + assert a_cshape == b_cshape and b_cshape == out_cshape and a_cshape is not None + if a.dtype == numpy.double: + fn = _np_helper.NPomp_dmul + elif a.dtype == numpy.complex128: + fn = _np_helper.NPomp_zmul + else: + return numpy.multiply(a, b, out=out) + fn(ctypes.c_size_t(a_cshape[0]), + ctypes.c_size_t(a_cshape[1]), + a.ctypes.data_as(ctypes.c_void_p), + ctypes.c_size_t(lda), + b.ctypes.data_as(ctypes.c_void_p), + ctypes.c_size_t(ldb), + out.ctypes.data_as(ctypes.c_void_p), + ctypes.c_size_t(ld_out)) + return out + def ndarray_pointer_2d(array): '''Return an array that contains the addresses of the first element in each row of the input 2d array. diff --git a/pyscf/lib/test/test_numpy_helper.py b/pyscf/lib/test/test_numpy_helper.py index dce15c50ba..fdb088dc53 100644 --- a/pyscf/lib/test/test_numpy_helper.py +++ b/pyscf/lib/test/test_numpy_helper.py @@ -254,6 +254,42 @@ def test_ndarray_pointer_2d(self): addr = lib.ndarray_pointer_2d(a) self.assertTrue(all(addr == a.ctypes.data + numpy.array([0, 24, 48]))) + def test_omatcopy(self): + a = numpy.random.random((5,5)) + b = numpy.empty_like(a) + lib.omatcopy(a, out=b) + self.assertTrue(numpy.all(a == b)) + a = numpy.random.random((403,410)).T + b = numpy.empty_like(a) + lib.omatcopy(a, out=b) + self.assertTrue(numpy.all(a == b)) + + def test_zeros(self): + a = lib.zeros((100,100), dtype=numpy.double) + self.assertTrue(numpy.all(a == 0)) + self.assertTrue(a.dtype == numpy.double) + a = lib.zeros((100,100), dtype=numpy.complex128) + self.assertTrue(numpy.all(a == 0)) + self.assertTrue(a.dtype == numpy.complex128) + a = lib.zeros((100,100), dtype=numpy.int32) + self.assertTrue(numpy.all(a == 0)) + self.assertTrue(a.dtype == numpy.int32) + + def test_entrywise_mul(self): + a = numpy.random.random((101,100)) + b = numpy.random.random((101,100)) + prod = lib.entrywise_mul(a, b) + self.assertTrue(numpy.allclose(prod, a * b)) + a = numpy.random.random((101,100)) + b = numpy.random.random((101,100)) + a = a + a*1j + b = b + b*1j + prod = lib.entrywise_mul(a, b) + self.assertTrue(numpy.allclose(prod, a * b)) + # inplace test + lib.entrywise_mul(a, b, out=b) + self.assertTrue(numpy.allclose(prod, b)) + if __name__ == "__main__": print("Full Tests for numpy_helper") unittest.main() From 8215a21c8ecebe3d57e3a9e06b68fcb8074f41cc Mon Sep 17 00:00:00 2001 From: Thijs Vogels Date: Mon, 25 Nov 2024 05:00:02 +0100 Subject: [PATCH 3/4] Fix biased implementation for the becke radi method (#2499) * Fix biased implementation for the becke radi method This fixed the Becke radial integration scheme as described in #2490. * Update test_grids.py * Update radi.py --------- Co-authored-by: Qiming Sun --- pyscf/dft/radi.py | 11 +++++++++-- pyscf/dft/test/test_grids.py | 2 +- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/pyscf/dft/radi.py b/pyscf/dft/radi.py index 718aea118c..43bfb7b71d 100644 --- a/pyscf/dft/radi.py +++ b/pyscf/dft/radi.py @@ -47,14 +47,21 @@ def murray(n, *args, **kwargs): raise RuntimeError('Not implemented') -# Gauss-Chebyshev of the first kind, and the transformed interval [0,\infty) def becke(n, charge, *args, **kwargs): '''Becke, JCP 88, 2547 (1988); DOI:10.1063/1.454033''' if charge == 1: rm = BRAGG_RADII[charge] else: rm = BRAGG_RADII[charge] * .5 - t, w = numpy.polynomial.chebyshev.chebgauss(n) + + # Points and weights for Gauss-Chebyshev quadrature points of the second kind + # The weights are adjusted to integrate a function on the interval [-1, 1] with uniform weighting + # instead of weighted by sqrt(1 - t^2) = sin(i pi / (n+1)). + i = numpy.arange(n) + 1 + t = numpy.cos(i * numpy.pi / (n + 1)) + w = numpy.pi / (n + 1) * numpy.sin(i * numpy.pi / (n + 1)) + + # Change of variables to map the domain to [0, inf) r = (1+t)/(1-t) * rm w *= 2/(1-t)**2 * rm return r, w diff --git a/pyscf/dft/test/test_grids.py b/pyscf/dft/test/test_grids.py index ba261f2f36..33dd9c02ec 100644 --- a/pyscf/dft/test/test_grids.py +++ b/pyscf/dft/test/test_grids.py @@ -91,7 +91,7 @@ def test_radi(self): grid.radi_method = radi.becke grid.build(with_non0tab=False) - self.assertAlmostEqual(lib.fp(grid.weights), 818061.875131255, 7) + self.assertAlmostEqual(lib.fp(grid.weights), 780.7183109298, 9) def test_prune(self): grid = gen_grid.Grids(h2o) From 9f8fe3021ec5eeaa66c0a5b14b5e5de32a2a35dc Mon Sep 17 00:00:00 2001 From: matthew-hennefarth Date: Sun, 24 Nov 2024 22:05:08 -0600 Subject: [PATCH 4/4] Verify Examples (#2379) * add verify_examples * make executable * minor fix relating to cwd * add doc string and usage example --- tools/verify_examples.py | 257 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 257 insertions(+) create mode 100755 tools/verify_examples.py diff --git a/tools/verify_examples.py b/tools/verify_examples.py new file mode 100755 index 0000000000..3f9d484c5e --- /dev/null +++ b/tools/verify_examples.py @@ -0,0 +1,257 @@ +#!/usr/bin/env python +# Copyright 2014-2024 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 +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Verify Examples +=============== + +Author: Matthew R. Hennefarth + +Script used to automatically run and verify PySCF example codes terminate +successfully. For any job that does not terminate normally, the stderr of the +example will be printed to the output. This script will exit with 0 only if all +examples terminate normally. + +Initially introduced in [PR 2379](https://github.com/pyscf/pyscf/pull/2379). + +Usage +------------- + +From the main pyscf repository directory, the tests can be run as +```sh +./tools/verify_examples.py examples +``` +This will run all example files (which can be very long). To run only a subset +of examples, provide instead a path to a subdirectory. For example, to run only +the example files in `pyscf/examples/gto` the command +```sh +./tools/verify_examples.py examples/gto +``` +It is also possible to run the examples in parallel using the `-j` or `--jobs` +flag (this is similar to make). As an example, to run the jobs in parallel over +4 threads, +```sh +./tools/verify_examples.py -j 8 +``` +Note that the environmental variable such as `OMP_NUM_THREADS` should be set to +an appropriate value such that number of jobs * OMP_NUM_THREADS does not exceed +the maximum number of cores on the computer. + +""" + +import os +import sys +import time +import subprocess +import argparse +import logging + +import multiprocessing as mp +from glob import glob +from enum import Enum + +logging.basicConfig(level=logging.DEBUG, format="%(message)s") + +logger = logging.getLogger() + + +class StdOutFilter(logging.Filter): + def filter(self, record): + return record.levelno < logging.ERROR + + +stdout_handler = logging.StreamHandler(sys.stdout) +stdout_handler.setLevel(logging.INFO) +stdout_handler.addFilter(StdOutFilter()) + +stderr_handler = logging.StreamHandler(sys.stderr) +stderr_handler.setLevel(logging.ERROR) + +logger.handlers = [] +logger.addHandler(stdout_handler) +logger.addHandler(stderr_handler) + + +class ANSIColors(Enum): + RESET = "\033[0m" + RED = "\033[31m" + GREEN = "\033[32m" + + +def colorize(text, color): + if sys.stdout.isatty(): + return f"\033[{color.value}{text}{ANSIColors.RESET.value}" + else: + return text + + +class Status(Enum): + OK = colorize("ok", ANSIColors.GREEN) + FAIL = colorize("FAILED", ANSIColors.RED) + + +def get_path(p): + if not os.path.isdir(p): + raise ValueError("Path does not point to directory") + + if os.path.basename(p) == "examples": + return p + + if os.path.isdir(os.path.join(p, "examples")): + return os.path.join(p, "examples") + + return p + + +class ExampleResults: + def __init__(self): + self.common_prefix = "" + self.failed_examples = [] + self.passed = 0 + self.failed = 0 + self.filtered = 0 + self.time = 0.0 + self.status = Status.OK + + +def run_example(progress, nexamples, example, failed_examples, common_prefix): + idx, lock = progress + + status = Status.OK + directory = os.path.dirname(example) + try: + subprocess.run( + ["python3", os.path.basename(example)], + cwd=directory, + capture_output=False, + stderr=subprocess.PIPE, + stdout=subprocess.DEVNULL, + check=True, + text=True, + ) + except subprocess.CalledProcessError as e: + status = Status.FAIL + failed_examples.append((example, e.stderr)) + + with lock: + idx.value += 1 + percent = int(100 * (idx.value) / nexamples) + + message = ( + f"[{percent:3}%]: {os.path.relpath(example, common_prefix)} ... {status.value}" + ) + logger.info(message) + + +def run_examples(example_path, num_threads): + examples = [ + y for x in os.walk(example_path) for y in glob(os.path.join(x[0], "*.py")) + ] + # remove symlinks? + # examples = list(set([os.path.realpath(e) for e in examples])) + + examples = sorted(examples, key=lambda e: e.split("/")) + + results = ExampleResults() + results.common_prefix = os.path.dirname(os.path.commonpath(examples)) + results.filtered = 0 + + with mp.Manager() as manager: + failed_examples = manager.list() + progress = (manager.Value("i", 0), manager.Lock()) + + logger.info("") + logger.info(f"running {len(examples)} examples") + tic = time.perf_counter() + with mp.Pool(num_threads) as pool: + pool.starmap( + run_example, + [ + ( + progress, + len(examples), + example, + failed_examples, + results.common_prefix, + ) + for example in examples + ], + ) + results.time = time.perf_counter() - tic + results.failed_examples = list(failed_examples) + + results.failed = len(results.failed_examples) + results.passed = len(examples) - results.failed + results.status = Status.FAIL if results.failed else Status.OK + + return results + + +def log_failures(results): + logger.info("") + logger.info("failures: ") + logger.info("") + + for e, msg in results.failed_examples: + logger.info(f"---- {os.path.relpath(e, results.common_prefix)} stderr ----") + logger.info(msg) + + logger.info("") + logger.info("failures:") + for e, _ in results.failed_examples: + logger.info(f" {os.path.relpath(e, results.common_prefix)}") + + +def main(): + parser = argparse.ArgumentParser(description="Verify pyscf examples") + parser.add_argument( + "path", + type=str, + default="examples", + help="Path to examples directory (default: ./)", + ) + parser.add_argument( + "-j", + "--jobs", + type=int, + default=1, + help="Number of parallel threads (default: 1)", + ) + args = parser.parse_args() + + example_path = get_path(args.path) + + results = run_examples(example_path, args.jobs) + + if results.status is Status.FAIL: + log_failures(results) + + logger.info("") + logger.info( + f"example results: {results.status.value}. {results.passed} passed; {results.failed} failed; {results.filtered} filtered out; finished in {results.time:.2f}s" + ) + logger.info("") + + if results.status is Status.OK: + sys.exit(0) + else: + logger.error( + f"{ANSIColors.RED.value}error{ANSIColors.RESET.value}: examples failed" + ) + sys.exit(1) + + +if __name__ == "__main__": + main()