error when running equilibrium on CaO-MgO-SiO2 system

[Skyinner]

cao_mgo_sio2.tdb.zip
code:

import matplotlib.pyplot as plt
from pycalphad import equilibrium
from pycalphad import Database, Model
import pycalphad.variables as v
db = Database('cao_mgo_sio2.tdb')
my_phases = list(db.phases.keys())
eq = equilibrium(db, ['CA', 'MG','SI','O', 'VA'], my_phases, {v.X('CA'): (0.03571428571428572),v.X('MG'): (0.03571428571428572),v.X('SI'): (0.28571428571428575) ,v.T: (1000), v.P: 101325})
print(eq)

error message:

ValueError                                Traceback (most recent call last)
<ipython-input-23-65b52dd3bc9c> in <module>
      1 my_phases = list(db.phases.keys())
----> 2 eq = equilibrium(db, ['CA', 'MG','SI','O'], my_phases, {v.X('CA'): (0.03571428571428572),v.X('MG'): (0.03571428571428572),v.X('SI'): (0.28571428571428575) ,v.T: (1000), v.P: 101325})
      3 print(eq)

~/anaconda3/lib/python3.8/site-packages/pycalphad/core/equilibrium.py in equilibrium(dbf, comps, phases, conditions, output, model, verbose, broadcast, calc_opts, to_xarray, scheduler, parameters, solver, callables, phase_records, **kwargs)
    243     output = sorted(output)
    244     if phase_records is None:
--> 245         models = instantiate_models(dbf, comps, active_phases, model=model, parameters=parameters)
    246         phase_records = build_phase_records(dbf, comps, active_phases, conds, models,
    247                                             output='GM', callables=callables,

~/anaconda3/lib/python3.8/site-packages/pycalphad/core/utils.py in instantiate_models(dbf, comps, phases, model, parameters, symbols_only)
    418         mod = models_defaultdict[name]
    419         if isinstance(mod, type):
--> 420             models_dict[name] = mod(dbf, comps, name, parameters=parameters)
    421         else:
    422             models_dict[name] = mod

~/anaconda3/lib/python3.8/site-packages/pycalphad/model.py in __init__(self, dbe, comps, phase_name, parameters)
    200 
    201         self.models = OrderedDict()
--> 202         self.build_phase(dbe)
    203 
    204         for name, value in self.models.items():

~/anaconda3/lib/python3.8/site-packages/pycalphad/model.py in build_phase(self, dbe)
    404         self.models.clear()
    405         for key, value in self.__class__.contributions:
--> 406             self.models[key] = S(getattr(self, value)(dbe))
    407 
    408     def _array_validity(self, constituent_array):

~/anaconda3/lib/python3.8/site-packages/pycalphad/model.py in reference_energy(self, dbe)
    640         phase = dbe.phases[self.phase_name]
    641         param_search = dbe.search
--> 642         pure_energy_term = self.redlich_kister_sum(phase, param_search,
    643                                                    pure_param_query)
    644         return pure_energy_term / self._site_ratio_normalization

~/anaconda3/lib/python3.8/site-packages/pycalphad/model.py in redlich_kister_sum(self, phase, param_search, param_query)
    606                     va_subls = [(v.Species('VA') in phase.constituents[idx]) for idx in range(len(phase.constituents))]
    607                     # The last index that contains a vacancy
--> 608                     va_subl_idx = (len(phase.constituents) - 1) - va_subls[::-1].index(True)
    609                     va_present = any((v.Species('VA') in c) for c in param['constituent_array'])
    610                     if va_present and (max(len(c) for c in param['constituent_array']) == 1):

ValueError: True is not in list

[bocklund]

Thank you, @Skyinner, for opening this issue with all the key details. I made a minor edit to your text to use backticks (`) instead of quotes (") to get the code formatting.

The error is caused by Model class assuming that an ionic liquid has a vacancy in the second sublattice in the code below. Vacancies are usually used to charge balance when there are positive cations in the first sublattice and no charged anions in the second sublattice. In this TDB, that occurs when the SiO2 neutral is used. Having vacancies present would break the pseudo-binary aspect of this TDB by allowing pure metallic liquids for species in the first sublattice to form.

I have seen in the Thermo-Calc example TCEX 17 that vacancy species can be removed from the second sublattice in pseudo-binary diagrams. Maybe making the below code path not add the vacancy term if vacancies are not present would still work.

Whoever works on this issue should carefully compare the energies from pycalphad to the energies from Thermo-Calc and a test should be added for this case.

pycalphad/pycalphad/model.py

Lines 599 to 619 in 62ad209

	# Special normalization rules for parameters apply under this model
	# If there are no anions present in the anion sublattice (only VA and neutral
	# species), then the energy has an additional Q*y(VA) term
	anions_present = any([m.species.charge < 0 for m in mixing_term.free_symbols])
	if not anions_present:
	pair_rule = {}
	# Cation site fractions must always appear with vacancy site fractions
	va_subls = [(v.Species('VA') in phase.constituents[idx]) for idx in range(len(phase.constituents))]
	# The last index that contains a vacancy
	va_subl_idx = (len(phase.constituents) - 1) - va_subls[::-1].index(True)
	va_present = any((v.Species('VA') in c) for c in param['constituent_array'])
	if va_present and (max(len(c) for c in param['constituent_array']) == 1):
	# No need to apply pair rule for VA-containing endmember
	pass
	elif va_subl_idx > -1:
	for sym in mixing_term.free_symbols:
	if sym.species.charge > 0:
	pair_rule[sym] = sym * v.SiteFraction(sym.phase_name, va_subl_idx, v.Species('VA'))
	mixing_term = mixing_term.xreplace(pair_rule)
	# This parameter is normalized differently due to the variable charge valence of vacancies
	mixing_term *= self.site_ratios[va_subl_idx]