Recent Posts

Pages: 1 2 3 4 5 6 7 8 9 10
alphaMELTS for MATLAB/Python / Re: Total failure of Matlab when using MELTS
« Last post by Paula on March 04, 2021, 12:09:01 pm »
Hi Simon,

So if you run your script as is then on the last step (i=14; index=15 in the _tbl.txt files) it adds a third cpx and a third spinel. The compositions are virtually identical to the second cpx and spinel respectively, so it really shouldn't have added them. If you try to add a break point within MATLAB, e.g. with "dbstop if error", it gives a memory error just like it did when you tried restarting the script. I suspect that is all linked but I'm not sure how.

The simplest way to get the script to run to i=1000 is to add the following:
path.engine.setSystemProperties("Limit number", "spinel 2");
path.engine.setSystemProperties("Limit number", "clinopyroxene 2");
It you limit spinel but not cpx, which is limited to 3 by default anyway, then it gets to i=599 and crashes when hornblende comes in. Likewise if you limit cpx to 2 but not spinel. In both cases there is additional weirdness. I'm not sure why it is adding the extra phases in the way it is but I have some ideas, which may also suggest a more general workaround.

There are several other things to try in addition, and one or more might make your runs go more smoothly:
  • Increase the temperature increment if you can. Unless you are trying to model the granite ternary minimum, something like 1 or 2oC should work pretty well. With 0.1oC I'm not that surprised you are hitting problems with the very small energy differences between the various pyroxenes etc.
  • Possibly remove the tiny amount of Cr2O3 from the bulk, as it overstabilizes spinel.
  • Start the run at a higher temperature. The terminal output shows that the algorithm is struggling on the first calculation, and I definitely saw one instance where it put metastable hornblende in and then dropped it on the next go.

Anyway, I will try to debug this behavior properly at some point. In the meantime, hopefully that gives you and other something to work with to make the code more reliable.

alphaMELTS for MATLAB/Python / Re: Total failure of Matlab when using MELTS
« Last post by Paula on March 03, 2021, 10:38:38 am »
Short answer is that it's getting stuck in a local minimum, which is why it fails. Not unrelated, it's trying to exsolve multiple instances of cpx, spinel etc. that are probably not real. I think ultimately that this so much adding and dropping of phases leads to the segmentation fault, though I haven't tried running it through the debugger or memory checker to find out exactly where yet. Debugging the mixed C and MATLAB in this way is possible in the Linux version of MATLAB but it is extremely slow...

In the meantime, I'm going to outline first how to go about tracking down problems in alphaMELTS for MATLAB/Python, and then some suggested strategies to work around the specific problem that Simon describes.

So first thing is always start MATLAB from the command line on Mac or Linux so that MELTS-related output goes to the terminal. On Windows a console will open automatically. For Python it will depend on which IDE you are using and whether you are on Windows or not. If you are having trouble seeing the output for an example that does work, like the tutorial, then let me know.

Most of the time you do not need to see this output, but if something does go wrong it can help with diagnosis. It is also useful for checking whether a particular setting has been picked up. If you set a breakpoint just after the setSystemProperties() line in the script then you should be able to see a "Processed line..." message in the terminal output.

The other trick is that, by default, the latest version of alphaMELTS for MATLAB/Python will write out alphaMELTS-style output files. (There are other output options but that is for a different thread.) If the MATLAB/Python version crashes then afterwards navigate in terminal to where the output files are and do "run-alphamelts.command -x". This will make the Phase_mass_tbl.txt file without doing any extra calculations, and that file can be very useful for working out what happened.

More soon...
This is one of those "not so much a bug as a feature" situations. It's deliberate that you might want, say, to equilibrate the model and then do "Supplemental Calculator" type calculations for actual samples to compare. The allowed phase names for calcPhaseProperties etc. are a regular expression:

(phase name in MELTS) + (optional digits, as output by calcEquilibriumState, but could also be user defined) + (optional '_' followed by alphanumeric string e.g. sample name).

The code doesn't have a built in way to tell if names like "spinel1" are the result of a previous equilibration calculation or a user-defined name for an instance of that phase. It would be possible to add that check (and I will), but really the plan was that getListProperty would be used on lists built up with addNodeBefore/After calls, and/or by joining different lists together. So just make sure to use addNodeAfter before each equilibration that you want to keep the output from. Consider that a bit like using menu option 4 in standalone alphaMELTS.

If you need to do a few extra calculations to get to where you need to be at the start of the run, then you can do those in place. Then use addNodeAfter and First.clearNode to clean up the list before getListProperty. This is a bit like menu option 3 in standalone alphaMELTS, except that at the moment it's your responsibility to discard the first node before the main calculation.

Note that in the Python version you should always use clearNode() and not removeNode(), as clearNode() tidies up the list after removing the node i.e. makes sure Prev, Next, First, Last is correct for all remaining nodes. removeNode() is a private method in MATLAB so this is not an issue in that version.

Likewise the copy. If you assign path2 = path1 then path1 and path2 point to the same object. This behavior is inherited from the matlab.mixin.Copyable class that MELTSdynamic builds on. Users going from MATLAB to Python often get tripped up in a similar way. So the intended way to start a new list is:

path2 = path1.copy;
        ... or ...
path2 = path1.copyAndKeepOutput;

Both copy the system settings and input parameters to the first node in a new list. The first one resets the output, the second one does exactly what it says on the tin. For an example of copyAndKeepOutput see tutorial.m.

Hope that makes sense!
Dear Paula,

Thank you very much for the detailed guidance. Following your advice, I found a little strange result. alphamelts 1.9 is used, and the steps are:
step 1: choose 1 (to read .melts file)
step2: LPUM.melts (input .melts file name)
step 3: 7 (impose initial entropy, enthalpy or volume)
step 4: 1 (to use current total enthalpy)
step 5: 3 (single/batch calculation)
step 6: 1 (initial superliquidus guess)
step 7: 7 (impose initial entropy, enthalpy or volume)
step 8: -116000.0 (new system total enthalpy J per 100 g)
step 9: 3 (single/batch calculation)
step 10: 0 (quit)

steps 1-7 are used to get current system total enthalpy, -115900.0 J per 100 g, for example. I need extract dQ=100.0 J per 100g from the system, so -116000.0 is input. These will give new temperature that is compatible with dQ. Setting file is:

ALPHAMELTS_DELTAT -0.1 (very small if dQ is small)

Accordingly, Log fO2 Path : IW and Log fO2 Delta: -1.0.

The above, if it runs just once, gives new temperature after dQ is removed. However, if the outcome above (new liquid composition, new T) is set as input of the next calculation, the steps 1-7 would give a new current system total enthalpy that is different from the result of previous enthalpy, that's strange. Say, H0 is initial system total enthalpy, H is new system total enthalpy:
calculation 1: H0=-115900, H=-116000
calculation 2: H0=-116037, H= -116137
calculation 3: H0=-116147, H=-116247
as you can see, H in calculation 1 is different from H0 in calculation, etc. Since at T~1800 C, system is always superliquidus, thus liquid composition stays constant. So, the only variable is temperature T. Why H in calculation is different from H0 in calculation 2? This is puzzling. Could you please give me some hints?

Hi Bhuvan,

It's Fs = En + 2Hd - 2Di, but that's probably just a typo. The reason you are getting negative Wo is that you are not including the Ca in the hedenbergite or other components.

For a pyroxene that lies in the quadrilateral, using independent end members En, Di, Hd:

total Ca = moles(Di) + moles(Hd)
total Fe = moles(Hd)
total Mg = 2*moles(En) + moles(Di)

I don't have Morimoto 1988 to hand but it looks like the Wo apex is Ca2Si2O6, rather than CaSiO3 like the mineral, which means:

molfrac(Wo) = ( molfrac(Di) + molfrac(Hd) ) / 2

molfrac(Fs) = ( molfrac(Hd) ) / 2

More generally for cpx using the MELTS set of end members:

molfrac(Wo) = ( sum of molfracs(Di,  Hd, alumino-buffonite,  buffonite, essenite) ) / 2*( sum of molfracs(En, Di, Hd) )

molfrac(Fs) = ( molfrac(Hd) ) / 2*( sum of molfracs(En, Di, Hd) )

Check out GeoPlotters ( for more information. It's an excellent resource.

Hi all,

I am trying to calculate the molefractions of Ferrosilite (Fs) and Wollastonite (Wo) from the given mole fractions of Diopside (Di), Hedenbergite (Hd) and Clinoenstatite (Cli-En) so that I can plot the compositional variation of pyroxenes during fractional crystallization in the pyroxene quadrilateral diagram of Morimoto, 1988.

Like you have mentioned, I have used the relation: Fs = Cli-En  + Hd - 2Di and
similarly, I have used the relation: 2Di = 2Wo + Cli-En to calculate the mole fraction of Wollastonite, but ended up with a negative value for Wo. As pyroxene quadrilateral is a ternary diagram, negative values do not plot within the ternary.

Is there a mistake in my approach? Any leads would be highly appreciated.
No worries Paula.  :)
Hi Simon,

No it's not properly implemented yet, sorry. I'm a bit snowed with other projects (and pandemic life!) at the moment, so not sure how soon I'll get back to working on this. Definitely intending to have it working by the Goldschmidt workshop, which will be virtual at the end of June. Will keep you posted if there's anything to test sooner.

Pages: 1 2 3 4 5 6 7 8 9 10