add pyscf interface and lcao convolutions for NLDF

.github/workflows/full_env.yml

@@ -15,6 +15,7 @@ dependencies:
   - numba<=0.58
   - numpy>=1.20.3,<1.25
   - pyscf>=2.1
+  - pytest
   - pyyaml
   - scikit-learn>=1.0.1
   - scipy>=1.7.1

.github/workflows/nocomp_env.yml

@@ -11,6 +11,7 @@ dependencies:
   - numba<=0.58
   - numpy>=1.20.3,<1.25
   - pyscf>=2.1
+  - pytest
   - pyyaml
   - scikit-learn>=1.0.1
   - scipy>=1.7.1

ciderpress/density.py

@@ -0,0 +1,145 @@
+#!/usr/bin/env python
+# CiderPress: Machine-learning based density functional theory calculations
+# Copyright (C) 2024 The President and Fellows of Harvard College
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>
+#
+# Author: Kyle Bystrom <kylebystrom@gmail.com>
+#
+
+import numpy
+import numpy as np
+
+LDA_FACTOR = -3.0 / 4.0 * (3.0 / np.pi) ** (1.0 / 3)
+GG_SMUL = 1.0
+GG_AMUL = 1.0
+GG_AMIN = 1.0 / 18
+CFC = (3.0 / 10) * (3 * np.pi**2) ** (2.0 / 3)
+DESC_VERSION_LIST = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k"]
+
+UEG_VECTORS = {
+    "b": np.array([1.0, 0.0, 1.0, 2.0, 2.0, 2.0, 0.0, 0.0]),
+    "d": np.array([1.0, 0.0, 1.0, 2.0, 2.0, 2.0, 0.0, 0.0]),
+}
+
+
+def ldax(n):
+    """
+    LDA exchange energy density
+    Args:
+        n: Density
+    """
+    return LDA_FACTOR * n ** (4.0 / 3)
+
+
+def ldaxp(n):
+    """
+    Fully spin-polarized LDA exchange energy density
+    Args:
+        n: Density
+    """
+    return 0.5 * ldax(2 * n)
+
+
+def lsda(nu, nd):
+    """
+    LSDA exchange energy density
+    Args:
+        nu: Spin-up Density
+        nd: Spin-down Density
+    """
+    return ldaxp(nu) + ldaxp(nd)
+
+
+def get_ldax_dens(n):
+    """
+    LDA exchange energy density
+    Args:
+        n: Density
+    """
+    return LDA_FACTOR * n ** (4.0 / 3)
+
+
+def get_ldax(n):
+    """
+    LDA Exchange energy per particle
+    Args:
+        n: Density
+    """
+    return LDA_FACTOR * n ** (1.0 / 3)
+
+
+def get_xed_from_y(y, rho):
+    """
+    Get the exchange energy density (n * epsilon_x)
+    from the exchange enhancement factor y
+    and density rho.
+    """
+    return rho * get_x(y, rho)
+
+
+def get_x(y, rho):
+    return (y + 1) * get_ldax(rho)
+
+
+def get_y_from_xed(xed, rho):
+    """
+    Get the exchange enhancement factor minus one.
+    """
+    return xed / (get_ldax_dens(rho) - 1e-12) - 1
+
+
+def get_gradient_magnitude(rho_data):
+    return np.linalg.norm(rho_data[1:4, :], axis=0)
+
+
+def get_normalized_grad(rho, mag_grad):
+    sprefac = 2 * (3 * np.pi * np.pi) ** (1.0 / 3)
+    n43 = rho ** (4.0 / 3)
+    s = mag_grad / (sprefac * n43 + 1e-16)
+    return s
+
+
+def get_single_orbital_tau(rho, mag_grad):
+    return mag_grad**2 / (8 * rho + 1e-16)
+
+
+def get_uniform_tau(rho):
+    return (3.0 / 10) * (3 * np.pi**2) ** (2.0 / 3) * rho ** (5.0 / 3)
+
+
+def get_regularized_tau(tau, tau_w, tau_unif):
+    alpha = (tau - tau_w) / (tau_unif + 1e-4)
+    return alpha**3 / (alpha**2 + 1e-3)
+
+
+def get_normalized_tau(tau, tau_w, tau_unif):
+    return (tau - tau_w) / (tau_unif + 1e-16)
+
+
+def get_dft_input(rho_data):
+    rho = rho_data[0, :]
+    mag_grad = get_gradient_magnitude(rho_data)
+    s = get_normalized_grad(rho, mag_grad)
+    tau_w = get_single_orbital_tau(rho, mag_grad)
+    tau_unif = get_uniform_tau(rho)
+    alpha = get_normalized_tau(rho_data[5], tau_w, tau_unif)
+    return rho, s, alpha, tau_w, tau_unif
+
+
+def get_dft_input_gga(rho_data):
+    rho = rho_data[0, :]
+    mag_grad = get_gradient_magnitude(rho_data)
+    s = get_normalized_grad(rho, mag_grad)
+    return rho, s

ciderpress/dft/baselines.py

@@ -0,0 +1,118 @@
+import numpy as np
+
+from ciderpress.density import LDA_FACTOR
+
+
+def nsp_rho_basline(X0T, e=None, dedx=None):
+    nspin, nfeat, nsamp = X0T.shape
+    if e is None:
+        e = np.zeros(nsamp)
+    if dedx is None:
+        dedx = np.zeros_like(X0T)
+    nspin = X0T.shape[0]
+    e[:] += X0T[:, 0].mean(0)
+    dedx[0, :] += 1.0 / nspin
+
+
+def _lda_x_helper(X0T, e, dedx):
+    rho = X0T[0]
+    e[:] += LDA_FACTOR * rho ** (4.0 / 3)
+    dedx[0] += 4.0 / 3 * LDA_FACTOR * rho ** (1.0 / 3)
+
+
+def _vi_x_damp_helper(X0T, e, dedx):
+    rho = X0T[0]
+    damp = 2.0 / (1.0 + 0.5 * X0T[3]) ** 2
+    ddamp = -1 * damp / (1.0 + 0.5 * X0T[3])
+    e[:] += LDA_FACTOR * damp * rho ** (4.0 / 3)
+    dedx[0] += 4.0 / 3 * LDA_FACTOR * damp * rho ** (1.0 / 3)
+    dedx[3] += LDA_FACTOR * ddamp * rho ** (4.0 / 3)
+
+
+def _pbe_x_helper(X0T, e, dedx):
+    rho = X0T[0]
+    s2 = X0T[1]
+    kappa = 0.804
+    mu = 0.2195149727645171
+    mk = mu / kappa
+    fac = 1.0 / (1 + mk * s2)
+    fx = 1 + kappa - kappa * fac
+    dfx = mu * fac * fac
+    e[:] += LDA_FACTOR * rho ** (4.0 / 3) * fx
+    dedx[0] += 4.0 / 3 * LDA_FACTOR * rho ** (1.0 / 3) * fx
+    dedx[1] += LDA_FACTOR * rho ** (4.0 / 3) * dfx
+
+
+def _chachiyo_x_helper(X0T, e, dedx):
+    rho = X0T[0]
+    s2 = X0T[1]
+    c = 4 * np.pi / 9
+    x = c * np.sqrt(s2)
+    dx = c / (2 * np.sqrt(s2))
+    chfx = (3 * x**2 + np.pi**2 * np.log(1 + x)) / (
+        (np.pi**2 + 3 * x) * np.log(1 + x)
+    )
+    dchfx = (
+        -3 * x**2 * (np.pi**2 + 3 * x)
+        + 3 * x * (1 + x) * (2 * np.pi**2 + 3 * x) * np.log(1 + x)
+        - 3 * np.pi**2 * (1 + x) * np.log(1 + x) ** 2
+    ) / ((1 + x) * (np.pi**2 + 3 * x) ** 2 * np.log(1 + x) ** 2)
+    dchfx *= dx
+    chfx[s2 < 1e-8] = 1 + 8 * s2[s2 < 1e-8] / 27
+    dchfx[s2 < 1e-8] = 8.0 / 27
+    e[:] += LDA_FACTOR * rho ** (4.0 / 3) * chfx
+    dedx[0] += 4.0 / 3 * LDA_FACTOR * rho ** (1.0 / 3) * chfx
+    dedx[1] += LDA_FACTOR * rho ** (4.0 / 3) * dchfx
+
+
+def _sl_x_helper(X0T, _helper):
+    nspin, nfeat, nsamp = X0T.shape
+    e = np.zeros(nsamp)
+    dedx = np.zeros_like(X0T)
+    nspin = X0T.shape[0]
+    for s in range(nspin):
+        _helper(X0T[s], e, dedx[s])
+    e[:] /= nspin
+    dedx[:] /= nspin
+    return e, dedx
+
+
+def zero_xc(X0T):
+    nspin, nfeat, nsamp = X0T.shape
+    e = np.zeros(nsamp)
+    dedx = np.zeros_like(X0T)
+    return e, dedx
+
+
+def one_xc(X0T):
+    nspin, nfeat, nsamp = X0T.shape
+    e = np.ones(nsamp)
+    dedx = np.zeros_like(X0T)
+    return e, dedx
+
+
+def lda_x(X0T):
+    return _sl_x_helper(X0T, _lda_x_helper)
+
+
+def nlda_x_damp(X0T):
+    return _sl_x_helper(X0T, _vi_x_damp_helper)
+
+
+def gga_x_chachiyo(X0T):
+    return _sl_x_helper(X0T, _chachiyo_x_helper)
+
+
+def gga_x_pbe(X0T):
+    return _sl_x_helper(X0T, _pbe_x_helper)
+
+
+BASELINE_CODES = {
+    "RHO": nsp_rho_basline,
+    "ZERO": zero_xc,
+    "ONE": one_xc,
+    "LDA_X": lda_x,
+    "NLDA_X_DAMP": nlda_x_damp,
+    "GGA_X_PBE": gga_x_pbe,
+    "GGA_X_CHACHIYO": gga_x_chachiyo,
+}

ciderpress/dft/debug_numint.py

@@ -2,8 +2,8 @@
 
 import numpy as np
 
+from ciderpress.dft.plans import get_cider_exponent
 from ciderpress.lib import load_library as load_cider_library
-from ciderpress.new_dft.plans import get_cider_exponent
 
 libcider = load_cider_library("libmcider")