While working on solver improvements I came across a few phase diagram points that kept failing even when the test suite would pass. These should probably be added to the suite.

db_alzn = Database('alzn_mey.tdb')
my_phases_alzn = ['LIQUID', 'FCC_A1', 'HCP_A3']
equilibrium(db_alzn, ['AL', 'ZN', 'VA'] , my_phases_alzn, {v.X('ZN'):0.04,
                                                       v.T: 400, v.P:101325}, verbose=True)
equilibrium(db_alzn, ['AL', 'ZN', 'VA'] , my_phases_alzn, {v.X('ZN'):0.7+1e-12,
                                                       v.T: 730, v.P:101325}, verbose=True)
equilibrium(db_alzn, ['AL', 'ZN', 'VA'] , my_phases_alzn, {v.X('ZN'):0.48,
                                                       v.T: 560, v.P:101325}, verbose=True)
equilibrium(db_alzn, ['AL', 'ZN', 'VA'] , my_phases_alzn, {v.X('ZN'):0.3,
                                                       v.T: 620, v.P:101325}, verbose=True)
equilibrium(db_alzn, ['AL', 'ZN', 'VA'] , my_phases_alzn, {v.X('ZN'):0.22,
                                                       v.T: 580, v.P:101325}, verbose=True)

Hello Richard, Shana Tova. I tried to calculate entropy of formation of a compound as function of composition. The calculation of entropy is done by numerical differentiation of Gibbs energy of formation. See code below. The result, as can be seen in the attached figure is noisy. Perhaps this can be corrected by calculating the entropy not by numerical differentiation. Is there a way to use symbolic algebra for this?

Thanks, Eli

from pycalphad import Database, equilibrium
import pycalphad.variables as v
from numpy import *
from pylab import *




db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')



def sb(T,phase_name,x_liquid,x_phase,db):
    dT=1.0

    f1m=dG(T-1*dT,phase_name,x_liquid,x_phase,db)
    f0=dG(T+0*dT,phase_name,x_liquid,x_phase,db)

    DG=f0-f1m
    s=DG/(1.0*dT)
    print(T,x_liquid,s)
    return s



def dG(T,phase_name,x_liquid,x_phase,db):
        data = equilibrium(db, ['AL', 'FE'], 'LIQUID', {v.X('FE'): x_liquid, v.T: T, v.P: 1e5},verbose=False)
        muFe=data['MU'].sel( component='FE').data[0].flatten()
        muAl=data['MU'].sel( component='AL').data[0].flatten()
        data = equilibrium(db, ['AL', 'FE','VA'], phase_name, {v.X('FE'): x_phase, v.T: T, v.P: 1e5},verbose=False)
        G=data['GM'].data[0].flatten()
        dG=(G-(1-x_phase)*muAl-(x_phase)*muFe)

        return dG


n=100.0
irange=arange(0,n+1,1)    
xrange=irange/n
xrange[0]=1.0e-4
xrange[-1]=.9999
T=1000*ones(size(irange))
S=zeros(size(irange))

for i in irange: 
    S[i]=sb(T[i],'AL13FE4',xrange[i],4/17.0,db)


plot(xrange,S)
xlabel('x(liquid,Fe)')
ylabel('Delta(S) formation of AL13FE4' )


show(block=False)

The reason is, e.g., GM_MIX, is only subtracting out the contribution from self.models['ref'], which only includes the chemical energy contribution to each end-member. For phases with magnetic, Einstein, or some custom end-member contributions, the result will not be correct.

Standard properties like GM, CPM, etc., appear to be fine.

Checklist

• [x] The documentation examples have been regenerated if the Jupyter notebooks in the examples/ have changed. To regenerate the documentation examples, run jupyter nbconvert --to rst --output-dir=docs/examples examples/*.ipynb from the top level directory)
• [x] If any dependencies have changed, the changes are reflected in the
  • [x] setup.py
  • [x] .travis.yml (deleted)
  • [x] appveyor.yml (deleted)
  • [x] conda_recipe/meta.yaml (deleted)
  • [x] environment-dev.yml (added)

By means of pycalphad, the curved surface of the chemical potentials calculated are discontinuous. As is shown in the figure1, this is the scatter diagram of the chemical potentials of component Ni of the Al-Ni-Cr ternary system at 1273K. In the figure1, some of the chemical potentials change dramatically in different compositions. Selecting a composition corresponding to the dramatically variational chemical potential from figure1, calculate the chemical potentials of different compositions around that composition. As is shown figure2, this is the diagram of the chemical potentials of component Ni of compositions of the area. How to solve the problem that the chemical potentials are discontinuous calculated by pycalphad?

Chemical potentials calculated by pycalphad are dis.docx

setuptools_scm is the PyPA blessed single source versioning solution. We are moving away from versioneer so we don't have to vendor it ourselves and becaue it depends on distutils which is deprecated in Python 3.10.

• Adds setuptools_scm as a runtime dependency so editable local installs can get the version dynamically
• Remove versioneer vendored module and shipped files from MANIFEST.in. We can now properly use setuptools.build_meta as our build system instead of setuptools.build_meta:__legacy__.
• Changes the version scheme to be strict PEP 440. For example, our current version as of a commit on this branch is 0.8.6.dev19+gf3765968. More version scheme details at https://github.com/pypa/setuptools_scm#default-versioning-scheme

Checklist

• [x] The documentation examples have been regenerated if the Jupyter notebooks in the examples/ have changed. To regenerate the documentation examples, run jupyter nbconvert --to rst --output-dir=docs/examples examples/*.ipynb from the top level directory)
• [x] If any dependencies have changed, the changes are reflected in the
  • [x] setup.py
  • [x] environment-dev.yml

I checked the generated sdist and verified that the version in PKG-INFO matches the the one seen by python setup.py --version. The wheels show the version correctly in the CI as well.

After downloading the pycalphad 0.7, using this scripts and B-Mg.tdb database, I got the following errors.

import matplotlib.pyplot as plt
from pycalphad import Database, binplot
import pycalphad.variables as v

db_b_mg = Database('B-Mg.tdb')
my_phases_b_mg = ['GAS', 'LIQUID', 'BETA_RHOMBO_B', 'HCP_A3', 'MGB2', 'MGB4', 'MGB7']
fig = plt.figure(figsize=(9,6))
binplot(db_b_mg, ['B', 'MG', 'VA'] , my_phases_b_mg, {v.X('B'):(0,1,0.05), v.T: (300, 4300, 20), v.P:101325},  ax=fig.gca())
plt.xlim((0,1))
plt.show() 

 runfile('C:/Users/student/Desktop/Pycalphad/pycalphad_binplot_Mg-B.py', wdir='C:/Users/student/Desktop/Pycalphad')
Traceback (most recent call last):

  File "<ipython-input-4-5a2237aca5f4>", line 1, in <module>
    runfile('C:/Users/student/Desktop/Pycalphad/pycalphad_binplot_Mg-B.py', wdir='C:/Users/student/Desktop/Pycalphad')

  File "C:\Users\student\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 880, in runfile
    execfile(filename, namespace)

  File "C:\Users\student\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 102, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/student/Desktop/Pycalphad/pycalphad_binplot_Mg-B.py", line 2, in <module>
    from pycalphad import Database, binplot

  File "C:\Users\student\Anaconda3\lib\site-packages\pycalphad\__init__.py", line 25, in <module>
    import pycalphad.variables as v

AttributeError: module 'pycalphad' has no attribute 'variables'

edited by @bocklund to reformat traceback into code block

Not sure if this is a newbie thing, a problem in my Cython module or in pycalphad: when installing pycalphad via setup.py, I get problems compiling the hyperplane.pyx file:

[1/1] Cythonizing pycalphad/core/hyperplane.pyx

Error compiling Cython file:
------------------------------------------------------------
...
        trial_simplices[i, :] = best_guess_simplex
    cdef double[::1,:,:] trial_matrix = np.empty((num_components, num_components, num_components), order='F')
    cdef double[::1,:] candidate_tieline = np.empty((num_components, num_components), order='F')
    cdef double[::1] candidate_energies = np.empty(num_components)
    cdef double[::1] candidate_potentials = np.empty(num_components)
    cdef bint[::1] bounding_indices = np.ones(num_components, dtype=np.int32)
             ^
------------------------------------------------------------

pycalphad/core/hyperplane.pyx:90:14: Invalid base type for memoryview slice: bint

Error compiling Cython file:
------------------------------------------------------------
...
        bounding_indices[...] = True
        # Should be exactly one candidate simplex
        candidate_simplex = trial_simplices[saved_trial, :]
        for i in range(candidate_simplex.shape[0]):
            idx = candidate_simplex[i]
            candidate_tieline[i, :] = compositions[idx]
                                                 ^
------------------------------------------------------------

pycalphad/core/hyperplane.pyx:127:50: Memoryview 'double[::1]' not conformable to memoryview 'double[:]'.
Traceback (most recent call last):
  File "setup.py", line 27, in <module>
    ext_modules=cythonize(['pycalphad/core/hyperplane.pyx', 'pycalphad/core/eqsolver.pyx']),
  File "/home/antonio_rocha/anaconda3/lib/python3.5/site-packages/Cython/Build/Dependencies.py", line 877, in cythonize
    cythonize_one(*args)
  File "/home/antonio_rocha/anaconda3/lib/python3.5/site-packages/Cython/Build/Dependencies.py", line 997, in cythonize_one
    raise CompileError(None, pyx_file)
Cython.Compiler.Errors.CompileError: pycalphad/core/hyperplane.pyx```

Can someone help me on this?

Hello Richard and Brandon, I still find problems with the equilibrium solver. See below some apparently simple calculations that give wrong results. The calculations are done on the Al-Fe system on the Al-rich side. At several temperatures, the calculation results in a wrong composition for the FCC_A1 phase. Here is the script (after the script I give the results): #!/usr/bin/env python3

-- coding: utf-8 --

""" Created on Sun Jun 18 13:45:49 2017

@author: ebrosh """

from pycalphad import Database, equilibrium import pycalphad.variables as v

db = Database('/home/ebrosh/Documents/ebroshWorks/tcworks/Al-Fe-solidification/alfe_sei.TDB')

print('conditions') conditions={v.X('FE'): .0028156, v.T: 931.385693359375, v.P: 1.0e5} print(conditions) data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium') print(data['Phase'].values.squeeze()) print('compositions of equilibrium phases') print(data['X'].values.squeeze())

print('conditions') conditions={v.X('FE'): .00460485, v.T: 929.933, v.P: 1.0e5} print(conditions) data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium') print(data['Phase'].values.squeeze()) print('compositions of equilibrium phases') print(data['X'].values.squeeze())

print('conditions') conditions={v.X('FE'): .0047534, v.T: 929.933, v.P: 1.0e5} print(conditions) data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium') print(data['Phase'].values.squeeze()) print('compositions of equilibrium phases') print(data['X'].values.squeeze())

print('conditions') conditions={v.X('FE'): .005477, v.T: 929.518, v.P: 1.0e5} print(conditions) data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium') print(data['Phase'].values.squeeze()) print('compositions of equilibrium phases') print(data['X'].values.squeeze())

this last equilibrium is correct

print('conditions') conditions={v.X('FE'): .005477, v.T: 929, v.P: 1.0e5} print(conditions) data = equilibrium(db, ['AL', 'FE','VA'],['LIQUID','FCC_A1'], conditions,verbose=False,calc_opts={'pdens': 2000})

print('phases in equilibrium') print(data['Phase'].values.squeeze()) print('compositions of equilibrium phases') print(data['X'].values.squeeze()) print('note that the last one is correct')

Results of running the script: conditions {P: 100000.0, T: 931.385693359375, X_FE: 0.0028156} phases in equilibrium ['FCC_A1' 'LIQUID'] compositions of equilibrium phases [[ **1.00000000e+00 2.66513395e-11]** [ 9.97001429e-01 2.99857143e-03]] conditions {P: 100000.0, T: 929.933, X_FE: 0.00460485} phases in equilibrium ['FCC_A1' 'LIQUID'] compositions of equilibrium phases [[ 1.00000000e+00 2.75875971e-11] [ 9.95001245e-01 4.99875457e-03]] conditions {P: 100000.0, T: 929.933, X_FE: 0.0047534} phases in equilibrium ['FCC_A1' 'LIQUID'] compositions of equilibrium phases [[ 1.00000000e+00 2.74369001e-11] [ 9.95001245e-01 4.99875456e-03]] conditions {P: 100000.0, T: 929.518, X_FE: 0.005477} phases in equilibrium ['FCC_A1' 'LIQUID'] compositions of equilibrium phases [[ 1.00000000e+00 2.65337850e-11] [ 9.94442085e-01 5.55791522e-03]] conditions {P: 100000.0, T: 929, X_FE: 0.005477} phases in equilibrium ['FCC_A1' 'LIQUID'] compositions of equilibrium phases [[ 9.99845081e-01 1.54919219e-04] [ 9.93613993e-01 6.38600654e-03]] note that the last one is correct

Can this be fixed ?

Best regards, Eli

Dear Richard, Dear Brandon,

while doing pip3 install pycalphad I have some troubles:

    pycalphad/core/eqsolver.cpp:1351:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _obj;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1352:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _grad;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1353:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _hess;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1354:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _internal_cons;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1355:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _internal_jac;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1356:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _internal_cons_hess;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1357:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _multiphase_cons;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1358:14: error: ‘LLVMDoubleVisitor’ in namespace ‘SymEngine’ does not name a type
       SymEngine::LLVMDoubleVisitor _multiphase_jac;
                  ^~~~~~~~~~~~~~~~~
    pycalphad/core/eqsolver.cpp:1359:15: error: ‘LLVMDoubleVisitor’ is not a member of ‘SymEngine’
       std::vector<SymEngine::LLVMDoubleVisitor>  _masses;
                   ^~~~~~~~~
    pycalphad/core/eqsolver.cpp:1359:15: error: ‘LLVMDoubleVisitor’ is not a member of ‘SymEngine’
    pycalphad/core/eqsolver.cpp:1359:43: error: template argument 1 is invalid
       std::vector<SymEngine::LLVMDoubleVisitor>  _masses;
                                               ^
    pycalphad/core/eqsolver.cpp:1359:43: error: template argument 2 is invalid
    pycalphad/core/eqsolver.cpp:1360:15: error: ‘LLVMDoubleVisitor’ is not a member of ‘SymEngine’
       std::vector<SymEngine::LLVMDoubleVisitor>  _massgrads;
                   ^~~~~~~~~
    pycalphad/core/eqsolver.cpp:1360:15: error: ‘LLVMDoubleVisitor’ is not a member of ‘SymEngine’
    pycalphad/core/eqsolver.cpp:1360:43: error: template argument 1 is invalid
       std::vector<SymEngine::LLVMDoubleVisitor>  _massgrads;
                                               ^
    pycalphad/core/eqsolver.cpp:1360:43: error: template argument 2 is invalid
    pycalphad/core/eqsolver.cpp:1361:15: error: ‘LLVMDoubleVisitor’ is not a member of ‘SymEngine’
       std::vector<SymEngine::LLVMDoubleVisitor>  _masshessians;
                   ^~~~~~~~~
    pycalphad/core/eqsolver.cpp:1361:15: error: ‘LLVMDoubleVisitor’ is not a member of ‘SymEngine’
    pycalphad/core/eqsolver.cpp:1361:43: error: template argument 1 is invalid
       std::vector<SymEngine::LLVMDoubleVisitor>  _masshessians;
                                               ^
    pycalphad/core/eqsolver.cpp:1361:43: error: template argument 2 is invalid
    error: command 'x86_64-linux-gnu-gcc' failed with exit status 1
    ----------------------------------------
ERROR: Command "/usr/bin/python3 -u -c 'import setuptools, tokenize;__file__='"'"'/tmp/pip-install-1_nh2rra/pycalphad/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-y_mi9r8j/install-record.txt --single-version-externally-managed --compile" failed with error code 1 in /tmp/pip-install-1_nh2rra/pycalphad/

Python is 3.5. NumPy, SciPy, Cython, Ipopt are present. SymEngine is of version 0.4.0. Is this a known problem?

Opening this PR for feedback.

Current status

I sketched out example implementations of MassFraction and Composition objects (the old Composition object was refactored to MoleFraction).

MassFraction objects are basically sister objects to MoleFraction, just using W instead of X.

Composition objects are convenience objects that would describe a complete composition for a system as would be used in a calculation. The key ideas of a Composition are

1. All Compositions are single point compositions. a. Any mixing, broadcasting, or paths between Composition objects are not Compositions themselves, but simply dicts mapping v.X or v.W to a float, tuple, or array, as in the equilibrium API. b. The Composition object could be a platform to create such dictionaries, but are fundamentally not those things.
2. Can be created from mass or mole fractions and easily convert between them
3. Dependent components can be changed
4. Equality convenience checks between two Compositions
5. Two compositions can be mixed (to create an new point Composition by Composition.mix) and linearly interpolated (to create a dict of {v.X/v.W: [list of values]} by Composition.interpolate_path)

The usual way I imagine a user would use this would be:

from pycalphad import Database, equilibrium, variables as v
dbf = Database('multicomponent.tdb')
comps = ['FE', 'CR', 'NI', 'MO']
phases = list(dbf.phases.keys())
ss316 = v.Composition(dbf, {v.W('CR'): 0.17, v.W('NI'): 0.12, v.W('MO'): 0.025}, 'FE')
conds = {v.P: 101325, v.T: (300, 2000, 20), v.N: 1, **ss316.mole_fractions}
eq = equilibrium(dbf, comps, phases, conds)

Some things that are still open:

1. A v.X/v.W of a non-pure element species will currently break things. Does it make sense to/should species support be implemented?
2. Compositions are currently instantiated by either passing a dictionary of masses or a Database (where the mass dictionary will be created from the element refdata). The main API I imagine that users will want to is Composition(dbf, {v.X('A'): ...}, 'B'). In principle Composition objects should be independent of any Database, besides a user might not want to load a Database (or might not have a database with the elements they want to manipulate, or might use fictitious elements). Should we provide a hardcoded mass dictionary for pure elements as the default and allow databases or custom mass dicts to be input as an optional argument?
3. Should the conditions argument in equilibrium allow for either {statevars + v.X} or {statevars + v.Composition} to be specified? What about {statevars + v.W} (i.e. convert to v.Composition and then to v.X internally)

calling pts = halton(4, 49, scramble=False) produces a zero in pts[48,3] which is wrong

the correct value is 0.00291545

the problem occurs in the mod_matrix calculation when calling np.floor(sum_matrix) for the problematic entry, sum_matrix is 0.9999999999999999 and floor brings it down to 0, while it should result in 1

adding an epsilon at float precision to sum_matrix does the job, but I have not checked whether it can cause a problem at another place

version: pycalphad 0.10.1

• Adds a new Workspace object and associated imperative API. This object now underlies the equilibrium function, which remains for backwards compatibility. A Workspace object can be mutated after creation. The Workspace.get() function accepts ComputableProperty objects as input and returns arrays; the Workspace.plot() function works similarly, but returns a matplotlib figure. Fixes #154
• Adds a new PhaseRecordFactory object and associated connections to PhaseRecord, which enables deferred compilation of symbolically-defined model properties until they are called in the PhaseRecord API. The PhaseRecord object also gains prop and prop_grad functions for computing properties and property gradients of arbitrary symbolic attributes of Model instances.
• Creates a "Computable Property Framework" (CPF) for calculating properties of the system, individual phases, or some combination of these. It also provides a framework for defining custom properties.
• Adds support for "dot" derivatives (constrained total derivatives) of thermodynamic properties (fixes #261 )
• Adds support for T0 (tzero) calculation
• Adds support for driving force calculation, including nucleation driving forces, via the IsolatedPhase and DormantPhase meta-properties. Meta-properties are objects that know how to create new properties (as defined by the CPF). In this case, these meta-properties know how to start sub-calculations . The current way of finding starting points for these calculations is not ideal and prone to failure (basically reuses CompositionSet objects from the workspace equilibrium calculation), but it works well in the typical case.
• Adds physical units support via pint. All ComputableProperty objects have "implementation units" and "display units." Implementation units are what are assumed by the underlying numerical computation (which is done without units, for performance). The display units are what the user gets when Workspace.get or Workspace.plot are called. ComputableProperty.__getitem__ allows the display units to be changed for individual calls to Workspace.get, Workspace.plot, or when setting Workspace.conditions.
• Also regarding physical units, a special unit conversion context is defined to allow conversions between per-mole properties and per-mass properties. There's room to add per-volume here as well, but that may not make it into this PR.
• Associated tests for these new features.

Documentation

• Adds several examples for the Computable Property API, including driving force, T0, and isolated energy surfaces. A demo of export-to-DataFrame is also shown using the new PandasRenderer.
• Other examples have been updated to use the new APIs. Examples using the legacy APIs continue to be supported, but have been moved to a new "Advanced Examples" section. These sections will continue to be expanded and reorganized to accommodate our evolving best practices.

• [x] If any dependencies have changed, the changes are reflected in the
  • [x] setup.py (runtime requirements)

Complex multi-sublattice phases with lots of parameters and using the magnetic ordering model can challenge pycalphad's Just-In-Time (JIT) compiler, especially for computation of second derivatives. In the worst case, the compiler will hang indefinitely and consume RAM until the process is killed.

The reason for this is the increase in algorithmic complexity that comes from having deep piecewise temperature branching in the Model object's representation of the Gibbs energy. However, a common case for the piecewise parameter description is that there is really only one nonzero "branch" for the entire temperature range. While TDB formally wraps every parameter in a piecewise, the Model object is free to discard trivial branches at build time. That is the approach used in the patch for this PR.

This can be done with no loss of important information, and in fact brings pycalphad more in line with how TC will "extrapolate" outside of temperature bounds for parameters. (Note, however, that the behavior is still not matched perfectly here, since pycalphad will still not extrapolate outside temperature bounds in cases where there is more than one nonzero branch.)

This PR includes a test for such a difficult case, where the Model object for the corresponding phase has a Gibbs energy Hessian that cannot be built by the develop JIT compiler. The patch is able to reduce the number of Piecewise nodes in the Gibbs energy's abstract syntax tree by more than half. For the sake of efficiency, instead of a full correctness test we only test that the number of nodes is reduced by half.

In addition, this PR includes a change to the point sampling algorithm in calculate. Currently the sampler (when fixed_grid is True) tries to add additional points between all pairwise combinations of endmembers. For certain multi-component, multi-sublattice phases, there can be thousands of endmembers and, thus, millions of endmember pairs. The proposed change detects this case; when there are more than 100,000 points to be added, the algorithm does not add any. All endmembers are still added, and this change does not affect the random sampling portion of the algorithm.

For certain pathological cases, this will reduce memory consumption by over 90% and resolve classes of memory errors for users attempting multi-component calculations with complex databases.

These are some minimal changes to allow for case sensitivity in DAT compound names (fixes #425). This doesn't address the fundamental issue that DAT and pycalphad think of constituents differently (#419) but it will hopefully enable a greater subset of DAT files to work without issue. In particular, we no longer convert the entire input file to uppercase before parsing. Instead, selective uppercasing is applied where it is "safer" to extract unambiguous names (though it is still not safe in general).

There is a test to make sure the database from the original issue parses, but this still needs a validation test to make sure the parsed database is correct.

When trying to parse in a DAT file which uses Cobalt 'Co' as system component and contains carbon monoxide 'CO' as gas phase, then a ValueError for using duplicate species name will be raised:

db = pycalphad.Database('Bugged_DAT.dat')

Traceback (most recent call last):
  File "C:\Users\timm\anaconda3\envs\factsage-data-extraction-3-8\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "C:\Users\timm\anaconda3\envs\factsage-data-extraction-3-8\lib\site-packages\pycalphad\io\database.py", line 115, in __new__
    return cls.from_file(fname, fmt=fmt)
  File "C:\Users\timm\anaconda3\envs\factsage-data-extraction-3-8\lib\site-packages\pycalphad\io\database.py", line 220, in from_file
    format_registry[fmt.lower()].read(dbf, fd)
  File "C:\Users\timm\anaconda3\envs\factsage-data-extraction-3-8\lib\site-packages\pycalphad\io\cs_dat.py", line 1196, in read_cs_dat
    parsed_phase.insert(dbf, header.pure_elements, header.gibbs_coefficient_idxs, header.excess_coefficient_idxs)
  File "C:\Users\timm\anaconda3\envs\factsage-data-extraction-3-8\lib\site-packages\pycalphad\io\cs_dat.py", line 504, in insert
    raise ValueError(f"A Species named {sp.name} (defined for phase {self.phase_name}) already exists in the database's species ({dbf.species}), but the constituents do not match.")
ValueError: A Species named CO (defined for phase GAS_IDEAL) already exists in the database's species ({Species('CO', 'CO1.0'), Species('N', 'N1.0'), Species('H', 'H1.0'), Species('O', 'O1.0'), Species('NI', 'NI1.0'), Species('FE', 'FE1.0'), Species('CU', 'CU1.0'), Species('B', 'B1.0'), Species('SC', 'SC1.0'), Species('C', 'C1.0'), Species('MN', 'MN1.0'), Species('K', 'K1.0'), Species('CL', 'CL1.0'), Species('S', 'S1.0'), Species('P', 'P1.0'), Species('MG', 'MG1.0'), Species('CA', 'CA1.0'), Species('SI', 'SI1.0')}), but the constituents do not match.

This bug has been found while trying to parse exported DAT files from FactSage 6.2. This has been tested under Python 3.9 and pycalphad 0.10.1. I've attached a zip file containing a minimal DAT file, a script for reproducing the error and a list of packages versions, in case that's needed.

script.zip