Convergence and com errors in Phreeqc python coupling

[process]

Hi all,
I am studying CaP phases precipitation in different pHs and temperatures. I created a nested loop for pH and temperature and also defined my other parameters. I have tried different ways to solve the problem but to no avail. Could someone help me out? I have pasted and attached the file as well.

   
# Import standard library modules first.
import os
# Then get third party modules.
from win32com.client import Dispatch
#import matplotlib.pyplot as plt
import numpy as np
from math import exp,sqrt,log10
import scipy.optimize as optimize
import xlsxwriter

def selected_array(db_path, input_string):
    """Load database via COM and run input string."""
    dbase = Dispatch('IPhreeqcCOM.Object')
    dbase.LoadDatabase(db_path)
    dbase.RunString(input_string)
    return dbase.GetSelectedOutputArray()

def phreecalc_pitz(input_string):
    """Get results from PHREEQC"""
    pitzer_result  = selected_array('pitzer.dat', input_string)
    return pitzer_result

def phreecalc_minteq(input_string):
    """Get results from PHREEQC"""
    result  = selected_array('minteq.v4.dat', input_string)
    return result

pH=np.linspace(2,13,7)
temp=list(reversed(np.linspace(25,40,5)))

P=np.zeros(len(pH));Alkalinity=np.zeros(len(pH))
Ca=np.zeros(len(pH));Mg=np.zeros(len(pH))
CaCO3=np.zeros(len(pH));CaH2PO4=np.zeros(len(pH));CaHPO4=np.zeros(len(pH))
CaPO4=np.zeros(len(pH));si_Amor_CalciumPhosphate=np.zeros(len(pH));si_Ca3PO42beta=np.zeros(len(pH))
si_CaHPO4=np.zeros(len(pH));si_Ca4HPO433H2O=np.zeros(len(pH));si_Hydroxylapatite=np.zeros(len(pH))
Monetite=np.zeros(len(pH));Struvite=np.zeros(len(pH));Ca3PO42beta=np.zeros(len(pH))
Amor_CalciumPhosphate=np.zeros(len(pH));Ca4HPO433H2O=np.zeros(len(pH));CaHPO4=np.zeros(len(pH))
CaHPO42H2O=np.zeros(len(pH));Calcite=np.zeros(len(pH));Hydroxylapatite=np.zeros(len(pH))




for i in range(len(pH-1)):
    print pH
    for i in range(len(temp)):
        print (pH,temp)
    input_string = """
PHASES
Struvite
    MgNH4PO4(H2O)6 = 6H2O + Mg+2 + NH4+ + PO4-3
    log_k     -13.26
    delta_h   22.6 kJ
Fix_H+
    H+ = H+
    log_k     0
MgSO4:7H2O
    MgSO4(H2O)7 = 7H2O + Mg+2 + SO4-2
    log_k     -2.14
    delta_h   11.81 kJ
Mg(OH)2
    Mg(OH)2 = Mg+2 + 2OH-
    log_k     16.84
    delta_h   -113.46 kJ
Amor_CalciumPhosphate
    Ca3(PO4)2H2O = 3Ca+2 + H2O + 2PO4-3
    log_k     -25.46
    delta_h   -440 kJ
Monetite
    CaHPO4 = Ca+2 + HPO4-2
    log_k     -6.81
CaCl2
    CaCl2 = Ca+2 + 2Cl-
    log_k     100
H2O
    H2O = H2O
    log_k     0
NH4+
    NH4+ = H+ + NH3
    log_k     100
HPO4-2
    HPO4-2 = H+ + PO4-3
    log_k     100
PO4-3
    PO4-3 = PO4-3
    log_k     100
NaOH
    NaOH=Na+ + OH-
    log_k     100
SOLUTION 1
    temp      %f
    pH        %f
    pe        4
    redox     pe
    units     mmol/l
    density   1
    C         5
    Ca        20
    Cl        9.3
    K         0.6
    Mg        1.15
    Na        8.7
    P         4
    S         1.15
    -water    1 # kg
EQUILIBRIUM_PHASES 1
    Amor_CalciumPhosphate 0 0
    Ca3(PO4)2(beta) 0 0
    Ca4H(PO4)3:3H2O 0 0
    CaHPO4    0 0
    CaHPO4:2H2O 0 0
    Calcite   0 0
    Fix_H+    -7.2 NaOH      1
    Hydroxylapatite 0 0
    Monetite  0 0
SELECTED_OUTPUT 1
    -file                 F:\Ben Gurion University\PhD Programme\PHREEQC\Aqueuos Chemistry Modeling Course\Python Phreeqc\CaP Ageing conditions_output.xls
    -reset                false
    -simulation           false
    -pH                   false
    -temperature          false
    -temperature          false
    -alkalinity           false
    -totals               P  Ca  Mg  Alkalinity
    -molalities           CaCO3  CaH2PO4+  CaHPO4  CaPO4-
                          H3PO4  PO4-3
    -activities           CaCO3  CaH2PO4+  CaHPO4  CaPO4-
                          PO4-3  MgCO3
    -equilibrium_phases   Amor_CalciumPhosphate  Ca3(PO4)2(beta)  Ca4H(PO4)3:3H2O  CaHPO4
                          CaHPO4:2H2O  Calcite  Hydroxylapatite  Monetite
                          Struvite
    -saturation_indices   Hydroxylapatite  Ca3(PO4)2(beta)  Ca4H(PO4)3:3H2O  CaHPO4
                          CaHPO4:2H2O  Amor_CalciumPhosphate  Monetite  Hydroxylapatite
                          Struvite
    -active               true
    -user_punch           false

         END"""%(pH,temp)
           
    sol=phreecalc_minteq(input_string)
   
   
    P= sol[2][0]; Alkalinity= sol[2][1]; Ca= sol[2][2]; Mg= sol[2][3]
    CaCO3= sol[2][4]; CaH2PO4= sol[2][5]; CaHPO4= sol[2][6]; CaPO4= sol[2][7]
    Amor_CalciumPhosphate= sol[2][33]; si_Ca3PO42beta= sol[2][29]
    CaHPO4= sol[2][31]; Ca4HPO433H2O= sol[2][30]; si_Hydroxylapatite= sol[2][28]
    Monetite= sol[2][34]; Struvite= sol[2][35]; Ca3PO42beta= sol[2][13]
    Amor_CalciumPhosphate= sol[2][11]; Ca4HPO433H2O= sol[2][15]; CaHPO4= sol[2][17]
    CaHPO42H2O= sol[2][19]; Calcite= sol[2][21];  Hydroxylapatite= sol[2][23]
    print pH
   
   
"""writing results to excell"""
# Create a workbook and add a worksheet.
workbook = xlsxwriter.Workbook('CaP phase ageing.xlsx')
worksheet = workbook.add_worksheet()

#write titles
worksheet.write(0, 0, 'pH')
worksheet.write(0, 1, 'P')
worksheet.write(0, 2, 'Alkalinity')
worksheet.write(0, 3, 'Ca')
worksheet.write(0, 4, 'Mg')
worksheet.write(0, 5, 'm_CaCO3')
worksheet.write(0, 6, 'm_CaH2PO4+')
worksheet.write(0, 7, 'm_CaHPO4')
worksheet.write(0, 8, 'm_CaPO4-')
worksheet.write(0, 9, 'si_Amor_CalciumPhosphate')
worksheet.write(0, 10, 'si_Ca3PO42beta')
worksheet.write(0, 11, 'si_CaHPO4')
worksheet.write(0, 12, 'si_Ca4H(PO4)3:3H2O')
worksheet.write(0, 13, 'si_Hydroxylapatite')
worksheet.write(0, 14, 'si_Monetite')
worksheet.write(0, 15, 'si_Struvite')
worksheet.write(0, 16, 'E_Ca3(PO4)2(beta)')
worksheet.write(0, 17, 'E_Amor_CalciumPhosphate')
worksheet.write(0, 18, 'E_Ca4H(PO4)3:3H2O')
worksheet.write(0, 19, 'E_CaHPO4')
worksheet.write(0, 20, 'E_CaHPO4:2H2O')
worksheet.write(0, 21, 'E_Calcite')
worksheet.write(0, 22, 'E_Hydroxylapatite')

# Start from the first cell. Rows and columns are zero indexed.
row = 1

# Iterate over the data and write it out row by row.
while row <= len(pH):
   
    worksheet.write(row, 0, pH[row-1])
   
    worksheet.write(row, 1, P[row-1])
   
    worksheet.write(row, 2, Alkalinity[row-1])
   
    worksheet.write(row, 3, Ca[row-1])
   
    worksheet.write(row, 4, Mg[row-1])
   
    worksheet.write(row, 5, CaCO3[row-1])
   
    worksheet.write(row, 6, CaH2PO4[row-1])
   
    worksheet.write(row, 7, CaHPO4[row-1])
   
    worksheet.write(row, 8, CaPO4[row-1])
   
    worksheet.write(row, 9, Amor_CalciumPhosphate[row-1])
   
    worksheet.write(row, 10, si_Ca3PO42beta[row-1])
   
    worksheet.write(row, 11, CaHPO4[row-1])
   
    worksheet.write(row, 12, Ca4HPO433H2O[row-1])
   
    worksheet.write(row, 13, si_Hydroxylapatite[row-1])
   
    worksheet.write(row, 14, Monetite[row-1])
   
    worksheet.write(row, 15, Struvite[row-1])
   
    worksheet.write(row, 16, Ca3PO42beta[row-1])
   
    worksheet.write(row, 17, Amor_CalciumPhosphate[row-1])
   
    worksheet.write(row, 18, Ca4HPO433H2O[row-1])
   
    worksheet.write(row, 19, CaHPO4[row-1])
   
    worksheet.write(row, 20, CaHPO42H2O[row-1])
   
    worksheet.write(row, 21, Calcite[row-1])
   
    worksheet.write(row, 22, Hydroxylapatite[row-1])
   
       
    row=row+1
   
workbook.close()   


|Errors are given below:

File "<COMObject IPhreeqcCOM.Object>", line 2, in RunString

com_error: (-2147352567, 'Exception occurred.', (0, u'IPhreeqcCOM.Object.4', u'ERROR:                    C has not converged. Total: 5.009604e-003\tCalculated: 1.083580e-004\tResidual: 4.901246e-003\n\nERROR:                   Ca has not converged. Total: 2.003842e-002\tCalculated: 1.981032e-002\tResidual: 2.280914e-004\n\nERROR:                   Cl has not converged. Total: 9.317863e-003\tCalculated: 0.000000e+000\tResidual: 9.317863e-003\n\nERROR:                    K has not converged. Total: 6.011525e-004\tCalculated: 4.726239e-007\tResidual: 6.006799e-004\n\nERROR:                   Mg has not converged. Total: 1.152209e-003\tCalculated: 1.139093e-003\tResidual: 1.311548e-005\n\nERROR:                   Na has not converged. Total: 8.716711e-003\tCalculated: 1.107553e-004\tResidual: 8.605956e-003\n\nERROR:                    P has not converged. Total: 4.007683e-003\tCalculated: 2.869977e-006\tResidual: 4.004813e-003\n\nERROR:                    S has not converged. Total: 1.152209e-003\tCalculated: 6.000000e+003\tResidual: -5.999999e+003\n\nERROR:                   Mu Ionic strength has not converged. \tResidual: 9.993255e+010\n\nERROR:               A(H2O) Activity of water has not converged. \tResidual: 1.189890e+002\n\nERROR: Model failed to converge for initial solution  1.\n', None, 0, -2147352567), None)

[dlparkhurst]

I think you have set pH to 25.

[process]

Hi David,

Below are my settings for pH and temp:
pH=np.linspace(2,13,7)
temp=list(reversed(np.linspace(25,40,5)))

Also, under solutions PH = %f .

I cannot seem to find where pH is set at 25.

[dlparkhurst]

I'm not that fluent in python, but it looks to me that the first argument is the temperature and the second is the pH in the following:

SOLUTION 1
    temp      %f
    pH        %f

But the first argument is pH and the second is temperature here:

END"""%(pH,temp)

You should print your input_string. If the string looks correct and the error persists, then include it in your next posting.

Note that the error is in an initial solution calculation, that is the error is not related to equilibrium phases, the reaction calculation, or selected output. I did not try every combination of pH and temperature, but I tried all of the extremes with no problem in the initial solution calculation. However, a pH of 25 gives the error that you show.

[process]

Hi David,

Thanks so much. I made the change at the END to be the same sequence at the beginning:
temp      %f
    pH        %f
END"""%(temp[t],pH)
I wasn't getting the convergence errors.
I know included the temperature to be true (argument) in the input string. However, I get 0.00 on the excel file. It does write the temperature I set in the solution. Below is the input string(attached is the file).

 '\n        SOLUTION 1\n    temp      0.000000\n    pH        14.000000\n    pe        4\n    redox     pe\n    units     mmol/l\n    density   1\n    C         5\n    Ca        20\n    Cl        9.3\n    K         0.6\n    Mg        1.15\n    Na        8.7\n    P         4\n    S         1.15\n    -water    1 # kg\n    PHASES\n    Struvite\n    MgNH4PO4(H2O)6 = 6H2O + Mg+2 + NH4+ + PO4-3\n    log_k     -13.26\n    delta_h   22.6 kJ\n    Fix_H+\n    H+ = H+\n    log_k     0\n    MgSO4:7H2O\n    MgSO4(H2O)7 = 7H2O + Mg+2 + SO4-2\n    log_k     -2.14\n    delta_h   11.81 kJ\n    Mg(OH)2\n    Mg(OH)2 = Mg+2 + 2OH-\n    log_k     16.84\n    delta_h   -113.46 kJ\n    Amor_CalciumPhosphate\n    Ca3(PO4)2H2O = 3Ca+2 + H2O + 2PO4-3\n    log_k     -25.46\n    delta_h   -440 kJ\n    Monetite\n    CaHPO4 = Ca+2 + HPO4-2\n    log_k     -6.81\n    CaCl2\n    CaCl2 = Ca+2 + 2Cl-\n    log_k     100\n    H2O\n    H2O = H2O\n    log_k     0\n    NH4+\n    NH4+ = H+ + NH3\n    log_k     100\n    HPO4-2\n    HPO4-2 = H+ + PO4-3\n    log_k     100\n    PO4-3\n    PO4-3 = PO4-3\n    log_k     100\n    EQUILIBRIUM_PHASES 1\n    Amor_CalciumPhosphate 0 0\n    Ca3(PO4)2(beta) 0 0\n    Ca4H(PO4)3:3H2O 0 0\n    CaHPO4    0 0\n    CaHPO4:2H2O 0 0\n    Calcite   0 0\n    Fix_H+    -7.2 NaOH      1\n    Hydroxylapatite 0 0\n    Monetite  0 0\n        SELECTED_OUTPUT 1\n    -file                 F:\\Ben Gurion University\\PhD Programme\\PHREEQC\\Aqueuos Chemistry Modeling Course\\Python Phreeqc\\CaP Ageing conditions.xls\n    -reset                false\n    -temperature          true\n    -totals               P  Alkalinity  Ca  Mg\n    -molalities           CaCO3  CaH2PO4+  CaHPO4  CaPO4-\n                          H3PO4  PO4-3\n    -equilibrium_phases   Amor_CalciumPhosphate  Ca3(PO4)2(beta)  Ca4H(PO4)3:3H2O  CaHPO4\n                          CaHPO4:2H2O  Calcite  Hydroxylapatite  Monetite\n                          Struvite\n    -saturation_indices   Hydroxylapatite  Ca3(PO4)2(beta)  Ca4H(PO4)3:3H2O  CaHPO4\n                          CaHPO4:2H2O  Amor_CalciumPhosphate  Monetite  Hydroxylapatite\n                          Struvite  Calcite\n    -active               true\n    \n         END

[dlparkhurst]

I don't understand your problem. You set the temperature to be 0.0 in the SOLUTION 1 definition ('\n        SOLUTION 1\n    temp      0.000000\n    pH        14.000000\n). Seems like you should get 0.0 in your Excel file.

There will be two calculations performed for each RunString, an initial solution calculation (row 1 of the output) and a reaction calculation where the solution reacts with EQUILIBRIUM_PHASES (row 2). You are setting varying temperatures and pH in your SOLUTION definition. The temperature you set in SOLUTION will be used by default in the reaction calculation. (Alternatively, you can set the temperature of the reaction calculation with REACTION_TEMPERATURE.) However, you are specifying a pH of 7.2 in the reaction calculation by specifying Fix_H+ -7.2. So, no matter the starting pH, the reaction pH is intended to be fixed. Is this what you wanted to do? It seems more likely that you wanted the pH of the reaction solution to range over a set of values; that is, varying the value in Fix_H+.

If you start at a pH of 14, you will not be able to attain a pH of 7.2 by adding or removing NaOH. You can only remove as much Na as is in the SOLUTION definition, which probably is not enough to get down to pH 7.2. You will need to add an acid. However, if you start at pH 2, you will indeed need to add a base. Assuming you always want pH 7.2, there is a trick shown in a footnote of example 8 to add either acid or base as needed.

Finally, a simplification would be to first define PHASES and SELECTED_OUTPUT in a string and run it with RunString. These would then remain defined for subsequent simulations. So, you could loop through definitions of SOLUTION and EQUILIBRIUM_PHASES. If the SOLUTION definition is constant, you could add it to the definition of PHASES and SELECTED_OUTPUT, in which case you only define it once. Then you would define each simulation for example as USE solution 1; REACTION_TEMPERATURE; ...;EQUILIBRIUM_PHASES 1; ... Fix_H+ x.x NaOH 10; ...