diff --git a/ref/genxc/genxc.h b/ref/genxc/genxc.h
new file mode 100644
index 0000000000000000000000000000000000000000..4187d040afa43ed87139239c9f12a6e34c5014bb
--- /dev/null
+++ b/ref/genxc/genxc.h
@@ -0,0 +1,711 @@
+
+double mzk (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {
+ return ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*(-(r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))
+)) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*((r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))
+));
+}
+
+double func_mf_der_r0 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {
+ return (r0 + r1)*(((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*(-(r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.18028119628842601*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, -0.66666666666666674)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ )) - 0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ 0
+ )
+ : (
+ 1.3333333333333333*pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 0.33333333333333326)*((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1 || p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 0
+ )
+ : (
+ (r0 - r1)/pow(r0 + r1, 2) - 1/(r0 + r1)
+ ))
+ ))
+)) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*((r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ 0.027190191964917291*pow(M_PI, -1.6666666666666665)*pow(r0, -3.6666666666666665)*s0*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))/pow(0.023448539988350059*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 1, 2) - 0.18028119628842601*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, -0.66666666666666674)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ )) - 0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ 0
+ )
+ : (
+ 1.3333333333333333*pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 0.33333333333333326)*((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1 || p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 0
+ )
+ : (
+ -(r0 - r1)/pow(r0 + r1, 2) + 1.0/(r0 + r1)
+ ))
+ ))
+))) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*(-(r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))
+)) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*((r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))
+));
+}
+
+double func_mf_der_r1 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {
+ return (r0 + r1)*(((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*(-(r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ 0.027190191964917291*pow(M_PI, -1.6666666666666665)*pow(r1, -3.6666666666666665)*s2*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))/pow(0.023448539988350059*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 1, 2) - 0.18028119628842601*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, -0.66666666666666674)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ )) - 0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ 0
+ )
+ : (
+ 1.3333333333333333*pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 0.33333333333333326)*((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1 || p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 0
+ )
+ : (
+ (r0 - r1)/pow(r0 + r1, 2) + 1.0/(r0 + r1)
+ ))
+ ))
+)) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*((r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.18028119628842601*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, -0.66666666666666674)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ )) - 0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ 0
+ )
+ : (
+ 1.3333333333333333*pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 0.33333333333333326)*((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1 || p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 0
+ )
+ : (
+ -(r0 - r1)/pow(r0 + r1, 2) - 1/(r0 + r1)
+ ))
+ ))
+))) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*(-(r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))
+)) + ((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*((r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.54084358886527806*pow(M_PI, -0.33333333333333331)*(1.804 - 0.6464160000000001/(0.018852626150633447*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 0.80400000000000005))*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))
+));
+}
+
+double func_mf_der_s0 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {
+ return (r0 + r1)*((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*((r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.010196321986843986*pow(M_PI, -1.6666666666666665)*pow(r0, -2.6666666666666665)*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ (r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))/pow(0.023448539988350059*pow(M_PI, -1.3333333333333333)*pow(r0, -2.6666666666666665)*s0 + 1, 2)
+));
+}
+
+double func_mf_der_s1 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {
+ return 0;
+}
+
+double func_mf_der_s2 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {
+ return (r0 + r1)*((p_a_dens_threshold >= 0.5*pow(r0 + r1, 1.0)*(-(r0 - r1)/(r0 + r1) + 1)) ? (
+ 0
+)
+: (
+ -0.010196321986843986*pow(M_PI, -1.6666666666666665)*pow(r1, -2.6666666666666665)*pow(r0 + r1, 0.33333333333333331)*((((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold <= p_a_zeta_threshold
+ )
+ : (
+ ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold >= 2 - p_a_zeta_threshold
+ )
+ : (
+ p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1
+ ))
+ ))) ? (
+ pow(p_a_zeta_threshold, 1.3333333333333333)
+ )
+ : (
+ pow(((p_a_zeta_threshold >= -(r0 - r1)/(r0 + r1) + 1) ? (
+ p_a_zeta_threshold - 1
+ )
+ : ((p_a_zeta_threshold >= (r0 - r1)/(r0 + r1) + 1) ? (
+ 1 - p_a_zeta_threshold
+ )
+ : (
+ -(r0 - r1)/(r0 + r1)
+ ))) + 1, 1.3333333333333333)
+ ))/pow(0.023448539988350059*pow(M_PI, -1.3333333333333333)*pow(r1, -2.6666666666666665)*s2 + 1, 2)
+));
+}
+
+
+void func_pol_gen(
+ const xc_func_type *p, int order, const double *my_rho, const double *my_sigma, double *izk,
+ double *ivrho, double *ivsigma, double *iv2rho2, double *iv2rhosigma, double *iv2sigma2,
+ double *iv3rho3, double *iv3rho2sigma, double *iv3rhosigma2, double *iv3sigma3,
+ double *iv4rho4, double *iv4rho3sigma, double *iv4rho2sigma2, double *iv4rhosigma3, double *iv4sigma4) {
+ const double x0 = my_rho[0] - my_rho[1];
+ const double x1 = my_rho[0] + my_rho[1];
+ const double x2 = 1.0/x1;
+ const double x3 = x0*x2;
+ const double x4 = x3 + 1;
+ const double x5 = 0.5*pow(x1, 1.0);
+ const double x6 = p->dens_threshold >= x4*x5;
+ const double x7 = pow(x1, 0.33333333333333331);
+ const double x8 = pow(p->zeta_threshold, 1.3333333333333333);
+ const double x9 = p->zeta_threshold >= x4;
+ const double x10 = p->zeta_threshold <= p->zeta_threshold;
+ const double x11 = -x3;
+ const double x12 = x11 + 1;
+ const double x13 = p->zeta_threshold >= x12;
+ const double x14 = -p->zeta_threshold;
+ const double x15 = p->zeta_threshold >= x14 + 2;
+ const double x16 = ((x9) ? (
+ x10
+)
+: (
+ ((x13) ? (
+ x15
+ )
+ : (
+ x9
+ ))
+));
+ const double x17 = p->zeta_threshold - 1;
+ const double x18 = x14 + 1;
+ const double x19 = ((x9) ? (
+ x17
+)
+: ((x13) ? (
+ x18
+)
+: (
+ x3
+))) + 1;
+ const double x20 = ((x16) ? (
+ x8
+)
+: (
+ pow(x19, 1.3333333333333333)
+));
+ const double x21 = x20*x7;
+ const double x22 = pow(M_PI, -1.3333333333333333);
+ const double x23 = 0.018852626150633447*x22;
+ const double x24 = pow(my_rho[0], -2.6666666666666665);
+ const double x25 = my_sigma[0]*x24;
+ const double x26 = 1.804 - 0.6464160000000001/(x23*x25 + 0.80400000000000005);
+ const double x27 = pow(M_PI, -0.33333333333333331);
+ const double x28 = 0.54084358886527806*x27;
+ const double x29 = x26*x28;
+ const double x30 = p->dens_threshold >= x12*x5;
+ const double x31 = pow(my_rho[1], -2.6666666666666665);
+ const double x32 = my_sigma[2]*x31;
+ const double x33 = 1.804 - 0.6464160000000001/(x23*x32 + 0.80400000000000005);
+ const double x34 = ((x13) ? (
+ x10
+)
+: (
+ ((x9) ? (
+ x15
+ )
+ : (
+ x13
+ ))
+));
+ const double x35 = ((x13) ? (
+ x17
+)
+: ((x9) ? (
+ x18
+)
+: (
+ x11
+))) + 1;
+ const double x36 = ((x34) ? (
+ x8
+)
+: (
+ pow(x35, 1.3333333333333333)
+));
+ const double x37 = x33*x36;
+ const double x38 = x28*x7;
+ const double x39 = ((x30) ? (
+ 0
+)
+: (
+ -x37*x38
+)) + ((x6) ? (
+ 0
+)
+: (
+ -x21*x29
+));
+ const double x40 = 0.18028119628842601*pow(x1, -0.66666666666666674)*x27;
+ const double x41 = -x37*x40;
+ const double x42 = x13 || x9;
+ const double x43 = -x2;
+ const double x44 = x0/pow(x1, 2);
+ const double x45 = 1.3333333333333333*pow(x35, 0.33333333333333326);
+ const double x46 = x33*x38;
+ const double x47 = -x20*x26*x40;
+ const double x48 = 0.023448539988350059*x22;
+ const double x49 = x21/pow(x25*x48 + 1, 2);
+ const double x50 = pow(M_PI, -1.6666666666666665);
+ const double x51 = 0.027190191964917291*x50;
+ const double x52 = -x44;
+ const double x53 = 1.3333333333333333*pow(x19, 0.33333333333333326);
+ const double x54 = x29*x7;
+ const double x55 = x36*x7/pow(x32*x48 + 1, 2);
+ const double x56 = 0.010196321986843986*x50;
+ izk[0] = x39;
+ ivrho[0] = x1*x39;
+ ivrho[1] = x1*(((x30) ? (
+ 0
+)
+: (
+ x41 - x46*((x34) ? (
+ 0
+ )
+ : (
+ x45*((x42) ? (
+ 0
+ )
+ : (
+ x43 + x44
+ ))
+ ))
+)) + ((x6) ? (
+ 0
+)
+: (
+ pow(my_rho[0], -3.6666666666666665)*my_sigma[0]*x49*x51 + x47 - x54*((x16) ? (
+ 0
+ )
+ : (
+ x53*((x42) ? (
+ 0
+ )
+ : (
+ x2 + x52
+ ))
+ ))
+))) + x39;
+ ivsigma[0] = x1*(((x30) ? (
+ 0
+)
+: (
+ pow(my_rho[1], -3.6666666666666665)*my_sigma[2]*x51*x55 + x41 - x46*((x34) ? (
+ 0
+ )
+ : (
+ x45*((x42) ? (
+ 0
+ )
+ : (
+ x2 + x44
+ ))
+ ))
+)) + ((x6) ? (
+ 0
+)
+: (
+ x47 - x54*((x16) ? (
+ 0
+ )
+ : (
+ x53*((x42) ? (
+ 0
+ )
+ : (
+ x43 + x52
+ ))
+ ))
+))) + x39;
+ ivsigma[1] = x1*((x6) ? (
+ 0
+)
+: (
+ -x24*x49*x56
+));
+ ivsigma[2] = 0;
+}
+
diff --git a/ref/genxc/gga_x_pbe_genxc_funcs.h b/ref/genxc/gga_x_pbe_genxc_funcs.h
index 2871438baaff36caa682840fd114a161a0494c1c..d933083463fc73083415006b6192f41975e2aee3 100644
--- a/ref/genxc/gga_x_pbe_genxc_funcs.h
+++ b/ref/genxc/gga_x_pbe_genxc_funcs.h
@@ -1,6 +1,8 @@
#define maple2c_order 4
#define MAPLE2C_FLAGS (XC_FLAGS_I_HAVE_EXC | XC_FLAGS_I_HAVE_VXC | XC_FLAGS_I_HAVE_FXC | XC_FLAGS_I_HAVE_KXC | XC_FLAGS_I_HAVE_LXC)
+#include "genxc.h"
+
void func_unpol(
const xc_func_type *p, int order, const double *my_rho, const double *my_sigma, double *izk,
double *ivrho, double *ivsigma, double *iv2rho2, double *iv2rhosigma, double *iv2sigma2,
@@ -40,7 +42,6 @@ void func_unpol(
#endif
-
const gga_x_pbe_params *params = p->params;
const double mu = MU_PBE;
const double kappa = 0.8040;
@@ -215,6 +216,17 @@ void func_unpol(
}
void func_pol(
+ const xc_func_type *p, int order, const double *my_rho, const double *my_sigma, double *izk,
+ double *ivrho, double *ivsigma, double *iv2rho2, double *iv2rhosigma, double *iv2sigma2,
+ double *iv3rho3, double *iv3rho2sigma, double *iv3rhosigma2, double *iv3sigma3,
+ double *iv4rho4, double *iv4rho3sigma, double *iv4rho2sigma2, double *iv4rhosigma3, double *iv4sigma4) {
+ func_pol_gen(p, order,my_rho,my_sigma,izk,
+ ivrho,ivsigma,iv2rho2,iv2rhosigma,iv2sigma2,
+ iv3rho3,iv3rho2sigma,iv3rhosigma2,iv3sigma3,
+ iv4rho4,iv4rho3sigma,iv4rho2sigma2, iv4rhosigma3,iv4sigma4);
+}
+
+void func_pol_ref(
const xc_func_type *p, int order, const double *my_rho, const double *my_sigma, double *izk,
double *ivrho, double *ivsigma, double *iv2rho2, double *iv2rhosigma, double *iv2sigma2,
double *iv3rho3, double *iv3rho2sigma, double *iv3rhosigma2, double *iv3sigma3,
@@ -316,7 +328,6 @@ void func_pol(
#endif
-
const gga_x_pbe_params *params = (gga_x_pbe_params * )(p->params);
const double mu = 0.249;
const double kappa = 0.8040;
@@ -374,7 +385,8 @@ void func_pol(
t77 = 0.1e1 + kappa * (0.1e1 - kappa / t72);
t81 = my_piecewise3(t53, 0, -0.3e1 / 0.8e1 * t5 * t62 * t77);
if(izk != NULL && (p->info->flags & XC_FLAGS_HAVE_EXC))
- izk[0] = t52 + t81;
+// izk[0] = t52 + t81;
+ izk[0] = mzk(my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold);
#ifndef XC_DONT_COMPILE_VXC
@@ -409,7 +421,8 @@ void func_pol(
t129 = t5 * t126 * t77 / 0.8e1;
t131 = my_piecewise3(t53, 0, -0.3e1 / 0.8e1 * t5 * t122 * t77 - t129);
if(ivrho != NULL && (p->info->flags & XC_FLAGS_HAVE_VXC))
- ivrho[0] = t52 + t81 + t6 * (t115 + t131);
+// ivrho[0] = t52 + t81 + t6 * (t115 + t131);
+ ivrho[0] = func_mf_der_r0(my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold);
t135 = my_piecewise5(t10, 0, t14, 0, -t7 - t84);
t138 = my_piecewise3(t20, 0, 0.4e1 / 0.3e1 * t23 * t135);
@@ -428,21 +441,25 @@ void func_pol(
t164 = t159 * t63 * t162;
t168 = my_piecewise3(t53, 0, -0.3e1 / 0.8e1 * t5 * t150 * t77 - t129 + t155 * t164 / 0.24e2);
if(ivrho != NULL && (p->info->flags & XC_FLAGS_HAVE_VXC))
- ivrho[1] = t52 + t81 + t6 * (t144 + t168);
+// ivrho[1] = t52 + t81 + t6 * (t144 + t168);
+ ivrho[1] = func_mf_der_r1(my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold);
t171 = t28 * t33;
t173 = t105 * t171 * t39;
t176 = my_piecewise3(t1, 0, -t102 * t173 / 0.64e2);
if(ivrho != NULL && (p->info->flags & XC_FLAGS_HAVE_VXC))
- ivsigma[0] = t6 * t176;
+// ivsigma[0] = t6 * t176;
+ ivsigma[0] = func_mf_der_s0(my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold);
if(ivrho != NULL && (p->info->flags & XC_FLAGS_HAVE_VXC))
- ivsigma[1] = 0.0e0;
+// ivsigma[1] = 0.0e0;
+ ivsigma[1] = func_mf_der_s1(my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold);
t178 = t158 * t171 * t68;
t181 = my_piecewise3(t53, 0, -t155 * t178 / 0.64e2);
if(ivrho != NULL && (p->info->flags & XC_FLAGS_HAVE_VXC))
- ivsigma[2] = t6 * t181;
+// ivsigma[2] = t6 * t181;
+ ivsigma[2] = func_mf_der_s2(my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold);
#ifndef XC_DONT_COMPILE_FXC
diff --git a/src/main.py b/src/main.py
index 215ee12909723a7b3244c84d22832d25a1de4c5e..46d138746efcbe1456c6117941537edf22350e30 100644
--- a/src/main.py
+++ b/src/main.py
@@ -1,27 +1,56 @@
import sympy as sp
# import math functions and constants directly to allow switching to implementations of other packages easily
-from sympy import exp, log, sqrt
-from sympy import pi
+from sympy import log, sqrt
# rs, z, t = sp.symbols('rs, z, t')
-rs, z, xt, xs0, xs1 = sp.symbols('rs, z, xt, xs0, xs1')
+# rs, z, xt, xs0, xs1 = sp.symbols('rs, z, xt, xs0, xs1')
-p_a_zeta_threshold = 1 # TODO: arbitrary value; sp.Symbol('p_a_zeta_threshold')
+p_a_zeta_threshold = sp.Symbol('p_a_zeta_threshold')
+p_a_dens_threshold = sp.Symbol('p_a_dens_threshold')
# p->zeta_threshold
# define my_piecewise3(c, x1, x2) ((c) ? (x1) : (x2))
def my_piecewise3(cond, left, right):
- return sp.Piecewise((left, cond), (right, True))
+ # print(cond)
+ pw = sp.Piecewise((left, cond), (right, True))
+ # print(pw)
+ # print(sp.simplify(pw))
+ # print()
+ # return sp.simplify(pw)
+ return pw
# define my_piecewise5(c1, x1, c2, x2, x3) ((c1) ? (x1) : ((c2) ? (x2) : (x3)))
+def my_piecewise5(c1, x1, c2, x2, x3):
+ return sp.Piecewise((x1, c1), (x2, c2), (x3, True))
# from https://gitlab.com/libxc/libxc/-/blob/master/maple/util.mpl
+X2S = 1 / (2 * (6 * sp.pi ** 2) ** sp.S(1 / 3))
+X2S_2D = 1 / (2 * (4 * sp.pi) ** sp.S(1 / 2))
+
+MU_GE = sp.S(10 / 81)
+X_FACTOR_C = 3 / 8 * (3 / sp.pi) ** (1 / 3) * 4 ** (2 / 3)
+X_FACTOR_2D_C = 8 / (3 * sp.sqrt(sp.pi))
+K_FACTOR_C = 3 / 10 * (6 * sp.pi ** 2) ** (2 / 3)
+
+# # 1D
+# DIMENSIONS = 1
+# RS_FACTOR = 1 / 2
+# # 2D
+# DIMENSIONS = 2
+# RS_FACTOR = 1 / sqrt(sp.pi)
+# LDA_X_FACTOR = -X_FACTOR_2D_C
+# 3D
+DIMENSIONS = 3
+RS_FACTOR = sp.S(3 / (4 * sp.pi)) ** sp.S(1 / 3)
+LDA_X_FACTOR = -X_FACTOR_C
+
+
def opz_pow_n(z, n):
return my_piecewise3(1 + z <= p_a_zeta_threshold, p_a_zeta_threshold ** n, (1 + z) ** n)
@@ -38,6 +67,46 @@ def tt(rs, z, xt):
return xt / (4 * 2 ** sp.S(1 / 3) * mphi(z) * sqrt(rs))
+def r_ws(n):
+ return RS_FACTOR / n ** sp.S(1 / DIMENSIONS)
+
+
+def n_total(rs):
+ return (RS_FACTOR / rs) ** DIMENSIONS
+
+
+def n_spin(rs, z):
+ return (1 + z) * n_total(rs) / 2
+
+
+# Screening for extreme values of zeta
+def z_thr(z):
+ return my_piecewise5(
+ 1 + z <= p_a_zeta_threshold, p_a_zeta_threshold - 1,
+ 1 - z <= p_a_zeta_threshold, 1 - p_a_zeta_threshold,
+ z)
+
+
+# Generate exchange and kinetic functionals from the expression for the
+# enhancement factor
+def lda_x_spin(rs, z):
+ return LDA_X_FACTOR * opz_pow_n(z, 1 + 1 / DIMENSIONS) * 2 ** (-1 - 1 / DIMENSIONS) * (RS_FACTOR / rs)
+
+
+def lda_k_spin(rs, z):
+ return K_FACTOR_C * opz_pow_n(z, 5 / 3) * 2 ** (-5 / 3) * (RS_FACTOR / rs) ** 2
+
+
+# Screen out small densities as well as the zero component from fully spin polarized densities due to terms of the form (1+z)^{-n}
+def screen_dens(rs, z):
+ return n_spin(rs, z) <= p_a_dens_threshold
+
+
+def gga_exchange(func, rs, z, xs0, xs1):
+ return (my_piecewise3(screen_dens(rs, z), 0, lda_x_spin(rs, z_thr(z)) * func(xs0))
+ + my_piecewise3(screen_dens(rs, -z), 0, lda_x_spin(rs, z_thr(-z)) * func(xs1)))
+
+
# from https://gitlab.com/libxc/libxc/-/blob/master/maple/lda_exc/lda_c_pw.mpl
# add None to honor Maple numbering starting with 1
@@ -71,45 +140,217 @@ def f_pw(rs, zeta):
# from https://gitlab.com/libxc/libxc/-/blob/master/maple/gga_exc/gga_c_pbe.mpl
-def mbeta(rs, t):
- return 0.06672455060314922
+# def mbeta(rs, t):
+# return 0.06672455060314922
+#
+#
+# mgamma = (1 - log(2)) / pi ** 2
+# BB = 1
+#
+#
+# def tp(rs, z, xt):
+# return tt(rs, z, xt)
+#
+#
+# # (* Equation (8) *)
+# def A(rs, z, t):
+# return mbeta(rs, t) / (mgamma * (exp(-f_pw(rs, z) / (mgamma * mphi(z) ** 3)) - 1))
+#
+#
+# # (* Equation (7) *) helpers
+# def f1(rs, z, t):
+# return t ** 2 + BB * A(rs, z, t) * t ** 4
+#
+#
+# def f2(rs, z, t):
+# return mbeta(rs, t) * f1(rs, z, t) / (mgamma * (1 + A(rs, z, t) * f1(rs, z, t)))
+#
+#
+# # (* Equation (7) *)
+# def fH(rs, z, t):
+# return mgamma * mphi(z) ** 3 * log(1 + f2(rs, z, t))
+#
+#
+# def f_pbe(rs, z, xt, xs0, xs1):
+# return f_pw(rs, z) + fH(rs, z, tp(rs, z, xt))
+#
+#
+# def f(rs, z, xt, xs0, xs1):
+# return f_pbe(rs, z, xt, xs0, xs1)
+
+
+# func = f(rs, z, xt, xs0, xs1)
+# print(func)
+# print(sp.ccode(func))
-mgamma = (1 - log(2)) / pi ** 2
-BB = 1
+# https://gitlab.com/libxc/libxc/-/blob/master/maple/gga_exc/gga_x_pbe.mpl
-tp = tt # tp := (rs, z, xt) -> tt(rs, z, xt) ?
+# (* standard PBE *)
+# $ifdef gga_x_pbe_params
+params_a_kappa = sp.S(0.8040)
+# params_a_mu = sp.S(0.2195149727645171)
+params_a_mu = sp.S(0.249) # TODO: why?
-# (* Equation (8) *)
-def A(rs, z, t):
- return mbeta(rs, t) / (mgamma * (exp(-f_pw(rs, z) / (mgamma * mphi(z) ** 3)) - 1))
+# # (* PBE_SOL *)
+# # $ifdef gga_x_pbe_sol_params
+# params_a_kappa = sp.S(0.8040)
+# params_a_mu = MU_GE
+#
+# # $ifdef gga_x_pbe_tca_params
+# params_a_kappa = sp.S(1.227)
+# params_a_mu = sp.S(0.2195149727645171)
-# (* Equation (7) *) helpers
-def f1(rs, z, t):
- return t ** 2 + BB * A(rs, z, t) * t ** 4
+def pbe_f0(s):
+ return 1 + params_a_kappa * (1 - params_a_kappa / (params_a_kappa + params_a_mu * s ** 2))
-def f2(rs, z, t):
- return mbeta(rs, t) * f1(rs, z, t) / (mgamma * (1 + A(rs, z, t) * f1(rs, z, t)))
+def pbe_f(x):
+ return pbe_f0(X2S * x)
-# (* Equation (7) *)
-def fH(rs, z, t):
- return mgamma * mphi(z) ** 3 * log(1 + f2(rs, z, t))
+def f(rs, z, xt, xs0, xs1):
+ return gga_exchange(pbe_f, rs, z, xs0, xs1)
-def f_pbe(rs, z, xt, xs0, xs1):
- return f_pw(rs, z) + fH(rs, z, tp(rs, z, xt))
+# https://gitlab.com/libxc/libxc/-/blob/master/scripts/maple2c.pl
+# zk is energy per unit particle
+def mzk(r0, r1, s0, s1, s2):
+ return f(r_ws(dens(r0, r1)), zeta(r0, r1), xt(r0, r1, s0, s1, s2), xs0(r0, r1, s0, s2), xs1(r0, r1, s0, s2))
-def f(rs, z, xt, xs0, xs1):
- return f_pbe(rs, z, xt, xs0, xs1)
+# mf is energy per unit volume
+def mf(r0, r1, s0, s1, s2):
+ return dens(r0, r1) * mzk(r0, r1, s0, s1, s2)
-func = f(rs, z, xt, xs0, xs1)
-print(func)
-print(sp.ccode(func))
-# https://gitlab.com/libxc/libxc/-/blob/master/maple/gga_exc/gga_x_pbe.mpl
+unpol = False
+if unpol:
+ def dens(r0, r1):
+ return r0
+
+
+ def zeta(r0, r1):
+ return 0
+
+
+ def xs0(r0, r1, sigma0, sigma2):
+ return sqrt(sigma0 / 4) / ((r0 / 2) ** (1 + sp.S(1 / DIMENSIONS)))
+
+
+ def xs1(r0, r1, sigma0, sigma2):
+ return sqrt(sigma0 / 4) / ((r0 / 2) ** (1 + sp.S(1 / DIMENSIONS)))
+
+
+ def xt(r0, r1, sigma0, sigma1, sigma2):
+ return sqrt(sigma0) / r0 ** (1 + sp.S(1 / DIMENSIONS))
+
+
+else: # pol
+ def dens(r0, r1):
+ return r0 + r1
+
+
+ def zeta(r0, r1):
+ return (r0 - r1) / (r0 + r1)
+
+
+ def xs0(r0, r1, sigma0, sigma2):
+ return sqrt(sigma0) / r0 ** (1 + sp.S(1 / DIMENSIONS))
+
+
+ def xs1(r0, r1, sigma0, sigma2):
+ return sqrt(sigma2) / r1 ** (1 + sp.S(1 / DIMENSIONS))
+
+
+ def xt(r0, r1, sigma0, sigma1, sigma2):
+ return sqrt(sigma0 + 2 * sigma1 + sigma2) / (r0 + r1) ** (1 + sp.S(1 / DIMENSIONS))
+
+# func = f(rs, z, xt, xs0, xs1)
+r0, r1, s0, s1, s2 = sp.symbols("r0, r1, s0, s1, s2")
+
+func_mzk = mzk(r0, r1, s0, s1, s2)
+func_mf = mf(r0, r1, s0, s1, s2)
+func_mf_der_r0 = func_mf.diff(r0)
+func_mf_der_r1 = func_mf.diff(r1)
+func_mf_der_s0 = func_mf.diff(s0)
+func_mf_der_s1 = func_mf.diff(s1)
+func_mf_der_s2 = func_mf.diff(s2)
+
+# func_mzk = sp.simplify(func_mzk)
+# func_mf = sp.simplify(func_mf)
+# func_mf_der_r0 = sp.simplify(func_mf_der_r0)
+# func_mf_der_r1 = sp.simplify(func_mf_der_r1)
+# func_mf_der_s0 = sp.simplify(func_mf_der_s0)
+# func_mf_der_s1 = sp.simplify(func_mf_der_s1)
+# func_mf_der_s2 = sp.simplify(func_mf_der_s2)
+
+unknowns = ["izk[0]", "ivrho[0]", "ivrho[1]", "ivsigma[0]", "ivsigma[1]", "ivsigma[2]"]
+updates = [func_mzk, func_mf, func_mf_der_r0, func_mf_der_r1, func_mf_der_s0, func_mf_der_s1, func_mf_der_s2]
+
+print(f"""
+double mzk (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {'{'}
+ return {sp.ccode(func_mzk)};
+{'}'}
+
+double func_mf_der_r0 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {'{'}
+ return {sp.ccode(func_mf_der_r0)};
+{'}'}
+
+double func_mf_der_r1 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {'{'}
+ return {sp.ccode(func_mf_der_r1)};
+{'}'}
+
+double func_mf_der_s0 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {'{'}
+ return {sp.ccode(func_mf_der_s0)};
+{'}'}
+
+double func_mf_der_s1 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {'{'}
+ return {sp.ccode(func_mf_der_s1)};
+{'}'}
+
+double func_mf_der_s2 (double r0, double r1, double s0, double s1, double s2, double p_a_zeta_threshold, double p_a_dens_threshold) {'{'}
+ return {sp.ccode(func_mf_der_s2)};
+{'}'}
+""")
+
+# my_rho[0], my_rho[1], my_sigma[0], my_sigma[1], my_sigma[2], p->zeta_threshold, p->dens_threshold
+var_mapping = {r0: sp.Symbol("my_rho[0]"), r1: sp.Symbol("my_rho[1]"),
+ s0: sp.Symbol("my_sigma[0]"), s1: sp.Symbol("my_sigma[1]"), s2: sp.Symbol("my_sigma[2]"),
+ p_a_zeta_threshold: sp.Symbol("p->zeta_threshold"), p_a_dens_threshold: sp.Symbol("p->dens_threshold")}
+
+# print(f"""
+# void func_pol_gen(
+# const xc_func_type *p, int order, const double *my_rho, const double *my_sigma, double *izk,
+# double *ivrho, double *ivsigma, double *iv2rho2, double *iv2rhosigma, double *iv2sigma2,
+# double *iv3rho3, double *iv3rho2sigma, double *iv3rhosigma2, double *iv3sigma3,
+# double *iv4rho4, double *iv4rho3sigma, double *iv4rho2sigma2, double *iv4rhosigma3, double *iv4sigma4) {'{'}
+# izk[0] = {sp.ccode(func_mzk.subs(var_mapping))};
+# ivrho[0] = {sp.ccode(func_mf_der_r0.subs(var_mapping))};
+# ivrho[1] = {sp.ccode(func_mf_der_r1.subs(var_mapping))};
+# ivsigma[0] = {sp.ccode(func_mf_der_s0.subs(var_mapping))};
+# ivsigma[1] = {sp.ccode(func_mf_der_s1.subs(var_mapping))};
+# ivsigma[2] = {sp.ccode(func_mf_der_s2.subs(var_mapping))};
+# {'}'}
+# """)
+
+replacements, reduced_updates = sp.cse(updates)
+
+print(f"""
+void func_pol_gen(
+ const xc_func_type *p, int order, const double *my_rho, const double *my_sigma, double *izk,
+ double *ivrho, double *ivsigma, double *iv2rho2, double *iv2rhosigma, double *iv2sigma2,
+ double *iv3rho3, double *iv3rho2sigma, double *iv3rhosigma2, double *iv3sigma3,
+ double *iv4rho4, double *iv4rho3sigma, double *iv4rho2sigma2, double *iv4rhosigma3, double *iv4sigma4) {'{'}""")
+
+for cs in replacements:
+ print(f" const double {sp.ccode(cs[0])} = {sp.ccode(cs[1].subs(var_mapping))};")
+
+for unknown, update in zip(unknowns, reduced_updates):
+ print(f" {unknown} = {sp.ccode(update.subs(var_mapping))};")
+
+print(f"""{'}'}
+""")