diff --git a/ciderpress/gpaw/atom_analysis_utils.py b/ciderpress/gpaw/atom_descriptor_utils.py similarity index 99% rename from ciderpress/gpaw/atom_analysis_utils.py rename to ciderpress/gpaw/atom_descriptor_utils.py index 93b1316..f5743dc 100644 --- a/ciderpress/gpaw/atom_analysis_utils.py +++ b/ciderpress/gpaw/atom_descriptor_utils.py @@ -389,7 +389,6 @@ def _get_paw_helper2( dydf_obg = dFdf_oag Nalpha = self.Nalpha - # TODO consts should be removed to make use of feat normalizers x_sig = np.zeros((nspin, nfeat, ngrid)) dxdf_oig = np.zeros((norb, nfeat, ngrid)) for s in range(nspin): @@ -406,7 +405,7 @@ def _get_paw_helper2( C_pg = self.C_aip[a, i_g].T pa_g = C_pg[0] + dq0_g * (C_pg[1] + dq0_g * (C_pg[2] + dq0_g * C_pg[3])) x_sig[s, i] += pa_g * y_sbg[s, a] - x_sig[s, i] *= ((self.consts[i, 1] + self.consts[-1, 1]) / 2) ** 1.5 + # x_sig[s, i] *= ((self.consts[i, 1] + self.consts[-1, 1]) / 2) ** 1.5 for o in range(norb): s = p_o[o][0] for i in range(nfeat): @@ -425,7 +424,7 @@ def _get_paw_helper2( dpa_g = 0.5 * (C_pg[1] + dq0_g * (2 * C_pg[2] + dq0_g * 3 * C_pg[3])) dxdf_oig[o, i] += dpa_g * dadf_g * y_sbg[s, a] dxdf_oig[o, i] += pa_g * dydf_obg[o, a] - dxdf_oig[o, i] *= ((self.consts[i, 1] + self.consts[-1, 1]) / 2) ** 1.5 + # dxdf_oig[o, i] *= ((self.consts[i, 1] + self.consts[-1, 1]) / 2) ** 1.5 return x_sig, dxdf_oig diff --git a/ciderpress/gpaw/calculator.py b/ciderpress/gpaw/calculator.py index a949e79..1740b14 100644 --- a/ciderpress/gpaw/calculator.py +++ b/ciderpress/gpaw/calculator.py @@ -40,7 +40,7 @@ def get_cider_functional( _force_nonlocal=False, ): """ - Initialize a CIDER surrogate hybrid exchange function of the form + Initialize a CIDER surrogate hybrid XC functional of the form E_xc = E_x^sl * (1 - xmix) + E_c^sl + xmix * E_x^CIDER where E_x^sl is given by the xkernel parameter, E_c^sl is given by the ckernel parameter, and E_x^CIDER is the ML exchange energy contains in mlfunc. @@ -50,9 +50,9 @@ def get_cider_functional( but correctness is not guaranteed because the nonlocal features will not be correct due to the lack of core electrons. - NOTE: If the mlfunc is determined to be semi-local, all the + NOTE: If the mlfunc is determined to be semilocal, all the internal settings are ignored, and a simpler, more efficient - class is returned to evaluate the semi-local functional. + class is returned to evaluate the semilocal functional. Args: mlfunc (MappedXC, str): An ML functional object or a str @@ -162,7 +162,7 @@ class is returned to evaluate the semi-local functional. msg = "Only implemented for b and d version, found {}" raise ValueError(msg.format(mlfunc.desc_version)) - const_list = _get_const_list(mlfunc) + const_list = get_const_list(mlfunc.settings.nldf_settings) nexp = 4 if use_paw: @@ -299,7 +299,7 @@ def cider_functional_from_dict(d): else: # xc_params should have Nalpha, lambd, encut. # For PAW, it should also have pasdw_ovlp_fit, pasdw_store_funcs. - const_list = _get_const_list(mlfunc) + const_list = get_const_list(mlfunc.settings.nldf_settings) xc = cls(cider_kernel, nexp, const_list, **(d["xc_params"])) else: xc = cls(LibXC(d["name"])) @@ -307,11 +307,13 @@ def cider_functional_from_dict(d): return xc -def _get_const_list(mlfunc): - settings = mlfunc.settings.nldf_settings +def get_const_list(settings): + """ + settings should be NLDFSettings object + """ thetap = 2 * np.array(settings.theta_params) vvmul = thetap[0] / (2 * settings.feat_params[1][0]) - fac_mul = thetap[2] if mlfunc.settings.sl_settings.level == "MGGA" else thetap[1] + fac_mul = thetap[2] if settings.sl_level == "MGGA" else thetap[1] consts = np.array([0.00, thetap[0], fac_mul, thetap[0] / 64]) / vvmul const_list = np.stack([0.5 * consts, 1.0 * consts, 2.0 * consts, consts * vvmul]) return const_list diff --git a/ciderpress/gpaw/cider_fft.py b/ciderpress/gpaw/cider_fft.py index e166418..948a788 100644 --- a/ciderpress/gpaw/cider_fft.py +++ b/ciderpress/gpaw/cider_fft.py @@ -1237,8 +1237,8 @@ def from_mlfunc( class CiderGGA(_CiderBase, GGA): def __init__(self, cider_kernel, nexp, consts, **kwargs): - if cider_kernel.mlfunc.desc_version != "d": - raise ValueError("Wrong mlfunc version") + if cider_kernel.mlfunc.settings.sl_settings.level != "GGA": + raise ValueError("CiderGGA only supports GGA functionals!") _CiderBase.__init__(self, cider_kernel, nexp, consts, **kwargs) GGA.__init__(self, LibXC("PBE"), stencil=2) @@ -1335,9 +1335,8 @@ def integrate(a1_g, a2_g=None): class CiderMGGA(_CiderBase, MGGA): def __init__(self, cider_kernel, nexp, consts, **kwargs): - - if cider_kernel.mlfunc.desc_version != "b": - raise ValueError("Wrong mlfunc version") + if cider_kernel.mlfunc.settings.sl_settings.level != "MGGA": + raise ValueError("CiderMGGA only supports MGGA functionals!") _CiderBase.__init__(self, cider_kernel, nexp, consts, **kwargs) MGGA.__init__(self, LibXC("PBE"), stencil=2) diff --git a/ciderpress/gpaw/analysis.py b/ciderpress/gpaw/descriptors.py similarity index 90% rename from ciderpress/gpaw/analysis.py rename to ciderpress/gpaw/descriptors.py index 4616ba5..fe5ceb4 100644 --- a/ciderpress/gpaw/analysis.py +++ b/ciderpress/gpaw/descriptors.py @@ -6,7 +6,14 @@ from gpaw.sphere.lebedev import Y_nL, weight_n from gpaw.xc.gga import calculate_sigma -from ciderpress.gpaw.atom_analysis_utils import ( +from ciderpress.dft.settings import ( + FeatureSettings, + FracLaplSettings, + NLDFSettings, + SDMXBaseSettings, + SemilocalSettings, +) +from ciderpress.gpaw.atom_descriptor_utils import ( calculate_paw_cider_features_p1, calculate_paw_cider_features_p2, calculate_paw_cider_features_p2_noderiv, @@ -15,6 +22,7 @@ get_features_with_sl_noderiv, get_features_with_sl_part, ) +from ciderpress.gpaw.calculator import get_const_list from ciderpress.gpaw.cider_paw import ( CiderGGA, CiderGGAHybridKernel, @@ -25,6 +33,154 @@ ) +class XCShell: + def __init__(self, settings): + if isinstance(settings, SemilocalSettings): + self.settings = FeatureSettings(sl_settings=settings) + elif isinstance(settings, NLDFSettings): + if settings.sl_level == "MGGA": + sl_settings = SemilocalSettings("npa") + else: + sl_settings = SemilocalSettings("np") + self.settings = FeatureSettings( + sl_settings=sl_settings, + nldf_settings=settings, + ) + elif isinstance(settings, FracLaplSettings): + self.settings = FeatureSettings( + sl_settings=SemilocalSettings("npa"), + nlof_settings=settings, + ) + elif isinstance(settings, SDMXBaseSettings): + self.settings = FeatureSettings( + sl_settings=SemilocalSettings("npa"), + sdmx_settings=settings, + ) + else: + raise ValueError("Unsupported settings") + + +def get_features( + calc, + settings, + p_i=None, + use_paw=True, + screen_dens=True, + **kwargs, + # version="b", + # a0=1.0, + # fac_mul=0.03125, + # vvmul=1.0, + # amin=0.015625, + # qmax=300, + # lambd=1.8, +): + """ + Compute grid weights, feature vectors, and (optionally) + feature vector derivatives with respect to orbital + occupation. + + Args: + calc: a converged GPAW calculation + p_i (optional): A list of (s, k, n) indexes for orbitals + for which to differentiate features with respect to + occupation + use_paw (bool, True): Whether to use PAW corrections. + version (str, 'b'): Descriptor version b, d, or l. + l stands for local. + + Returns: + if p_i is None: + weights, features + else: + weights, features, [list of feature derivatives] + """ + + if isinstance(settings, str): + if settings != "l": + raise ValueError("Settings must be Settings object or letter l") + kcls = CiderMGGAHybridKernel + cls = FeatSLPAW if use_paw else FeatSL + # TODO this part is messy/unncessary, should refactor to avoid + const_list = np.ones((4, 4)) + cider_kernel = kcls( + XCShell(SemilocalSettings("nst")), 0, "GGA_X_PBE", "GGA_C_PBE" + ) + nexp = len(const_list) + elif isinstance(settings, SemilocalSettings): + if settings.level == "MGGA": + kcls = CiderMGGAHybridKernel + cls = FeatSLPAW if use_paw else FeatSL + else: + kcls = CiderGGAHybridKernel + cls = FeatSLPAW if use_paw else FeatSL + # TODO this part is messy/unncessary, should refactor to avoid + const_list = np.ones((4, 4)) + cider_kernel = kcls(None, 0, "GGA_X_PBE", "GGA_C_PBE") + nexp = len(const_list) + elif isinstance(settings, NLDFSettings): + if settings.sl_level == "MGGA": + kcls = CiderMGGAHybridKernel + cls = FeatMGGAPAW if use_paw else FeatMGGA + else: + kcls = CiderGGAHybridKernel + cls = FeatGGAPAW if use_paw else FeatGGA + const_list = get_const_list(settings) + cider_kernel = kcls(XCShell(settings), 0, "GGA_X_PBE", "GGA_C_PBE") + nexp = len(const_list) + elif isinstance(settings, FracLaplSettings): + raise NotImplementedError( + "Fractional Laplacian-based orbital descriptors have not yet been implemented for GPAW." + ) + elif isinstance(settings, SDMXBaseSettings): + raise NotImplementedError( + "SDMX descriptors have not yet been implemented for GPAW." + ) + else: + raise ValueError("Unsupported settings") + # TOD clean this up a bit + kwargs["xmix"] = 0.0 + kwargs["encut"] = kwargs.get("qmax") or kwargs.get("encut") or 300 + kwargs["lambd"] = kwargs.get("lambd") or 1.8 + empty_xc = cls( + cider_kernel, + nexp, + const_list, + pasdw_ovlp_fit=True, + pasdw_store_funcs=False, + **kwargs, + ) + empty_xc.set_grid_descriptor(calc.density.finegd) + empty_xc.initialize(calc.density, calc.hamiltonian, calc.wfs) + empty_xc.set_positions(calc.spos_ac) + if calc.density.nt_sg is None: + calc.density.interpolate_pseudo_density() + if empty_xc.orbital_dependent: + calc.converge_wave_functions() + if p_i is None: + return get_features_and_weights(calc, empty_xc, screen_dens=screen_dens) + else: + if use_paw: + DD_aop = get_empty_atomic_matrix(calc.density, len(p_i)) + calculate_single_orbital_atomic_density_matrix( + calc.wfs, DD_aop, p_i, calc.density.ncomponents + ) + else: + DD_aop = None + drhodf_ixR = get_drho_df(calc, p_i) + drhodf_ixg = interpolate_drhodf( + empty_xc.gd, empty_xc.distribute_and_interpolate, drhodf_ixR + ) + return get_features_and_weights_deriv( + calc, + empty_xc, + p_i, + drhodf_ixg, + DD_aop=DD_aop, + screen_dens=screen_dens, + ) + + def calc_fixed(bd, fref_sqn, f_qn): if bd.nbands == fref_sqn.shape[2]: nspins = fref_sqn.shape[0] @@ -55,6 +211,7 @@ def calculate_single_orbital_atomic_density_matrix(self, D_aop, p_o, nspin): """ Get the atomic density matrices for the band given by index tuple p + Args: self: GPAW Wavefunctions object D_asp: stores the single-orbital density matrix @@ -84,12 +241,12 @@ def get_empty_atomic_matrix(dens, norb): class FixedGOccupationNumbers(FixedOccupationNumbers): """ Same as FixedOccupationNumbers but with a perturb_occupation_ - functionn for testing incremental occupation numbers. + function for testing incremental occupation numbers. """ def perturb_occupation_(self, p, delta): """ - Change the occupation number of a band at a k-point + Change the occupation number of a band at every k-point Args: p: (s, k, n) or (s, n). k is ignored delta: amount to change occupation number @@ -580,11 +737,8 @@ def calculate_6d_integral_deriv( if n_g is not None: dfeatdf = np.zeros([nexp - 1] + list(n_g.shape)) for i in range(nexp - 1): - # TODO this should be changed to use the feature normalizers - const = ((self.consts[i, 1] + self.consts[-1, 1]) / 2) ** 1.5 - const = const + self.consts[i, 0] / (self.consts[i, 0] + const) for ind, a in enumerate(self.alphas): - dfeatdf[i, :] += const * p_iag[i, a] * self.rbuf_ag[a] + dfeatdf[i, :] += p_iag[i, a] * self.rbuf_ag[a] self.timer.start("6d comm fwd") if n_g is not None: @@ -615,8 +769,8 @@ def get_features_on_grid_deriv(self, p_j, drhodf_jxg, DD_aop=None): Derivatives of the density (0), gradient (1,2,3), and kinetic energy density (4) with respect to orbital occupation number. - include_density: whether to return density or just - feature values + DD_aop (dict of numpy array): derivatives of atomic density + matrices with respect to orbital occupation. Returns: @@ -923,98 +1077,3 @@ def interpolate_drhodf(gd, interp, drhodf_ixR): for x in range(drhodf_ixR.shape[1]): interp(drhodf_ixR[i, x], drhodf_ixg[i, x]) return drhodf_ixg - - -def get_features( - calc, - p_i=None, - use_paw=True, - version="b", - a0=1.0, - fac_mul=0.03125, - vvmul=1.0, - amin=0.015625, - qmax=300, - lambd=1.8, - screen_dens=True, -): - """ - Compute grid weights, feature vectors, and (optionally) - feature vector derivatives with respect to orbital - occupation. - - Args: - calc: a converged GPAW calculation - p_i (optional): A list of (s, k, n) indexes for orbitals - for which to differentiate features with respect to - occupation - use_paw (bool, True): Whether to use PAW corrections. - version (str, 'b'): Descriptor version b, d, or l. - l stands for local. - - Returns: - if p_i is None: - weights, features - else: - weights, features, [list of feature derivatives] - """ - from ciderpress.descriptors import XCShell - - mlfunc = XCShell(version if version != "l" else "b", a0, fac_mul, amin, vvmul) - if version == "b": - kcls = CiderMGGAHybridKernel - cls = FeatMGGAPAW if use_paw else FeatMGGA - elif version == "d": - kcls = CiderGGAHybridKernel - cls = FeatGGAPAW if use_paw else FeatGGA - elif version == "l": - kcls = CiderMGGAHybridKernel - cls = FeatSLPAW if use_paw else FeatSL - else: - raise ValueError("Version not supported.") - cider_kernel = kcls(mlfunc, 0, "GGA_X_PBE", "GGA_C_PBE") - consts = np.array([0.00, mlfunc.a0, mlfunc.fac_mul, mlfunc.amin]) - const_list = np.stack( - [0.5 * consts, 1.0 * consts, 2.0 * consts, consts * mlfunc.vvmul] - ) - nexp = 4 - empty_xc = cls( - cider_kernel, - nexp, - const_list, - Nalpha=None, - lambd=lambd, - encut=qmax, - xmix=0.0, - pasdw_ovlp_fit=True, - pasdw_store_funcs=False, - ) - empty_xc.set_grid_descriptor(calc.density.finegd) - empty_xc.initialize(calc.density, calc.hamiltonian, calc.wfs) - empty_xc.set_positions(calc.spos_ac) - if calc.density.nt_sg is None: - calc.density.interpolate_pseudo_density() - if empty_xc.orbital_dependent: - calc.converge_wave_functions() - if p_i is None: - return get_features_and_weights(calc, empty_xc, screen_dens=screen_dens) - else: - if use_paw: - DD_aop = get_empty_atomic_matrix(calc.density, len(p_i)) - calculate_single_orbital_atomic_density_matrix( - calc.wfs, DD_aop, p_i, calc.density.ncomponents - ) - else: - DD_aop = None - drhodf_ixR = get_drho_df(calc, p_i) - drhodf_ixg = interpolate_drhodf( - empty_xc.gd, empty_xc.distribute_and_interpolate, drhodf_ixR - ) - return get_features_and_weights_deriv( - calc, - empty_xc, - p_i, - drhodf_ixg, - DD_aop=DD_aop, - screen_dens=screen_dens, - ) diff --git a/ciderpress/gpaw/paw_cider_kernel.py b/ciderpress/gpaw/paw_cider_kernel.py index 6015754..4c6d750 100644 --- a/ciderpress/gpaw/paw_cider_kernel.py +++ b/ciderpress/gpaw/paw_cider_kernel.py @@ -187,7 +187,9 @@ def get_paw_atom_feat(self, n_sg, sigma_xg, y_sbg, F_sbg, ae, include_density): n_sg, sigma_xg, y_sbg, F_sbg, ae=ae ) if include_density: - from ciderpress.gpaw.atom_analysis_utils import get_features_with_sl_noderiv + from ciderpress.gpaw.atom_descriptor_utils import ( + get_features_with_sl_noderiv, + ) feat_sig = get_features_with_sl_noderiv(feat_sig, n_sg, sigma_xg, None) # _feat_sig = feat_sig @@ -353,7 +355,9 @@ def get_paw_atom_feat(self, n_sg, sigma_xg, y_sbg, F_sbg, ae, include_density): n_sg, tau_xg, y_sbg, F_sbg, ae=ae ) if include_density: - from ciderpress.gpaw.atom_analysis_utils import get_features_with_sl_noderiv + from ciderpress.gpaw.atom_descriptor_utils import ( + get_features_with_sl_noderiv, + ) feat_sig = get_features_with_sl_noderiv(feat_sig, n_sg, sigma_xg, tau_sg) # _feat_sig = feat_sig diff --git a/ciderpress/gpaw/tests/test_descriptors.py b/ciderpress/gpaw/tests/test_descriptors.py new file mode 100644 index 0000000..156960d --- /dev/null +++ b/ciderpress/gpaw/tests/test_descriptors.py @@ -0,0 +1,498 @@ +import unittest + +import numpy as np +from ase.build import bulk +from ase.dft.bandgap import bandgap +from ase.parallel import parprint +from ase.units import Ha +from gpaw import GPAW, PW, Mixer +from gpaw.mpi import world +from gpaw.xc import XC +from numpy.testing import assert_almost_equal + +from ciderpress.gpaw.calculator import get_cider_functional +from ciderpress.gpaw.descriptors import ( + get_drho_df, + get_features, + get_homo_lumo_fd_, + interpolate_drhodf, + run_constant_occ_calculation_, +) +from ciderpress.gpaw.xc_tools import non_self_consistent_eigenvalues as nscfeig + + +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, + ) + + +def get_static_xc_difference(calc, xc): + if isinstance(xc, (str, dict)): + xc = XC(xc) + xc.set_grid_descriptor(calc.density.finegd) + xc.initialize(calc.density, calc.hamiltonian, calc.wfs) + xc.set_positions(calc.spos_ac) + return calc.hamiltonian.get_xc_difference(xc, calc.density) * Ha + + +# TODO add unit tests for band, kpt, and domain parallelization + + +def _perturb_calc_density_(calc, p, delta): + calc.wfs.occupations.perturb_occupation_(p, delta, div_by_wt=True) + calc.wfs.calculate_occupation_numbers() + calc.wfs.calculate_atomic_density_matrices(calc.density.D_asp) + calc.density.calculate_pseudo_density(calc.wfs) + calc.density.interpolate_pseudo_density() + + +def run_fd_deriv_test(xc, use_pp=False, spinpol=False): + 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)}, + occupations={"name": "fermi-dirac", "width": 0.0}, + spinpol=spinpol, + setups="sg15" if use_pp else "paw", + txt="si.txt", + ) + + etot = si.get_potential_energy() + gap, vbm, cbm, dev, dec, pvbm, pcbm = get_homo_lumo_fd_(si.calc, delta=1e-6) + parprint(dev - vbm, dec - cbm, dev, dec) + assert_almost_equal(dev, vbm, 5) + assert_almost_equal(dec, cbm, 5) + return etot + + +def run_nscf_eigval_test(xc0, xc1, spinpol=False, use_pp=False, safe=True): + if use_pp is False: + from ciderpress.gpaw.interp_paw import DiffGGA + + if isinstance(xc0, str): + xc0 = DiffGGA(XC(xc0).kernel) + if isinstance(xc1, str): + xc1 = DiffGGA(XC(xc1).kernel) + + k = 3 + si = bulk("Si") + si.calc = GPAW( + mode=PW(250), + mixer=Mixer(0.7, 5, 50.0), + xc=xc0, + kpts=(k, k, k), + convergence={"energy": 1e-8, "density": 1e-10}, + parallel={"domain": min(2, world.size)}, + occupations={"name": "fermi-dirac", "width": 0.0}, + spinpol=spinpol, + setups="sg15" if use_pp else "paw", + txt="si.txt", + ) + delta = 1e-5 + si.get_potential_energy() + gap, p_vbm, p_cbm = bandgap(si.calc) + run_constant_occ_calculation_(si.calc) + + if safe: + # CIDER changes the xc_corrections and uses different grids, + # so we need to check this. + ediff00 = si.calc.get_xc_difference(xc0) / Ha + _perturb_calc_density_(si.calc, p_vbm, -delta) + ediff10 = get_static_xc_difference(si.calc, xc0) / Ha + _perturb_calc_density_(si.calc, p_vbm, delta) + _perturb_calc_density_(si.calc, p_cbm, delta) + ediff20 = get_static_xc_difference(si.calc, xc0) / Ha + _perturb_calc_density_(si.calc, p_cbm, -delta) + ediff01 = si.calc.get_xc_difference(xc1) / Ha + eig_vbm, ei_vbm, en_vbm = nscfeig( + si.calc, xc1, n1=p_vbm[2], n2=p_vbm[2] + 1, kpt_indices=[p_vbm[1]] + ) + eig_cbm, ei_cbm, en_cbm = nscfeig( + si.calc, xc1, n1=p_cbm[2], n2=p_cbm[2] + 1, kpt_indices=[p_cbm[1]] + ) + eigdiff_vbm = (en_vbm - ei_vbm)[p_vbm[0], 0, 0] / Ha + eigdiff_cbm = (en_cbm - ei_cbm)[p_cbm[0], 0, 0] / Ha + _perturb_calc_density_(si.calc, p_vbm, -delta) + ediff11 = get_static_xc_difference(si.calc, xc1) / Ha + _perturb_calc_density_(si.calc, p_vbm, delta) + _perturb_calc_density_(si.calc, p_cbm, delta) + ediff21 = get_static_xc_difference(si.calc, xc1) / Ha + _perturb_calc_density_(si.calc, p_cbm, -delta) + ediff0 = ediff01 - ediff00 + ediff1 = ediff11 - ediff10 + ediff2 = ediff21 - ediff20 + else: + ediff0 = si.calc.get_xc_difference(xc1) / Ha + eig_vbm, ei_vbm, en_vbm = nscfeig( + si.calc, xc1, n1=p_vbm[2], n2=p_vbm[2] + 1, kpt_indices=[p_vbm[1]] + ) + eig_cbm, ei_cbm, en_cbm = nscfeig( + si.calc, xc1, n1=p_cbm[2], n2=p_cbm[2] + 1, kpt_indices=[p_cbm[1]] + ) + eigdiff_vbm = (en_vbm - ei_vbm)[p_vbm[0], 0, 0] / Ha + eigdiff_cbm = (en_cbm - ei_cbm)[p_cbm[0], 0, 0] / Ha + _perturb_calc_density_(si.calc, p_vbm, -delta) + ediff1 = ( + get_static_xc_difference(si.calc, xc1) + - get_static_xc_difference(si.calc, xc0) + ) / Ha + _perturb_calc_density_(si.calc, p_vbm, delta) + _perturb_calc_density_(si.calc, p_cbm, delta) + ediff2 = ( + get_static_xc_difference(si.calc, xc1) + - get_static_xc_difference(si.calc, xc0) + ) / Ha + fd_eigdiff_vbm = (ediff0 - ediff1) / delta + fd_eigdiff_cbm = (ediff2 - ediff0) / delta + parprint(eigdiff_cbm - eigdiff_vbm, fd_eigdiff_cbm - fd_eigdiff_vbm) + parprint(eigdiff_cbm, eigdiff_vbm, fd_eigdiff_cbm, fd_eigdiff_vbm) + assert_almost_equal(eigdiff_vbm, fd_eigdiff_vbm, 5) + assert_almost_equal(eigdiff_cbm, fd_eigdiff_cbm, 5) + + +def run_drho_test(spinpol=False, use_pp=False): + k = 3 + si = bulk("Si") + si.calc = GPAW( + mode=PW(250), + mixer=Mixer(0.7, 5, 50.0), + xc="PBE", + kpts=(k, k, k), + convergence={"energy": 1e-8, "density": 1e-10}, + parallel={"domain": min(2, world.size)}, + occupations={"name": "fermi-dirac", "width": 0.0}, + spinpol=spinpol, + setups="sg15" if use_pp else "paw", + # symmetry='off', + txt="si.txt", + ) + delta = 1e-4 + si.get_potential_energy() + gap, p_vbm, p_cbm = bandgap(si.calc) + run_constant_occ_calculation_(si.calc) + p_i = [p_vbm, p_cbm] + drhodf_ixR = get_drho_df(si.calc, p_i) + nt0_sG = si.calc.density.nt_sG.copy() + nt0_sg = si.calc.density.nt_sg.copy() + + _perturb_calc_density_(si.calc, p_vbm, -delta) + + nt1_sG = si.calc.density.nt_sG.copy() + nt1_sg = si.calc.density.nt_sg.copy() + dnt_sG = (nt1_sG - nt0_sG) / delta * -1 + dnt_sg = (nt1_sg - nt0_sg) / delta * -1 + diff = drhodf_ixR[0, 0] - dnt_sG + idiff = si.calc.hamiltonian.xc.gd.integrate(diff * nt0_sG) + assert_almost_equal(idiff, 0) + + drhodf_ixg = interpolate_drhodf( + si.calc.hamiltonian.xc.gd, + si.calc.density.distribute_and_interpolate, + drhodf_ixR, + ) + ind = 0 + diff = drhodf_ixg[ind, 0] - dnt_sg + idiff = si.calc.hamiltonian.xc.gd.integrate(diff * nt0_sg) + assert_almost_equal(idiff, 0) + + +def run_nl_feature_test(xc, use_pp=False, spinpol=False, baseline="PBE"): + nspin = 2 if spinpol else 1 + k = 3 + si = bulk("Si") + if use_pp is False: + from ciderpress.gpaw.interp_paw import DiffGGA + + baseline = DiffGGA(XC(baseline).kernel) + + si.calc = GPAW( + mode=PW(250), + mixer=Mixer(0.7, 5, 50.0), + xc=xc, + kpts=(k, k, k), + convergence={"energy": 1e-8, "density": 1e-10}, + parallel={"domain": min(2, world.size)}, + occupations={"name": "fermi-dirac", "width": 0.0}, + spinpol=spinpol, + setups="sg15" if use_pp else "paw", + txt="si.txt", + ) + + mlfunc = xc.cider_kernel.mlfunc + kwargs = dict( + settings=mlfunc.settings.nldf_settings, + use_paw=not use_pp, + screen_dens=False, + qmax=xc.encut, + lambd=xc.lambd, + ) + + si.get_potential_energy() + + calc = si.calc + gap, p_vbm, p_cbm = bandgap(si.calc) + delta = 1e-6 + run_constant_occ_calculation_(si.calc) + + ediff = si.calc.get_xc_difference(baseline) / Ha + xmix = xc.cider_kernel.xmix + feat0, wt0 = get_features(si.calc, **kwargs) + feat, dfeat_j, wt = get_features(si.calc, p_i=[p_vbm, p_cbm], **kwargs) + if use_pp: + assert feat.shape[-1] == np.prod(xc.gd.get_size_of_global_array()) + assert feat.shape[-1] == dfeat_j.shape[-1] + assert feat.shape[-1] == wt.size + assert_almost_equal(feat, feat0) + assert_almost_equal(wt, wt0) + desc = feat.copy() + etst = 0 + X0TN = mlfunc.settings.normalizers.get_normalized_feature_vector(desc) + eps, deps = mlfunc(X0TN) + etst = np.sum(eps * wt) + # TODO why is the precision so low here? (err ~ 10^-6) + # TODO uncomment, fix spinpol + parprint(xmix * etst, -ediff) + assert_almost_equal(xmix * etst, -ediff, 6) + + eig_vbm, ei_vbm, en_vbm = nscfeig( + si.calc, baseline, n1=p_vbm[2], n2=p_vbm[2] + 1, kpt_indices=[p_vbm[1]] + ) + eig_cbm, ei_cbm, en_cbm = nscfeig( + si.calc, baseline, n1=p_cbm[2], n2=p_cbm[2] + 1, kpt_indices=[p_cbm[1]] + ) + eigdiff_vbm = (en_vbm - ei_vbm)[p_vbm[0], 0, 0] / Ha + eigdiff_cbm = (en_cbm - ei_cbm)[p_cbm[0], 0, 0] / Ha + eigdiff_gap = eigdiff_vbm - eigdiff_cbm + + _perturb_calc_density_(calc, p_vbm, -delta) + feat1, wt1 = get_features(si.calc, **kwargs) + deriv1 = (feat0 - feat1) / delta + _perturb_calc_density_(calc, p_vbm, delta) + _perturb_calc_density_(calc, p_cbm, delta) + feat2, wt2 = get_features(si.calc, **kwargs) + _perturb_calc_density_(calc, p_cbm, -delta) + deriv2 = (feat2 - feat0) / delta + derivs = [deriv1, deriv2] + # TODO features don't match but give same energy differences. + # Not sure why this is. Definitely related to symmetry, but + # sum(dfeat * rho) should still be consistent with FD and is not. + # for i in range(feat0[0].shape[0]): + # for j, p in enumerate([p_vbm, p_cbm]): + # s = p[0] + # rho = feat0[s, 0] + # assert_almost_equal( + # np.sum(dfeat_j[j, i] * rho), np.sum(derivs[j][s, i] * rho), 5 + # ) + + etst = 0 + etst2 = 0 + p_j = [p_vbm, p_cbm] + pwt = [-1, 1] + X0TN = mlfunc.settings.normalizers.get_normalized_feature_vector(desc) + eps, deps = mlfunc(X0TN) + deps = mlfunc.settings.normalizers.get_derivative_wrt_unnormed_features(desc, deps) + for j in range(2): + p = p_j[j] + s = p[0] + dfeat_tmp = dfeat_j[j] + dfeat_tmp2 = derivs[j][s] + deps_tmp = np.einsum("xg,xg->g", deps[s], dfeat_tmp) + etst += pwt[j] * np.sum(deps_tmp * wt) + deps_tmp = np.einsum("xg,xg->g", deps[s], dfeat_tmp2) + etst2 += pwt[j] * np.sum(deps_tmp * wt) + # TODO why is this not consistent depending on symmetry? + parprint(xmix * etst, xmix * etst2, eigdiff_gap) + assert_almost_equal(xmix * etst, eigdiff_gap, 5) + assert_almost_equal(xmix * etst2, xmix * etst, 4) + assert_almost_equal(xmix * etst2, eigdiff_gap, 4) + + +def run_sl_feature_test(use_pp=False, spinpol=False): + # TODO precision is poor for Si. Why is this? + k = 3 + si = bulk("Ge") + from gpaw.xc.libxc import LibXC + + from ciderpress.gpaw.interp_paw import DiffGGA, DiffMGGA + + if use_pp: + xc0name = "PBE" + xc1name = "PBEsol" + xc0 = xc0name + xc1 = xc1name + else: + xc0name = "MGGA_X_SCAN+MGGA_C_SCAN" + xc1name = "GGA_X_PBE+GGA_C_PBE" + xc0 = DiffMGGA(LibXC(xc0name)) + xc1 = DiffGGA(LibXC(xc1name)) + + si.calc = GPAW( + mode=PW(250), + mixer=Mixer(0.7, 5, 50.0), + xc=xc0, + kpts=(k, k, k), + convergence={"energy": 1e-8, "density": 1e-10}, + parallel={"domain": min(2, world.size)}, + occupations={"name": "fermi-dirac", "width": 0.0}, + spinpol=spinpol, + setups="sg15" if use_pp else "paw", + txt="si.txt", + ) + + kwargs = dict( + settings="l", + use_paw=not use_pp, + screen_dens=False, + ) + + si.get_potential_energy() + gap, p_vbm, p_cbm = bandgap(si.calc) + # run_constant_occ_calculation_(si.calc) + # etot = si.get_potential_energy() / Ha + + ediff = si.calc.get_xc_difference(xc1) / Ha + feat0, wt0 = get_features(si.calc, **kwargs) + feat, dfeat_j, wt = get_features(si.calc, p_i=[p_vbm, p_cbm], **kwargs) + + from pyscf.dft.numint import NumInt + + assert_almost_equal(feat, feat0) + assert_almost_equal(wt, wt0) + rho = feat.copy() + if not spinpol: + rho = rho.squeeze(0) + + ni = NumInt() + if use_pp: + exc0, vxc0 = ni.eval_xc_eff(xc0name, rho, 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] + # exc0, vxc0 = eval_xc('PBE', rho, 1 if spinpol else 0)[:2] + # exc1, vxc1 = eval_xc('PBE', rho, 1 if spinpol else 0)[:2] + exc_tot_0 = np.sum(exc0 * feat[:, 0].sum(0) * wt) + exc_tot_1 = np.sum(exc1 * feat[:, 0].sum(0) * wt) + assert_almost_equal(exc_tot_1 - exc_tot_0, ediff, 4) + + eig_vbm, ei_vbm, en_vbm = nscfeig( + si.calc, xc1, n1=p_vbm[2], n2=p_vbm[2] + 1, kpt_indices=[p_vbm[1]] + ) + eig_cbm, ei_cbm, en_cbm = nscfeig( + si.calc, xc1, n1=p_cbm[2], n2=p_cbm[2] + 1, kpt_indices=[p_cbm[1]] + ) + eigdiff_vbm = (en_vbm - ei_vbm)[p_vbm[0], 0, 0] / Ha + eigdiff_cbm = (en_cbm - ei_cbm)[p_cbm[0], 0, 0] / Ha + + dde_list = [] + for p, dfeat in zip([p_vbm, p_cbm], dfeat_j): + if spinpol: + s = p[0] + vtmp0 = vxc0[s] + vtmp1 = vxc1[s] + else: + vtmp0 = vxc0 + vtmp1 = vxc1 + de0 = np.sum(dfeat[: len(vtmp0)] * vtmp0 * wt) + de1 = np.sum(dfeat[: len(vtmp1)] * vtmp1 * wt) + dde = de1 - de0 + dde_list.append(dde) + + parprint(dde_list, eigdiff_vbm, eigdiff_cbm) + assert_almost_equal(dde_list[0], eigdiff_vbm, 4) + assert_almost_equal(dde_list[1], eigdiff_cbm, 4) + assert_almost_equal(dde_list[1] - dde_list[0], eigdiff_cbm - eigdiff_vbm, 4) + + +class TestDescriptors(unittest.TestCase): + def test_occ_deriv_rho(self): + run_drho_test(spinpol=False, use_pp=True) + run_drho_test(spinpol=False, use_pp=False) + + def test_eigval(self): + for use_pp in [True, False]: + xc = get_xc( + "functionals/CIDER23X_NL_GGA.yaml", use_paw=not use_pp, force_nl=True + ) + run_nscf_eigval_test( + "PBE", xc, spinpol=False, use_pp=use_pp, safe=not use_pp + ) + run_nscf_eigval_test( + "PBE", xc, spinpol=True, use_pp=use_pp, safe=not use_pp + ) + xc = get_xc( + "functionals/CIDER23X_NL_MGGA_DTR.yaml", + use_paw=not use_pp, + force_nl=True, + ) + run_nscf_eigval_test( + "PBE", xc, spinpol=False, use_pp=use_pp, safe=not use_pp + ) + run_nscf_eigval_test( + "PBE", xc, spinpol=True, use_pp=use_pp, safe=not use_pp + ) + + def test_sl_features(self): + # run_sl_feature_test(spinpol=False, use_pp=True) + run_sl_feature_test(spinpol=False) + run_sl_feature_test(spinpol=True) + + def test_nl_features(self): + SPINPOL = False + for use_pp in [True, False]: + xc = get_xc( + "functionals/CIDER23X_NL_GGA.yaml", use_paw=not use_pp, force_nl=True + ) + baseline = "0.75_GGA_X_PBE+1.00_GGA_C_PBE" + run_nl_feature_test(xc, use_pp=use_pp, spinpol=SPINPOL, baseline=baseline) + + baseline = "0.75_GGA_X_PBE+1.00_GGA_C_PBE" + xc = get_xc("functionals/CIDER23X_NL_MGGA_DTR.yaml", use_paw=not use_pp) + run_nl_feature_test(xc, spinpol=SPINPOL, use_pp=use_pp, baseline=baseline) + + baseline = "0.75_GGA_X_PBE+1.00_GGA_C_PBE" + xc = get_xc("functionals/CIDER23X_NL_MGGA_DTR.yaml", use_paw=not use_pp) + run_nl_feature_test(xc, spinpol=True, use_pp=use_pp, baseline=baseline) + + def test_eigval_vs_fd(self): + run_fd_deriv_test("PBE") + + run_fd_deriv_test("MGGA_X_R2SCAN+MGGA_C_R2SCAN") + + xc = get_xc("functionals/CIDER23X_NL_GGA.yaml") + run_fd_deriv_test(xc) + + xc = get_xc("functionals/CIDER23X_NL_MGGA_DTR.yaml") + run_fd_deriv_test(xc) + + xc = get_xc("functionals/CIDER23X_NL_MGGA_DTR.yaml") + run_fd_deriv_test(xc, spinpol=True) + + xc = get_xc("functionals/CIDER23X_SL_GGA.yaml") + run_fd_deriv_test(xc) + + xc = get_xc("functionals/CIDER23X_SL_MGGA.yaml") + run_fd_deriv_test(xc) + + xc = get_xc("functionals/CIDER23X_NL_GGA.yaml", use_paw=False) + run_fd_deriv_test(xc, use_pp=True) + + xc = get_xc("functionals/CIDER23X_NL_MGGA.yaml", use_paw=False) + run_fd_deriv_test(xc, use_pp=True) + + +if __name__ == "__main__": + unittest.main()