Running 'stack_mosaic.R' [178s]
  Running 'stack_solubility.R' [163s]
  Running 'tinytest.R' [27s]
Running the tests in 'tests/stack_mosaic.R' failed.
Complete output:
  > # Test variations on stacked mosaic diagrams culminating with stack_mosaic() wrapper function
  > # 20220722 v1 adapted from Mosaic Stacking 2 from multi-metal.Rmd
  > 
  > library(CHNOSZ)
  CHNOSZ version 2.1.0 (2024-02-11)
  reset: creating "thermo" object
  OBIGT: loading default database with 1940 aqueous, 3486 total species
  > 
  > # Use low resolution for routine package checking
  > res <- 20
  > ## Increase resolution for nicer-looking PDF (the one kept in the test directory)
  > #res <- 200
  > # Define system
  > pH <- c(0, 14, res)
  > O2 <- c(-48, -33, res)
  > T <- 200
  > logmS <- -2
  > m_NaCl <- 0.1
  > logm_aq <- -6 # for both Fe- and Cu-bearing aq species
  > # Define basis species
  > S.aq <- c("H2S", "HS-", "HSO4-", "SO4-2")
  > # Define minerals
  > Fe.cr <- c("pyrite", "pyrrhotite", "magnetite", "hematite")
  > Fe.abbrv <- c("Py", "Po", "Mag", "Hem")
  > FeCu.cr <- c("chalcopyrite", "bornite")
  > Cu.cr <- c("copper", "cuprite", "tenorite", "chalcocite", "covellite")
  > FeCu.abbrv <- c("Ccp", "Bn", "Cu", "Cpr", "Tnr", "Cct", "Cv")
  > # Define aqueous species
  > iFe.aq <- retrieve("Fe", c("S", "O", "H", "Cl"), "aq")
  > Fe.aq <- info(iFe.aq)$name
  check.GHS: calculated ΔG°f of HFeO2(aq) differs by -4600 cal mol-1 from database value
  check.GHS: calculated ΔG°f of FeO2-(aq) differs by -4641 cal mol-1 from database value
  > iCu.aq <- retrieve("Cu", c("S", "O", "H", "Cl"), "aq")
  > Cu.aq <- info(iCu.aq)$name
  check.EOS: calculated Cp° of Cu+2(aq) differs by 3.62 cal K-1 mol-1 from database value
  check.EOS: calculated V° of CuCl4-2(aq) differs by -1.07 cm3 mol-1 from database value
  check.EOS: calculated V° of CuO2-2(aq) differs by -1.09 cm3 mol-1 from database value
  check.EOS: calculated Cp° of CuHS(aq) differs by 7.35 cal K-1 mol-1 from database value
  check.EOS: calculated V° of CuHS(aq) differs by -16.2 cm3 mol-1 from database value
  > nacl <- NaCl(T = T, P = "Psat", m_tot = m_NaCl)
  > 
  > setup <- function() {
  +   reset()
  +   # Setup basis species
  +   basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+", "Cl-"))
  +   basis("H2S", logmS)
  +   basis("Cl-", log10(nacl$m_Cl))
  + }
  > 
  > ref_Cu_craq <- function() {
  +   # Load Cu-bearing minerals
  +   species(Cu.cr)
  +   # Add aqueous species 20210220
  +   species(iCu.aq, logm_aq, add = TRUE)
  + 
  +   mCu <- mosaic(list(S.aq), pH = pH, O2 = O2, T = T, IS = nacl$IS)
  +   diagram(mCu$A.species)
  + }
  > 
  > ref_FeCu_cr <- function() {
  +   # Load Fe-bearing minerals
  +   species(Fe.cr)
  +   mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
  +   dFe <- diagram(mFe$A.species, lwd = 0, names = FALSE, plot.it = FALSE)
  + 
  +   # Load Cu-bearing minerals
  +   species(c(FeCu.cr, Cu.cr))
  +   # Mosaic with all Fe species as basis species
  +   mFeCu <- mosaic(list(S.aq, Fe.cr), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant))
  + 
  +   diagram(mFeCu$A.species)
  + }
  > 
  > setup_FeCu <- function() {
  +   # Load Fe-bearing minerals
  +   species(Fe.cr)
  +   # Add aqueous species 20210220
  +   species(iFe.aq, logm_aq, add = TRUE)
  +   mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
  +   dFe <- diagram(mFe$A.species, lwd = 0, names = FALSE, plot.it = FALSE)
  + 
  +   # Load Cu-bearing minerals
  +   species(c(FeCu.cr, Cu.cr))
  +   # Add aqueous species 20210220
  +   species(iCu.aq, logm_aq, add = TRUE)
  +   
  +   # Return Fe diagram for stacked mosaic calculations
  +   dFe
  + }
  > 
  > 
  > test_FeCu_old <- function() {
  +   dFe <- setup_FeCu()
  + 
  +   # NOTE: limitation in mosaic() when using both solid and aqueous basis species (i.e. c(Fe.cr, Fe.aq)):
  +   #   Only the activity of the first-defined basis species is used for all basis species.
  +   #   The first-defined basis species is pyrite (logact = 0), but logact < 0 for the aq Fe species.
  +   # Workaround: Adjust standard Gibbs energies of aq Fe species to virtually change their activities
  +   DG_J <- convert(-logm_aq, "G", T = convert(T, "K"))
  +   # We should use calories here because the database values are in calories 20220604
  +   stopifnot(all(info(iFe.aq)$E_units == "cal"))
  +   DG_cal <- convert(DG_J, "cal")
  +   G.orig <- info(iFe.aq)$G
  +   G.new <- G.orig + DG_cal
  +   mod.OBIGT(iFe.aq, G = G.new)
  + 
  +   # Mosaic with all Fe species as basis species
  +   mFeCu <- mosaic(list(S.aq, c(Fe.cr, Fe.aq)), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant))
  + 
  +   diagram(mFeCu$A.species)
  + }
  > 
  > test_FeCu_new <- function() {
  +   dFe <- setup_FeCu()
  + 
  +   # Mosaic with all Fe species as basis species
  +   # Use loga_aq argument to control the activity of aqueous species in mosaic calculation 20220722
  +   # c(NA, logm_aq) means to use:
  +   #   basis()'s value for logact of aqueous S species
  +   #   logm_aq for logact of aqueous Fe species
  +   mFeCu <- mosaic(list(S.aq, c(Fe.cr, Fe.aq)), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant), loga_aq = c(NA, logm_aq))
  + 
  +   diagram(mFeCu$A.species)
  + }
  > 
  > test_FeCu_stack <- function() {
  +   # Aqueous S species
  +   bases <- S.aq
  +   # Fe-bearing minerals
  +   species1 <- Fe.cr
  +   # Cu-bearing and FeCu-bearing minerals
  +   species2 <- Cu.cr
  +   species12 <- FeCu.cr
  + 
  +   # Fe-bearing aqueous species
  +   species1 <- c(species1, Fe.aq)
  +   # Cu-bearing aqueous species
  +   species2 <- c(species2, Cu.aq)
  + 
  +   sm <- stack_mosaic(bases, species1, species2, species12, pH = pH, O2 = O2, T = T, IS = nacl$IS, loga_aq = logm_aq, plot.it = FALSE)
  +   diagram(sm[[2]])
  + }
  > 
  > ## Setup plot
  > pdf("stack_mosaic.pdf", width = 8, height = 9)
  > mat <- matrix(c(1,6,2, 6,3,6, 4,6,5), nrow = 3, byrow = TRUE)
  > layout(mat, widths = c(2, 3, 2), heights = c(2, 3, 2, 3))
  > par(mar = c(4, 4, 3, 1))
  > 
  > ## Run tests
  > # This first one is a simple mosaic diagram (no stacking)
  > setup()
  reset: resetting "thermo" object
  OBIGT: loading default database with 1940 aqueous, 3486 total species
       Cl Cu Fe H O S  Z ispecies    logact state
  Cu+   0  1  0 0 0 0  1      801  0.000000    aq
  FeS2  0  0  1 0 0 2  0     2082  0.000000    cr
  H2S   0  0  0 2 0 1  0       65 -2.000000    aq
  O2    0  0  0 0 2 0  0     2679  0.000000   gas
  H2O   0  0  0 2 1 0  0        1  0.000000   liq
  H+    0  0  0 1 0 0  1        3  0.000000    aq
  Cl-   1  0  0 0 0 0 -1       29 -1.022401    aq
  > ref_Cu_craq()
  info.character: found H2S(aq); also available in gas
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 33 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, Cu+2, CuCl3-2, CuCl+, CuCl3-, CuCl4-2, CuOH+, HCuO2-, CuO2-2, Cu+, Cu(OH)2-, Cu(HS)2-, CuCl2-, HS-, HSO4-, SO4-2 (Bdot equation)
  nonideal: calculations for H2S, CuCl2, CuO, CuOH, CuHS, CuCl, H2S (Setchenow equation)
  subcrt: 3 polymorphs for chalcocite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of species for all 4 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of Cu+ in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  > title("Cu\ncr+aq", font.main = 1)
  > 
  > setup()
  reset: resetting "thermo" object
  OBIGT: loading default database with 1940 aqueous, 3486 total species
       Cl Cu Fe H O S  Z ispecies    logact state
  Cu+   0  1  0 0 0 0  1      801  0.000000    aq
  FeS2  0  0  1 0 0 2  0     2082  0.000000    cr
  H2S   0  0  0 2 0 1  0       65 -2.000000    aq
  O2    0  0  0 0 2 0  0     2679  0.000000   gas
  H2O   0  0  0 2 1 0  0        1  0.000000   liq
  H+    0  0  0 1 0 0  1        3  0.000000    aq
  Cl-   1  0  0 0 0 0 -1       29 -1.022401    aq
  > ref_FeCu_cr()
  info.character: found H2S(aq); also available in gas
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 15 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, HS-, HSO4-, SO4-2 (Bdot equation)
  nonideal: calculations for H2S, H2S (Setchenow equation)
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of species for all 4 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of FeS2 in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  info.character: found H2S(aq); also available in gas
  info.character: found pyrrhotite(cr) with 2 polymorphic transitions
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 22 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, HS-, HSO4-, SO4-2 (Bdot equation)
  nonideal: calculations for H2S, H2S (Setchenow equation)
  subcrt: 3 polymorphs for chalcopyrite ... polymorph 1 is stable
  subcrt: 3 polymorphs for bornite ... polymorph 1 is stable
  subcrt: 3 polymorphs for chalcocite ... polymorph 2 is stable
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of basis species group 2: pyrite pyrrhotite magnetite hematite
  mosaic: calculating affinities of species for all 16 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of Cu+ in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  > title("Fe-Cu\ncr", font.main = 1)
  > 
  > T1result <- "(TODO: looks ok)"
  > 
  > setup()
  reset: resetting "thermo" object
  OBIGT: loading default database with 1940 aqueous, 3486 total species
       Cl Cu Fe H O S  Z ispecies    logact state
  Cu+   0  1  0 0 0 0  1      801  0.000000    aq
  FeS2  0  0  1 0 0 2  0     2082  0.000000    cr
  H2S   0  0  0 2 0 1  0       65 -2.000000    aq
  O2    0  0  0 0 2 0  0     2679  0.000000   gas
  H2O   0  0  0 2 1 0  0        1  0.000000   liq
  H+    0  0  0 1 0 0  1        3  0.000000    aq
  Cl-   1  0  0 0 0 0 -1       29 -1.022401    aq
  > d_FeCu_old <- test_FeCu_old()
  info.character: found H2S(aq); also available in gas
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 28 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, FeCl+, FeCl+2, Fe+2, Fe+3, FeOH+, HFeO2-, FeOH+2, FeO+, FeO2-, HS-, HSO4-, SO4-2 (Bdot equation)
  nonideal: calculations for H2S, FeCl2, FeO, HFeO2, FeSO4, H2S (Setchenow equation)
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of species for all 4 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of FeS2 in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  check.GHS: calculated ΔG°f of HFeO2(aq) differs by -4600 cal mol-1 from database value
  check.GHS: calculated ΔG°f of FeO2-(aq) differs by -4641 cal mol-1 from database value
  check.GHS: calculated ΔG°f of HFeO2(aq) differs by -4600 cal mol-1 from database value
  check.GHS: calculated ΔG°f of FeO2-(aq) differs by -4641 cal mol-1 from database value
  mod.OBIGT: updated FeCl+(aq)
  mod.OBIGT: updated FeCl2(aq)
  mod.OBIGT: updated FeCl+2(aq)
  mod.OBIGT: updated Fe+2(aq)
  mod.OBIGT: updated Fe+3(aq)
  mod.OBIGT: updated FeOH+(aq)
  mod.OBIGT: updated FeO(aq)
  mod.OBIGT: updated HFeO2-(aq)
  mod.OBIGT: updated FeOH+2(aq)
  mod.OBIGT: updated FeO+(aq)
  mod.OBIGT: updated HFeO2(aq)
  mod.OBIGT: updated FeO2-(aq)
  mod.OBIGT: updated FeSO4(aq)
  info.character: found H2S(aq); also available in gas
  info.character: found pyrrhotite(cr) with 2 polymorphic transitions
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 43 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, Cu+2, CuCl3-2, CuCl+, CuCl3-, CuCl4-2, CuOH+, HCuO2-, CuO2-2, Cu+, Cu(OH)2-, Cu(HS)2-, CuCl2-, HS-, HSO4-, SO4-2, Fe+2, HFeO2-, FeO2- (Bdot equation)
  nonideal: calculations for H2S, CuCl2, CuO, CuOH, CuHS, CuCl, H2S, FeSO4 (Setchenow equation)
  subcrt: 3 polymorphs for chalcopyrite ... polymorph 1 is stable
  subcrt: 3 polymorphs for bornite ... polymorph 1 is stable
  subcrt: 3 polymorphs for chalcocite ... polymorph 2 is stable
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of basis species group 2: pyrite pyrrhotite magnetite hematite Fe+2 HFeO2- FeO2- FeSO4
  mosaic: calculating affinities of species for all 32 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of Cu+ in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  > title("Fe-Cu\ncr+aq", font.main = 1)
  > 
  > T2result <- "(TODO: looks ok)"
  > T3result <- "(TODO: looks different)"
  > 
  > setup()
  reset: resetting "thermo" object
  OBIGT: loading default database with 1940 aqueous, 3486 total species
       Cl Cu Fe H O S  Z ispecies    logact state
  Cu+   0  1  0 0 0 0  1      801  0.000000    aq
  FeS2  0  0  1 0 0 2  0     2082  0.000000    cr
  H2S   0  0  0 2 0 1  0       65 -2.000000    aq
  O2    0  0  0 0 2 0  0     2679  0.000000   gas
  H2O   0  0  0 2 1 0  0        1  0.000000   liq
  H+    0  0  0 1 0 0  1        3  0.000000    aq
  Cl-   1  0  0 0 0 0 -1       29 -1.022401    aq
  > d_FeCu_new <- test_FeCu_new()
  info.character: found H2S(aq); also available in gas
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 28 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, FeCl+, FeCl+2, Fe+2, Fe+3, FeOH+, HFeO2-, FeOH+2, FeO+, FeO2-, HS-, HSO4-, SO4-2 (Bdot equation)
  nonideal: calculations for H2S, FeCl2, FeO, HFeO2, FeSO4, H2S (Setchenow equation)
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of species for all 4 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of FeS2 in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  info.character: found H2S(aq); also available in gas
  info.character: found pyrrhotite(cr) with 2 polymorphic transitions
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 43 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, Cu+2, CuCl3-2, CuCl+, CuCl3-, CuCl4-2, CuOH+, HCuO2-, CuO2-2, Cu+, Cu(OH)2-, Cu(HS)2-, CuCl2-, HS-, HSO4-, SO4-2, Fe+2, HFeO2-, FeO2- (Bdot equation)
  nonideal: calculations for H2S, CuCl2, CuO, CuOH, CuHS, CuCl, H2S, FeSO4 (Setchenow equation)
  subcrt: 3 polymorphs for chalcopyrite ... polymorph 1 is stable
  subcrt: 3 polymorphs for bornite ... polymorph 1 is stable
  subcrt: 3 polymorphs for chalcocite ... polymorph 2 is stable
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of basis species group 2: pyrite pyrrhotite magnetite hematite Fe+2 HFeO2- FeO2- FeSO4
  mosaic: calculating affinities of species for all 32 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of Cu+ in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  > title("Fe-Cu\ncr+aq", font.main = 1, xpd = NA)
  > 
  > if(identical(d_FeCu_old$predominant, d_FeCu_new$predominant)) T4result <- "(OK)" else T4result <- "(FAILED)"
  > 
  > setup()
  reset: resetting "thermo" object
  OBIGT: loading default database with 1940 aqueous, 3486 total species
       Cl Cu Fe H O S  Z ispecies    logact state
  Cu+   0  1  0 0 0 0  1      801  0.000000    aq
  FeS2  0  0  1 0 0 2  0     2082  0.000000    cr
  H2S   0  0  0 2 0 1  0       65 -2.000000    aq
  O2    0  0  0 0 2 0  0     2679  0.000000   gas
  H2O   0  0  0 2 1 0  0        1  0.000000   liq
  H+    0  0  0 1 0 0  1        3  0.000000    aq
  Cl-   1  0  0 0 0 0 -1       29 -1.022401    aq
  > d_FeCu_stack <- test_FeCu_stack()
  info.character: found H2S(aq); also available in gas
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 28 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, FeCl+, FeCl+2, Fe+2, Fe+3, FeOH+, HFeO2-, FeOH+2, FeO+, FeO2-, HS-, HSO4-, SO4-2 (Bdot equation)
  nonideal: calculations for H2S, FeCl2, FeO, HFeO2, FeSO4, H2S (Setchenow equation)
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of species for all 4 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of FeS2 in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  info.character: found H2S(aq); also available in gas
  info.character: found pyrrhotite(cr) with 2 polymorphic transitions
  affinity: temperature is 200 ºC
  affinity: pressure is Psat
  affinity: ionic strength is 0.0949726913845165
  affinity: variable 1 is pH at 20 values from 0 to 14
  affinity: variable 2 is log10(f_O2) at 20 values from -48 to -33
  subcrt: 43 species at 200 ºC and 15.54 bar (wet) [energy units: J]
  nonideal: calculations for Cu+, Cl-, Cu+2, CuCl3-2, CuCl+, CuCl3-, CuCl4-2, CuOH+, HCuO2-, CuO2-2, Cu+, Cu(OH)2-, Cu(HS)2-, CuCl2-, HS-, HSO4-, SO4-2, Fe+2, HFeO2-, FeO2- (Bdot equation)
  nonideal: calculations for H2S, CuCl2, CuO, CuOH, CuHS, CuCl, H2S, FeSO4 (Setchenow equation)
  subcrt: 3 polymorphs for chalcopyrite ... polymorph 1 is stable
  subcrt: 3 polymorphs for bornite ... polymorph 1 is stable
  subcrt: 3 polymorphs for chalcocite ... polymorph 2 is stable
  subcrt: 3 polymorphs for pyrrhotite ... polymorph 2 is stable
  mosaic: calculating affinities of basis species group 1: H2S HS- HSO4- SO4-2
  mosaic: calculating affinities of basis species group 2: pyrite pyrrhotite magnetite hematite Fe+2 HFeO2- FeO2- FeSO4
  mosaic: calculating affinities of species for all 32 combinations of the basis species
  balance: on moles of H2S in formation reactions
  equilibrate: n.balance is 1 1 1 1
  equilibrate: loga.balance is -2
  equilibrate: using boltzmann method
  balance: on moles of Cu+ in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  balance: on moles of Cu+ in formation reactions
  diagram: using maximum affinity method for 2-D diagram
  > title("Fe-Cu\ncr+aq", font.main = 1, xpd = NA)
  > 
  > if(identical(d_FeCu_new$predominant, d_FeCu_stack$predominant)) T5result <- "(OK)" else T5result <- "(FAILED)"
  > 
  > ## Add arrows and text
  > plot.new()
  > text(0.5, 0.99, "mosaic()\nWITHOUT 'loga_aq'", font = 2, xpd = NA)
  > 
  > arrows(0.35, 0.9, 0.65, 0.9)
  > text(0.5, 0.93, paste0("Test 1: Cu mineral reactions unaffected\n", T1result))
  > 
  > arrows(0.1, 0.72, 0.25, 0.65)
  > text(0.05, 0.62, paste0("Test 2: aq reactions unaffected\n", T2result), xpd = NA)
  > 
  > arrows(0.9, 0.72, 0.75, 0.65)
  > text(0.9, 0.62, paste0("Test 3: cr reactions unaffected\n", T3result), xpd = NA)
  > 
  > abline(h = 0.25, lty = 2, col = 8)
  > text(0.4, 0.23, "mosaic()\nWITH 'loga_aq'", font = 2)
  > lines(c(0.5, 0.5), c(0.15, 0.25), lty = 2, col = 8)
  > text(0.6, 0.23, "stack_mosaic()", font = 2)
  > 
  > arrows(0.25, 0.4, 0.1, 0.3)
  > text(0.05, 0.35, paste0("Test 4: identical diagram\n", T4result), xpd = NA)
  > 
  > arrows(0.35, 0.1, 0.65, 0.1)
  > text(0.5, 0.05, paste0("Test 5: identical diagram\n", T5result))
  > 
  > dev.off()
  null device 
            1 
  > 
  > proc.time()
     user  system elapsed 
    13.42    0.62   14.01 
Running the tests in 'tests/stack_solubility.R' failed.
Complete output:
  > # Make stacked mosaic diagrams and solubility contours
  > # 20220723 jmd
  > 
  > # - Is solubility() function basically correct?
  > #   Check: Solubility contour (right column for each stack) coincides with mineral-aqueous stability boundary
  > #   NOTE: Contour and aqueous species are both set to logarithm of activity = -6
  > 
  > # - Are solubilities of bimetallic minerals (i.e. chalcopyrite) correct?
  > #   Check: ???
  > #   TODO: Is there a reference for the correct result?
  > 
  > # - Is an Fe(first)-Cu stack the same as a Cu(first)-Fe stack?
  > #   Check: Chalcopyrite field has consistent extents in Stack 2 (Diagram 4) and Stack 4 (Diagram 8)
  > #   TODO: Chalcopyrite field doesn't change, but assemblages with other minerals are different
  > 
  > # - Are solubilities correct for different stoichiometries?
  > #   Check: Change chalcopyrite (Cu:Fe = 1:1) to bornite (Cu:Fe = 5:1)
  > #   TODO: Solubility contours for bornite are incorrect
  > 
  > # NOTE: Among these diagrams, I think that Diagrams 4 and 8 are likely to be
  > # the most accurate depictions of the solubilities of Cu and Fe, respectively
  > # (Bornite and chalcopyrite should both be present, after the solubility() calculations for bornite are fixed.)
  > 
  > # This is a long test: don't run it for routine package checking
  > if(FALSE) {
  + 
  + library(CHNOSZ)
  + reset()
  + 
  + pdf("stack_solubility.pdf", width = 6, height = 9)
  + par(mfrow = c(4, 2))
  + 
  + # Define system
  + res <- 200
  + pH <- c(0, 14, res)
  + O2 <- c(-48, -33, res)
  + T <- 200
  + logmS <- -2
  + m_NaCl <- 0.1
  + logm_aq <- -6 # for both Fe- and Cu-bearing aq species
  + # Define basis species
  + S.aq <- c("H2S", "HS-", "HSO4-", "SO4-2")
  + # Define minerals
  + Fe.cr <- c("pyrite", "pyrrhotite", "magnetite", "hematite")
  + FeCu.cr <- c("chalcopyrite", "bornite")
  + Cu.cr <- c("copper", "cuprite", "tenorite", "chalcocite", "covellite")
  + # Define aqueous species
  + iFe.aq <- retrieve("Fe", c("S", "O", "H", "Cl"), "aq")
  + Fe.aq <- info(iFe.aq)$name
  + iCu.aq <- retrieve("Cu", c("S", "O", "H", "Cl"), "aq")
  + Cu.aq <- info(iCu.aq)$name
  + # Apply NaCl concentration
  + nacl <- NaCl(T = T, P = "Psat", m_tot = m_NaCl)
  + 
  + ### Setup basis species for Fe-Cu stacks
  + reset()
  + basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+", "Cl-"))
  + basis("H2S", logmS)
  + basis("Cl-", log10(nacl$m_Cl))
  + 
  + ## Diagram 1: Only Fe-bearing minerals
  + species(Fe.cr)
  + # Mosaic with S species as basis species
  + mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
  + dFe <- diagram(mFe$A.species, col = 2)
  + 
  + ## Diagram 2: Cu- (and FeCu-) bearing minerals and aqueous species
  + species(c("chalcopyrite", Cu.cr))
  + species(Cu.aq, -6, add = TRUE)
  + mFeCu <- mosaic(list(S.aq, Fe.cr), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant))
  + diagram(mFeCu$A.species)
  + 
  + ## Diagram 2a: Overlay single solubility contour
  + species(c("chalcopyrite", Cu.cr))
  + sout1 <- solubility(iCu.aq, bases = list(S.aq, Fe.cr), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant))
  + diagram(sout1, add = TRUE, col = 4, lwd = 2, levels = -6)
  + 
  + mtext("Stack 1: Fe minerals only -> Cu(Fe) minerals and Cu aqueous", adj = 1.1, line = 1.1)
  + 
  + ## Diagram 3: Fe-bearing minerals and aqueous species
  + species(Fe.cr)
  + species(Fe.aq, -6, add = TRUE)
  + mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
  + dFe <- diagram(mFe$A.species, col = 2)
  + 
  + ## Diagram 4: Cu- (and FeCu-) bearing minerals and aqueous species
  + species(c("chalcopyrite", Cu.cr))
  + species(Cu.aq, -6, add = TRUE)
  + mFeCu <- mosaic(list(S.aq, c(Fe.cr, Fe.aq)), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant), loga_aq = c(NA, logm_aq))
  + diagram(mFeCu$A.species)
  + 
  + ## Diagram 4a: Overlay single solubility contour
  + species(c("chalcopyrite", Cu.cr))
  + sout2 <- solubility(iCu.aq, bases = list(S.aq, c(Fe.cr, Fe.aq)), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dFe$predominant), loga_aq = c(NA, logm_aq))
  + diagram(sout2, add = TRUE, col = 4, lwd = 2, levels = -6)
  + 
  + mtext("Stack 2: Fe minerals and aqueous -> Cu(Fe) minerals and Cu aqueous", adj = 1, line = 1.1)
  + 
  + ### Setup basis species for Cu-Fe stacks
  + reset()
  + basis(c("Fe+2", "copper", "H2S", "oxygen", "H2O", "H+", "Cl-"))
  + basis("H2S", logmS)
  + basis("Cl-", log10(nacl$m_Cl))
  + 
  + ## Diagram 5: Only Cu-bearing minerals
  + species(Cu.cr)
  + mCu <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
  + dCu <- diagram(mCu$A.species, col = 4)
  + 
  + ## Diagram 6: Fe- (and CuFe-) bearing minerals and aqueous species
  + species(c("chalcopyrite", Fe.cr))
  + species(Fe.aq, -6, add = TRUE)
  + mCuFe <- mosaic(list(S.aq, Cu.cr), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dCu$predominant))
  + diagram(mCuFe$A.species)
  + 
  + ## Diagram 6a: Overlay single solubility contour
  + species(c("chalcopyrite", Fe.cr))
  + sout3 <- solubility(iFe.aq, bases = list(S.aq, Cu.cr), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dCu$predominant))
  + diagram(sout3, add = TRUE, col = 2, lwd = 2, levels = -6)
  + 
  + mtext("Stack 3: Cu minerals only -> Fe(Cu) minerals and Fe aqueous", adj = 1.1, line = 1.1)
  + 
  + ## Diagram 7: Cu-bearing minerals and aqueous species
  + species(Cu.cr)
  + species(Cu.aq, -6, add = TRUE)
  + mCu <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
  + dCu <- diagram(mCu$A.species, col = 4)
  + 
  + ## Diagram 8: Fe- (and CuFe-) bearing minerals and aqueous species
  + species(c("chalcopyrite", Fe.cr))
  + species(Fe.aq, -6, add = TRUE)
  + mCuFe <- mosaic(list(S.aq, c(Cu.cr, Cu.aq)), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dCu$predominant), loga_aq = c(NA, logm_aq))
  + diagram(mCuFe$A.species)
  + 
  + ## Diagram 8a: Overlay single solubility contour
  + species(c("chalcopyrite", Fe.cr))
  + sout4 <- solubility(iFe.aq, bases = list(S.aq, c(Cu.cr, Cu.aq)), pH = pH, O2 = O2, T = T, IS = nacl$IS, stable = list(NULL, dCu$predominant), loga_aq = c(NA, logm_aq))
  + diagram(sout4, add = TRUE, col = 2, lwd = 2, levels = -6)
  + 
  + mtext("Stack 4: Cu minerals and aqueous -> Fe(Cu) minerals and Fe aqueous", adj = 1, line = 1.1)
  + 
  + dev.off()
  + 
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.23    0.06    0.26