IsotoPy — Version 0.45

Several scripts that allow the calculation of the equilibrium fractionation factor and other parameters related to isotopic fractionation in experiments and from statistical mechanics theory.

The scripts are accessible through starting the Graphical User Interface from the IsotoPy script file in the parent directory. The IsotoPy.py file is based on several scripts for embedding matplotlib into wxPython and packages listed below.

Errors/bugs are entirely my own and should be commented upon for improvement.

Tjerk Heijboer 2008-11-07

Installation:

Before being able to use this program, you need to have a working installation of: python (2.5/2.6), numpy / scipy (1.04 / 0.60), matplotlib (0.81) and wx (2.8.4.0). Python can be downloaded from www.python.org, numpy/scipy can be downloaded from www.scipy.org, matplotlib can be downloaded from www.matplotlib.sourceforge.net. Wx for python should be automatically installed with a functioning python installation. IsotoPy can be downloaded from https://sourceforge.net/projects/isotopy/.

Alternatively for Windows users, you can download all packages at once by going to:

www.enthought.com/products/epdacademic.php?ver=Win32Intel.

Note the software has been developed on Windows 2000 with python 2.5, but has also been tested on Windows XP, Windows Vista, Mac OSX 10.4 and Ubuntu 8.04, Linux (python 2.6), so it may work on other platforms as well.

If the packages are installed and functioning the Isotopy.py file can be used on any platform by simply typing:

cd go/to/correct/directory/

python Isotopy.py

in a terminal/dos-shell or alternatively on unix/linux type platforms by making the python file executable:

cd go/to/correct/directory/

chmod 755 Isotopy.py

./Isotopy.py

make sure that the .py files have the correct line endings for your system! One can also double-click/Right click with the mouse on the Isotopy.py file and it should open the first window (Fig. 1 below).

Information on how to do a calculation using the IsotoPy program:

Fig.1 Intro screen on ubuntu 9.04

Start the program and go to “Calculation type” in the “File” menu.

1. Diffusion Calculation:

Fig. 2 Select a type of calculation you want to do: 9 types. 1) diffusion is to calculate/find the equilibrium fractionation factor and the diffusion coefficient by applying a Crank-Nicholson finite difference approach to the diffusion equation from experimental data on isotopic exchange. 2) Using the Kieffer/Polyakov approach to calculate theoretical values of the equilibrium fractionation factor. 3) Doing the same as 1. but using pseudo first order rate equation. 4) A comparison of the Northrop and Clayton approach. 5) Plot the activation energies verses temperature, with optional calculation of the uncertainties in the determination of Ea. 6) Doing the same as 1. but not by diffusion but by a simple incremental recrystallization. 7) Same as 1 but for complete equilibration during recrystallization. 8) Same as 7 but for multiple isotopes of the same elements, more than 2. 9) Synthesis using approximately the same method as in 7. but for a chemical reaction in which dissolution-precipitation was an important mechanism.

Fig. 3 First tab in the diffusion calculation mode: Left hand side are numerical parameters and yes/no switches, giving different output from the program. More explanation of these parameters is included when letting the mouse rest for several seconds on one name. The right hand side is a text field which will be filled when opening an input file such as here:

Fig. 4 First tab with experimental input data. Experimental data can be found in the “inputfiles” folder with the same name as the mode of the program. Note the experimental data is in python scripting language, e.g. it is possible to do everything that is possible in python in this text field, such as if-elif-else statements, small calculations such as done in this figure for the Vratio variable. Series of numbers are stored mostly in the “tuple” type, which are comma delimited series of numbers or strings. To make numbers floating point decimals one can add a “.” after an integer. The numerical parameters are read first and can be reset by typing a certain value for that parameter in the “experimental data” text field. The names in the left hand side bar are the same names as the actual variables.

Fig. 5) Run a calculation with the numerical parameters and the experimental data as shown. The second possibility is to stop a running calculation using “stop”. You can only run one calculation at the same time. This will generate the following tabs on screen:

Fig. 6 Output tab with some information about the calculation. This tab is always generated for each type of calculation. It is refreshed after each subsequent calculation as well.

Fig. 7 Figures of the numerical experiment (a) concentration profile through the mineral (blue lines), change in isotopic composition of the mineral at the boundary (red dot), change in isotopic composition of the fluid in equilibrium with the boundary (green diamond) (b) change in total isotopic compositions with time (c) Change in fraction of exchange (F) with time (d) Northrop and Clayton diagram, here not calculated (e) contour diagram, here not calculated (f) log(1-F) d(log(1-F)/dt), giving information on the rate constant and the order of the reaction. Note that there are some numerical glitches with the Crank-Nicholson scheme, because of the initial starting condition of the outer boundary is set to the value in equilibrium with the fluid, which causes a discontinuity in the grid. The final result however is not affected by this. After a calculation the graphics can be automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

Fig. 8 (a) normalized change in apparent equilibrium fractionation factor with time using the approach of Criss (1999), (b) same but not normalized, (c) change in rate constant with time using a pseudo first order rate equation. (d) the chi-squared ellipse near the minimized values of D and ae. After a calculation the graphics will be automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

Numbered output files (contour grid, minimization path and concentration profiles) are written to different text files in a directory created in the same directory in which the input file is situated. The new directory has the same name as the input file with extension removed. Also post-script (.ps) or scaled vector graphics (.svg) files of the figures shown here are written there. These files can be read by programs such as Inkscape (free software) or Adobe Illustrator.

2. Incremental/equilibrated recrystallization:

Fig. 9 a,b,c: Same as Fig. 2

Fig. 10 Same as Fig. 3 but then with numerical parameters and experimental data necessary to do the incremental recrystallization calculation. Note that isotopic and mass data are nearly the same as in Fig. 3. Standard deviations of the isotopic data are used in the minimization routine to minimize to the smallest difference of the chi-squared value. E.g. measured initial and final values may not be the values with the lowest chi-squared values. (also: lines starting with “#” are comment lines)

Fig. 11 Output figures from the incremental recrystallization routine. (a) Change in isotopic composition with extent of the reaction (ksi) (b) change in fraction of exchange (F) with extent of the reaction, (c) change in final fractionation factor (af) with extent of the reaction (ksi) and (d) contour plot of the log(chi squared) values for a grid of (ksi and equilibrium fractionation factor (ae) values). No figures are plotted on the second tab. After a calculation the graphics will be automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

Fig. 12 Results of the equilibrium recrystallization routine, where all plotted parameters in the sub figures are the same as in Fig. 11. After a calculation the graphics will be automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

Fig 13. Input screen for the multi-isotope calculation. Note this calculation is still under development and cannot do a minimization yet. Note that for ae and the delta values there are two values, for both minor isotopes for one experimental run. This is necessary because the ae's for these two isotopes are not a-priori known unless there is a dependency between the two. These data are an example of data for 1 experimental run. e.g. for oxygen the d018 and d017 values.

Fig. 14 Example of output figures similar to figures 11 and 12 apart from the contour diagram in (d)

4. Synthesis

Fig. 15 same as Fig. 2

Fig. 16 Input screen for numerical simulation of synthesis or mineral replacement by a dissolution-precipitation mechanism. The “nu” parameter is the amount of the atom in the formula unit of the phase. the “num” parameter is the stoichiometry of the reaction.

Fig. 17 Results of the synthesis routine giving similar results as figures 11,12 and 14. for a given synthesis experiment. Note “ksi” is defined slightly different here. (a) change in isotopic composition of the phases with change in extent of the reaction. (red dashed = final mineral, blue thick = total mineral, green = fluid) (b) change in fraction of exchange with change in extent of the reaction (ksi) with as numbers of the lines the mole fraction of mineral present in the experiment. After a calculation the graphics are automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

Fig. 18 (a) change in isotopic composition of the phases with change in extent of the reaction and different mole fractions of initial mineral present. (red dashed = final mineral, blue thick = total mineral, green = fluid) (b) change in fraction of exchange with change in extent of the reaction (ksi) with as numbers of the lines the mole fraction of mineral present in the experiment.

4. Kieffer/Polyakov-type statistical mechanical calculations:

Fig. 19 Same as Fig. 2

Fig. 20 Input screen for Kieffer/Polyakov type theoretical statistical mechanical calculation of reduced partition function ratios/Thermodynamics from data on vibrational frequencies.

Fig. 21 (a) Vibrational model for lower frequencies (cm^-1), (b) Comparison how well the thermal pressure and potential pressure match for different temperatures (10⁹ Pa). (c) Cp function, comparison with thermodynamic data. (d) Calculated isothermal bulk modulus, (e) Entropy change with temperature, (f) Internal energy change with temperature, (g) Helmholtz free energy change with temperature.

Fig. 22 (a) Change in molar volume compared to volume at P=P0 and T=298.15 °C, comparison of experimental data and calculation based on grüneisen parameters. (b) change in equilibrium fractionation factor with temperature. (c) similar to (b) but different units. (d) change in reduced partition function ratio with temperature. (e) contour diagram of change in reduced partition function ration with temperature and pressure. After a calculation the graphics will be automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

4. Comparison of Northrop and Clayton approach, Criss et al. 1987, and Ne equation, Heijboer et al, in prep.

Fig. 23 Same as Fig. 2

Fig. 24. Input tab. with similar variables as before.

Fig. 25 (a) equilibrium fractionation factor (ae) versus fraction of exchange (F), when comparing the equations of Criss 1999 and Northrop and Clayton. (b) equilibrium fractionation factor when using equation … versus the rate-constant times the duration k(t-t₀). (c) same as (b) when using the fluid. After a calculation the graphics will be automatically saved to a directory with the same name as the input file from the 'inputfiles' directory.

4. Comparison of Northrop and Clayton approach, Criss et al. 1987, and Ne equation, Heijboer et al, in prep. Using the new_NC calculation

Differences between calculation 6. and 7. are that the equations used in 6. are transformed or solved for the final d-value allowing a minimization procedure similar to the diffusion model and the recrystallization method to be used. The figures drawn and the input data are basically the same as for the recrystallization and synthesis procedures apart from the use of “kt”, the rate constant times the duration of the experiment as an unknown value (similar to the extent and D). These unknowns have in common that they all relate to the amount of reaction or the speed by which the reaction proceeds.

Fig. 26 The same as Fig. 2

Fig. 27 Input screen of the NC_New calculation. Note the inclusion of an extra error_parameter which will be reset during the minimization.

Fig. 28 (a) shows the change in total isotopic composition of the mineral and the fluid with change in duration for a given calculated value of the rate constant k. (b) Contour diagram near the minimized values of ae and kt showing in this case a linear minimum.

4. Plotting of D₀ and Activation energy

Plotting activation energy and D0 from diffusion constants derived by the diffusion calculations in calculation type 1.

Fig. 29 Same as Fig. 2

Fig. 30 Input data, temperature and diffusion coefficients calculated from calculation type 1.

Fig. 31 Calculated values of log D0 and Ea or Q, the activation energy with the specific error.

Fig. 32 calculated D0 and Activation energies in Arrhenius plot of log(D) vs. 1000/T. Experiments of Vennemann and O’Neil (1996).

Miscellaneous options that are possible for all or most calculations

1. Graphics format change

Fig. 33 Change of graphics output formats

Fig. 34 Choose between .ps, .jpg, .svg, .png graphic types to export the images.

2. Keep the figure during multiple runs

Fig. 35 keep figure is an on-off switch to be able to plot multiple experiments in one figure for the same type of experiment.

Fig. 36 The switch can result in a figure like this...

3. Calculate the data from multiple experiments at the same time

Fig. 37 Select multiple files in the open file menu and press open.

Fig. 38 This results in a table being produced from these different input files and a spreadsheet (the second tab) Data can be easily and comprehensively manipulated.

Fig. 39 Spreadsheet of the experimental data. Note that all values for ae and some other parameters are copied to fill in for the number of runs that were done. Only the first value is actually used. Data can be simply saved to a .csv file for manipulation in a spreadsheet program.

Detailed explanation:

The program is divided into the main program (Isotopy.py), which is just a graphical user interface allowing the visualization of the calculations. This program, dynamically determines which of the subdirectories in the IsotoPyx_y directories is a python package (e.g. the directory termed Diffusion, which contains an __init__.py file) and then looks for the main.py file, the pars.py file and the help_strings.py file.

The main.py file contains the calculation code or (in the case of the kieffer routine) the call to the calculations, the pars.py file of a certain directory contains numerical parameters (e.g. whether to plot a diagram or not) and the help_strings.py file contains the help to these parameters that are set and used in the graphical interface (Fig. 28). Each .py file can be opened and analysed/changed to the users needs. In addition to the files, that should be located in a sub-directory, the IsotoPy.py file can open an input file which in default cases should be located inside a directory at the same level as the main.py file (e.g. /IsotoPyx_y/Diffusion/inputfiles/ Input.py). These .py files contain the variables pertaining to the experiment from which a calculation is done (e.g. data for a chlorite-water exchange experiment) The opened file can be changed by typing in the text field and saved. The changed code is then used to run a particular calculation. Numbered output files are generated in a subdirectory of the input files, which in the case of the diffusion routine stores most of the produced data, together with the graphics. In the case of the kieffer routine an additional package of python files is located inside the /IsotoPyx_y/routine/thermo/ sub directory which do all the detailed calculations.

Fig. 40 The program has the following flow scheme.

Example of a pars.py file:

adjustment = True,False

min_adjustment = True,False

ae_min_max_step = 0.99,1.01, 1000.0,\

1.5 , 2.0, 1000.0,\

0.9 , 1.1, 1000.0,\

1.7 , 2.5, 1000.0,\

1.2 , 1.6, 1000.0,\

0.95,1.05, 1000.0,\

0.95,1.00, 1000.0

atom = "oxygen","hydrogen","carbon","XX"

and its corresponding help_strings.py file

str_adjustment = "adjustment: ","adjustment",True

tt_str_adjustment = "set True or False whether to adjust the initial measured d values of the fluid"

str_min_adjustment = "min adjustment: ","min_adjustment",False

tt_str_min_adjustment = "whether to adjust the final isotopic composition of the mineral"

str_atom = "atom: ","atom","hydrogen"

tt_str_atom = "set which atom is used for the isotopic ratios"

str_aelist = "aelist: ","ae_min_max_step",0

tt_str_aelist = """tuple containing the min,max and the number of steps of the eq. fractionation factor in the contour diagram """

For each parameter (for example adjustment) in pars.py there is a string parameter such as str_adjustment in Help_Strings.py that gives the text placed on the left side of the screen, the name of the parameter and the default value. The tt_str_ parameter gives information about what the parameter stores. This information is given when hovering with the mouse pointer over the parameter in the graphical interface.

When adding or changing the parameters in pars.py a str_ and tt_str_ parameter must be given in the help_strings.py file so that each parameter can be correctly shown in the graphical interface.

A calculation module can be simply added to the package by creating a directory containing the following:

/__init__.py

/main.py

/pars_directoryname.py

/help_strings.py

/inputfiles/

the main.py file should at least contain the following class definition:

class Main:

def __init__(self):

self.fig1name ="not used"

self.fig2name ="not used"

def main(self):

return 0