! this file demonstrates the use of vertex and frendly to calculate
! a silica chemical potential (u(SiO2)) vs. fluid composition (YCO2)
! Schreinemakers diagram for the system CaO-SiO2-Al2O3-SiO2-H2O-CO2-O2
! at constant pressure, temperature, and oxygen chemical potential.
! this example differs from the previous ones in that chemical
! potentials are used as independent variables, a difficulty in
! such calculations is that usually users have no idea what range
! of chemical potential values are appropriate for a calculation.
! In this case there are two possible approaches, either one
! calculates an initial phase diagram for a very wide range of
! chemical potentials, and then examines and refines interesting
! regions of the resulting diagram; or one calculates limits on
! the chemical potentials based on saturation values and activity
! models. Frendly is useful in the second approach is applied here.
! For this problem, which is relevant to skarn formation pressure
! and temperature are assumed (T=863 K, P=2000 bar).
! The first
! task is to estimate the range of u(SiO2) of interest. The upper
! limit can be taken as the saturation value of quartz, and the
! lower limit as the Gibb free energy of quartz with some reduced
! activity, say 0.01 in the present case (alternatively a "low"
! u(SiO2) could be obtained by taking the G for the reaction
! En - Fo at 863 and 2000 bar).
! Note that Vertex will not recognize saturation of a Mobile component
! (in fact this would be inconsistent with its algorithm), and so
! it is the users responsibility to make sure that reasonable limits
! on chemical potentials are used.
! u(O2) must also be input in the present problem, and here it
! has been taken arbitrarily to be that corresponding to
! log f(O2) = -10 (very oxidizing conditions). More generally
! i recommend using frendly to calculate a limit on u(O2) from
! phase equilibria (e.g., QFM), or better still to create a
! buffer phase which gives u(O2) consistent with some equilibria
! such as QFM as a function of, at least, P and T (see documentation
! section 1.4, example 1.4).
! Note that fluid equations of state programmed in rlib.f do not
! take account of the effects of u(O2) on fluid speciation; thus,
! unless you add the equations yourself, results will only be
! valid under (oxidizing) conditions in which the fluid can be
! assumed to be essentially a binary mixture of H2O and CO2.
! here the user first runs frendly to get u(SiO2) and u(O2), i hope
! most of this is fairly self-explanatory.
buzzard{jamie}16: frendly
Enter the name of the thermo data file, e.g.,"hp90ver.dat"
the name must be < 15 characters and left justified,.
hp90ver.dat
Hi! i am a user freindly program.
what is your name?
! the program starts by asking you some irritating questions (feel
! free to get rid of them).
jamie
Hi jamie !, i like you and i hope we have fun!
can you say "therm-o-dy-nam-ics" (y/n)?
y
You can use frendly for the following options:
1) calculate equilibrium coordinates for a reaction.
2) calculate thermodynamic properties for a phase or reaction relative to
the reference state.
3) calculate the change in thermodynamic properties from one p-t-x
condition to another.
4) create new entries in the thermodynamic data file for frendly or vertex.
5) quit.
If you select options 1-3 you may also modify thermodynamic data, the
modified data can then be stored as a new entry in the thermodynamic data
file.
Enter the number identifying the option (above) of your choice:
2
Do you want a list of the phases in the database (y/n)?
n
! as with most of the Perplex programs a blank answer to a
! a y/n question is read as no (but be careful, this isn't
! true all the time, please tell me if you find examples).
Do you want to calculate thermodynamic properties for a reaction (y/n)?
n
Enter the name of the phase that you want to calculate thermodynamic properties for:
Q
Enter the activity of: Q (enter 1.0 for all fluid species)
1
standard state properties follow:
g(j/mole) s(j/mole k) v(j/mole bar)
-856400. 41.5000 2.26900
heat capacity (j/mole k) function:
a b*t c/t**2 d*t**2 e/t**1/2
97.9000 -0.335000E-02 -636200. 0. -774.000
f/t g/t**3
0. 0.
volumetric (j/mole bar) function:
b2*(tr-t) b3*(exp(-t/300)-exp(tr/300)
0.800000E-04 0.
b4*(p-pr) b5*(exp(-p/35.d3)-exp(-pr/35.d3)
-0.592000E-05 0.
super function for g (j/mole):
-848473. 706.425 2.24515 97.9000 -0.167500E-02
-318100. 0. -3096.00 0. 0.
0.800000E-04 0. -0.296000E-05 0. 0.
0.
now jamie there is one thing i want you to remember and that is that the
values for g and h that i calculate for you are "apparent" free energies and
enthalpies. jamie if you dont know what that means, read helgeson et al. 1978.
jamie now i am going to prompt you for the conditions at which i should
calculate thermodynamic properties, when you want to stop just enter zeroes.
ok, here we go, its fun time jamie !!
Enter p(bars) and t(k) (zeroes to quit):
2000 863
At 863.000 k and 2000.00 bar:
! the value of G here is the saturation value for Q, i.e., disregarding
! kinetic problems, the upper limit for u(SiO2) at 863 K and 2000 bar.
G(kj) = -894.718
H(kj) = -805.492
A(kj) = -899.323
U(kj) = -810.097
S(j/k) = 103.391
V(j/bar)= 2.30235
Cp(j/k) = 67.8066
Do you want to modify or output the thermodynamic parameters of a phase (y/n)?
! to get the lower limit on u(SiO2), the user now must lower the
! activity of quartz
y
Do you only want to ouput data (y/n)?
n
Enter the number of the phase you want to modify:
1) Q
1
Thermodynamic properties are calculated from the reference state constants
g, s, v, and an activity coefficient together with the isobaric heat capacity function:
cp(pr,t)=a+b*t+c/(t*t)+d*t*t+e/t**(1/2)+f/t+g/t**3
and the volumetric function:
v(p,t)=v(pr,tr)-b2*(t-tr)+b3*(exp(-t/300)-exp(tr/300))+b4*(p-pr)
+b5*(exp(-p/35000)-exp(-pr/35000))+b6*(p-pr)**2+b7*(t-tr)**2
Enter the number of the parameter you wish to modify:
1) standard free energy g (j) 2) standard entropy s (j/k)
3) standard volume v (j/bar) 4) heat capacity coefficient a
5) heat capacity coefficient b 6) heat capacity coefficient c
7) heat capacity coefficient d 8) heat capacity coefficient e
9) heat capacity coefficient f 10) heat capacity coefficient g
11) volumetric coefficient b2 12) volumetric coefficient b3
13) volumetric coefficient b4 14) volumetric coefficient b5
15) volumetric coefficient b6 16) volumetric coefficient b7
17) activity 18) reaction coefficient
and enter a zero when you are finished
17
Enter a name (<8 characters left justified) to distinguish the modified
version of Q
WARNING: if you intend to store the modified data in the data file this
name must be unique.
q
The old value for the activity of q was 1.00000
Enter the new value:
0.01
Enter the number of the parameter you wish to modify:
1) standard free energy g (j) 2) standard entropy s (j/k)
3) standard volume v (j/bar) 4) heat capacity coefficient a
5) heat capacity coefficient b 6) heat capacity coefficient c
7) heat capacity coefficient d 8) heat capacity coefficient e
9) heat capacity coefficient f 10) heat capacity coefficient g
11) volumetric coefficient b2 12) volumetric coefficient b3
13) volumetric coefficient b4 14) volumetric coefficient b5
15) volumetric coefficient b6 16) volumetric coefficient b7
17) activity 18) reaction coefficient
and enter a zero when you are finished
0
Do you want to store q as a permanent entry
in the thermodynamic data file (hp90ver.dat ) (y/n)?
n
Do you want to modify the properties of another phase (y/n)?
n
! the modified activity will be apllied herein, i.e., the user
! can now calculate a lower limit for u(SiO2):
Enter p(bars) and t(k) (zeroes to quit):
2000 863
At 863.000 k and 2000.00 bar:
G(kj) = -927.751
H(kj) = -805.481
A(kj) = -932.378
U(kj) = -810.108
S(j/k) = 141.680
V(j/bar)= 2.31345
Cp(j/k) = 67.8066
Do you want to modify or output the thermodynamic parameters of a phase (y/n)?
n
Enter p(bars) and t(k) (zeroes to quit):
0 0
You can use frendly for the following options:
1) calculate equilibrium coordinates for a reaction.
2) calculate thermodynamic properties for a phase or reaction relative to
the reference state.
3) calculate the change in thermodynamic properties from one p-t-x
condition to another.
4) create new entries in the thermodynamic data file for frendly or vertex.
5) quit.
If you select options 1-3 you may also modify thermodynamic data, the
modified data can then be stored as a new entry in the thermodynamic data
file.
! now the user calculates the oxygen chemical potential corresponding
! to log fO2 = -10.
Enter the number identifying the option (above) of your choice:
2
Do you want a list of the phases in the database (y/n)?
Do you want to calculate thermodynamic properties for a reaction (y/n)?
Enter the name of the phase that you want to calculate thermodynamic properties for:
O2
Enter the activity of: O2 (enter 1.0 for all fluid species)
! note that here i have first calculated u(O2) for a fugacity
! (= activity in this case) of log fO2 = -20. this value is subsequently
! modified.
1.d-20
standard state properties follow:
g(j/mole) s(j/mole k) v(j/mole bar)
-14.9084 205.200 0.
heat capacity (j/mole k) function:
a b*t c/t**2 d*t**2 e/t**1/2
48.3000 -0.691000E-03 499200. 0. -420.700
f/t g/t**3
0. 0.
volumetric (j/mole bar) function:
b2*(tr-t) b3*(exp(-t/300)-exp(tr/300)
0. 0.
b4*(p-pr) b5*(exp(-p/35.d3)-exp(-pr/35.d3)
0. 0.
super function for g (j/mole):
62998.3 164.009 0. 48.3000 -0.345500E-03
249600. 0. -1682.80 0. 0.
0. 0. 0. 0. 0.
0.
now jamie there is one thing i want you to remember and that is that the
values for g and h that i calculate for you are "apparent" free energies and
enthalpies. jamie if you dont know what that means, read helgeson et al. 1978.
jamie now i am going to prompt you for the conditions at which i should
calculate thermodynamic properties, when you want to stop just enter zeroes.
ok, here we go, its fun time jamie !!
Enter p(bars) and t(k) (zeroes to quit):
2000 863
At 863.000 k and 2000.00 bar:
G(kj) = -457.160
H(kj) = 79.1279
A(kj) = -457.160
U(kj) = 79.1279
S(j/k) = 621.422
V(j/bar)= 0.
Cp(j/k) = 34.0530
Do you want to modify or output the thermodynamic parameters of a phase (y/n)?
y
Do you only want to ouput data (y/n)?
Enter the number of the phase you want to modify:
1) O2
1
Thermodynamic properties are calculated from the reference state
constants g, s, v, and an activity coefficient together with the
isobaric heat capacity function:
cp(pr,t)=a+b*t+c/(t*t)+d*t*t+e/t**(1/2)+f/t+g/t**3
and the volumetric function:
v(p,t)=v(pr,tr)-b2*(t-tr)+b3*(exp(-t/300)-exp(tr/300))+b4*(p-pr)
+b5*(exp(-p/35000)-exp(-pr/35000))+b6*(p-pr)**2+b7*(t-tr)**2
Enter the number of the parameter you wish to modify:
1) standard free energy g (j) 2) standard entropy s (j/k)
3) standard volume v (j/bar) 4) heat capacity coefficient a
5) heat capacity coefficient b 6) heat capacity coefficient c
7) heat capacity coefficient d 8) heat capacity coefficient e
9) heat capacity coefficient f 10) heat capacity coefficient g
11) volumetric coefficient b2 12) volumetric coefficient b3
13) volumetric coefficient b4 14) volumetric coefficient b5
15) volumetric coefficient b6 16) volumetric coefficient b7
17) activity 18) reaction coefficient
and enter a zero when you are finished
17
Enter a name (<8 characters left justified) to distinguish the modified
version of O2 . WARNING: if you intend to store the modified
data in the data file this name must be unique.
o2
The old value for the activity of o2 was 0.100000E-19
Enter the new value:
1.d-10
Enter the number of the parameter you wish to modify:
1) standard free energy g (j) 2) standard entropy s (j/k)
3) standard volume v (j/bar) 4) heat capacity coefficient a
5) heat capacity coefficient b 6) heat capacity coefficient c
7) heat capacity coefficient d 8) heat capacity coefficient e
9) heat capacity coefficient f 10) heat capacity coefficient g
11) volumetric coefficient b2 12) volumetric coefficient b3
13) volumetric coefficient b4 14) volumetric coefficient b5
15) volumetric coefficient b6 16) volumetric coefficient b7
17) activity 18) reaction coefficient
and enter a zero when you are finished
0
Do you want to store o2 as a permanent entry
in the thermodynamic data file (hp90ver.dat ) (y/n)?
Do you want to modify the properties of another phase (y/n)?
Enter p(bars) and t(k) (zeroes to quit):
2000 863
At 863.000 k and 2000.00 bar:
G(kj) = -291.941
H(kj) = 79.1279
A(kj) = -291.941
U(kj) = 79.1279
S(j/k) = 429.976
V(j/bar)= 0.
Cp(j/k) = 34.0525
Do you want to modify or output the thermodynamic parameters of a phase (y/n)?
n
Enter p(bars) and t(k) (zeroes to quit):
0 0
You can use frendly for the following options:
1) calculate equilibrium coordinates for a reaction.
2) calculate thermodynamic properties for a phase or reaction relative to
the reference state.
3) calculate the change in thermodynamic properties from one p-t-x
condition to another.
4) create new entries in the thermodynamic data file for frendly or vertex.
5) quit.
If you select options 1-3 you may also modify thermodynamic data, the
modified data can then be stored as a new entry in the thermodynamic data
file.
Enter the number identifying the option (above) of your choice:
5
Have a nice day jamie !
! The user is now armed with a range of values for u(SiO2) (roughly
! -925000 to -890000 j/mol-SiO2) and a value for u(O2) (roughly
! -300000 j/mol-O2).
! The next step is to run build to prepare the input for Vertex:
! NOTE: this calculation is slightly oversimplified in that only
! two mineral solutions, grandite and epidote, are considered.
! However the result is more easily understood. This is also a
! problem in which you may want to try calculating the pseudo-
! univariant equilibria as well, because they don't result in
! a completely illegible diagram (The interpretaion of pseudo-
! univariant equilibria is discussed in Connolly and Kerrick
! 1987 CALPHAD, Connolly 1990 AJS, and Connolly and Trommsdorff
! 1991 CMP)
! An interesting point demonstrated here is that increasing a
! a chemical potential has a similar effect to increasing the
! thermal potential, hence the qualitative similarity of the
! calculated u(SiO2)-YCO2 diagram with the more common T-YCO2
! diagram.
buzzard{jamie}17:
NO is the default answer to all Y/N prompts
Enter the name of the input file for Vertex (i.e. the file
OUTPUT by this program) < 15 characters and left justified:
in6.dat
The file in6.dat exists, do you really want to overwrite it (y/n)?
y
Enter the thermo data file name (e.g. hp90ver.dat),
< 15 characters, left justified:
hp90ver.dat
Do you want to generate a print file (Y/N)?
y
Do you want to generate a graphics file (Y/N)?
y
Enter a name for the print file, <15 characters, left justified:
print6.out
Enter a name for the graphics file, <15 characters, left justified:
plot6.out
Specify the kind of phase diagram calculation:
0 - for a Composition diagram
1 - for a Schreinemakers-type diagram
3 - for a Mixed-variable diagram
1
Specify the reliability level [1-5, default is 5]:
1 - gives lowest efficiency, highest reliability
5 - gives highest efficiency, lowest reliability
High values increase the probability that a curve
may be only partially determined or skipped entirely.
5
Print dependent potentials for chemographies?
Answer no if you do not know what this means.
n
The data base components are:
NA2O MGO AL2O3 SIO2 K2O CAO TIO2 MNO FEO O H2O CO2
Do you want to redefine them (Y/N)?
n
Do you want to do calculations with a saturated phase (Y/N)?
The phase is: FLUID
Its components can be: H2O CO2
Its compositional variable is called: X(CO2)
y
Enter the number of components in the FLUID
(1 or 2 for COH buffered fluids):
2
Do you want to do calculations with saturated components (Y/N)?
n
Do you want to treat the chemical potential of a
component as an INDEPENDENT variable (Y/N)?
y
NOTE: The chemical potential of a component, i.e. a "mobile component,"
is symbolized, here and in VERTEX, by an uppercase U followed by the component
name enclosed in parentheses, e.g. the chemical potential of A is written U(A).
Select "mobile" components from the following set:
NA2O MGO AL2O3 SIO2 K2O CAO TIO2 MNO FEO O
How many "mobile" components (maximum 2)?
2
Enter component names, left justified, one per line:
SIO2
O
Select remaining components from the following set:
NA2O MGO AL2O3 K2O CAO TIO2 MNO FEO
How many thermodynamic components (minimum 2, maximum 7)?
3
Enter component names, left justified, one per line:
AL2O3
CAO
FEO
working...
Do you want to exclude phases from your calculations (Y/N)?
n
Select the x-axis variable:
1 - P(bars)
2 - T(K)
3 - X(CO2)
4 - U(SIO2 )
5 - U(O )
3
Enter the minimum and maximum values, respectively, for: X(CO2)
0 .8
Select the y-axis variable:
2 - T(K)
3 - P(bars)
4 - U(SIO2 )
5 - U(O )
4
Enter the minimum and maximum values, respectively, for: U(SIO2 )
! be careful with negative numbers, build won't stop you from
! entering min > max, though vertex will.
-925000
-890000
Calculate sections as a function of a third variable (Y/N)?
n
Specify the sectioning value for: P(bars)
2000
Specify the sectioning value for: T(K)
863
Specify the sectioning value for: U(O )
-300000
Select the equation of state for the FLUID phase
1 - MRK (DeSantis et al 1974)
2 - HSMRK (Kerrick and Jacobs 1981)
3 - Hybrid MRK-HSMRK
4 - Saxena and Fei 1987, pseudo-virial expansion
5 - Bottinga and Richet 1982, RK
6 - Holland and Powell 1990, CORK
7 - Hybrid Haar/HSMRK
8 - Graphite buffered COH-MRK fluid
6
Do you want to treat solution phases (Y/N)?
y
Enter the name of the file which contains the solution models you want to
use (e.g.solut.dat), left justified, < 15 characters:
solut.dat
**warning ver025** missing all endmembers for GlTr(i)
**warning ver025** missing all endmembers for Tschermaki
**warning ver026** only one endmember for Pl(i)
**warning ver026** only one endmember for MaPa(i)
**warning ver026** only one endmember for Pl(hB)
**warning ver026** only one endmember for Pl(h)
**warning ver026** only one endmember for DiCats(i)
**warning ver025** missing all endmembers for Ab(h)
**warning ver025** missing all endmembers for Kf(h)
**warning ver025** missing all endmembers for Ab
**warning ver025** missing all endmembers for Kf
**warning ver025** missing all endmembers for Kspar(h)
**warning ver025** missing all endmembers for Kspar
**warning ver025** missing all endmembers for Sanidine
**warning ver025** missing all endmembers for Omph(i)
**warning ver026** only one endmember for Chl(i)
**warning ver025** missing all endmembers for Bio(i)
**warning ver026** only one endmember for St(i)
**warning ver026** only one endmember for Carp(i)
**warning ver026** only one endmember for Crd(i)
Fluid is not a possible solution phase because of endmember H2O
**warning ver026** only one endmember for Cumm(i)
**warning ver026** only one endmember for Anth(i)
**warning ver026** only one endmember for Tr(i)
**warning ver025** missing all endmembers for Mica
**warning ver025** missing all endmembers for Pa
**warning ver025** missing all endmembers for Mu
** warning ver111 ** the following endmembers are missing for solution gralad
fefe
it will be treated as a simpler solution between the remaining endmembers
**warning ver026** only one endmember for EnFs(i)
**warning ver026** only one endmember for HeDi(i)
**warning ver026** only one endmember for Talc(i)
**warning ver026** only one endmember for Ol(i)
**warning ver026** only one endmember for Dol(i)
**warning ver026** only one endmember for Sp(J&R)
**warning ver026** only one endmember for Sp(G&S)
**warning ver025** missing all endmembers for Neph(F&B)
**warning ver026** only one endmember for Sp(i)
** warning ver111 ** the following endmembers are missing for solution Gt(i-c)
py
it will be treated as a simpler solution between the remaining endmembers
** warning ver111 ** the following endmembers are missing for solution Gt(i-b)
py
it will be treated as a simpler solution between the remaining endmembers
**warning ver025** missing all endmembers for tr-ed
**warning ver025** missing all endmembers for mgal-amph
**warning ver026** only one endmember for Ctd(i)
**warning ver026** only one endmember for mg-fe-hb
** warning ver111 ** the following endmembers are missing for solution GrPyAlSp(B
py
it will be treated as a simpler solution between the remaining endmembers
In solution GrPyAlSp(B the following endmembers have fixed compositions:
spess X = 0.100000
To relax or change these constraints you must modify solut.dat .
** warning ver111 ** the following endmembers are missing for solution GrPyAl(B)
py
it will be treated as a simpler solution between the remaining endmembers
**warning ver025** missing all endmembers for JdDi(G1)
**warning ver025** missing all endmembers for JdDi(G2)
**warning ver025** missing all endmembers for JdDi(W?)
Select phases from the following list, enter the
names one at a time and left justified, ENTER A BLANK WHEN YOU ARE doNE.
gralad GrAd(E&W) GrAd(i) EpCz(i) h-EpCz(i) Gt(i-c)
Gt(i-b) GrPyAlSp(B GrPyAl(B)
EpCz(i)
GrAd(i)
Calculate high variance phase fields (Y/N)?
! here the user is requesting that pseudounivariant equilibria are
! calculated, you may want to simplify the calculation by answering
! no here (alternatively you can suppress the plotting of the
! pseudounivariant curves in PSVDRAW as done later in this example)
y
Enter a one-line title for your calculation:
sample 6
buzzard{jamie}18: vertex
Enter the name of the computational option file (i.e., the
file created with BUILD) the name must be < 15 characters
and left justified.
in6.dat
Reading thermodynamic data from file: hp90ver.dat
Writing print output to file: print.out6
Writing plot output to file: plot.out6
Reading solution models from file: solut.dat
Initial number of divariant fields to be tested is: 21
Testing initial assemblage # 1 # of new divariant assemblages is: 0
finished with equilibrium ( 1).
finished with equilibrium ( 2).
finished with equilibrium ( 1).
finished with equilibrium ( 3).
finished with equilibrium ( 4).
.
.
blah blah blah
.
.
Testing new divariant assemblage # 235 number remaining to be tested is: 0
! the user now runs psvdraw twice, the first time suppressing the
! pseudounivariant curves and reaction labels, and the second time
! showing the pseudounivariant curves and suppressing invariant
! point labels.
buzzard{jamie}25: psvdraw
Enter plot file name
plot.out6
PostScript will be written to file:
plot.out6.ps
Do you want to modify the default plot (y/n)?
y
modify x-y limits (y/n)?
fit curves with B-splines (y/n)?
y
suppress curve labels (y/n)?
y
suppress point labels (y/n)?
n
suppress high variance curves (y/n)?
n
restrict phase fields (y/n)?
buzzard{jamie}26: lp plot.out6.ps
buzzard{jamie}27: psvdraw
Enter plot file name
plot.out6
PostScript will be written to file:
plot.out6.ps
Do you want to modify the default plot (y/n)?
y
modify x-y limits (y/n)?
n
fit curves with B-splines (y/n)?
y
suppress curve labels (y/n)?
y
suppress point labels (y/n)?
n
suppress high variance curves (y/n)?
y
restrict phase fields (y/n)?
buzzard{jamie}28: lp plot.out6.ps