Fix internal DOF, now we're perfectly matching release version

espei/error_functions/non_equilibrium_thermochemical_error.py

@@ -18,6 +18,8 @@
 from pycalphad import Workspace
 from pycalphad.property_framework import IsolatedPhase
 from pycalphad.property_framework.metaproperties import find_first_compset
+from pycalphad.core.solver import Solver, SolverResult
+
 
 from espei.datasets import Dataset
 from espei.core_utils import ravel_conditions, get_prop_data, filter_temperatures
@@ -30,6 +32,45 @@
 
 _log = logging.getLogger(__name__)
 
+class NoSolveSolver(Solver):
+    def solve(self, composition_sets, conditions):
+        """
+        Initialize a solver, but don't actually do anything - i.e. return the SolverResult as if the calculation converged.
+
+        """
+        spec = self.get_system_spec(composition_sets, conditions)
+        self._fix_state_variables_in_compsets(composition_sets, conditions)
+        state = spec.get_new_state(composition_sets)
+        # DO NOT: converged = spec.run_loop(state, 1000)
+
+        if self.remove_metastable:
+            phase_idx = 0
+            compsets_to_remove = []
+            for compset in composition_sets:
+                # Mark unstable phases for removal
+                if compset.NP <= 0.0 and not compset.fixed:
+                    compsets_to_remove.append(int(phase_idx))
+                phase_idx += 1
+            # Watch removal order here, as the indices of composition_sets are changing!
+            for idx in reversed(compsets_to_remove):
+                del composition_sets[idx]
+
+        phase_amt = [compset.NP for compset in composition_sets]
+
+        x = composition_sets[0].dof
+        state_variables = composition_sets[0].phase_record.state_variables
+        num_statevars = len(state_variables)
+        for compset in composition_sets[1:]:
+            x = np.r_[x, compset.dof[num_statevars:]]
+        x = np.r_[x, phase_amt]
+        chemical_potentials = np.array(state.chemical_potentials)
+
+        if self.verbose:
+            print('Chemical Potentials', chemical_potentials)
+            print(np.asarray(x))
+        return SolverResult(converged=True, x=x, chemical_potentials=chemical_potentials)
+
+
 # TODO: make into an object similar to how ZPF data works?
 FixedConfigurationCalculationData = NewType("FixedConfigurationCalculationData", Dict[str, Any])
 
@@ -287,7 +328,8 @@ def get_thermochemical_data(dbf, comps, phases, datasets, model=None, weight_dic
                 model_cls = model.get(phase_name, Model)
                 mod = model_cls(dbf, comps, phase_name, parameters=symbols_to_fit)
                 if prop.endswith('_FORM'):
-                    output = ''.join(prop.split('_')[:-1])
+                    output = ''.join(prop.split('_')[:-1])+"R"
+                    # print("SHIFTING", output)
                     mod.shift_reference_state(ref_states, dbf, contrib_mods={e: symengine.S.Zero for e in exclusion})
                 else:
                     output = prop
@@ -318,7 +360,6 @@ def get_thermochemical_data(dbf, comps, phases, datasets, model=None, weight_dic
     return all_data_dicts
 
 
-
 def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfigurationCalculationData, parameters):
     species = calc_data['species']
     phase_name = calc_data['phase_name']
@@ -343,20 +384,16 @@ def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfig
         for sv_key, sv_val in str_statevar_dict.items():
             cond_dict.update({sv_key: sv_val[index]})
 
-        # here we build the internal DOF as a dictionary to use as conditions
+        # Build internal DOF as if they were used in conditions
         dof = {}
         for site_frac in range(len(constituent_list)):
             comp = constituent_list[site_frac]
             occupancy = calc_data['calculate_dict']['points'][index,site_frac]
             sublattice = sublattice_list[site_frac]
             dof.update({v.Y(phase_name,sublattice,comp): occupancy})
 
-        # TODO: convergence seems like it can get messed up if including DOF as conditions.
-        # We temporarily disable that, but be aware that this will allow the minimizer to minimizer internal DOF, so it will only truly work for phases with (# internal DOF) = (# composition conditions)
-        # wks = Workspace(database=dbf, components=species, phases=phase_name, models=models, conditions={**cond_dict, **dof})
-        # TODO: If tests are passing with this, that's bad. We need a phase with more internal DOF (think Laves 2SL)
-        #  -- I believe the test_equilibrium_thermochemcial_error_species fails due to this, but more investigation needed
         # TODO: active_pure_elements should be replaced with wks.components when wks.components no longer includes phase constituent Species
+        # Build composition conditions, probably not necessary given that we don't actually solve anything, but still useful in terms of derivatives probably.
         active_pure_elements = [list(x.constituents.keys()) for x in species]
         active_pure_elements = sorted(set(el.upper() for constituents in active_pure_elements for el in constituents) - {"VA"})
         ind_comps = len(active_pure_elements) - 1
@@ -369,17 +406,16 @@ def compute_fixed_configuration_property_differences(dbf, calc_data: FixedConfig
         # We workaround this by replacing the phase_record_factory models with ours, but this is definitely a hack we'd like to avoid.
         wks = Workspace(database=dbf, components=species, phases=[phase_name], conditions={**cond_dict}, models=models, phase_record_factory=phase_records, parameters=parameters)
         wks.phase_record_factory.models = models
-        # Sometimes in miscibility gaps the solver has trouble converging for
-        # IsolatedPhase if the first compset it finds is at the edge of
-        # composition space. Here, since we know the degrees of freedom, we
-        # initialize the compset for isolated phase correctly in the first
-        # place. This is kind of a hack until I figure out full DOF support.
+        # We then get a composition set and we use a special "NoSolveSolver" to
+        # ensure that we don't change from the data-specified DOF.
         compset = find_first_compset(phase_name, wks)
         new_sitefracs = np.array([sf for _, sf in sorted(dof.items(), key=lambda y: (y[0].phase_name, y[0].sublattice_index, y[0].species.name))])
-        new_sitefracs[new_sitefracs < 1e-14] = 1e-14  # fix zeros, otherwise we get ZeroDivisionErrors in the minimizer.
-        new_statevars = np.array(compset.dof[:len(compset.phase_record.state_variables)])  # not actually changed
-        compset.update(new_sitefracs, compset.NP, new_statevars)
-        results = wks.get(IsolatedPhase(compset, wks=wks)(output))
+        new_statevars = np.array(compset.dof[:len(compset.phase_record.state_variables)])  # no updates expected
+        compset.update(new_sitefracs, 1.0, new_statevars)
+        iso_phase = IsolatedPhase(compset, wks=wks)
+        iso_phase.solver = NoSolveSolver()
+        wks.solver = NoSolveSolver()
+        results = wks.get(iso_phase(output))
         sample_differences = results - sample_values[index]
         differences.append(sample_differences)
     return differences

tests/test_error_functions.py

@@ -100,6 +100,19 @@ def test_fixed_configuration_residual_function(datasets_db):
     assert np.isclose(likelihood, -14.28729, rtol=1e-6)
 
 
+def test_fixed_configuration_residual_with_internal_degrees_of_freedom(datasets_db):
+    """Unstable endmembers in phases that have internal degrees of freedom should retain fixed internal DOF"""
+    dbf = Database(CU_MG_TDB)
+    datasets_db.insert(CU_MG_HM_FORM_LAVES_C15_ANTISITE)
+
+    residual_func = FixedConfigurationPropertyResidual(dbf, datasets_db, phase_models=None, symbols_to_fit=[])
+
+    # Expect that the database has the same energy as the dataset (the former was fit to the latter)
+    residuals, weights = residual_func.get_residuals(np.asarray([]))
+    assert len(residuals) == len(weights)
+    assert np.allclose(residuals, [0.0])
+
+
 def test_get_thermochemical_data_filters_configurations_when_all_configurations_are_invalid(datasets_db):
     datasets_db.insert(CU_MG_HM_MIX_CUMG2_ALL_INVALID)  # No valid configurations
 

tests/testing_data.py

@@ -33,9 +33,9 @@
 $ Generated by brandon (pycalphad 0.5.2.post1)
 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
 
-ELEMENT CU BLANK 0 0 0 !
-ELEMENT MG BLANK 0 0 0 !
-ELEMENT VA BLANK 0 0 0 !
+ELEMENT CU FCC_A1 0 0 0 !
+ELEMENT MG HCP_A3 0 0 0 !
+ELEMENT VA VACUUM 0 0 0 !
 
 FUNCTION GFCCCU 298.15 GHSERCU#; 3200.0 N !
 FUNCTION GFCCMG 298.15 -0.9*T + GHSERMG# + 2600; 3000.0 N !
@@ -342,6 +342,28 @@
 }
 
 
+CU_MG_HM_FORM_LAVES_C15_ANTISITE = yaml.load("""{
+  "components": ["CU", "MG", "VA"],
+  "phases": ["LAVES_C15"],
+  "solver": {
+    "sublattice_site_ratios": [2, 1],
+    "sublattice_configurations": [["MG", "CU"]],
+    "mode": "manual"
+  },
+  "conditions": {
+    "P": 101325,
+    "T": [298.15],
+  },
+
+  "output": "HM_FORM",
+    "values":   [[[34720.0]]],
+  "reference": "FAKE DATA",
+  "comment": "Cu2:Mg is the stable endmember, Mg2:Cu is the antisite endmember that should have high energy."
+}
+""", Loader=YAML_LOADER)
+
+
+
 CU_MG_HM_MIX_CUMG2_ANTISITE = yaml.load("""{
   "components": ["CU", "MG", "VA"],
   "phases": ["CUMG2"],