Bumped to 0.3.2 - Added Pourbaix()

DESCRIPTION

 Package: RedModRphree
 Title: Programmable Interface to the PHREEQC Geochemical Solver 
-Version: 0.3.1
+Version: 0.3.2
 Authors@R: c(person(given = "Marco",
                     family = "De Lucia",
                     email = "delucia@gfz-potsdam.de",

NAMESPACE

@@ -20,11 +20,14 @@ export(FindAllSpeciesNames)
 export(FindAllTotNames)
 export(FindLogK)
 export(FindPhase)
+export(FormSelectedOutput)
 export(FormulaFromBal)
+export(FormulaToExpression)
 export(Matplot)
 export(MatplotSingle)
 export(ParseFormula)
 export(PlotModsInSample)
+export(Pourbaix)
 export(RGetPhases)
 export(RPhreeExt)
 export(RPhreeFile)
@@ -32,7 +35,7 @@ export(RPhreeWriteInp)
 export(RPinfo)
 export(RReadOut)
 export(RReadOutKin)
-export(RReadPhases)
+export(RReadOutVal)
 export(ReactTranspBalanceEq)
 export(ReactTranspBalanceKin)
 export(RecomposeState)
@@ -49,6 +52,7 @@ export(pH2Act)
 export(splitMultiFix)
 export(splitMultiKin)
 export(stopmsg)
+export(strrev)
 import(doParallel)
 import(foreach)
 import(graphics)

R/Rphree_Files.R

 ### File-related utility functions for RedModRphree
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
-### Time-stamp: "Last modified 2021-04-12 12:39:15 delucia"
+### Time-stamp: "Last modified 2021-04-28 16:52:34 delucia"
 
 
 ##' Reads a normal PHREEQC input file and prepares it for
@@ -88,12 +88,12 @@ RReadOut <- function(out)
 {
     if (length(out)==1) {## is it a buffer or a file?
         ## note that empty lines are removed!
-        cat(paste("RReadOut:: opening file",out,"\n"))
+        msg("opening file",out)
         tot <- RPhreeFile(out, is.db=FALSE, tabs=TRUE) 
     } else {
         ## remove empty lines and tabs as in RPhreeFile, store
         ## everything in "tot"
-        cat("RReadOut:: scanning the buffer... \n")
+        msg("scanning the buffer...")
         tot <- sub(' +$', '', out)  ## remove spaces at the end of the lines
         tot <- sub("#.*$","", tot)  ## remove everything after the #'s
         tot <- gsub('\t', '  ', tot)  ## substitute tabs with spaces
@@ -113,7 +113,7 @@ RReadOut <- function(out)
     ## This is unique (only for calculated solution)
     endsim <- grep('End of simulation.',tot,fixed=TRUE)
 
-    cat(paste("RReadOut::",ntot," simulation in the given output"))
+    msg(ntot," simulations in the given output")
     
     ## "phase assemblage" is unique, but could not be there
     has_pphases <- TRUE
@@ -145,7 +145,7 @@ RReadOut <- function(out)
         error <- grep("^ERROR",tot[ solutions[n] : endsim[n] ])
         if (length(error) > 0) {
             res[[n]] <- "error"
-            cat(" ERROR!!\n")
+            msg(" ERROR in simulation",n,", skipping")
             next
         }
 
@@ -239,37 +239,34 @@ RReadOut <- function(out)
         
     }  ## end of loop over simulations
 
-    cat(" OK\n")
+    msg(" Finished, bye!")
     
     return(res)
     
 }
 
-##' Kinetics simulation: reads a phreeqc output file and forms an
-##' output list - as if the calculation was made through
-##' Rphree.
-##'
-##' 
+##' Kinetics simulation: reads a phreeqc output file or buffer and
+##' forms an output list.
 ##' @title RReadOutKin, import the output file of a kinetic simulation
 ##'     into R
-##' @param out The PHREEQC output file.
+##' @param out The PHREEQC output file or buffer
 ##' @param strip logical. If TRUE, the "ListInfo" element will be
-##' appended to the list.
+##'     appended to the list
 ##' @param verbose logical. If TRUE more output is given to the
-##' console (for debugging).
-##' @return  An output list as per Rphree call.
+##'     console (for debugging)
+##' @return An output list
 ##' @author MDL
 ##' @export
 RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
 {
     if (length(out)==1) {## is it a buffer or a file?
         ## note that empty lines are removed!
-        cat(paste("RReadOutKin:: opening file",out,"\n"))
+        msg("opening file",out)
         tot <- RPhreeFile(out, is.db=FALSE, tabs=TRUE) 
     } else {
         ## remove empty lines and tabs as in RPhreeFile, store
         ## everything in "tot"
-        cat("RReadOutKin:: scanning the buffer... \n")
+        msg("scanning the buffer...")
         tot <- sub(' +$', '', out)  ## remove spaces at the end of the lines
         tot <- sub("#.*$","", tot)  ## remove everything after the #'s
         tot <- gsub('\t', '  ', tot)  ## substitute tabs with spaces
@@ -281,7 +278,8 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
     if (length(toremove)>0) {
         toremove <- sort(c(toremove, toremove+1))
         tot <- tot[-toremove]
-        if (verbose) cat("RReadOutKin:: Appears to be PHREEQC version 3\n")
+        if (verbose)
+            msg("Appears to be PHREEQC version 3")
     }
     
     times <- grep("Time step:", tot, fixed=TRUE)
@@ -290,7 +288,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
     years <- round(as.numeric(sub('\ .*$','',gsub('.*:\ ','',tot[times]))),1)
     if (verbose) {
        cat(out,":\n")
-       cat(paste("RReadOutKin:: Found ",ntot,"time steps, the last at time",years[ntot],"\n"))
+       msg("Found ",ntot,"time steps, the last at time",years[ntot])
     }
     endkin  <- grep('-Phase assemblage-',tot,fixed=TRUE)
     endpure <- grep('-Solution composition-',tot,fixed=TRUE)[-1]
@@ -314,7 +312,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
     ## loop over all steps
     for (n in seq_along(times)) {
     	if (verbose)
-	    cat(paste("RReadOutKin:: Reading", n, "solution, time ",years[n],"\n"))
+	    msg("Reading", n, "solution, time ",years[n])
         ## find the last solution
         start <- times[n]+2
 
@@ -327,7 +325,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
         
         colnames(kin) <- c("moles","delta")
         if (verbose)
-            cat(paste("RReadOutKin:: Read kinetic block ", n, "\n"))
+            msg("Read kinetic block ", n)
         
         ## next block in output file is EQUILIBRIUM_PHASES
         startpure <- endkin[n]+3
@@ -338,7 +336,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
         
         colnames(pure) <- c("moles","delta")
         if (verbose)
-            cat(paste("RReadOutKin:: Read pphases block ", n, "of length",npure+1," \n"))
+            msg("Read pphases block ", n, "of length",npure+1)
         
         ## now solutes  
         startcomp <- endpure[n] + 2
@@ -350,7 +348,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
         
         colnames(comp) <- c("molal","moles")
         if (verbose)
-            cat(paste("RReadOutKin:: Read total solutes block ", n, "\n"))
+            msg("Read total solutes block ", n)
 
         ## desc: pH, pe, ecc
         block <- tot[(endcomp[n]+1):(enddesc[n]-1)]
@@ -394,7 +392,7 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
         colnames(species) <- c("molal","act")
         
         if (verbose)
-            cat(paste("RReadOutKin:: Read speciation ", n, "\n"))
+            msg("Read speciation ", n)
 
         ## Now saturation indexes
         block <- tot[(endspec[n] + 2) :( endsim[n]-1)]
@@ -419,3 +417,198 @@ RReadOutKin <- function(out, strip=TRUE, verbose=FALSE)
     return(res)
 }
 
+
+##' Reads a phreeqc output file or buffer and forms a results list.
+##' This version looks for a precise valence of an element. It is used
+##' primarily by \code{Pourbaix} to track the species and the SIs
+##' which must be retained for the diagram.
+##' @title RReadOutVal, parse a PHREEQC output file/buffer in search
+##'     of specific valence state of an element
+##' @param out The PHREEQC output file or buffer
+##' @param valence the specified element valence, e.g., Fe(2)
+##' @return An output list, as if the simulation would have being run
+##'     through Rphree (the same blocks and the same names are
+##'     returned)
+##' @author MDL
+##' @export
+RReadOutVal <- function(out, valence)
+{
+    if (length(out)==1) {## is it a buffer or a file?
+        ## note that empty lines are removed!
+        msg("opening file",out)
+        tot <- RPhreeFile(out, is.db=FALSE, tabs=TRUE) 
+    } else {
+        ## remove empty lines and tabs as in RPhreeFile, store
+        ## everything in "tot"
+        msg("scanning the buffer...")
+        tot <- sub(' +$', '', out)  ## remove spaces at the end of the lines
+        tot <- sub("#.*$","", tot)  ## remove everything after the #'s
+        tot <- gsub('\t', '  ', tot)  ## substitute tabs with spaces
+        tot <- tot[tot!=""]        ## remove empty elements
+    }
+
+    ## some extra lines in PHREEQC version 3 we need to take care of
+    toremove <- grep("^\\*\\*For", tot)
+    if (length(toremove)>0) {
+        toremove <- sort(c(toremove, toremove+1))
+        tot <- tot[-toremove]
+    }
+   
+    solutions <- grep("Beginning of initial solution calculations", tot)
+    ntot <- length(solutions)
+
+    ## This is unique (only for calculated solution)
+    endsim <- grep('End of simulation.',tot,fixed=TRUE)
+
+    msg(ntot," simulations in the buffer/file")
+
+    ## "phase assemblage" is unique, but could not be there
+    has_pphases <- TRUE
+    endkin  <- grep('-Phase assemblage-',tot,fixed=TRUE)
+    if (length(endkin)==0) {
+        ## This means we don't have any pure phase in the simulations,
+        ## and that the only solutions to read are the initial ones!
+        has_pphases <- FALSE
+        initials <- seq(1,ntot)
+    } else {
+        ## this index stores the lines marking initial solutions
+        initials <- -(seq(1,ntot)*2 -1)
+       
+    }
+    ## for these, don't consider the "initial solution" instead
+    endpure <- grep('-Solution composition-',tot,fixed=TRUE)[initials]
+    endcomp <- grep('-Description of solution-',tot,fixed=TRUE)[initials]
+    enddesc <- grep('-Distribution of species-',tot,fixed=TRUE)[initials]
+    endspec <- grep('-Saturation indices-',tot,fixed=TRUE)[initials]
+
+    ## create the container
+    res <- vector(mode="list",length=ntot)
+
+    ## loop over all steps
+    for (n in seq_along(solutions)) {
+        
+        ## check if there is a fatal error
+        error <- grep("^ERROR",tot[ solutions[n] : endsim[n] ])
+        if (length(error) > 0) {
+            res[[n]] <- "error"
+            stopmsg(" ERROR!!")
+            next
+        }
+
+        ## next block in output file is the "optional"
+        ## EQUILIBRIUM_PHASES
+        if (has_pphases) {
+            startpure <- endkin[n]+3
+            npure <- endpure[n] - startpure - 1
+        
+            pure_conn <- textConnection( tot[startpure:(startpure+npure)])
+            pure <- utils::read.table( pure_conn, row.names=1, fill=TRUE, as.is=TRUE)
+            close(pure_conn)
+            indreactants <- which(pure[,2]=="reactant")
+            if (length(indreactants) > 0)
+                {
+                    for (ir in indreactants) {
+                        rownames(pure)[ir-1] <- paste(rownames(pure)[ir-1], rownames(pure)[ir],sep="_")
+                        pure[(ir-1),c(4,5,6)] <- pure[ir,c(3,4,5)]
+                    }
+                    pure <- pure[-indreactants,]
+                }
+            pure <- pure[,c(5,6)]
+            
+            colnames(pure) <- c("moles","delta")
+            ##    if (verbose)
+            ##      cat(paste(":: Read pphases block ", n, "of length",npure+1," \n"))
+        } else {
+            pure <- NA
+        }
+        
+        ## now total solutes  
+        startcomp <- endpure[n] + 2
+        ncomp <- endcomp[n]-startcomp
+        comp_conn <- textConnection(tot[startcomp:(startcomp+ncomp-1)])
+        comp <- utils::read.table(comp_conn,row.names=1,fill=TRUE)
+        close(comp_conn)
+        colnames(comp) <- c("molal","moles")
+        ##    if (verbose)
+        ##      cat(paste(":: Read total solutes block ", n, "\n"))
+        
+        ## desc: pH, pe, ecc
+        block <- tot[(endcomp[n]+1):(enddesc[n]-1)]
+
+        pH <- as.numeric(unlist(strsplit(grep('^pH',block, value=TRUE)," +"))[3])
+        pe <- as.numeric(unlist(strsplit(grep('^pe',block, value=TRUE)," +"))[3])
+        temp <- as.numeric(unlist(strsplit(grep('^Temperature ',block,value=TRUE)," = "))[2])
+        water <- as.numeric(unlist(strsplit(grep('Mass of water',block,fixed=TRUE,value=TRUE)," = "))[2])
+        WaterActivity <- as.numeric(unlist(strsplit(grep('Activity of water',block,fixed=TRUE,value=TRUE)," = "))[2])
+        IonStr <- as.numeric(unlist(strsplit(grep('Ionic strength',block,fixed=TRUE,value=TRUE)," = "))[2])
+        TotAlk <- as.numeric(unlist(strsplit(grep('Total alkalinity',block,fixed=TRUE,value=TRUE)," = "))[2])
+        Tot_C  <- as.numeric(unlist(strsplit(grep('Total carbon',block,fixed=TRUE,value=TRUE)," = "))[2])
+        Tot_CO2  <- as.numeric(unlist(strsplit(grep('Total CO2',block,fixed=TRUE,value=TRUE)," = "))[2])
+        ElBal <- as.numeric(unlist(strsplit(grep('Electrical balance',block,fixed=TRUE,value=TRUE)," = "))[2])
+        Per_Error <- as.numeric(unlist(strsplit(grep('Percent error',block,fixed=TRUE,value=TRUE)," = "))[2])
+        Iter <- as.numeric(unlist(strsplit(grep('^Iterations',block,value=TRUE)," = "))[2])
+        Tot_H <- as.numeric(unlist(strsplit(grep('Total H',block,fixed=TRUE,value=TRUE)," = "))[2])
+        Tot_O <- as.numeric(unlist(strsplit(grep('Total O',block,fixed=TRUE,value=TRUE)," = "))[2])
+                block <- tot[(endcomp[n]+1):(enddesc[n]-1)]
+
+        desc <- data.frame(val=c(pH=pH,pe=pe,temp=temp,water=water,WaterActivity=WaterActivity,
+                               Tot_Alk=TotAlk, IonicStr=IonStr, ElectrBal=ElBal,
+                               Per_Error=Per_Error, iterations=Iter, Tot_CO2=Tot_CO2, Tot_H=Tot_H, Tot_O=Tot_O))
+
+
+        ## next block in output is the speciation
+        if (missing(valence)) {
+            block <- tot[(enddesc[n] + 4) :( endspec[n] - 1)]
+            ## find out the short lines
+            shorts <- nchar(block) 
+            excl <- which(shorts < mean(shorts))
+            block <- block[-excl]
+            
+            block_conn <- textConnection(block)
+            ## Now we can read.table it
+            tmp <- utils::read.table( block_conn, fill=TRUE, as.is=TRUE)[,c(1,2,3)]
+            close(block_conn)
+            ## remove duplicated
+            rnames <- tmp$V1
+            dup <- which(!duplicated(rnames))
+            species <- tmp[dup,c(2,3)]
+            rownames(species) <- rnames[dup]
+            colnames(species) <- c("molal","act")
+        } else {
+            ## we only select species of the given valence state
+            block <- tot[(enddesc[n] + 4) :( endspec[n] - 1)]
+            ## find out the short lines
+            shorts <- nchar(block) 
+            excl <- which(shorts < mean(shorts))
+            ind_val <- grep(valence, block[excl], fixed=TRUE)
+            nb <- block[ seq(excl[ind_val]+1, excl[ind_val+1]-1)]
+
+            block_conn <- textConnection(nb)
+            ## Now we can read.table it
+            tmp <- utils::read.table( block_conn, fill=TRUE, as.is=TRUE)[,c(1,2,3)]
+            close(block_conn)
+            ## remove duplicated
+            rnames <- tmp$V1
+            dup <- which(!duplicated(rnames))
+            species <- tmp[dup,c(2,3)]
+            rownames(species) <- rnames[dup]
+            colnames(species) <- c("molal","act")
+        }
+
+        ## Now saturation indexes
+        block <- tot[(endspec[n] + 2) :( endsim[n] - 4)]
+        block_conn <- textConnection(block)
+        SI <- utils::read.table(block_conn, fill=TRUE, row.names=1, as.is=TRUE)
+        close(block_conn)
+        names(SI) <- c("SI","IAP","logK","formula")
+        
+        ## finally pack all together
+        res[[n]] <- list(desc=desc,pphases=pure,tot=comp, SI=SI, species=species)
+        
+    }  ## end of loop over simulations
+
+    msg(" OK")
+    
+    return(res)
+    
+}

R/Rphree_Pourbaix.R

+### Licence: LGPL version 2.1
+## Time-stamp: "Last modified 2021-04-28 17:04:27 delucia"
+
+##' @title Revert a string
+##' @param x character vector. All element of the vector will be
+##'     reverted
+##' @return the reverted character vector
+##' @export
+##' @author MDL
+strrev <- function(x) {
+    res <- x
+    for (i in seq_along(x)) {
+        nc <- nchar(x[i])
+        res[i] <- paste0(substring(x[i], nc:1, nc:1), collapse = "")
+    }
+    res
+}
+
+##' @title Converts a PHREEQC formula into R expression for nice
+##'     annotations
+##' @param ff character vector of formulas which will be translated
+##' @return character vector of valid R expressions
+##' @author MDL
+##' @export
+FormulaToExpression <- function(ff) {
+### transform PHREEQC formulas into R's expressions
+
+    add.sup   <- FALSE
+    add.colon <- FALSE
+
+    
+    ## Body is everything which precedes a "-", a "+" or a ":"
+    body <- sub("[\\:\\+\\-].*$", "", ff)
+    ##body <- sub("[+-:].*$", "", ff)
+
+    ## first we take the exponents for the charge
+    has.sup <- grep("[\\+\\-]", ff)
+    if (length(has.sup)>0) {
+        add.sup <- TRUE
+        super <- sub("(^.*)([\\+\\-])", "\\2", ff[has.sup])
+        super[nchar(super)==2] <- strrev(super[nchar(super)==2])
+        esup <- paste0("^{",super,"phantom()}")
+    }
+
+    ## treat the ":"
+    has.colon <- grep(":", ff, fixed=TRUE)
+    if (length(has.colon)>0) {
+        add.colon  <- TRUE
+        aftercolon <- sub("^.*\\:", ":", ff[has.colon])
+        after <- sub("(\\:[0-9]*)H2O", "\\1*H[2]*O", aftercolon)
+        after <- sub(":", "*':'*", after)
+    }
+
+    has.sub <- grep("[0-9]+", body)
+    body[has.sub] <- gsub("([0-9]+)", "[\\1]", body[has.sub])
+
+    body <- gsub("\\]([A-Z])", "]*\\1" , body)
+    formula <- body
+    if (add.sup)
+        formula[has.sup]   <- paste0(formula[has.sup], esup)
+    if (add.colon)
+        formula[has.colon] <- paste0(formula[has.colon],after)
+    return(formula)
+}
+
+##' This function takes the string output of a PHREEQC calculation and
+##' forms the "SELECTED_OUTPUT" for further simulations including all
+##' totals and the activities and saturation indeces of all species
+##' and minerals which contain the optional \emph{element}
+##'
+##' @title Form SELECTED_OUTPUT from a PHREEQC string output
+##' @param sim the string output of a previous PHREEQC simulation
+##' @param element optional, the single element which must be present
+##'     in all species and minerals included in the SELECTED_OUTPUT
+##' @return A character vector containing the SELECTED_OUTPUT block
+##' @author MDL
+##' @export
+FormSelectedOutput <- function(sim, element) {
+    AllSpecies <- rownames(sim$species)
+    AllTot     <- rownames(sim$tot)
+    AllSI      <- rownames(sim$SI)
+    
+    if (!missing(element)) {
+        ## first we grab the formula from the SI tab
+        mins   <- sim$SI$formula
+        indmin <- grep(element, mins)
+        AllSI  <- AllSI[indmin]
+        indsp  <-  grep(element, AllSpecies) 
+        AllSpecies  <- AllSpecies[indsp]
+    }
+    SI      <- paste("    -saturation_indices", paste(AllSI, collapse=" "))
+    Totals  <- paste("    -totals", paste(AllTot, collapse=" "))
+    Species <- paste("    -activities", paste(AllSpecies, collapse=" "))
+    ## make sure there is no "H2O"
+    Species <- sub(" H2O ", " ", Species, fixed=TRUE)
+    SelOut <- c("SELECTED_OUTPUT",
+                "    -high_precision  true",
+                "    -reset   false",
+                "    -simulation  true",
+                "    -pH      true",
+                "    -pe      true",
+                "    -water   true",
+                "    -charge_balance true",
+                "    -ionic_strength true",
+                "    -percent_error  true",
+                Totals, Species, SI)
+
+    return(SelOut)
+}
+
+
+##' This function uses PHREEQC to calculate and plot a Pourbaix
+##' diagram for the specified system at a given T and P. The
+##' geochemical system is specified as a single PHREEQC input
+##' specifying one "SOLUTION" while the actual values of pe and pH at
+##' which the predominance shall be calculated must be specified by
+##' the user at input. To try and prevent numerical instabilities,
+##' which are quite frequent, the function only actually performs
+##' PHREEQC calculations inside the water stability region.
+##'
+##' @title Pourbaix (pe/pH) diagrams using PHREEQC
+##' @param element character of length 1. The element for which the
+##'     Pourbaix diagram is being calculated
+##' @param pe numeric vector, grid values for redox potential at which
+##'     predominance is calculated
+##' @param pH numeric vector, grid values for pH at which predominance
+##'     is calculated
+##' @param Tc temperature in Celsius
+##' @param Patm pressure in atmosphere
+##' @param base character vector. The template PHREEQC input script,
+##'     mostly containing only one SOLUTION block
+##' @param first Optional character vector containing i.e. own PHASES
+##'     definitions which comes on top of the base input
+##' @param db PHREEQC chemical database
+##' @param aqonly logical, defaults to FALSE. If TRUE, only aqueous
+##'     species are considered
+##' @param suppress optional character vector. Species and minerals
+##'     contained here are not considered in the diagram
+##' @param plot logical, defaults to TRUE. Do you want to plot the
+##'     diagram?
+##' @param ann.title Optional string. If specified, it will replace
+##'     main title of the picture
+##' @param ann.sub Optional string. If specified, it will replace
+##'     subtitle of the picture (usually this specifies the
+##'     concentrations)
+##' @param colors Optional character vector containing the colors. If
+##'     missing, grDevices::terrain.colors is used internally for
+##'     plotting
+##' @param procs how many CPUs should be used for the computations?
+##'     Defaults to 1. Parallelization is achieved through DoParallel
+##'     and foreach, and should work on both *nix and Windows using
+##'     fork and sockets respectively
+##' @param ... further graphic parameters passed to the \code{image}
+##'     for plotting
+##' @return invisibly returns a list containing all calculated
+##'     simulations, input script for PHREEQC simulations, formulas,
+##'     names, etc
+##' @author MDL
+##' @export
+##' @examples
+##' base <- c("SOLUTION 1",
+##'           "units mol/kgw",
+##'           "temp 25.0",
+##'           "pressure 1",
+##'           "-water 1",
+##'           "pH 7",
+##'           "pe 4.0",
+##'           "Na 1",
+##'           "C 0.0125",
+##'           "Ca 0.025",
+##'           "Mg 0.025",
+##'           "Cl 2.002 charge",
+##'           "END")
+##' par(mfrow=c(1,2))
+##' a <- Pourbaix(element="Ca", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25,
+##'               Patm=1, base=base, first="", ann.title="Custom title", db=llnl.dat)
+##' b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25,
+##'               Patm=1, base=base, first="", ann.sub="Custom subtitle", db=llnl.dat)
+##'
+##' ## Now we restrict the speciation to a fixed valence state of the element,
+##' ## here Fe(2) and Cu(2), only aqueous species:
+##' base <- c("SOLUTION 1",
+##'           "units mol/kgw",
+##'           "temp 25.0",
+##'           "pressure 1",
+##'           "-water 1",
+##'           "pH 7",
+##'           "pe 4.0",
+##'           "Ca        0.025",
+##'           "Cu        0.001",
+##'           "Cl        0.05122 charge",
+##'           "Fe        3.0E-04",
+##'           "END")
+##' par(mfrow=c(1,2))
+##' a <- Pourbaix(element="Fe(2)", pe=seq(-10,12, length=51), pH=seq(2,10, length=51),
+##'               Patm=1, base=base, aqonly=TRUE,  db=llnl.dat)
+##' b <- Pourbaix(element="Cu(2)", pe=seq(-10,12, length=51), pH=seq(2,10, length=51),
+##'               Patm=1, base=base, aqonly=TRUE,  db=llnl.dat)
+Pourbaix <- function(element, pe=seq(-5,8), pH=seq(6,12), Tc=25, Patm=1, base, first, db,
+                     aqonly=FALSE, suppress, ann.title, ann.sub, colors, procs=1, plot=TRUE,
+                     ...)
+{  
+    phreeqc::phrLoadDatabaseString(db)
+    phreeqc::phrSetOutputStringsOn(TRUE)
+
+    ## treat "first" as a really optional argument
+    if (missing(first))
+        first <- ""
+    ## testrun of base script to get the output
+    aa  <- phreeqc::phrRunString(c(first, base))
+    con <- phreeqc::phrGetOutputStrings()
+
+    has_valence <- FALSE
+    if (missing(element)) {