Skip to content

Commit

Permalink
clean up gpaw descriptors and unit test a bit
Browse files Browse the repository at this point in the history
  • Loading branch information
@kylebystrom
kylebystrom committed Oct 20, 2024
1 parent 1270bc8 commit 84e88b3
Show file tree
Hide file tree
Showing 2 changed files with 13 additions and 48 deletions.
53 changes: 9 additions & 44 deletions ciderpress/gpaw/fast_descriptors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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":
Expand All @@ -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:])
Expand Down
8 changes: 4 additions & 4 deletions ciderpress/gpaw/tests/test_descriptors.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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))
Expand Down Expand Up @@ -555,15 +555,15 @@ 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]
exc1, vxc1 = ni.eval_xc_eff(xc1name, rho[..., :4, :], xctype="GGA")[:2]
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]]
Expand Down

0 comments on commit 84e88b3

Please sign in to comment.