condopacpint.f

1. (changes results, small) Checked by KC (Nov 10, 2026) The mixture's derivatives are multiplied by $\log K$.

(line 68↗; line 72↗) QODC = CONDL * (WtH1*QODCX + …) multiplies the weighted derivatives by CONDL = $\log K$ itself. The derivative of the harmonic sum $1/K = \sum_k w_k/K_k$ is $\partial\log K/\partial\log\rho = \sum_k (w_k K/K_k)\,\partial\log K_k/\partial\log\rho$, with no factor $\log K$: For pure hydrogen ($w_\mathrm{H} = 1$) the code gives $\log K$ times the right value. The commented-out older lines (lines 69–70↗, lines 73–74↗) weight each derivative by its CONDX, as the fix below does, but also multiply by CONDL. What follows: QODC and QOTC are wrong wherever Potekhin's opacity is computed ($\log\rho \ge -6$), and GETOPAC folds them into the opacity's derivatives QOD and QOT with the weight $\kappa_\mathrm{rad}/(\kappa_\mathrm{rad} + \kappa_\mathrm{cond})$ (getopac.f, lines 219–222↗), so the error counts where conduction matters (degenerate cores). Only with LcondOpacP (off by default, set in all 51 decks of the test suite and recommended by the user guide).

A possible fix:

condopacpint.f

@@ lines 65–71 @@
 c           CONDL = LOG10( (X + 1D0) / 2D0 ) -
 c     *             LOG10( X*CONDX + .5D0*Y*CONDY + .5D0*z*CONDZ)

-           QODC = CONDL * (WtH1*QODCX + WtHe4*QODCY + WtOx*QODCZ)
+C 11/09/26 KC via Claude
+C code changed to take the derivative of the harmonic sum: each
+C composition weighted by its share of 1/K, WT*CONDX*K, K = 10**CONDL.
+           QODC = 10.0D0**CONDL*(WtH1*CONDX*QODCX + WtHe4*CONDY*QODCY
+     *            + WtOx*CONDZ*QODCZ)
 c           QODC = 2D0 * CONDL / ( 1D0 + X ) *
 c     *          (X*CONDX*QODCX +.5D0*Y*CONDY*QODCY +.5D0*Z*COBDZ*QODCZ)

How the code itself shows it. The routine's own lines fix what QODC must be. CondOpacP, which it calls, defines its outputs as $\log_{10}$ of the thermal conductivity $K$ and plain log derivatives (condopacp.f, lines 22–24↗: "DRK - d log kappa / d log rho," "DTK - d log kappa / d log T"). lines 76–77↗ build the opacity as OCL = -3.5194D0+3.0D0*TL-DL-CONDL, that is $\kappa_\mathrm{cond} = 16\sigma T^3/(3\rho K)$, and its density derivative as -1.0D0-QODC, so QODC must be $\partial$CONDL$/\partial\log\rho$: For pure, fully ionized hydrogen that is hydrogen's DRK itself, while line 68↗ gives CONDL times it. The same routine's branch below $\log\rho = -6$ uses DRK exactly so, with no factor CONDL, both as the derivative and as the slope it extrapolates along (line 84↗, lines 97–99↗). YREC's other conductive opacity, Hubbard and Lampe's in GETOPAC, returns QODC = 0.3901D0*TL - 0.2772D0*DL - 3.5853D0 (getopac.f, lines 206–211↗), exactly the density derivative of its own formula, and GETOPAC combines either option's QODC with the radiative one in the same way (getopac.f, lines 219–222↗); both options must therefore deliver the same quantity, and only this branch does not. The factor is not small: CONDL, $\log_{10}K$, runs from about 1.1 to 19.3 in the table condall06.d. The combined derivatives reach the iteration's matrix (coefft.f, lines 263–264↗) and, in rotating runs, the meridional circulation through $\partial\ln\chi/\partial\ln T = 3 - \partial\ln\kappa/\partial\ln T$ (setupv.f, lines 431–433↗), each weighted by $\kappa_\mathrm{rad}/(\kappa_\mathrm{rad}+\kappa_\mathrm{cond})$, which goes to 1 where conduction dominates: The error counts most in degenerate cores, which no run with the fix reached. The fix above is the derivative of the harmonic sum the routine builds at line 64↗; the commented-out older lines (lines 69–70↗) already weight each element by its CONDX.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): 34 runs, 26 of YREC's cases (from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning) and 8 made-up ones (each one of YREC's cases with a few settings changed), of which one stopped in its first second. The step ran in all 33 others, from 0.06 million to 74 million times a run, at $\log T$ = 3 to 8.40 and $\log\rho$ = −6 to 7.45, and the derivative differed from the right one by more than 0.001 in nearly every call (74,434,834 of 74,434,834 in testsuite/Test_m9p0feh+0p0GS98_base_TAHB). With the fix, testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS moves by up to $1.4\times10^{-5}$ dex in $\log L$ and $7.1\times10^{-6}$ dex in $\log T_c$, against $1\times10^{-9}$ and $8\times10^{-11}$ for the null change used as a noise floor (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$): real, but small (259 models in both, end age 5.23952 Gyr instead of 5.23984). On the 3 $M_\odot$ giant testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2, the change is within the noise: $\log L$ moves by up to 0.0037 dex and $\log T_c$ by 0.0015, against 0.0076 and 0.0008 for the null change (1,847 models instead of 1,853). Not run with the fix: a degenerate core (testsuite/Test_m1p0feh+0p0GS98_base_HeIgnite), where conduction matters most. QOD and QOT enter the solver's derivatives (TPGRAD, COEFFT), which fits so small a change. Intent: The commented-out older lines also multiply by CONDL.

2. (changes results, small) Checked by KC (Nov 10, 2026) The temperature derivative is taken from the density derivatives.

(line 72↗) QOTC is computed from QODCX, QODCY, QODCZ, not from QOTCX, QOTCY, QOTCZ (which the calls at lines 56–58↗ return and nothing reads), so $\partial\log\kappa/\partial\log T$ is $3 - \partial\log K/\partial\log\rho$ (and multiplied by $\log K$, the factor of suspected bug 1). The commented-out older lines lines 73–74↗ used QOTCX, QOTCY, QOTCZ. Same path as suspected bug 1: Potekhin's opacity with LcondOpacP (off by default, set in all 51 decks of the test suite), at $\log\rho \ge -6$.

A possible fix (in the form of the fix of suspected bug 1, each composition weighted by its share of $1/K$, with no factor $\log K$):

condopacpint.f

@@ lines 69–75 @@
 c           QODC = 2D0 * CONDL / ( 1D0 + X ) *
 c     *          (X*CONDX*QODCX +.5D0*Y*CONDY*QODCY +.5D0*Z*COBDZ*QODCZ)

-            QOTC = CONDL * (WtH1*QODCX + WtHe4*QODCY + WtOx*QODCZ)
+C 11/09/26 KC via Claude
+C code changed to take the temperature derivatives QOTCX, QOTCY and
+C QOTCZ, weighted as the density derivatives are.
+            QOTC = 10.0D0**CONDL*(WtH1*CONDX*QOTCX + WtHe4*CONDY*QOTCY
+     *             + WtOx*CONDZ*QOTCZ)
 c           QOTC = 2D0 * CONDL / ( 1D0 + X ) *
 c     *          (X*CONDX*QOTCX +.5D0*Y*CONDY*QOTCY +.5D0*Z*CONDZ*QOTCZ)

How the code itself shows it. QOTC must be $\partial$CONDL$/\partial\log T$: line 76↗ builds the opacity as OCL = -3.5194D0+3.0D0*TL-DL-CONDL and line 78↗ its temperature derivative as 3.0D0-QOTC. CondOpacP returns that derivative for each element as its last argument, "DTK - d log kappa / d log T" (condopacp.f, line 24↗), and the calls at lines 56–58↗ store it in QOTCX, QOTCY, and QOTCZ, which nothing reads; line 72↗ uses the density derivatives QODCX, QODCY, and QODCZ instead, the same expression as QODC on line 68↗. The commented-out older line it replaced (lines 73–74↗) reads QOTCX, QOTCY, and QOTCZ. YREC's other conductive opacity, Hubbard and Lampe's in GETOPAC, returns QOTC = 0.3901D0*DL - 0.6438D0*TL + 5.1324D0 (getopac.f, lines 206–211↗), the temperature derivative of its own formula, which GETOPAC combines with the radiative one as it does this routine's (getopac.f, lines 219–222↗).

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The step ran as for suspected bug 1 (in 33 of the 34 runs, from 0.06 million to 74 million times a run), and the temperature and density derivatives differed by more than 0.001 in nearly every call (74,434,834 of 74,434,834 in testsuite/Test_m9p0feh+0p0GS98_base_TAHB). With this bug's fix alone (the temperature derivatives, still times CONDL; measurement only), testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS moves by up to $5.0\times10^{-6}$ dex in $\log L$ and $6.0\times10^{-6}$ dex in $\log T_c$, against $1\times10^{-9}$ and $8\times10^{-11}$ for the null change: real, but small (259 models in both, end age 5.24007 Gyr instead of 5.23984). On the 3 $M_\odot$ giant the change is within the noise: $\log L$ moves by up to 0.0050 dex and $\log T_c$ by 0.00066, against 0.0076 and 0.0008 for the null change. Not run with the fix: a degenerate core. Intent: The commented-out older lines lines 73–74↗ took QOTCX, QOTCY, and QOTCZ.

3. (changes results, small) Checked by KC (Nov 10, 2026) A sign in the extrapolation below $\log\rho = -6$.

(line 92↗; line 64↗) In step 3, CONDL = LOG10(…), while in step 2 CONDL = -LOG10(…), with the same sum inside. Since $\sum_k w_k/K_k = 1/K$, step 2 gives $\log K$ but step 3 gives $-\log K$, so the extrapolated $\log\kappa_\mathrm{cond}$ (line 96↗) is off by $2\log K$. Probably a missing minus sign. The path: LcondOpacP, $\log\rho < -6$ and $\log R \ge 0$, so only for $\log T < 4$ (cool surface layers); GETOPAC uses the result whenever LCONDO is true (getopac.f, line 197↗). How much the total opacity changes there was not traced.

A possible fix:

condopacpint.f

@@ lines 89–95 @@
            CONDY = 10.0D0**(-CONDYL)
            CONDZ = 10.0D0**(-CONDZL)

-           CONDL = LOG10(WtH1*CONDX + WtHe4*CONDY + WtOx*CONDZ)
+C 11/09/26 KC via Claude
+C code changed to give CONDL the sign it has above, log K.
+           CONDL = -LOG10(WtH1*CONDX + WtHe4*CONDY + WtOx*CONDZ)
 c           CONDL = LOG10( (X + 1D0) / 2D0 ) -
 c     *            LOG10( X*CONDX + .5D0*Y*CONDY + .5D0*Z*CONDZ)

How the code itself shows it. The two steps compute the same quantity: lines 88–90↗ repeat lines 60–62↗, CONDX = 10.0D0**(-CONDXL) and its twins, so each CONDX is $1/K_k$ and the sum inside both logarithms is $1/K$. Step 2 takes CONDL = -LOG10(…) (line 64↗), and so does the older formula commented out under each step, CONDL = LOG10( (X + 1D0) / 2D0 ) - followed by LOG10( X*CONDX + .5D0*Y*CONDY + .5D0*Z*CONDZ) (lines 93–94↗; lines 65–66↗ the same, with z): a minus before the logarithm of the sum. Step 3's opacity, OCL = -3.5194D0+3.0D0*TL+6.0D0-CONDL (line 96↗), is step 2's OCL = -3.5194D0+3.0D0*TL-DL-CONDL (line 76↗) at DL = −6, which needs the same CONDL, $\log_{10}K$; with the opposite sign the two opacities do not meet at $\log\rho = -6$, where the extrapolation starts. The fix gives line 92↗ the sign of line 64↗ and of both commented formulas.

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The extrapolation ran in 21 of the 33 runs, from 41,369 to 3.6 million times a run, at $\log T$ = 3 to 4.00 and $\log\rho$ = −8.99 to −6, mostly the low-mass and solar models; CONDL had the wrong sign in every call. With the fix, the 0.04 $M_\odot$ brown dwarf examples/run_from_dbl_to_zams/Test_m0040_feh0_noGS_ZAMS moves by up to $2.4\times10^{-4}$ dex in $\log L$ and $3.6\times10^{-5}$ dex in $\log T_\mathrm{eff}$, against $2.5\times10^{-5}$ and $1.2\times10^{-5}$ for the null change: about ten times the noise in $\log L$, so real, but small (748 models in both, end age 8,692 Gyr instead of 8,681, against 8,687 for the null change). On testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS the change is $9\times10^{-9}$ dex in $\log L$, negligible. Intent: none found; the commented-out older line (lines 93–94↗) subtracts the log of the sum, as line 64↗ does, so only line 92↗ has the other sign.

4. (harmless in the runs) The step of the extrapolation is reversed.

(line 99↗) OCL = OCL + QODC*(-6.0D0-DL) adds the slope times $(-6 - \log\rho)$, where a linear extrapolation from $\log\rho = -6$ adds the slope times $(\log\rho + 6)$; the extrapolated opacity moves the wrong way as the density falls. Same path as suspected bug 3: LcondOpacP, $\log\rho < -6$ and $\log R \ge 0$, so only for $\log T < 4$ (cool surface layers).

A possible fix:

condopacpint.f

@@ lines 96–102 @@
            OCL = -3.5194D0+3.0D0*TL+6.0D0-CONDL
            QODC = -1.0D0-QODC
            QOTC = 3.0D0-QOTC
-           OCL = OCL + QODC*(-6.0D0-DL)
+C 11/09/26 KC via Claude
+C code changed to extrapolate from log rho = -6 down to DL.
+           OCL = OCL + QODC*(DL+6.0D0)
            OC = 10.0D0**OCL
            LCONDO=.TRUE.
          END IF

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): The step ran as for suspected bug 3 (in 21 of the 33 runs, from 41,369 to 3.6 million times a run), and the extrapolation went the wrong way by more than 0.01 dex in nearly every call (3,629,479 of 3,630,172 in testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS). With the fix, the brown dwarf examples/run_from_dbl_to_zams/Test_m0040_feh0_noGS_ZAMS moves by up to $9.5\times10^{-5}$ dex in $\log L$ and $7.6\times10^{-6}$ dex in $\log T_c$, against $2.5\times10^{-5}$ and $1.5\times10^{-6}$ for the null change: four to five times the noise, so small, and not clearly real (748 models in both). The track of testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS is identical with the fix. The extrapolated opacity is wrong in nearly every call, silently, but even the brown dwarf, which takes the path most often (3.2 million times), shows no change clearly above the noise. Intent: none found.