Skip to content

VdW features implementation #21

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
msamytawfik wants to merge 3 commits into
base: main
Choose a base branch
from
Open

VdW features implementation #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5e54e78
Implemented vdW features for CIDER. Pending tests.
msamytawfik May 1, 2025
74c326b
Fixed some bugs with the c code of the vdw features
msamytawfik May 1, 2025
8d1d6ac
Fixed some issues with vdw features; now they should be evaluating co…
msamytawfik May 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 30 additions & 6 deletions ciderpress/dft/plans.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,8 @@
"se_ar2": 1,
"se_a2r4": 2,
"se_erf_rinv": 3,
"rinv2": 4, # Placeholder for H_i term (1 / (a_i r^2 + 1))
"rinv4": 5, # Placeholder for G_i/H_i term (1 / (a_0 r^2 + 1))^2
}
VI_ID_MAP = {
"se": 0,
Expand All @@ -62,7 +64,7 @@
"se_grad": 7,
}


#TODO Apr25: this is where -1 default is used
def _get_ovlp_fit_interpolation_coefficients(
plan, arg_g, i=-1, local=True, vbuf=None, dbuf=None
):
Expand All @@ -83,10 +85,32 @@ def _get_ovlp_fit_interpolation_coefficients(
p = plan.empty_coefs(ngrids, local=local, buf=vbuf)
dp = plan.empty_coefs(ngrids, local=local, buf=dbuf)
if i == -1:
feat_id = VJ_ID_MAP["se"]
if not hasattr(plan.nldf_settings, 'vdw_param'):
feat_id = VJ_ID_MAP["se"]
elif not plan.nldf_settings.vdw_param:
feat_id = VJ_ID_MAP["se"]
else:
feat_id = VJ_ID_MAP["rinv4"]
# if plan.nldf_settings.feat_spec_list[i] == "rinv2_rinv4":
# feat_id = VJ_ID_MAP["rinv4"]
# elif plan.nldf_settings.feat_spec_list[i] == "se_rinv4":
# feat_id = VJ_ID_MAP["rinv4"]
# else:
# feat_id = VJ_ID_MAP["se"]
else:
assert 0 <= i < plan.nldf_settings.num_feat_param_sets
feat_id = VJ_ID_MAP[plan.nldf_settings.feat_spec_list[i]]
if not hasattr(plan.nldf_settings, 'vdw_param'):
feat_id = VJ_ID_MAP[plan.nldf_settings.feat_spec_list[i]]
elif not plan.nldf_settings.vdw_param:
feat_id = VJ_ID_MAP[plan.nldf_settings.feat_spec_list[i]]
else:
if plan.nldf_settings.feat_spec_list[i] == "rinv2_rinv4":
feat_id = VJ_ID_MAP["rinv2"]
elif plan.nldf_settings.feat_spec_list[i] == "se_rinv4":
feat_id = VJ_ID_MAP["se"]
else:
feat_id = VJ_ID_MAP[plan.nldf_settings.feat_spec_list[i]]
print(f"HERE: {feat_id}")
if feat_id == VJ_ID_MAP["se_erf_rinv"]:
extra_args = plan._get_extra_args(i)
num_extra_args = len(extra_args)
Expand Down Expand Up @@ -117,7 +141,7 @@ def _get_ovlp_fit_interpolation_coefficients(
ctypes.c_int(len(alphas)),
ctypes.c_int(feat_id),
extra_args,
)
) #TODO Apr25: c functions call
return p, dp


Expand Down Expand Up @@ -1506,7 +1530,7 @@ def get_drhodf_tuple(self, rho_data, drhodf_data, with_spin=False):
is_mgga=self.nldf_settings.sl_level == "MGGA",
)

def eval_feat_exp(self, rho_tuple, i=-1):
def eval_feat_exp(self, rho_tuple, i=-1): #TODO Apr25: check this
"""
Evaluate the exponents that determine the length-scale
of feature i at each grid point. i can take a range
Expand Down Expand Up @@ -1896,7 +1920,7 @@ def get_interpolation_coefficients(self, arg_g, i=-1, vbuf=None, dbuf=None):
ctypes.c_int(self.nalpha),
)
return p, dp
else:
else: #TODO Apr25: this is where coefficients are obtained
return self._get_interpolation_coefficients(
arg_g, i=i, vbuf=vbuf, dbuf=dbuf
)
Expand Down
41 changes: 36 additions & 5 deletions ciderpress/dft/settings.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
from abc import ABC, abstractmethod

import numpy as np
from scipy.special import gamma as gamma_func
from scipy.special import gamma as gamma_func, erfc

from ciderpress.dft.feat_normalizer import (
ConstantNormalizer,
Expand Down Expand Up @@ -953,8 +953,9 @@ def get_reasonable_normalizer(self):

se_rvec: squared-exponential times vector (r'-r)
"""

ALLOWED_J_SPECS = ["se", "se_ar2", "se_a2r4", "se_erf_rinv"]
#TODO: separate rinv4, rinv2
ALLOWED_J_SPECS = ["se", "se_ar2", "se_a2r4", "se_erf_rinv",
"se_rinv4", "rinv2_rinv4"] # Added vdW placeholders
"""
Allowed spec strings for version j features.
Version k features have the same allowed spec strings
Expand All @@ -967,6 +968,15 @@ def get_reasonable_normalizer(self):

se_erf_rinv: squared-exponential * 1/r with short-range
erf damping

se_rinv4: squared-exponential * (1 / (a r^2 + 1))^2 for G_i feature

rinv2_rinv4: (1 / (a r^2 + 1)) * (1 / (a r^2 + 1))^2 for H_i feature

rinv2: (1 / (a r^2 + 1))

rinv4: (1 / (a r^2 + 1))^2

"""
ALLOWED_K_SPECS = ALLOWED_J_SPECS
"""
Expand Down Expand Up @@ -1008,6 +1018,8 @@ def get_reasonable_normalizer(self):
"se_grad": 1,
"se_rvec": -1,
"grad_rho": 4,
"se_rinv4": 0,
"rinv2_rinv4": 0,
}
"""
Uniform-scaling powers (USPs) descibe how features scale as the
Expand Down Expand Up @@ -1274,7 +1286,7 @@ def get_reasonable_normalizer(self):
raise NotImplementedError
return norms


#TODO: feat_specs, theta_spec default, -1 is the theta params when call_coeffs (passed integer)
class NLDFSettingsVJ(NLDFSettings):

version = "j"
Expand All @@ -1285,7 +1297,7 @@ def __init__(
theta_params,
rho_mult,
feat_specs,
feat_params,
feat_params, vdw_param = False, #TODO Apr25: should add optional argument for vdW feat_params
):
"""
Initialize NLDFSettingsVJ
Expand Down Expand Up @@ -1316,6 +1328,7 @@ def __init__(
super(NLDFSettingsVJ, self).__init__(sl_level, theta_params, rho_mult)
self.feat_params = feat_params
self.feat_specs = feat_specs
self.vdw_param = vdw_param
self._check_specs(self.feat_specs, ALLOWED_J_SPECS)
if len(self.feat_params) != len(self.feat_specs):
raise ValueError("specs and params must have same length")
Expand Down Expand Up @@ -1369,6 +1382,24 @@ def ueg_vector(self, rho=1.0):
elif spec == "se_erf_rinv":
expnt3 = expnt2 * params[-1]
integral *= np.sqrt(expnt / (expnt + expnt3))
elif spec == "rinv2_rinv4":
a_i = expnt2
a_0 = _get_ueg_expnt(a0t, t0t, rho)
if a_i < 1e-12 or a_0 < 1e-12:
integral = 0.0
else:
integral = np.pi**2 / ((np.sqrt(a_i / a_0) + 1)**2 * a_0**1.5)
elif spec == "se_rinv4":
a_i = expnt2
a_0 = _get_ueg_expnt(a0t, t0t, rho)
if a_i < 1e-12 or a_0 < 1e-12:
integral = 0.0
else:
ratio = a_i / a_0
sqrt_ratio = np.sqrt(ratio)
term1 = (2 * a_i + a_0) * np.exp(ratio) * erfc(sqrt_ratio) / (a_0 ** 2.5)
term2 = 2.0 * sqrt_ratio / (np.sqrt(np.pi) * a_0 ** 2)
integral = np.pi * np.pi * (term1 - term2)
else:
raise ValueError
ueg_feats.append(rho * rho_mult * integral)
Expand Down
76 changes: 76 additions & 0 deletions ciderpress/lib/mod_cider/cider_coefs.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,9 @@
#include <stdio.h>
#include <stdlib.h>

const double SQRT_PI = 1.77245385090551602729;
const double FOUR_PI = 4.0 * M_PI;

#define CIDER_FEAT_R0_GAUSSIAN 0
#define FILL_CIDER_R0_GAUSSIAN(ind) \
tmp = 1.0 / (exp_g[g] + alphas[a]); \
Expand Down Expand Up @@ -51,6 +54,67 @@
dp[ind] = -1.5 * p[ind] * tmp; \
tmp2 = tmp2 * tmp; \
dp[ind] -= 0.5 * p[ind] * tmp2 * tmp2 * extra_args[0] * alphas[a];
// Kernel: g(a, r) = 1 / (a r^2 + 1)
#define CIDER_FEAT_RINV2_GAUSSIAN 4
#define FILL_CIDER_RINV2_GAUSSIAN(ind) \
double aval = exp_g[g]; \
double alpha = alphas[a]; \
double aval_safe = fmax(aval, 1e-20); \
double alpha_safe = fmax(alpha, 1e-20); \
double sqrt_alpha = sqrt(alpha_safe); \
double sqrt_aval = sqrt(aval_safe); \
double a_pow_1_5 = aval_safe * sqrt_aval; \
double a_pow_2 = aval_safe * aval_safe; \
double a_pow_2_5 = a_pow_2 * sqrt_aval; \
double a_pow_3 = a_pow_2 * aval_safe; \
double a_pow_3_5 = a_pow_3 * sqrt_aval; \
double ratio = alpha_safe / aval_safe; \
double sqrt_ratio = sqrt(ratio); \
double exp_ratio = exp(ratio); \
double erfc_term = erfc(sqrt_ratio); \
double exp_erfc_term = exp_ratio * erfc_term; \
double term1_p = SQRT_PI / (2.0 * aval_safe * sqrt_alpha); \
double term2_p = (M_PI * exp_erfc_term) / (2.0 * a_pow_1_5); \
p[ind] = FOUR_PI * (term1_p - term2_p); \
double term1_dp = (M_PI * alpha_safe * exp_erfc_term) / (2.0 * a_pow_3_5); \
double term2_dp = (3.0 * M_PI * exp_erfc_term) / (4.0 * a_pow_2_5); \
double term3_dp = (SQRT_PI * sqrt_alpha) / (2.0 * a_pow_3); \
double term4_dp = SQRT_PI / (2.0 * a_pow_2 * sqrt_alpha); \
dp[ind] = FOUR_PI * (term1_dp + term2_dp - term3_dp - term4_dp);
// Kernel: g(a, r) = (1 / (a r^2 + 1))^2
#define CIDER_FEAT_RINV4_GAUSSIAN 5
#define FILL_CIDER_RINV4_GAUSSIAN(ind) \
double aval = exp_g[g]; \
double alpha = alphas[a]; \
double aval_safe = fmax(aval, 1e-20); \
double alpha_safe = fmax(alpha, 1e-20); \
double sqrt_alpha = sqrt(alpha_safe); \
double sqrt_aval = sqrt(aval_safe); \
double ratio = alpha_safe / aval_safe; \
double sqrt_ratio = sqrt(ratio); \
double exp_ratio = exp(ratio); \
double erfc_term = erfc(sqrt_ratio); \
double exp_erfc_term = exp_ratio * erfc_term; \
double a_plus_2alpha = aval_safe + 2.0 * alpha_safe; \
double a_pow_0_5 = sqrt_aval; \
double a_pow_1_5 = aval_safe * a_pow_0_5; \
double a_pow_2 = aval_safe * aval_safe; \
double a_pow_2_5 = a_pow_2 * a_pow_0_5; \
double a_pow_3_5 = a_pow_2_5 * aval_safe; \
double a_pow_3 = a_pow_2 * aval_safe; \
double term1_p_num = M_PI * exp_erfc_term * a_plus_2alpha; \
double term2_p_num = 2.0 * SQRT_PI * sqrt_aval * sqrt_alpha; \
double common_den_p = 4.0 * a_pow_2_5; \
p[ind] = FOUR_PI * (term1_p_num - term2_p_num) / common_den_p; \
double term1_d1 = SQRT_PI * sqrt_alpha * a_plus_2alpha / a_pow_1_5; \
double term2_d1 = M_PI * alpha_safe * exp_erfc_term * a_plus_2alpha / a_pow_2; \
double term3_d1 = SQRT_PI * sqrt_alpha / a_pow_0_5; \
double term4_d1 = M_PI * exp_erfc_term; \
double num_d1 = term1_d1 - term2_d1 - term3_d1 + term4_d1; \
double den_d1 = common_den_p; \
double num_d2 = term1_p_num - term2_p_num; \
double den_d2 = 8.0 * a_pow_3_5; \
dp[ind] = FOUR_PI * ( (num_d1 / den_d1) - (5.0 * num_d2 / den_d2) );

#define CIDER_GQ_LOOP(FEATNAME) \
for (g = 0; g < ngrids; g++) { \
Expand Down Expand Up @@ -100,6 +164,12 @@ void cider_coefs_gto_gq(double *p_ga, double *dp_ga, double *exp_g,
case CIDER_FEAT_ERF_GAUSSIAN:
#pragma omp for
CIDER_GQ_LOOP(ERF_GAUSSIAN);
case CIDER_FEAT_RINV2_GAUSSIAN:
#pragma omp for
CIDER_GQ_LOOP(RINV2_GAUSSIAN);
case CIDER_FEAT_RINV4_GAUSSIAN:
#pragma omp for
CIDER_GQ_LOOP(RINV4_GAUSSIAN);
default:
printf("INTERNAL CIDER ERROR\n");
}
Expand Down Expand Up @@ -131,6 +201,12 @@ void cider_coefs_gto_qg(double *p_ag, double *dp_ag, double *exp_g,
case CIDER_FEAT_ERF_GAUSSIAN:
#pragma omp for
CIDER_QG_LOOP(ERF_GAUSSIAN);
case CIDER_FEAT_RINV2_GAUSSIAN:
#pragma omp for
CIDER_QG_LOOP(RINV2_GAUSSIAN);
case CIDER_FEAT_RINV4_GAUSSIAN:
#pragma omp for
CIDER_QG_LOOP(RINV4_GAUSSIAN);
default:
printf("INTERNAL CIDER ERROR\n");
}
Expand Down
2 changes: 1 addition & 1 deletion examples/pyscf/simple_calc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
xkernel="GGA_X_PBE",
ckernel="GGA_C_PBE",
# exact exchange mixing parameter
xmix=0.25,
xmix=1.0,
)
ks = ks.density_fit()
ks.with_df.auxbasis = "def2-universal-jfit"
Expand Down
Loading