From 84e88b33c816414f044bd7b32e63c031ad2d856a Mon Sep 17 00:00:00 2001 From: kylebystrom Date: Sat, 19 Oct 2024 22:51:33 -0400 Subject: [PATCH] clean up gpaw descriptors and unit test a bit --- ciderpress/gpaw/fast_descriptors.py | 53 ++++------------------- ciderpress/gpaw/tests/test_descriptors.py | 8 ++-- 2 files changed, 13 insertions(+), 48 deletions(-) diff --git a/ciderpress/gpaw/fast_descriptors.py b/ciderpress/gpaw/fast_descriptors.py index c5b4d3c..8bf186c 100644 --- a/ciderpress/gpaw/fast_descriptors.py +++ b/ciderpress/gpaw/fast_descriptors.py @@ -24,7 +24,6 @@ from gpaw.occupations import FixedOccupationNumbers, OccupationNumberCalculator from gpaw.pw.density import ReciprocalSpaceDensity from gpaw.sphere.lebedev import Y_nL, weight_n -from gpaw.xc.gga import calculate_sigma from ciderpress.dft.plans import SemilocalPlan from ciderpress.dft.settings import ( @@ -111,21 +110,15 @@ def get_descriptors( raise ValueError("Settings must be Settings object or letter l") kcls = CiderMGGAHybridKernel cls = FeatSLMGGAPAW if use_paw else FeatSLMGGA - # TODO this part is messy/unncessary, should refactor to avoid cider_kernel = kcls(XCShell(RhoVectorSettings()), 0, "GGA_X_PBE", "GGA_C_PBE") elif isinstance(settings, SemilocalSettings): if settings.level == "MGGA": kcls = CiderMGGAHybridKernel cls = FeatSLMGGAPAW if use_paw else FeatSLMGGA - mode = "npa" else: kcls = CiderGGAHybridKernel cls = FeatSLGGAPAW if use_paw else FeatSLGGA - mode = "np" - # TODO this part is messy/unncessary, should refactor to avoid - cider_kernel = kcls( - XCShell(SemilocalSettings(mode)), 0, "GGA_X_PBE", "GGA_C_PBE" - ) + cider_kernel = kcls(XCShell(settings), 0, "GGA_X_PBE", "GGA_C_PBE") elif isinstance(settings, NLDFSettings): if settings.sl_level == "MGGA": kcls = CiderMGGAHybridKernel @@ -886,24 +879,9 @@ def _collect_feat(self, feat_xg): def get_features_on_grid(self): nt_sg = self._get_pseudo_density() - self.nspin = nt_sg.shape[0] - _, gradn_svg = calculate_sigma(self.gd, self.grad_v, nt_sg) + taut_sg = self._get_taut(nt_sg)[0] + rho_sxg = self._get_cider_inputs(nt_sg, taut_sg)[1] settings = self.cider_kernel.mlfunc.settings.sl_settings - if settings.level == "MGGA": - taut_sG = self.wfs.calculate_kinetic_energy_density() - if taut_sG is None: - taut_sG = self.wfs.gd.zeros(len(nt_sg)) - taut_sg = np.empty_like(nt_sg) - for taut_G, taut_g in zip(taut_sG, taut_sg): - if self.has_paw: - taut_G += 1.0 / self.wfs.nspins * self.tauct_G - self.distribute_and_interpolate(taut_G, taut_g) - tau_sg = taut_sg - rho_sxg = np.concatenate( - [nt_sg[:, np.newaxis], gradn_svg, tau_sg[:, np.newaxis]], axis=1 - ) - else: - rho_sxg = np.concatenate([nt_sg[:, np.newaxis], gradn_svg], axis=1) rho_sxg = rho_sxg.view() rho_sxg.shape = rho_sxg.shape[:2] + (-1,) if settings.mode != "l": @@ -915,29 +893,16 @@ def get_features_on_grid(self): def get_features_on_grid_deriv(self, p_j, drhodf_jxg, DD_aop=None): nt_sg = self._get_pseudo_density() - self.nspin = nt_sg.shape[0] - _, gradn_svg = calculate_sigma(self.gd, self.grad_v, nt_sg) - settings = self.cider_kernel.mlfunc.settings.sl_settings - if settings.level == "MGGA": - taut_sG = self.wfs.calculate_kinetic_energy_density() - if taut_sG is None: - taut_sG = self.wfs.gd.zeros(len(nt_sg)) - taut_sg = np.empty_like(nt_sg) - for taut_G, taut_g in zip(taut_sG, taut_sg): - if self.has_paw: - taut_G += 1.0 / self.wfs.nspins * self.tauct_G - self.distribute_and_interpolate(taut_G, taut_g) - tau_sg = taut_sg - rho_sxg = np.concatenate( - [nt_sg[:, np.newaxis], gradn_svg, tau_sg[:, np.newaxis]], axis=1 - ) - else: - rho_sxg = np.concatenate([nt_sg[:, np.newaxis], gradn_svg], axis=1) - drhodf_jxg = np.ascontiguousarray(drhodf_jxg[:, :4]) + taut_sg = self._get_taut(nt_sg)[0] + self.timer.start("Reorganize density") + rho_sxg = self._get_cider_inputs(nt_sg, taut_sg)[1] + drhodf_jxg = self._distribute_to_cider_grid(drhodf_jxg) + self.timer.stop() rho_sxg = rho_sxg.view() drhodf_jxg = drhodf_jxg.view() rho_sxg.shape = rho_sxg.shape[:2] + (-1,) drhodf_jxg.shape = drhodf_jxg.shape[:2] + (-1,) + settings = self.cider_kernel.mlfunc.settings.sl_settings if settings.mode != "l": plan = SemilocalPlan(settings, self.nspin) dfeat_jig = np.empty((len(p_j), settings.nfeat) + rho_sxg.shape[2:]) diff --git a/ciderpress/gpaw/tests/test_descriptors.py b/ciderpress/gpaw/tests/test_descriptors.py index 2f47706..4e2dc94 100644 --- a/ciderpress/gpaw/tests/test_descriptors.py +++ b/ciderpress/gpaw/tests/test_descriptors.py @@ -503,7 +503,7 @@ def run_nl_feature_test(xc, use_pp=False, spinpol=False, baseline="PBE"): def run_sl_feature_test(use_pp=False, spinpol=False): # TODO precision is poor for Si. Why is this? k = 3 - si = bulk("Ge") + si = bulk("Si") from gpaw.xc.libxc import LibXC if use_pp: @@ -512,7 +512,7 @@ def run_sl_feature_test(use_pp=False, spinpol=False): xc0 = xc0name xc1 = xc1name else: - xc0name = "MGGA_X_SCAN+MGGA_C_SCAN" + xc0name = "MGGA_X_R2SCAN+MGGA_C_R2SCAN" xc1name = "GGA_X_PBE+GGA_C_PBE" xc0 = DiffMGGA(LibXC(xc0name)) xc1 = DiffGGA(LibXC(xc1name)) @@ -555,7 +555,7 @@ def run_sl_feature_test(use_pp=False, spinpol=False): ni = NumInt() if use_pp: - exc0, vxc0 = ni.eval_xc_eff(xc0name, rho, xctype="GGA")[:2] + exc0, vxc0 = ni.eval_xc_eff(xc0name, rho[..., :4, :], xctype="GGA")[:2] exc1, vxc1 = ni.eval_xc_eff(xc1name, rho[..., :4, :], xctype="GGA")[:2] else: exc0, vxc0 = ni.eval_xc_eff(xc0name, rho, xctype="MGGA")[:2] @@ -563,7 +563,7 @@ def run_sl_feature_test(use_pp=False, spinpol=False): exc_tot_0 = np.sum(exc0 * feat[:, 0].sum(0) * wt) exc_tot_1 = np.sum(exc1 * feat[:, 0].sum(0) * wt) parprint(exc_tot_1 - exc_tot_0, ediff) - assert_almost_equal(exc_tot_1 - exc_tot_0, ediff, 7) + assert_almost_equal(exc_tot_1 - exc_tot_0, ediff, 5) eig_vbm, ei_vbm, en_vbm = nscfeig( si.calc, xc1, n1=p_vbm[2], n2=p_vbm[2] + 1, kpt_indices=[p_vbm[1]]