Using the basic Module

This chapter explains how to use ChemApp for Python’s basic module. Since this module is a direct port of the ChemApp API into Python, we only provide material to highlight the differences between the ChemApp for Python basic interface, and the original ChemApp API. For further information, please refer to the ChemApp documentation provided by GTT-Technologies.

API Changes in Python

The basic module is intended to be a copy of the original ChemApp API, which is available in languages such as Fortran and C, to Python. Python, however, is quite different in nature to these languages, and using its language features holds many benefits. It is able to help you write less code, make fewer mistakes, and get more done quickly. For these reasons we changed the original ChemApp API to be a bit more Pythonic.

Indexing: 0-based instead of 1-based

ChemApp uses names and indexes to identify system components, phases, phase constituents, sub-lattices, sub-lattice constituents, etc. The names are simply text strings like gas_ideal and SiO2. The indexes are integer numbers.

ChemApp uses the number 1 as index for the first element in a list. We call this one-based or 1-based indexing. Most modern programming environments like C, C++, Java, Go and Python use zero as index of the first element, which is referred to as (you guessed it) zero-based or 0-based indexing.

We decided to use zero-based indexing in ChemApp for Python, since that is what Python uses. This allows us to write less code, and make fewer mistakes. If we remained with ChemApp’s 1-based indexing, you would have to write a loop as follows:

for phase_index in range(1, chemapp.tqnop()+1):
  print(chemapp.tqgnp(phase_index))

The range(1, chemapp.tqnop()+1) bit gives us two opportunities to make a mistake. Now that we have chosen 0-based indexing, the loop can look like this:

for phase_index in range(chemapp.tqnop()):
  print(chemapp.tqgnp(phase_index))

You no longer needed to manage the gap between the Python world and the ChemApp world. ChemApp for Python does that internally.

Error Handling

The following code snippet shows how to process error information with the original ChemApp API in the C programming language:

int main()
{
  LI noerr, nphase;
  int i;

  // initialise chemapp
  tqini(&noerr);

  // invoke an error by trying to get the number of phases before
  // reading a thermochemical data-file
  tqnop(&nphase, &noerr);

  // check for error, and process it
  if (noerr) {
    printf(" ChemApp reports error no. %li\n", noerr);

    // retrieve the error message
    tqerr(TQERRMSG,&noerr);

    // display error message
    for(i=0;i<3;i++)
      printf("%s\n",TQERRMSG[i]);
  }

  return 0;
}

All the tq* functions in the original ChemApp API take an integer called noerr as the last parameter. ChemApp returns a non-zero value into this variable to indicate that an error occurred. This means that every call to ChemApp must be followed by an if statement to check the value of noerr. This is cumbersome, and results in many extra lines of code.

In ChemApp for Python, ChemApp’s error-handling mechanism is hidden. If an error occurs, a chemapp.core.ChemAppError exception is thrown. This means you will always be notified when an error occurs, without having to write extra code. This is what the preceeding block of C code looks like in Python:

from chemapp import basic as chemapp

chemapp.tqini()
nphase = chemapp.tqnop()

Executing this code results in the following output from Python:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "chemapp/basic.pyx", line 906, in chemapp.basic.tqnop
chemapp.core.ChemAppError: A thermodynamic data-file must be read first.

If you want to acquire the ChemApp error code (this can be very informative), it is possible to obtain it like this:

from chemapp import basic as chemapp
from chemapp.core import ChemAppError

chemapp.tqini()
try:
   ca.tqnop()
except ChemAppError as ca_err:
   print(ca_err.errno)

which will result in the following output:

104

You can, of course, also use Python’s exception handling infrastructure to catch and handle the error yourself.

from chemapp import basic as chemapp
from chemapp.core import ChemAppError

try:
    chemapp.tqini()
    nphase = chemapp.tqnop()
except ChemAppError:
    print('A ChemApp error occurred!')

Function Return Values

The original ChemApp API passes results back to the user through variables provided as parameters to the tq* function calls. Here is a C example for the tqnosc function that determines the number of system components in the loaded data file. The variable nscom is used to store the result.

int main()
{
    LI noerr, nscom, i;
    char fname[13], scname[TQSTRLEN];

    tqini(&noerr);  // initialise chemapp.

    strcpy(fname,"cosi.dat");  // specify the file name.
    tqopna(fname,10,&noerr);  // open "cosi.dat" for reading.
    tqrfil(&noerr);  // read the file.
    tqclos(10,&noerr);  // close the file.

    // get the number of system components.
    tqnosc(&nscom, &noerr);

    printf("Number of system component in file %s: %li\n", fname, nscom);

    return 0;
}

ChemApp for Python passes results back by function return values, and not function parameters, which is more intuitive. Here is the equivalent Python code:

from chemapp import basic as chemapp

chemapp.tqini()  # initialise chemapp.

fname = "cosi.dat"  # specify the file name.
chemapp.tqopna(fname, 10)  # open "cosi.dat" for reading.
chemapp.tqrfil()  # read the file.
chemapp.tqclos(10)  # close the file.

# get the number of system components.
nscom = chemapp.tqnosc()

print(f'Number of system components in file {fname}: {nscom}')

For functions that return more than one result, such as tqstpc, which obtains the stoichiometry and molar mass of a phase constituent, ChemApp for Python returns the results as a tuple that can easily be assigned to multiple variables when calling the function. The following example demonstrates this.

from chemapp import basic as chemapp

chemapp.tqini()  # initialise chemapp.

fname = "cosi.dat"  # specify the file name.
chemapp.tqopna(fname, 10)  # open "cosi.dat" for reading.
chemapp.tqrfil()  # read the file.
chemapp.tqclos(10)  # close the file.

# get stoichiometry and molar mass of constituent '2' in phase '1'.
stoich, mm = chemapp.tqstpc(1, 2)

Function Parameter Default Values

Some of the ChemApp tq* functions are very often used with the same set of parameters. It therefore makes sense to assign default values to these parameters to reduce the amount of code that you need to write. C and Fortran do not make provision for this, but Python does.

Some functions in ChemApp for Python basic therefore have default values for their arguments. The default values were chosen as the values that you are most likely to use. Now you can omit the parameter values altogether if they don’t differ from the defaults. This results in less code, and less clutter.

Let’s have a look at the tqconf function as an example. tqconf currently supports only one option E, which is used to for semi-automatic charge balance correction. This is usually done for the entire thermochemical system, but you can also specify the electron system component for a single phase. To do the correction for the entire system, you do the following in C:

tqconf("E", 0, 0, 0, &noerr);

This means that three of the function parameters, valuea, valueb, and valuec are not really used, but you still have to type them. ChemApp for Python assigns a default value of -1 to each parameter. (Remember that ChemApp for Python is zero-based, and ChemApp is one-based.) The ChemApp for Python equivalent of the function call is therefore:

chemapp.tqconf(ConfigurationOption.E)

Other functions that have default parameters values include tqce, tqcel, tqcen, tqcenl, tqgetr, and tqgdat.

Enumerations Instead of Strings

In the original ChemApp API, hard-coded strings are passed as option parameters to many tq* functions. These strings are listed in the user documentation and entered manually into the code. An example of this is seen in the tqgio function that takes an input-output option as first parameter. The possible values are FILE, ERROR, LIST, and LANGUAGE. In the C code snippet below, the option, FILE, is hard-coded into the tqgio function call.

int main()
{
  LI noerr, iovalue;

  tqini(&noerr);
  tqgio("FILE", &iovalue, &noerr);

  return 0;

}

Hardcoding such values can lead to errors, and these may be difficult to identify and fix later. Hard-coding is seen as bad coding practice and should be avoided if possible.

ChemApp for Python uses enumerations to overcome this problem. All the option lists used in ChemApp were turned into Python enumerations. You do not need to look up the option values any more if you are working with an IDE like PyCharm, for example. It will indicate the available options as you type code. Below is the Python equivalent to the C code above.

from chemapp.core import IoOption
from chemapp import basic as chemapp

chemapp.tqini()
iovalue = chemapp.tqgio(IoOption.FILE)

The enumerations used in ChemApp for Python basic, together with their associated functions are presented in the core module chapter.

Values Returned by tqsetc

In the original ChemApp API a call to tqsetc returns a number to identify the equilibrium condition that is set. For system-level conditions it is always the same value, but conditions related to incoming amounts, chemical potentials and relative activities have unique numbers.

In most cases the returned condition number is not needed in the calculation being done, and can therefore be ignored. In ChemApp for Python, you can simply not assigned returned value to a variable. Here is a Fortran example.

! declare variables
     INTEGER NUMCOM, ISN, IPB

! get system component indexes
     CALL TQINSC('Pb ', IPB, NOERR)
     CALL TQINSC('Sn ', ISN, NOERR)

! set the incoming amounts for the Sn and Pb
     CALL TQSETC('IA ', 0, ISN, 0.2D0, NUMCON, NOERR)
     CALL TQSETC('IA ', 0, IPB, 0.8D0, NUMCON, NOERR)

Here is a Python example in which we do want to use the returned condition numbers.

# get system component indexes
ipb = chemapp.tqinsc('Pb')
isn = chemapp.tqinsc('Sn')

# set the incoming amounts for the Sn and Pb
condition_sn = chemapp.tqsetc(ConditionVariable.IA, -1, isn, 0.2)
condition_pb = chemapp.tqsetc(ConditionVariable.IA, -1, ipb, 0.8)

Finally, in the next example we simply ignore the condition number.

# get system component indexes
ipb = chemapp.tqinsc('Pb')
isn = chemapp.tqinsc('Sn')

# set the incoming amounts for the Sn and Pb
chemapp.tqsetc(ConditionVariable.IA, -1, isn, 0.2)
chemapp.tqsetc(ConditionVariable.IA, -1, ipb, 0.8)

This is another example of how ChemApp for Python and Python help you to get more done with less fuss.

Worked Examples

This section presents the same set of worked examples that are included in the ChemApp documentation on the GTT-Technologies website. The examples demonstrate how to perform the same tasks, but with ChemApp for Python basic in Python rather than with ChemApp in C or Fortran.

Print Data File Contents

"""
Print data-file contents.
"""

# Import the chemapp light package.
from chemapp import basic as ca

# Determine the data file path.
from pathlib import Path
data_dir = Path(ca.__file__).parents[0] / 'data'

# Initialise the number of phase constituents which cannot be used as
# incoming species to 0.
nnperm = 0

# Initialise ChemApp for Python.
ca.tqini()

fname = str(data_dir / 'cosi.dat')
# Suggestion: for this program to be flexible, add code that asks for
# the name of the data-file instead.

# Open the thermochemical data-file cosi.dat (system C-O-Si) C for
# reading.
ca.tqopna(fname, 10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

print(f'\nContents of data-file "{fname}"')


# **********************************************************************
# Information on system components.
# **********************************************************************

# Get the number of system components.
nscom = ca.tqnosc()
print(f'Number of system components: {nscom:>2}')
# The :>2 in {nscom:>2} is to right align the text by a width of 2
# characters.

# Print the names of all system components.
print('List of system components:')

for i in range(nscom):
    scname = ca.tqgnsc(i)
    print(f'{i:>5} : {scname}')


# **********************************************************************
# Information on phases.
# **********************************************************************

# Get the number of phases.
nphase = ca.tqnop()
print(f'\nNumber of phases: {nphase:>2}')

# Print the names of all phases, the associated model, and the number of
# phase constituents, if it is a mixture phase.
print('List of phases:')

for i in range(nphase):
    # Get the phase name.
    pname = ca.tqgnp(i)

    # Get the model name.
    model = ca.tqmodl(i).name

    # Get its number of phase constituents.
    npcon = ca.tqnopc(i)

    # If this phase has more than 1 phase constituent, it is a mixture
    # phase...
    if npcon > 1:
        print(f'{i:>5} : {pname:<25} {model:<8} {npcon} constituents')
    # ... otherwise it is a stoichiometric compound.
    else:
        print(f'{i:>5} : {pname:<25} {model:<8}')


# **********************************************************************
# Information on phase constituents of mixture phases.
# **********************************************************************

# Print information on phase constituents.
print('\nList of phase constituents of mixture phases:')

for i in range(nphase):
    # Get its number of phase constituents.
    npcon = ca.tqnopc(i)

    # If this phase has more than 1 phase constituent, it is a mixture
    # phase, and we want to print the list of phase constituents.
    if npcon > 1:
        pname = ca.tqgnp(i)

        print(f'List of phase constituents of phase {pname}')

        # Get and print the name of each phase constituent of the
        # current phase.
        for j in range(npcon):
            pcname = ca.tqgnpc(i, j)

            # Check whether the phase constituent may be used as incoming
            # species. If not, increase the total number of such phase
            # constituents by 1 and print a mark.
            inpcis = ca.tqpcis(i, j)
            nnperm = 0
            if not inpcis:
                nnperm += 1
                print(f'{j:>5} : {pcname:<5} X')
            else:
                print(f'{j:>5} : {pcname:<5}')

    # If there have been phase constituents which cannot be used as
    # incoming species, print a note.
    if nnperm > 0:
        print(f'Note: {nnperm} phase constituent(s) marked with "X"'
              f' is/are not permitted as incoming species.')

Equilibrium Composition of Mixture Phases

"""
Equilibrium composition of mixture phases.
"""

# Import the ChemApp for Python light edition and chemapp's core module.
from chemapp.core import (ConditionVariable, TargetVariable,
                            ResultVariable)
from chemapp import basic as ca

# Determine the data file path.
from pathlib import Path
data_dir = Path(ca.__file__).parents[0] / 'data'


# Initialise ChemApp for Python.
ca.tqini()

# Open the thermochemical data-file cosi.dat (system C-O-Si) for reading.
ca.tqopna(str(data_dir / 'cosi.dat'), 10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

# For the present example, the input composition consists of 2 mol of
# CO/GAS/ and 3 mol of quartz.
ico = ca.tqinpc('CO', 0)
ca.tqsetc(ConditionVariable.IA, 0, ico, 2.0)

isio2 = ca.tqinp('SiO2(quartz)')
ca.tqsetc(ConditionVariable.IA, isio2, 0, 3.0)

# The temperature at which the equilibrium is to be calculated is set to
# 1500 K.
ca.tqsetc(ConditionVariable.T, -1, -1, 1500.0)

# Calculate the equilibrium and print a ChemSage result table.
ca.tqcel(TargetVariable.NoTarget, -1, -1, [0.0, 0.0])


# First, get the equilibrium amount of the gas phase, which is the only
# mixture phase in the system. Note that since we did not change any
# units, the default amount unit (mol) is used.
agas = ca.tqgetr(ResultVariable.A, 0, -1)

print(f'\nEquilibrium amount of gas phase in mol: {agas:.5}.')

# To determine the mole fraction of the gas phase, relative to the whole
# system, we first need to get the amount contained in the system.
asys = ca.tqgetr(ResultVariable.A, -1, -1)

print(f'\nMole fraction of the gas phase: {agas/asys:.5}.')


# The equilibrium amounts of all phase constituents of the gas phase can
# be obtained by a call to tqgetr which passes an array together with the
# proper index/indexp combination.
values = ca.tqgetr(ResultVariable.A, 0, -2)

# Print the equilibrium amount of each phase constituent.
npcon = ca.tqnopc(0)
print(f'\nEquilibrium amounts of gas phase constituents.')

for i in range(npcon):
    name = ca.tqgnpc(0, i)
    print(f'{i:>5} : {name:<25} {values[i]:^.5E} mol')


# The mole fractions of the phase constituents can be calculated using
# the values just retrieved.

# Print the mole fractions of each phase constituent.
print(f'\nMole fractions of gas phase constituents.')
for i in range(npcon):
    name = ca.tqgnpc(0, i)
    print(f'{i:>5} : {name:<25} {values[i]/agas:.5E}')

# Information relating the system components (in the present case C, O,
# and Si) to the phases can also be obtained.

# Get the mole fraction of the system components in the gas phase.
nscom = ca.tqnosc()
print(f'\nMole fractions of system components in the gas phase.')
for i in range(nscom):
    name = ca.tqgnsc(i)
    value = ca.tqgetr(ResultVariable.XP, 0, i)
    print(f'{i:>5} : {name:<25} {value:.5E}')


# To get the fraction of a system component in the gas phase, we have to
# make use of the information contained in the stoichiometry matrix.
isi = ca.tqinsc('Si')
fracsi = 0.0

# Get the equilibrium amount of each gas phase constituent, then use the
# stoichiometry matrix returned by tqstpc to determine how much Si is
# contained in the constituent, and sum it all up.
for i in range(npcon):
    stoi, wmass = ca.tqstpc(0, i)
    value = ca.tqgetr(ResultVariable.A, 0, i)
    fracsi += stoi[isi] * value

# Get the total amount of Si in the system.
sitot = ca.tqgetr(ResultVariable.A, -1, isi)

print(f'\nMole fraction of system component Si in the gas phase: '
      f'{fracsi / sitot:.5E}')


# To demonstrate how ChemApp for Python can be used to check the relationship
# between partial and integral values of extensive properties, a
# different data-file for a system that contains non-ideal mixture
# phases is used.

# Open a new thermochemical data-file (pbsn.dat, system Pb-Sn) for
# reading.
ca.tqopna(str(data_dir / 'pbsn.dat'), 10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

# Find out the system component index numbers for Pb and Sn, as well as
# the phase index number for the liquid phase:
ipb = ca.tqinsc('Pb')
isn = ca.tqinsc('Sn')
iliq = ca.tqinp('LIQUID')

# Using the index numbers just determined, we can enter our incoming
# amounts, in this case an alloy consisting of 0.7 mol Sn and 0.3 mol
# Pb:
xsn = 0.7
ca.tqsetc(ConditionVariable.IA, -1, isn, xsn)
ca.tqsetc(ConditionVariable.IA, -1, ipb, 1.0 - xsn)

# Set the temperature to 600 K. Under these conditions, the liquid phase
# is the only stable phase in this system (see Figure).
ca.tqsetc(ConditionVariable.T, -1, -1, 600.0)

# Calculate the equilibrium.
ca.tqce(TargetVariable.NoTarget, -1, -1, [0.0, 0.0])

# Get the integral enthalpy of the liquid phase.
hint1 = ca.tqgetr(ResultVariable.H, iliq, -1)

# Get the partial molar enthalpies of Sn/LIQUID/ and Pb/LIQUID/.
hsn = ca.tqgetr(ResultVariable.HM, iliq, isn)
hpb = ca.tqgetr(ResultVariable.HM, iliq, ipb)

# Calculate the integral enthalpy of the liquid phase from the two
# partial values.
hint2 = xsn * hsn + (1.0 - xsn) * hpb

# Print both values.
print(f'\nhint1 = {hint1:>10.3f} J')
print(f'hint2 = {hint2:>10.3f} J')


# As expected, the two values are the same.

Phase Equilibria and Phase Target Calculations

"""
Phase equilibria and phase target calculations.
"""

# Note that this program contains phase target calculations,
# which cannot be performed with the 'light' version of ChemApp for Python.

from chemapp.core import(TargetVariable, ConditionVariable,
                           ResultVariable)
from chemapp import basic as ca

# Determine the data file path.
from pathlib import Path
data_dir = Path(ca.__file__).parents[0] / 'data'

# Initialise ChemApp for Python.
ca.tqini()

# Open the thermochemical data-file pbsn.dat (system Pb-Sn) for reading.
ca.tqopna(str(data_dir / 'pbsn.dat'),10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

# To find out about the phases contained in the data-file, print out the
# names of all phases and the associated models.
nphase = ca.tqnop()

for i in range (nphase):
    name = ca.tqgnp(i)
    model = ca.tqmodl(i).name
    print(f'{i:>5} : {name:<25} {model}')

# As can be seen from the short table printed above, all phases present
# in the data-file are mixture phases, since the model identifier
# 'PURE' does not appear anywhere. Also, the phase BCT_A5 is entered
# twice into the data-file. This is required because of the inherent
# metastable miscibility gap in the system.

# **********************************************************************
# Formation phase target calculation.
# Item 1 in the Pb-Sn phase diagram.
# **********************************************************************

# The first calculation will be the determination of the eutectic
# temperature. To do this, we need to define an alloy composition
# which, upon solidification, would undergo a eutectic transformation.
# An alloy consisting of 50% Sn and 50% Pb is a good candidate. We
# will tell ChemApp for Python to use an alloy of this composition, and ask for
# the temperature (which thus is the target variable) at which the
# liquid phase first becomes stable (which thus is the target).
# Applying the 'formation phase target liquid' thus means that ChemApp for Python
# will determine the point where the liquid phase is stable (activity
# is unity) and its equilibrium amount is the one selected by the
# programmer (in our case zero).

# We will use 'global conditions', because we are not interested in
# calculating extensive property (e.g. 'heat') balances.

# We will use the system components to input our incoming amounts. Note
# that it does not matter whether we input our incoming amounts as system
# components, orthrough any of the mixture phase constituents or
# stoichiometric compounds, since we are not interested in the extensive
# property balances.

# Find out the system component index numbers for Pb and Sn:
ipb = ca.tqinsc('Pb')
isn = ca.tqinsc('Sn')

# Using the index numbers just determined, we can enter our incoming
# amounts:

# Entering 0.5 mol Pb.
ca.tqsetc(ConditionVariable.IA, -1, ipb, 0.50)

# Entering 0.5 mol Sn.
ca.tqsetc(ConditionVariable.IA, -1, isn, 0.50)

# Check if we are working with the 'light' version.
# If we do, omit the following phase target calculation(s).
islite =  ca.tqlite()

if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n '
          f'*** since they are not possible with the\n '
          f'*** ''light'' version of ChemApp for Python.')

else:
    # The 'liquid' phase is the target phase. From the first output we
    # see that it has the index number 0, which can also be determined
    # by a call of tqinp:
    i_liq = ca.tqinp('LIQUID')

    # Using tqsetc, we define the _target_: formation phase target
    # 'liquid', zero amount.
    ca.tqsetc(ConditionVariable.A, i_liq, -1, 0.0)

    # The _target variable_ is passed to tqce/tqcel: The target variable
    # used is the temperature, with 300 K as the first guess. This
    # estimate does not need to be very accurate.
    vals = [300.0, 0.0]
    ca.tqce(TargetVariable.T, -1, -1, vals)

    # Retrieve the calculated temperature
    value = ca.tqgetr(ResultVariable.T, -1, -1)

    print(f'\nThe eutectic temperature is {value:<5.2f} K')

# **********************************************************************
# Phase equilibrium calculation (verification of the 'lever rule').
# Item 2 in the Pb-Sn phase diagram.
# **********************************************************************

# Now, for a temperature below the eutectic one, calculate the mole
# fraction of Sn in both eutectic phases (FCC and BCT_A5#1).

# The alloy composition is left unchanged (50% Sn, 50% Pb), but now the
# temperature is set to 400 K.
ca.tqsetc(ConditionVariable.T, -1, -1, 400.0)

# Calculate the equilibrium.
vals = [0.0, 0.0]
ca.tqce(TargetVariable.NoTarget, -1, -1, vals)

# To check the lever rule for this case, and compare it to the graphical
# representation in the Pb-Sn phase diagram, we need the mole fraction
# of Sn in the BCT_A5#1 and FCC phases (fsnbct and fsnfcc), the amounts
# of BCC and FCC phase (abct and afcc), plus the mole fraction of Sn in
# the alloy, which is 0.5.

# Get the equilibrium amount of phase BCT_A5#1 in mol.
i_bct = ca.tqinp('BCT_A5#1')
abct = ca.tqgetr(ResultVariable.A, i_bct, -1)

# Get the mole fraction of Sn in BCT_A5#1.
fsnbct = ca.tqgetr(ResultVariable.XP, i_bct, isn)

# Get the equilibrium amount of phase FCC in mol.
i_fcc = ca.tqinp('FCC')
afcc = ca.tqgetr(ResultVariable.A, i_fcc, -1)

# Get the mole fraction of Sn in FCC.
fsnfcc = ca.tqgetr(ResultVariable.XP, i_fcc, isn)

# For the lever rule condition to be satisfied, the following equation
# has to hold: (0.5 - fsnfcc) * afcc - (fsnbct - 0.5) * abct = eps with
# eps being zero within the numerical precision.
eps = (0.50 - fsnfcc) * afcc - (fsnbct - 0.50) * abct
print(f'\nThe value of eps is {eps:.5E}.')

# As you can see, the lever rule is sufficiently satisfied.

# Check if we are working with the 'light' version. If we do, omit the
# following phase target calculation(s).
islite = ca.tqlite()

if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n '
          f'*** since they are not possible with the\n '
          f'*** ''light'' version of ChemApp for Python.')

else:

    # *****************************************************************
    # Precipitation target calculation, with temperature as target.
    # variable Item 3 in the Pb-Sn phase diagram.
    # *****************************************************************

    # The next calculation is again a target calculation. This time, a
    # 'precipitation phase target' is used to determine a point on the
    # liquidus curve, the alloy composition being the same as before.
    # Using a precipitation phase target tells ChemApp for Python to search for the
    # condition under which a second phase becomes stable.

    # The subsequent call to tqsetc tells ChemApp for Python to set the following
    # _target_: The liquid phase is stable, and a second phase (which
    # does not need to be specified) has unit activity.
    ca.tqsetc(ConditionVariable.A, i_liq, -1, -1.0)

    # Use the temperature as _target variable_ and calculate the
    # equilibrium.
    vals = [1000.0, 0.0]
    ca.tqce(TargetVariable.T, -1, -1, vals)

    # Retrieve the calculated liquidus temperature.
    value = ca.tqgetr(ResultVariable.T, -1, -1)
    print(f'\nThe liquidus temperature is {value:<5.2f} k')

# **********************************************************************
# Precipitation target calculation, with composition as target variable.
# Item 4 in the Pb-Sn phase diagram.
# **********************************************************************

# A target variable which is especially useful in T-x diagrams with
# rather steep phase boundaries is the composition of an alloy. The
# following example determines the point of the FCC phase boundary
# at 500 K where LIQUID is formed (or precipitates).

# First, for convenience, remove _all_ previously set conditions. This is
# especially useful since we have to remove the conditions describing
# the fixed composition of the alloy (set toward the beginning of the
# program) anyway, as the composition is now the target variable and
# thus supposed to vary.
ca.tqremc(-2)

# Set the temperature to 500 K.
ca.tqsetc(ConditionVariable.T, -1, -1, 500.0)

# Set the target: The FCC phase is stable, and a second phase has unit
# activity.
ca.tqsetc(ConditionVariable.A, i_fcc, -1, -1.0)

# If the composition is the target variable and the total input amount
# should always be 1 mol, we need to vary two compositions symmetrically.
# In the present case, we tell ChemApp for Python to vary the input amount of
# Sn/FCC/ between 0 and 1 to find the target, and to vary Pb/FCC/
# symmetrically between 1 and 0.

# The upper and lower limits of the composition target variables are
# passed via the array vals to tqce. First we pass the limits for Sn/
# FCC/.
vals = [0.0, 1.0]

# We need to tell ChemApp for Python that we are not finished with our input, since
# we still need to specify the symmetrical limits for Pb/FCC/ later.
# This is done by calling tqce with the option 'IA0', which tells ChemApp for Python
# not to perform an equilibrium calculation, but to expect more incoming
# amounts.
ca.tqce(TargetVariable.IA0, i_fcc, isn, vals)

# Enter the symmetrical limits for Pb/FCC/.
vals[0] = 1.0 - vals[0]
vals[1] = 1.0 - vals[1]

# With the next call to tqce, the information ChemApp for Python needs to perform a
# composition target is complete, so we call ChemApp for Python with the option
# 'IA', upon which the equilibrium is calculated.

# ChemApp for Python has now varied the composition of the alloy until the desired
# point of the phase boundary has been found. The desired value (the mole
# fraction of Sn in the FCC phase) is retrieved.
fsnfcc = ca.tqgetr(ResultVariable.XP, i_fcc, isn)

print(f'\nThe mole fraction of Sn is {fsnfcc:.5}')

Simple Process Using Streams

"""
A simple process using streams.
"""

# Note that this program contains phase target calculations, which
# cannot be performed with the 'light' version of ChemApp for Python.

from chemapp.core import(
    ResultVariable, ConditionVariable, StreamVariable, Quantity,
    AmountUnit, TemperatureUnit, TargetVariable)
from chemapp import basic as ca

# import packages for plotting graphs.
import matplotlib.pyplot as plt
import pandas as pd

# Determine the data file path.
from pathlib import Path
data_dir = Path(ca.__file__).parents[0] / 'data'

# The following function calculates the sum of all equilibrium amounts
# of the condensed phases. It starts with phase #2, since the gas phase
# is #0. The index number of the last phase to be considered is passed
# as an argument to the function and is retrieved in the main program
# with a call to tqnop.

def getacp(lastph):
    getacp = 0.0
    # Add up the equilibrium amounts of all condensed phases.
    for i in range(lastph):
        amount = ca.tqgetr(ResultVariable.A, i, 0)
        getacp = getacp + amount
    return amount

# Main program.
# Initialise ChemApp for Python.
ca.tqini()

# Open data-file for reading. In this case, a data-file containing the
# thermochemical data for the system carbon-nitrogen-oxygen is selected.
ca.tqopna(str(Path(data_dir) / 'cno.dat'), 10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

# Since we will perform an extensive property target calculation later,
# we will use streams. This way we can assign an initial state regarding
# temperature and pressure to the incoming substances, which will be
# taken into account when ChemApp for Python calculates the extensive property
# (e.g. heat) balance.

# Since we are going to use degrees Celsius as temperature unit
# throughout the calculations, we will change the unit now.
ca.tqcsu(Quantity.Temperature, TemperatureUnit.C)

# According to the process scheme shown above, define the first stream.
# It is the incoming stream of mixed CO2 and CO at 200 C and 1 bar,
# which we will call 'CO2CO'. The amount of CO2 and CO in this stream
# is set to 0.5 and 2.5 mol, respectively.
TP = [200.0, 1.0]
ca.tqsttp('CO2CO', TP[0], TP[1])

# Set the constituent amounts of CO2 and CO. To do this, get the index
# number of these constituents of the gas phase, which is phase number
# 0.
ico2 = ca.tqinpc('CO2', 0)
ico = ca.tqinpc('CO', 0)

# Use the index numbers just obtained to set the constituent amount of
# CO2 in the stream 'CO2CO' to 0.5 mol.
ca.tqstca('CO2CO', 0, ico2, 0.5)

# Similarly, set the amount of CO in the same stream to 2.5 mol.
ca.tqstca('CO2CO', 0, ico, 2.5)

# Define the second stream ('O2N2'), which is air supplied at room
# temperature and 1 bar for the combustion. The air is assumed to
# consist of O2 and N2 at a ratio of 1:4. The amount of O2 is
# initially set to 0.5 mol (thus N2 = 2.0 mol), and will later
# be varied.
TP = [25.0, 1.0]
ca.tqsttp('O2N2', TP[0], TP[1])

# Get the index numbers of O2 and N2.
io2 = ca.tqinpc('O2', 0)
in2 = ca.tqinpc('N2', 0)

# Set the amount of O2.
ca.tqstca('O2N2', 0, io2, 0.5)

# Set the amount of N2.
ca.tqstca('O2N2', 0, in2, 2.0)

# Set the temperature of the reactor at which the equilibrium reaction
# (i.e. the combustion) will take place. It is initially set to 600 C
# at first, but will be varied later.
ca.tqstec(ConditionVariable.T, -1, 600.0)

# In the result table to be printed out subsequently, the amount unit
# for the condensed phases formed should be 'gram', while the input
# was given in 'mol'.
ca.tqcsu(Quantity.Amount, AmountUnit.gram)

# Calculate the equilibrium and print out the result table for a
# quick overview on the equilibrium state.
vals = [0.0, 0.0]
ca.tqcel(TargetVariable.NoTarget, -1, -1, vals)

# As can be seen from the above table, both the gas phase and carbon
# are stable under the conditions chosen. It can also be seen that
# N2, CO2, and CO are the only majority gas phase species. All
# other gas phase species are produced in negligible amounts only.

# The values for the amounts of condensed phases and the gas phase
# can be retrieved using tqgetr. Since the gas phase, if present in
# the data-file, has always the phase index number 0, this means
# that we have to add up all the amounts of the phases from phase
# #2 to the last one in order to retrieve the total amount of
# condensed phases.

# Find out how many phases we have in the data-file, and use this
# number in all subsequent calls to getacp.
nphase = ca.tqnop()

# Get the total equilibrium amount of the condensed phases, using
# the function getacp declared above.
acond = getacp(nphase)

print(f'Amount of condensed material: {acond:.4f} gram.')

# Retrieve the amount of gas phase. Note: The current amount unit
# is still 'gram'. If the amount of the gas phase is desired in
# 'mol' instead, tqcsu would need to be called again.
agas = ca.tqgetr(ResultVariable.A, -1, -1)

print(f'Amount of gas phase: {agas:.2f} gram.')

# Retrieve values for some extensive properties.  Since we used
# streams for this example, these extensive property values are
# calculated relative to the incoming streams, when the whole
# system is selected.

# Get the enthalpy change associated with this reaction.
result = ca.tqgetr(ResultVariable.H, -1, -1)
print(f'Delta H: {result:.4E} J')

# The negative sign of this value indicates that energy is released
# by the system.

# Get the change in volume of the process, which is the volume of
# the equilibrium gas phase minus the sum of the volumes of the
# incoming streams.
result =  ca.tqgetr(ResultVariable.V, -1, -1)
print(f'Delta V: {result:.2f} dm^3')

# Note: using the option 'V' with tqgetr normally would retrieve
# the _total_ volume of the system. There is one exception, in
# which case it retrieves the _change_ in volume, and that is the
# case when streams are used _and_ the whole system (i.e. tqgetr
# is called using '-1,-1', as above) is selected. To retrieve the
# _total_ volume instead in the present case, use the option 'VT':
result = ca.tqgetr(ResultVariable.VT, -1, -1)
print(f'Total V: {result:.2f} dm^3')

# Compare the two values for the volume with the corresponding
# ones in the output table produced earlier. The volume printed
# at the top of such an output table is always the _total_ volume,
# whereas the line specifying the extensive property balance toward
# the end of the table shows the _change_ in volume.

# tqstxp can be used to retrieve extensive property values for the
# incoming streams. Get the volume of the incoming CO2/CO stream:
result = ca.tqstxp('CO2CO', StreamVariable.V)
print(f'\nV(CO2CO): {result:.2f} dm^3')

# Get the value of the incoming air stream:
result = ca.tqstxp('O2N2', StreamVariable.V)
print(f'V(O2N2): {result:.3f} dm^3')

# For users of ChemApp versions prior to 3.2.0: Please note that
# the values returned by tqstxp always er to the _last_ calculated
# equilibrium state. Even if one uses tqstxp to retrieve properties
# of _incoming_ streams, one thus needs to perform an equilibrium
# calculation before the correct values are obtained. If tqstxp
# is called before the first equilibrium calculation has been
# performed using tqce or tqcel, VAL will be zero.

# As of version 3.2.0 of ChemApp, tqstxp can also be called
# _before_ the equilibrium calculation.

# As can be seen from the above results, this process, under the
# conditions we have chosen, results in solids to be formed in the
# reactor. We will now further investigate this systems by looking
# at how the amount of condensed material changes when
# a) the reaction temperature is varied, and
# b) the amount of air supplied is varied.

# *********************************************************************
# Varying the reaction temperature.
# *********************************************************************

# Since several individual equilibria will be calculated, the results
# (temperature and amount of condensed phases) will be stored to a
# file for further processing.
file_11 = open('caf103_1.rst', 'w')

# Vary the temperature between 200 C and 1000 C in steps of 5 C.
for j in range(200, 1001, 5):

    # Use the variable temp to hold the current temperature.
    temp = j
    ca.tqstec(ConditionVariable.T, -1, temp)

    # Calculate the equilibrium. In this case, we only need the value
    # for the amount of condensed phases formed, therefore tqce is
    # used and not tqcel.
    ca.tqce(TargetVariable.NoTarget, -1, -1, [0.0, 0.0])

    # Get the amount of condensed phases.
    acond = getacp(nphase)

    # Write the temperature and the amount of condensed phases to the
    # result file.
    file_11.write(f'{temp},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure.

# Read the created caf103_1.rst file, and name the first column
# 'temp' and name the second 'acond'.
data_1 = pd.read_csv('caf103_1.rst', names = ['temp', 'acond'])
plt.axis([200, 1000, 0, 10])
plt.ylabel('Amount/gram')
plt.xlabel('Temperature/C')
plt.plot(data_1['temp'], data_1['acond'], 'r-')
plt.show()

# *********************************************************************
# Varying the amount of air supplied.
# *********************************************************************

# Instead of the temperature, the amount of incoming air will be varied
# now. The temperature will be kept constant at the initial value of
# 600 C. To plot the results later on, we will again store them in a
# separate file.
file_11 = open('caf103_2.rst', 'w')

# Set the temperature for the equilibrium calculation (600 C).
ca.tqstec(ConditionVariable.T, -1, 600.0)

# Calculate 201 equilibria...
for j in range(0, 201):

    # ... with the O2 amount varying between 0 and 1 mol...
    iao2 = j * 0.005
    # ... and the N2 amount between 0 and 4 mol.
    ian2 = iao2 * 4.0

    # To make sure ChemApp for Python interprets these values as 'mol', set the
    # amount unit accordingly.
    ca.tqcsu(Quantity.Amount, AmountUnit.mol)

    # Set the amount of O2 in the air stream.
    ca.tqstca('O2N2', 0, io2, iao2)

    # Set the amount of N2 in the air stream.
    ca.tqstca('O2N2', 0, in2, ian2)

    # Calculate the equilibrium.
    vals = [0.0, 0.0]
    ca.tqce(TargetVariable.NoTarget, -1, -1, vals)

    # Set the amount unit to 'gram'.
    ca.tqcsu(Quantity.Amount, AmountUnit.gram)

    # Get the amount of condensed phases.
    acond = getacp(nphase)

    # Write the amounts of O2 and condensed phases to the result file.
    file_11.write(f'{iao2},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure:
data_2 = pd.read_csv('caf103_2.rst', names = ['iao2', 'acond'])
plt.axis([0, 1, 0, 12])
plt.ylabel('Amount/gram')
plt.xlabel('n(O2)/mol')
plt.plot(data_2['iao2'], data_2['acond'], 'r-')
plt.show()

# *********************************************************************
# Varying the amount of air supplied _and_ the temperature.
# *********************************************************************

# To investigate the combined influence of temperature and amount of air
# supplied on the resulting amount of condensed phases, both parameters
# are varied in nested loops. The results are again stored in a file for
# further processing.
file_11 = open('caf103_3.rst', 'w')

# The outer loop will vary the temperature between 200 C and 1000 C.
for j in range(200, 1001, 20):

    temp = j
    ca.tqstec(ConditionVariable.T , -1, temp)

    # The inner loop varies the amount of incoming air as above.
    for k in range(0, 201, 5):
        iao2 = k * 0.005
        ian2 = iao2 * 4.0

        ca.tqcsu(Quantity.Amount, AmountUnit.mol)
        # Set the incoming amount of O2.
        ca.tqstca('O2N2', 0, io2, iao2)

        # Set the incoming amount of N2.
        ca.tqstca('O2N2', 0, in2, ian2)
        ca.tqce(TargetVariable.NoTarget, -1, -1, [0.0, 0.0])
        ca.tqcsu(Quantity.Amount, AmountUnit.gram)

        # Get the amount of condensed phases
        acond = getacp(nphase)

        # Write the temperature and the amounts of O2 and condensed
        # phases to the result file.
        file_11.write(f'{temp},{iao2},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure:

data_3 = pd.read_csv('caf103_3.rst', names = ['temp', 'iao2', 'acond'])
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')
ax.set_xlim3d(200, 1000)
ax.set_ylim3d(0, 1)
ax.set_zlim3d(-10, 20)
ax.set_xlabel('Temperature/C')
ax.set_ylabel('n(O2)/mol')
ax.set_zlabel('Amount/gram')

x = data_3['temp']
y = data_3['iao2']
z = data_3['acond']
ax.scatter(x, y, z, color = 'red')
plt.show()

# Check if we are working with the 'light' version. If we do, omit
# the following phase target calculation(s).
islite = ca.tqlite()
if islite == True:
    print(f'*** Phase target calculations have been omitted here, "\n"'
          f'*** since they are not possible with the "\n"'
          f'*** ''light''version of ChemApp for Python.')

else:

    # *****************************************************************
    # Varying the amount of air supplied _and_ the temperature,
    # searching for the explicit condition where the amount of condensed
    # phase is zero.
    # *****************************************************************

    # Assuming that one is interested in the above process _not_
    # producing any condensed phases (i.e. complete combustion), it would
    # be helpful to have a diagram which corresponds to the bottom plane
    # of the above figure: A graph which, as a function of the incoming
    # amount of air, shows the temperature at which the amount of
    # condensed phases formed is zero.

    # For the calculation procedure with ChemApp for Python, this results in the
    # following steps:
    # 1) We are looking for the temperature at which condensed phases
    # _just_ start to form. In terms of ChemApp for Python, we thus have to define
    # a precipitation target with respect to the gas phase.
    # 2) For each value of the amount of air supplied, define the 'O2N2'
    # stream with the proper amounts of O2 and N2.
    # 3) ChemApp for Python is then supposed to vary the temperature until the
    # target (a phase precipitates from the gas phase) is met.
    # Temperature is thus the target variable. This is specified as a
    # parameter to tqce, together with an estimated value.
    # 4) The equilibrium is calculated for each of the amounts of air
    # supplied, then the equilibrium temperature is retrieved and stored
    # in the result file.
    file_11 = open('caf103_4.rst', 'w')

    # Reset the amount unit to 'mol'.
    ca.tqcsu(Quantity.Amount, AmountUnit.mol)

    # Define the target as a precipitation target ('-1') for the gas
    # phase (phase index 0).
    ca.tqstec(ConditionVariable.A, 0, -1.0)

    # The initial estimate for the temperature is passed via the first
    # element of vals to subsequent calls to tqce.
    vals = [600.0, 0.0]

    # Vary the amount of supplied air.
    for j in range(0, 200):
        iao2 = j * 0.005
        ian2 = iao2 * 4.0

        # Set the incoming amount of O2.
        ca.tqstca('O2N2', 0, io2, iao2)

        # Set the incoming amount of N2.
        ca.tqstca('O2N2 ', 0, in2, ian2)

        # Define the target variable (temperature), with the initial
        # estimate in the first element of vals, and calculate the
        # equilibrium.
        ca.tqce(TargetVariable.T, -1, -1, vals)

        # Retrieve the calculated equilibrium temperature.
        temp =  ca.tqgetr(ResultVariable.T, -1, -1)

        # Store the amount of O2 supplied and the calculated
        # equilibrium temperature to the result file.
        file_11.write(f'{iao2},{temp}\n')

    # Close data-file.
    file_11.close()

    # When plotted, the resulting graph looks like the following Figure.
    data_4 = pd.read_csv('caf103_4.rst', names = ['iao2', 'temp'])
    plt.axis([0, 1, 500, 800])
    plt.ylabel('Temperature/C')
    plt.xlabel('n(O2)/mol')
    plt.plot(data_4['iao2'], data_4['temp'], 'r-')
    plt.show()

One-dimensional Phase Mapping

"""
One-dimensional phase mapping calculations.

Note that this program contains phase target calculations that cannot
be performed with the 'light' version of ChemApp for Python.

The function bsort implements a simple bubble sort algorithm. It is used
to sort the array of temperature values for which phase equilibria are to
be calculated.
"""
import sys
from chemapp.core import (Status,
                        ConditionVariable, TargetVariable,
                        ResultVariable, PhaseMapVariable)
from chemapp import basic as ca

# Determine the data file path.
from pathlib import Path
data_dir = Path(ca.__file__).parents[0] / 'data'

# import packages for plotting graphs.
import matplotlib.pyplot as plt
import pandas as pd

# Sort the first nused elements of tarray, which has a total number of
# npts elements, in ascending order.
def bsort(npts, tarray, nused):
    for i in range(nused):
        for j in range(nused):
            if tarray[i] > tarray[j]:
                temp = tarray[i]
                tarray[i] = tarray[j]
                tarray[j] = temp
                result = tarray
    return result


# The function psp prints the value of the 'search variable'
# (temperature, pressure, or composition), the names of the
# stable phases, and their amounts. It is used after every
# call to tqmap. If the second parameter to psp is zero, we
# only need to print the first one.
def psp(value1, value2):
    eps = 1E-8

    if value2 != 0.0:

        print(f'{value1:>10.3f}{value2:>10.3f}')
    else:
        print(f'{value1:>10.2f}')

    # Find out how many phases we have in the system.
    nphase = ca.tqnop()
    # For every phase, check whether its activity is unity. Note that for
    # numerical reasons, ChemApp for Python will sometimes return an activity of
    # less than unity for a stable phase. To still be able to determine
    # whether a phase is stable, one should be on the safe side if the
    # phase's activity is compared against 1.0 - 1E-8, as in the
    # if-clause below. If the activity is larger than this value, the
    # phase can be considered to be stable.
    for i in range(nphase):
        act = ca.tqgetr(ResultVariable.AC, i, -1)

        if act >= 1.0 - eps:

            # Check if the phase has the status ENTERED, otherwise (i.e.
            # when the status of the phase is DORMANT) it might have an
            # activity > 1 and still have an equilibrium amount of zero:
            pstat = ca.tqgsp(i)
            if pstat == Status.ENTERED:

                # If it is stable _and_ its status is ENTERED, get the
                # name of that phase...
                pname = ca.tqgnp(i)

                # ...and its equilibrium amount...
                amt = ca.tqgetr(ResultVariable.A, i, -1)

                # ... and print both:
                print(f'{pname:>35} {amt:>13.2f}')


asize = 100
eps = 1E-3
nused = 0

# Initialise ChemApp for Python.
ca.tqini()

# Open the thermochemical data-file pbsn.dat (system Pb-Sn) for reading.
ca.tqopna(str(data_dir / 'pbsn.dat'), 10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

# Check if we are working with the 'light' version.
# If we do, omit the following phase target calculation(s).
islite = ca.tqlite()
if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n '
          f'*** since they are not possible with the\n '
          f'*** ''light'' version of ChemApp for Python.')
    sys.exit(0)

# The data-file for system Pb-Sn contains a second copy of the
# BCT_A5 phase, because of the inherent metastable miscibility gap.
# Since this miscibility gap will stay metastable under the
# conditions used in this example, we will eliminate the second
# copy of that phase, which means it will be disregarded
# completely in the following calculations.
bct2 = ca.tqinp('BCT_A5#2')
ca.tqcsp(bct2, Status.ELIMINATED)

# Determine the index numbers for the system components Pb and Sn.
ipb = ca.tqinsc('Pb')
isn = ca.tqinsc('Sn')

# As indicated in Figure (The Pb-Sn phase diagram), label "1",
# the first phase map should be done with the temperature as the
# search variable, for X(Sn) = 0.2.

# Set the incoming amounts for the Sn and Pb.
ca.tqsetc(ConditionVariable.IA, -1, isn, 0.2)
ca.tqsetc(ConditionVariable.IA, -1, ipb, 0.8)

# vals is passed to tqmap; this list of 2 float variables contains
# the lower and upper values of the search variable. Note that the
# first two calls to tqmap _always_ calculate the equilibria at
# vals[0] and vals[1], and thus do _not_ constitute phase boundaries.
vals = [400.0, 600.0]
tarray = []

# The purpose of the array tarray is to hold the temperature values
# at which the equilibria will be calculated later. This set of
# temperature values consists of,
# 1) a number of evenly spaced temperature intervals,
# 2) and the temperatures at which tqmap finds phase boundaries.

# Part of the array is first filled with evenly spaced temperature
# intervals. As upper and lower values we use the same values as we
# did for the initialisation of vals[0] and vals[1] above. The
# variable nused keeps track of how many items of the array are
# actually used:

for itemp in range(400, 601, 10):
    tarray.append(itemp)
    tarray[nused] = itemp
    nused += 1

# resno keeps track of the number of results we get from tqmap. This
# allows us to ignore the first two temperatures tqmap reports. As
# indicated above, the first two calls to tqmap do _not_ give us
# information about a phase boundary, but on the state at the lower
# and upper end of the search interval (i.e. at vals[0] and vals[1]).
# Since we have these temperatures in our temperature array tarray
# already (see the for loop above), we only need to pay attention
# when we retrieve the third and any subsequent result with tqmap.
resno = 0

# First call to tqmap. Note the option, PhaseMapVariable.TF, where
# 'T' indicates that temperature is the search variable, and 'F'
# ('first') that it is the first call to tqmap for the current
# phase map. The variable returned through icont with every call to
# tqmap will allow us to decide whether more calls to tqmap are
# necessary. Since this is only the first call to tqmap, we know
# anyway that we have to continue.
icont = ca.tqmap(PhaseMapVariable.TF, -1, -1, vals)
resno += 1

# Get the temperature. Since it is the first call to tqmap, we know
# that this will be the value of vals[0].
value1 = ca.tqgetr(ResultVariable.T, -1, -1)

# Print a header for the table to come.
print(f"\n{'Temperature/K':<26}{'Stable phases':<18}{'Amount/mol'}")

# Call psp (print stable phases), passing the value for the
# temperature we just retrieved. psp will print a formatted
# entry for the result table.
psp(value1, 0.0)

# tqmap is called again. Note the option, PhaseMapVariable.TN, where
# 'T' indicates that temperature is the search variable, and 'N'
# ('next') that it is _not_ the first call to tqmap for the current
# phase map.
icont = ca.tqmap(PhaseMapVariable.TN, -1, -1, vals)
resno += 1

# Get the temperature.
value1 = ca.tqgetr(ResultVariable.T, -1, -1)

psp(value1, 0)

while icont:
    # The variable icont returned by tqmap indicates whether we need
    # to make any more calls to tqmap, i.e. whether there are any
    # more phase boundaries within the search interval we specified.
    # If icont is positive, we need to continue.
    icont = ca.tqmap(PhaseMapVariable.TN, 0, 0, vals)
    value1 = ca.tqgetr(ResultVariable.T, -1, -1)
    resno += 1

    # If this is the second call to tqmap, we know that the
    # temperature retrieved here is vals[1], and also that tarray
    # already contains it. If, on the other hand, we are already
    # retrieving our third or any subsequent result (checked
    # through resno), we know it is a phase boundary.
    if resno > 2.0:
        # What follows here is necessary because we would like to
        # see the results of this one-dimensional phase map in a
        # graph later.

        # If we were just interested in the _position_ of this phase
        # boundary, and thus the _value_ of the search variable (the
        # temperature in this case), it would be enough to record
        # this piece of information. However, since we would like to
        # see a graphical representation of the phase map, we might
        # not get the desired results if we just add the temperature
        # found to tarray. Instead of calculating the equilibrium at
        # the temperature just retrieved, we better calculate it at
        # _two_ temperatures: One slightly below, and one slightly
        # above the temperature found. Thus these two temperatures
        # are added to tarray instead of the exact temperature for
        # the phase boundary (eps has been defined using a float at
        # the beginning of the program):
        nused += 1
        tarray.append(value1 - eps)
        nused += 1
        tarray.append(value1 + eps)

        # Print the entry for the result table. For this purpose, we
        # need the _exact_ temperature at the phase boundary, as
        # determined by tqmap:
        psp(value1, 0.0)

# When we arrive at this point of the program icont is zero, which
# tells us that there are no more phase boundaries left within the
# search interval we specified.

# The array tarray now contains all the temperatures at which we
# want to calculate the equilibria for our graphical phase map:
# the evenly spaced temperatures plus the two temperatures for
# every phase boundary: One slightly below and one slightly about
# the temperature of the original phase boundary.

# The only problem is that these temperatures are not in ascending
# order, since all the temperatures for the phase boundaries have
# simply been added to the end of the list of evenly spaced
# temperatures.

# This is revealed by looking at the first and last 10 values of
# tarray:
for itemp in range(nused):
    if itemp < 10 or itemp > nused -10:
        print(f'TARRAY({itemp:>2})= {tarray[itemp]:>8.5f}')

# To have the values for the graph displayed properly later, we thus
# need to sort this array using the function bsort:
bsort(asize, tarray, nused)

# We are now ready to calculate the equilibria for all the
# temperatures in tarray. These equilibrium calculations now
# being performed are _not_ target calculations of any kind
# any more, that's why we reset _all_ previous conditions
# including all targets and target variables.
ca.tqremc(-2)

# The following equilibria are to be calculated for the same alloy
# concentration of course, so we use the same incoming amounts again.
ca.tqsetc(ConditionVariable.IA, 0, isn, 0.2)
ca.tqsetc(ConditionVariable.IA, 0, ipb, 0.8)

# The graphical representation of the one-dimensional phase map will
# be done using two x/y graphs: one displaying the phase amounts vs.
# the temperature, the other one the phase activities vs. the
# temperature. We will write two output data-files for this purpose:
# caf104_1.rst will contain the phase amounts, whereas caf104_2.rst
# will contain the phase activities.
file_11 = open('caf104_1.rst', 'w')

file_12 = open('caf104_2.rst', 'w')

# For every temperature in tarray, calculate the equilibrium:
for itemp in range(nused):
    ca.tqsetc(ConditionVariable.T, 0, 0, tarray[itemp])
    # Note again that these are regular equilibrium calculations, no
    # targets are involved:
    ca.tqce(TargetVariable.NoTarget, 0, 0, vals)

    # Retrieve the equilibrium amounts of _all_ phases with a single
    # call to tqgetr. Note that amts is an array of appropriate size.
    amts_liquid = ca.tqgetr(ResultVariable.A, 0, -1)
    amts_bct2 = ca.tqgetr(ResultVariable.A, 1, -1)
    amts_fcca1 = ca.tqgetr(ResultVariable.A, 3, -1)

    # Write the temperature plus the phase amounts to the first
    # output data-file:
    file_11.write(f'{tarray[itemp]},{amts_liquid},{amts_bct2},'
                  f'{amts_fcca1}\n')
    # Retrieve the equilibrium activities of _all_ phases with a
    # single call to tqgetr. Note that scts is also an array of
    # appropriate size.
    acts_liquid = ca.tqgetr(ResultVariable.AC, 0, -1)
    acts_bct2 = ca.tqgetr(ResultVariable.AC, 1, -1)
    acts_fcca1 = ca.tqgetr(ResultVariable.AC, 3, -1)

    # Write the temperature plus the phase activities to the second
    # output data-file:
    file_12.write(f'{tarray[itemp]},{acts_liquid},{acts_bct2},'
                  f'{acts_fcca1}\n')
    # We have now calculated the equilibria for all temperature
    # values in tarray and can close the two output data-files.
file_11.close()
file_12.close()

data_1 = pd.read_csv('caf104_1.rst', names=['tarray', 'amts_liquid',
                     'amts_bct2', 'amts_fcca1'])
plt.axis([400, 600, 0, 1])
plt.ylabel('Amount/mol')
plt.xlabel('Temperature/K')
plt.plot(data_1['tarray'], data_1['amts_liquid'], 'r-+',
         linewidth=0.5, markersize=5, label='LIQUID')
plt.plot(data_1['tarray'], data_1['amts_bct2'], 'g-x',
         linewidth=0.5, markersize=5, label='BCT_A5#1')
plt.plot(data_1['tarray'], data_1['amts_fcca1'], 'm-s',
         linewidth=0.5, markersize=5, markerfacecolor='none',
         label='FCC_A1')
plt.legend(loc='best')
plt.show()

data_2 = pd.read_csv('caf104_2.rst', names=['tarray', 'acts_liquid',
                     'acts_bct2', 'acts_fcca1'])
plt.axis([400, 600, 0.55, 1])
plt.ylabel('Activity')
plt.xlabel('Temperature/K')
plt.plot(data_2['tarray'], data_2['acts_liquid'], 'r-+',
         linewidth = 0.5, markersize = 5, label = 'LIQUID')
plt.plot(data_2['tarray'], data_2['acts_bct2'], 'g-x',
         linewidth = 0.5, markersize = 5, label = 'BCT_A5#1')
plt.plot(data_2['tarray'], data_2['acts_fcca1'], 'm-s',
         linewidth = 0.5, markersize = 5, markerfacecolor='none',
         label = 'FCC_A1')
plt.legend(loc="best")
plt.show()

# When plotted, the two graphs, displaying the equilibrium phase
# amounts and phase activities vs. temperature look like Figures
# 3 and 4.

# For the next one-dimensional phase mapping calculations we will
# only use tabular output. First, remove all previously set
# conditions.
ca.tqremc(-2)

# We will now calculate another one-dimensional phase map in the same
# system, again with temperature as the search variable. The new
# alloy composition will cause tqmap to find the eutectic.
ca.tqsetc(ConditionVariable.IA, 0, isn, 0.5)
ca.tqsetc(ConditionVariable.IA, 0, ipb, 0.5)

# The temperature search interval used is the same:
vals = [400.0, 600.0]

# First call to tqmap (note again the 'F' in the option parameter)
icont = ca.tqmap(PhaseMapVariable.TF, 0, 0, vals)

# Retrieve the temperature, which we already know is value1.
value1 = ca.tqgetr(ResultVariable.T, -1, -1)

# Print a header for the table to come.
print(f'\n{"Temperature/K":<26}{"Stable phases":<18}{"Amount/mol"}')

# Call psp to print the information on the stable phases:
psp(value1, 0)

# tqmap is called with 'N' in the option parameter:
icont = ca.tqmap(PhaseMapVariable.TN, 0, 0, vals)

# Get the temperature...
value1 = ca.tqgetr(ResultVariable.T, -1, -1)

# ... and print the entry for the result table:
psp(value1, 0)
while icont:
    if icont > 0:
        icont = ca.tqmap(PhaseMapVariable.TN, 0, 0, vals)
        value1 = ca.tqgetr(ResultVariable.T, -1, -1)
        psp(value1, 0)

# Note that the above table lists a 3-phase equilibria now, which
# for the current system means we found the eutectic temperature
# (compare this with the results for the same equilibrium determined
# in worked example 3: Phase equilibria and phase target
# calculations).


# So far we have calculated one-dimensional phase maps for a
# temperature/composition phase diagram, in our case keeping the
# composition (and implicitly the pressure) constant and using the
# temperature as the search variable. The next phase map we will
# calculate at a constant temperature, this time varying the
# pressure. For this purpose, we will use a different thermochemical
# data-file for the system Fe2SiO4-Mg2SiO4 which includes
# pressure-dependent Gibbs-energy data.

# Open the thermochemical data-file femgsio4.dat (system Fe2SiO4-
# Mg2SiO4) for reading.
ca.tqopna(str(data_dir / 'femgsio4.dat'), 10)

# Read data-file.
ca.tqrfil()

# Close data-file.
ca.tqclos(10)

# Determine the index numbers for the stoichiometric compounds
# Fe2SiO4 and Mg2SiO4. These two phases will be used to define
# the incoming amounts for the next calculations. Note that they
# are included in the thermochemical data-file only for this
# purpose (to allow for a somewhat easier definition of incoming
# amounts), since the Gibbs energy functions for these two
# stoichiometric compounds are only dummies.
ipfesi = ca.tqinp('Fe2SiO4')
ipmgsi = ca.tqinp('Mg2SiO4')

# Select a composition of 2 mol-% Fe2SiO4, the rest is Mg2SiO4:
ca.tqsetc(ConditionVariable.IA, ipfesi, 0, 0.02)
ca.tqsetc(ConditionVariable.IA, ipmgsi, 0, 0.98)

# Select the range for the search variable, which is the pressure in
# the current case. We will search for any phase boundaries between
# 120 kbar and 150 kbar (see also the line marked '1' in Figure 2).
vals = [120000.0, 150000.0]

# The temperature is set to 1000 K.
ca.tqsetc(ConditionVariable.T, 0, 0, 1000.0)

# First call to tqmap. Note the option, PhaseMapVariable.PF, where
# 'P' indicates that pressure is the search variable, and 'F' that
# it is the first call to tqmap for the current phase map.
icont = ca.tqmap(PhaseMapVariable.PF, 0, 0, vals)

# Get the pressure. Since it is the first call to tqmap, we _know_
# that this will be the value of value1.
value1 = ca.tqgetr(ResultVariable.P, -1, -1)

# Print a header for the table to come.
print(f'\n{"Pressure/bar":<24}{"Stable phases":<16}{"Amount/mol"}')

# Call psp (print stable phases), passing the value for the pressure
# we just retrieved. psp will print a formatted entry for the result
# table.
psp(value1, 0.0)

# tqmap is called again. Note the option, PhaseMapVariable.PN, where
# 'P' indicates that temperature is the search variable, and 'N'
# that it is _not_ the first call to tqmap for the current phase
# map.
icont = ca.tqmap(PhaseMapVariable.PN, 0, 0, vals)

# Get the pressure.
value1 = ca.tqgetr(ResultVariable.P, -1, -1)

# Print the entry for the result table.
psp(value1, 0.0)

# If the value of icont is positive, we need to call ca.tqmap
# again...
while icont:
    icont = ca.tqmap(PhaseMapVariable.PN, 0, 0, vals)
    value1 = ca.tqgetr(ResultVariable.P, -1, -1)
    psp(value1, 0)

# The result table is complete now. Compare the phase boundaries
# found with the phase diagram in Figure 2. Note that the
# total amount in the system is always calculated to be 2 mol.
# This is due to the fact that the phase constituents of all
# mixture phases in this data-file (Spinell, Olivine, and
# \DF-phase) have been scaled with the factor 0.5 in order to
# obtain the correct ideal entropy between Fe and Mg.

# The last one-dimensional phase map will use composition as the
# search variable, keeping both temperature and pressure constant.
# As with other phase target calculations which use composition as
# the target variable, the input in this case is a bit more
# complicated (see also worked example 3: Phase equilibria and
# phase target calculations). This is due to the fact that when
# one is interested in comparing the result for this type of
# calculation with graphical phase diagram information, it is
# necessary to keep the total substance amount in the system
# constant (i.e. use concentrations rather than absolute
# amounts).

# First, reset all conditions and targets again:
ca.tqremc(-1)

# Set the pressure to 135 Kbar (see Figure 2):
ca.tqsetc(ConditionVariable.P, 0, 0, 135000.0)

# As explained above, we need to make more than one call to tqmap to
# tell ChemApp for Python about our choice of the search variable ranges. In the
# present case, we are working with a (pseudo-)binary system, so no
# more than two calls are necessary.

# With the first call, we define the composition range for Fe2SiO4 as
# having a minimum of 0 mol and a maximum of 0.2 mol, compare again
# to Figure 2. The option, PhaseMapVariable.IA0, is the same
# one as used in phase target calculations with composition as
# target variable (see again worked example 3: Phase equilibria and
# phase target calculations), and makes sure that tqmap does _not_
# actually calculate anything at this point ('0'), but only accepts
# the values passed as information on the incoming amounts ('IA').
vals = [0.0, 0.2]

icont = ca.tqmap(PhaseMapVariable.IA0, ipfesi, 0, vals)

# Now we need to tell ChemApp for Python to vary the composition of Mg2SiO4
# accordingly, keeping the total substance amount in the system
# constant. For the range of Fe2SiO4 selected above, this means
# setting vals[0] to 1 mol and vals[1] to 0.8 mol:
vals = [1.0, 0.8]

# In the subsequent call to tqmap, the option, PhaseMapVariable.IAF,
# is used, meaning that again information on incoming amounts is
# passed ('IA'), but this time ChemApp for Python is supposed to start
# calculation and determine the first result ('F'):
icont = ca.tqmap(PhaseMapVariable.IAF, ipmgsi, 0, vals)

# Retrieve the composition for the first result.
value1 = ca.tqgetr(ResultVariable.IA, ipfesi, -1)
value2 = ca.tqgetr(ResultVariable.IA, ipmgsi, -1)

# Print a header for the table to come.
print(f'\n{"X(Fe2SiO4)"}{"X(Mg2SiO4)":>15}{"Stable phases":>15}'
      f'{"Amount/mol":>15}')

# Pass both concentrations to psp to print the table entry.
psp(value1, value2)

# tqmap is called again. Note the option, PhaseMapVariable.IAN,
# where 'IA' indicates that composition is the search variable,
# and 'N' that it is _not_ the first call to TQMAP for the
# current phase map.
icont = ca.tqmap(PhaseMapVariable.IAN, 0, 0, vals)

# Retrieve the composition...
value1 = ca.tqgetr(ResultVariable.IA, ipfesi, 0)
value2 = ca.tqgetr(ResultVariable.IA, ipmgsi, 0)

# ... and add an entry to the output table.
psp(value1, value2)

# If the value of icont is positive, we need to ca.tqmap
# again...
while icont:
    icont = ca.tqmap(PhaseMapVariable.IAN, 0, 0, vals)
    value1 = ca.tqgetr(ResultVariable.IA, ipfesi, 0)
    value2 = ca.tqgetr(ResultVariable.IA, ipmgsi, 0)
    psp(value1, value2)

# When icont does not have a positive value any more, we know that we
# have found all possible phase boundaries within the search interval
# specified.