add more gpaw tests and move env yml files to avoid them being treate…

…d as ci

.github/workflows/mpi_tests.yml

@@ -26,7 +26,7 @@ jobs:
       uses: mamba-org/setup-micromamba@v1
       with:
         micromamba-version: '1.5.6-0' # any version from https://github.com/mamba-org/micromamba-releases
-        environment-file: .github/workflows/full_env.yml
+        environment-file: .github/full_env.yml
         environment-name: run_tests
         init-shell: >-
           bash

.github/workflows/tests.yml

@@ -26,7 +26,7 @@ jobs:
       uses: mamba-org/setup-micromamba@v1
       with:
         micromamba-version: '1.5.6-0' # any version from https://github.com/mamba-org/micromamba-releases
-        environment-file: .github/workflows/nocomp_env.yml
+        environment-file: .github/nocomp_env.yml
         environment-name: run_tests
         init-shell: >-
           bash

ciderpress/gpaw/tests/test_force.py

@@ -0,0 +1,107 @@
+import unittest
+
+from ase import Atoms
+from gpaw import GPAW, PW, FermiDirac
+
+from ciderpress.gpaw.cider_paw import get_cider_functional
+from ciderpress.gpaw.tests import equal
+
+
+def numeric_force(atoms, a, i, d=0.001, get_xc=None):
+    """Compute numeric force on atom with index a, Cartesian component i,
+    with finite step of size d
+    """
+    p0 = atoms.get_positions()
+    p = p0.copy()
+    p[a, i] += d
+    atoms.set_positions(p, apply_constraint=False)
+    # if get_xc is not None:
+    #    atoms.calc.set(xc=get_xc())
+    eplus = atoms.get_potential_energy()
+    p[a, i] -= 2 * d
+    atoms.set_positions(p, apply_constraint=False)
+    # if get_xc is not None:
+    #    atoms.calc.set(xc=get_xc())
+    eminus = atoms.get_potential_energy()
+    atoms.set_positions(p0, apply_constraint=False)
+    # if get_xc is not None:
+    #    atoms.calc.set(xc=get_xc())
+    return (eminus - eplus) / (2 * d)
+
+
+def test_cider_forces(functional, get_xc=None):
+    a = 5.45
+    bulk = Atoms(
+        symbols="Si8",
+        positions=[
+            (0, 0, 0.1 / a),
+            (0, 0.5, 0.5),
+            (0.5, 0, 0.5),
+            (0.5, 0.5, 0),
+            (0.25, 0.25, 0.25),
+            (0.25, 0.75, 0.75),
+            (0.75, 0.25, 0.75),
+            (0.75, 0.75, 0.25),
+        ],
+        pbc=True,
+    )
+    bulk.set_cell((a, a, a), scale_atoms=True)
+    if get_xc is not None:
+        functional = get_xc()
+    calc = GPAW(
+        h=0.15,
+        mode=PW(520),
+        # xc='PBE',
+        xc=functional,
+        nbands="150%",
+        occupations=FermiDirac(width=0.01),
+        # poissonsolver=PoissonSolver('fd', nn='M', relax='J'),
+        kpts=(4, 4, 4),
+        convergence={"energy": 1e-7},
+        parallel={"augment_grids": True},
+    )
+    bulk.calc = calc
+    f1 = bulk.get_forces()[0, 2]
+    bulk.get_potential_energy()
+
+    f2 = numeric_force(bulk, 0, 2, 0.001, get_xc=get_xc)
+    print((f1, f2, f1 - f2))
+    equal(f1, f2, 0.005)
+
+
+class TestForce(unittest.TestCase):
+    def test_sl_gga(self):
+        test_cider_forces(get_cider_functional("functionals/CIDER23X_SL_GGA.yaml"))
+
+    def test_nl_gga(self):
+        test_cider_forces(get_cider_functional("functionals/CIDER23X_SL_MGGA.yaml"))
+
+    def test_nl_gga(self):
+        def get_xc():
+            return get_cider_functional(
+                "functionals/CIDER23X_NL_GGA.yaml",
+                qmax=300,
+                lambd=1.8,
+                xmix=0.25,
+                pasdw_ovlp_fit=True,
+                pasdw_store_funcs=True,
+            )
+
+        test_cider_forces(get_xc())
+
+    def test_nl_mgga(self):
+        def get_xc():
+            return get_cider_functional(
+                "functionals/CIDER23X_NL_MGGA_DTR.yaml",
+                qmax=300,
+                lambd=1.8,
+                xmix=0.25,
+                pasdw_ovlp_fit=False,
+                pasdw_store_funcs=False,
+            )
+
+        test_cider_forces(get_xc())
+
+
+if __name__ == "__main__":
+    unittest.main()

ciderpress/gpaw/tests/test_read_write.py

@@ -0,0 +1,111 @@
+import unittest
+
+from ase.build import bulk
+from ase.parallel import parprint
+from gpaw import GPAW, PW, Mixer, restart
+from gpaw.mpi import world
+from gpaw.xc.libxc import LibXC
+from numpy.testing import assert_almost_equal, assert_equal
+
+from ciderpress.gpaw.cider_paw import CiderGPAW, get_cider_functional
+from ciderpress.gpaw.interp_paw import DiffGGA, DiffMGGA
+
+
+def test_load_write(xc, use_pp=False, is_cider=False, is_nl=False):
+    k = 3
+    si = bulk("Si")
+    kwargs = dict(
+        mode=PW(250),
+        mixer=Mixer(0.7, 5, 50.0),
+        xc=xc,
+        kpts=(k, k, k),
+        convergence={"energy": 1e-8},
+        parallel={"domain": min(2, world.size)},
+        occupations={"name": "fermi-dirac", "width": 0.0},
+        setups="sg15" if use_pp else "paw",
+        txt="si.txt",
+    )
+    si.calc = CiderGPAW(**kwargs)
+
+    e0 = si.get_potential_energy()
+    calc = si.calc
+
+    calc.write("_tmp.gpw")
+    si1, calc1 = restart("_tmp.gpw", Class=CiderGPAW)
+    if is_cider and is_nl:
+        xc0 = calc.hamiltonian.xc
+        xc1 = calc1.hamiltonian.xc
+        assert_equal(xc0.encut, xc1.encut)
+        assert_equal(xc0.Nalpha, xc1.Nalpha)
+        assert_equal(xc0.lambd, xc1.lambd)
+    calc1.set(mode=PW(250))
+    e1 = si1.get_potential_energy()
+    assert_almost_equal(e0, e1)
+    if is_cider and is_nl:
+        assert_almost_equal(xc0.alphas, xc1.alphas)
+    calc.set(mode=PW(320))
+    calc1.set(mode=PW(320))
+    e2 = calc.get_potential_energy()
+    e3 = calc1.get_potential_energy()
+    assert_almost_equal(e2, e3)
+    if is_cider:
+        new_kwargs = {k: v for k, v in kwargs.items()}
+        new_kwargs["xc"] = si.calc.hamiltonian.xc.todict()
+        new_kwargs["mlfunc_data"] = si.calc.hamiltonian.xc.get_mlfunc_data()
+        si.calc = CiderGPAW(**new_kwargs)
+        e4 = si.get_potential_energy()
+        assert_almost_equal(e0, e4)
+    # Test equivalence to normal GPAW
+    si.calc = GPAW(**kwargs)
+    e5 = si.get_potential_energy()
+    assert_almost_equal(e0, e5)
+
+
+def get_xc(fname, use_paw=True, force_nl=False):
+    return get_cider_functional(
+        fname,
+        qmax=300,
+        lambd=1.8,
+        xmix=0.25,
+        pasdw_ovlp_fit=True,
+        pasdw_store_funcs=False,
+        use_paw=use_paw,
+        _force_nonlocal=force_nl,
+    )
+
+
+class TestReadWrite(unittest.TestCase):
+    def test_nl_ml(self):
+        for fname in [
+            "functionals/CIDER23X_NL_GGA.yaml",
+            "functionals/CIDER23X_NL_MGGA.yaml",
+        ]:
+            for use_paw in [False, True]:
+                xc = get_xc(fname, use_paw)
+                parprint("TEST", fname, use_paw)
+                test_load_write(xc, use_pp=not use_paw, is_cider=True, is_nl=True)
+
+    def test_sl_ml(self):
+        for fname in [
+            "functionals/CIDER23X_SL_GGA.yaml",
+            "functionals/CIDER23X_SL_MGGA.yaml",
+        ]:
+            xc = get_xc(fname)
+            parprint("TEST", fname)
+            test_load_write(xc, is_cider=True)
+
+    def test_sl_ne(self):
+        for xc in [
+            DiffGGA(LibXC("PBE")),
+            DiffMGGA(LibXC("MGGA_X_R2SCAN+MGGA_C_R2SCAN")),
+        ]:
+            parprint("TEST", xc.kernel.name)
+            test_load_write(xc)
+
+        for xc in ["PBE", "MGGA_X_R2SCAN+MGGA_C_R2SCAN"]:
+            parprint("TEST", xc)
+            test_load_write(xc)
+
+
+if __name__ == "__main__":
+    unittest.main()

ciderpress/gpaw/tests/test_si_stress.py

@@ -0,0 +1,109 @@
+import unittest
+
+import numpy as np
+from ase.build import bulk
+from ase.parallel import parprint
+from gpaw import GPAW, PW, Mixer
+from gpaw.mpi import world
+
+from ciderpress.gpaw.cider_paw import get_cider_functional
+
+
+def test_pw_si_stress(xc, use_pp=False, s_numerical=None):
+    # This is based on a test from GPAW
+    k = 3
+    si = bulk("Si")
+    si.calc = GPAW(
+        mode=PW(250),
+        mixer=Mixer(0.7, 5, 50.0),
+        xc=xc,
+        kpts=(k, k, k),
+        convergence={"energy": 1e-8},
+        parallel={"domain": min(2, world.size)},
+        setups="sg15" if use_pp else "paw",
+        txt="si.txt",
+    )
+
+    si.set_cell(
+        np.dot(si.cell, [[1.02, 0, 0.03], [0, 0.99, -0.02], [0.2, -0.01, 1.03]]),
+        scale_atoms=True,
+    )
+
+    etot = si.get_potential_energy()
+    print(etot)
+
+    # Trigger nasty bug (fixed in !486):
+    si.calc.wfs.pt.blocksize = si.calc.wfs.pd.maxmyng - 1
+
+    s_analytical = si.get_stress()
+    parprint(s_analytical)
+    # TEST_CIDER_GGA Numerical
+    # [-0.00261187 -0.03790705 -0.03193711 -0.0209582   0.13427714  0.00928778]
+    # TEST_CIDER_MGGA Numerical
+    # [-0.00681636 -0.04026119 -0.03689781 -0.02227667  0.14441494  0.00907815]
+    if s_numerical is None:
+        s_numerical = si.calc.calculate_numerical_stress(si, 1e-5)
+    s_err = s_numerical - s_analytical
+
+    parprint("Analytical stress:\n", s_analytical)
+    parprint("Numerical stress:\n", s_numerical)
+    parprint("Error in stress:\n", s_err)
+    assert np.all(abs(s_err) < 1e-4)
+
+
+def get_xc(fname, use_paw=True):
+    return get_cider_functional(
+        fname,
+        qmax=120,
+        lambd=1.8,
+        xmix=0.25,
+        pasdw_ovlp_fit=True,
+        pasdw_store_funcs=True,
+        use_paw=use_paw,
+    )
+
+
+class TestStress(unittest.TestCase):
+    def test_nl_gga(self):
+        xc = get_xc("functionals/CIDER23X_NL_GGA.yaml")
+        s_numerical = np.array(
+            [-0.00261187, -0.03790705, -0.03193711, -0.0209582, 0.13427714, 0.00928778]
+        )
+        test_pw_si_stress(xc, s_numerical=s_numerical)
+
+    def test_nl_mgga(self):
+        xc = get_xc("functionals/CIDER23X_NL_MGGA_DTR.yaml")
+        s_numerical = np.array(
+            [-0.00681636, -0.04026119, -0.03689781, -0.02227667, 0.14441494, 0.00907815]
+        )
+        test_pw_si_stress(xc, s_numerical=s_numerical)
+
+    def test_sl_gga(self):
+        xc = get_xc("functionals/CIDER23X_SL_GGA.yaml")
+        # s_numerical = np.array([-0.00261187, -0.03790705, -0.03193711,
+        #                        -0.0209582,   0.13427714,  0.00928778])
+        test_pw_si_stress(xc, s_numerical=None)
+
+    def test_sl_mgga(self):
+        xc = get_xc("functionals/CIDER23X_SL_MGGA.yaml")
+        # s_numerical = np.array([-0.00681636, -0.04026119, -0.03689781,
+        #                        -0.02227667,  0.14441494,  0.00907815])
+        test_pw_si_stress(xc, s_numerical=None)
+
+    def test_nl_gga_pp(self):
+        xc = get_xc("functionals/CIDER23X_NL_GGA.yaml", use_paw=False)
+        s_numerical = np.array(
+            [0.00205983, -0.03604186, -0.02808641, -0.02021089, 0.1333823, 0.00980205]
+        )
+        test_pw_si_stress(xc, use_pp=True, s_numerical=s_numerical)
+
+    def test_nl_mgga_pp(self):
+        with self.assertRaises(NotImplementedError):
+            xc = get_xc("functionals/CIDER23X_NL_MGGA.yaml", use_paw=False)
+            # s_numerical = np.array([-0.00681636, -0.04026119, -0.03689781,
+            #                        -0.02227667,  0.14441494,  0.00907815])
+            test_pw_si_stress(xc, use_pp=True, s_numerical=None)
+
+
+if __name__ == "__main__":
+    unittest.main()

ciderpress/gpaw/tests/test_sph_harm.py

@@ -0,0 +1,32 @@
+import unittest
+
+import numpy as np
+from numpy.testing import assert_almost_equal
+from scipy.special import lpmn
+
+from ciderpress.dft.futil import fast_sph_harm as fsh
+
+
+class TestSphHarm(unittest.TestCase):
+    def test_sph_harm():
+        x = np.linspace(-2, 2, 401)
+        assert 0 in x
+        assert -1 in x
+        assert 1 in x
+        m = 8
+        n = 8
+        N = x.size
+        pm, pd = fsh.lpmn_vec(m, n, x)
+        assert pm.shape == (m + 1, n + 1, N)
+        assert pd.shape == (m + 1, n + 1, N)
+        pm_ref, pd_ref = np.zeros_like(pm), np.zeros_like(pd)
+        for i in range(N):
+            pm_ref[:, :, i], pd_ref[:, :, i] = lpmn(m, n, x[i])
+        diff = np.abs(pd_ref - pd)
+        diff[np.isnan(diff)] = 0.0
+        assert_almost_equal(pm, pm_ref, 14)
+        assert_almost_equal(pd, pd_ref, 14)
+
+
+if __name__ == "__main__":
+    unittest.main()