I would like to do this as well !
I have a work around. Just use fprintf to write the PHREEQC input file and then use ! to run PHREEQC and inport the results.
this makes a simple pe pH diagram for hydrous ferric oxide. I modified a local version of the PHREEQC database to have constant pH and pe "phases". So this won't run as I've written it for you I don't think but it gives you an idea of how to call PHREEQC from matlab.
I find it is pretty fast. Code also work on Windows but much much slower. I use iPHREEQC when I am on my windows computer.
anyway, here is the code example ...
% make pH pe diagram for iron salts in water
function makepHpe_simpleFesystem
figure(1); clf
pH=2:1:12; pe=-15:2:22; FeT=1e-5;
% only perform calculations if within the stability field of water
for i=1:size(pH,2)
for j=1:size(pe,2)
HFO(i,j)=NaN;
if pe(j)>=-pH(i)
if pe(j)<=21.6-pH(i)
HFO(i,j)=returnFespeciation(FeT,pH(i),pe(j));
end
end
end
end
surf(pH,pe,HFO'./FeT)
set(gca,'Clim',[0 1])
set(gca,'linewidth',2,'fontsize',10)
axis([min(pH) max(pH) min(pe) max(pe) 0 1.1]); view([0 90]);
hold on
plot(pH,-pH,'k','linewidth',2)
plot(pH,21.6-pH,'k','linewidth',2)
title('HFO PHREEQC database')
xlabel('pH'); ylabel('pe')
colorbar
end
function [HFO]=returnFespeciation(FeT,pH_val,pe_val)
temp_val=25;
fileID=fopen('runphreeqc.txt','w');
fprintf(fileID,'SOLUTION 1 Pure water\n');
fclose(fileID);
fileID=fopen('runphreeqc.txt','a');
pHtxt=['pH ',num2str(pH_val),'\n']; fprintf(fileID,pHtxt); % the initial pH condition
petxt=['pe ',num2str(pe_val),'\n']; fprintf(fileID,petxt); % the initial pe condition
temptxt=['temp ',num2str(temp_val),'\n']; fprintf(fileID,temptxt); % the temperature
fprintf(fileID,'-units mol/L\n'); % the unit of the input; usually mol/L is used
FeTtxt=['Fe(+3) ' num2str(FeT),'\n']; fprintf(fileID,FeTtxt); % total Fe
ClTtxt=['Cl ' num2str(FeT*3),'\n']; fprintf(fileID,ClTtxt); % total Cl (assume adding FeCl3)
fprintf(fileID,'EQUILIBRIUM_PHASES 1\n');
fprintf(fileID,' Fe(OH)3(a) 0.0 10 precipitate_only\n'); % Fe(OH)3 will precipitate out when the saturaiton index is greater than 0
fixpHtxt=[' pH_Fix -',num2str(pH_val),' HCl 10.0\n']; fprintf(fileID,fixpHtxt);
fprintf(fileID,'-force_equality true\n');
fixpetxt=[' pe_Fix ',num2str(-1*pe_val),' O2 \n']; fprintf(fileID,fixpetxt);
fprintf(fileID,'-force_equality true\n');
fprintf(fileID,'SELECTED_OUTPUT\n');
fprintf(fileID,'-file selected.out\n');
fprintf(fileID,'-selected_out true\n');
fprintf(fileID,'-user_punch true\n');
fprintf(fileID,'-reset false\n');
fprintf(fileID,'-simulation false\n');
fprintf(fileID,'-state false\n');
fprintf(fileID,'-distance false\n');
fprintf(fileID,'-time false\n');
fprintf(fileID,'-step false\n');
fprintf(fileID,'-ph true\n');
fprintf(fileID,'-pe true\n');
fprintf(fileID,'-reaction false\n');
fprintf(fileID,'-temperature false\n');
fprintf(fileID,'-alkalinity false\n');
fprintf(fileID,'-ionic_strength false\n');
fprintf(fileID,'-water false\n');
fprintf(fileID,'-charge_balance false\n');
fprintf(fileID,'-percent_error false\n');
fprintf(fileID,'-equilibrium_phases Fe(OH)3(a)');
fclose(fileID);
%in LINUX
!phreeqc runphreeqc.txt out.txt PHREEQC.dat
%can see output of that one - so track error messages
%str=['!phreeqc runphreeqc.txt out.txt PHREEQC.dat']; evalc(str); % so no screen output
str=['!tail -n +2 "selected.out" > "selected.tmp" && mv "selected.tmp" "selected.out"']; evalc(str); load selected.out; % way faster only works in linux
%t= importdata('selected.out', '\t', 1); selected=t.data; % works in linux or windows but slower
selected=selected(2,:);
pH=selected(1)
pe=selected(2)
HFO=selected(4)
%check pH was actually fixed. it is a little off without pH fixed phase
pHtst=pH-pH_val
petst=pe-pe_val
end
