CHNOSZ via Rpad

Web interface to activity diagrams

setup environment

choose example

        # keep track of how many times we've run
        if(!exists("REnvSetupNo")) assign("REnvSetupNo",1) else assign("REnvSetupNo",REnvSetupNo+1)
        cat(paste("WWW: REnvSetup (",REnvSetupNo,")\n"))
        # if we are being run more than once that means 
        # the Rpad localversion is in use and the variables 
        # stay around between page reloads.
        assign("SelectBasisGo",TRUE)
        assign("SelectSpeciesGo",TRUE)
        # load required libraries
        library(CHNOSZ)
        # set options
        options(width=120)
        #if(REnvSetupNo==1) {
          # hide the output of initial example setup
          sink("setup.txt")
          # some default defaults
          assign("data.keys","HDN+78,WEP+82,SH88,SHS89")
          #"HDN+78,WEP+82,SH88,SHS89,RH95,Gur96,SSW+97,SSW01,TS01,DLH06,LH06a,LH06b"
          assign("temperature",25); assign("pressure",1)
          assign("default.logact",numeric())
          assign("default.slogact",numeric())
          assign("default.coeff",numeric())
          assign("protein.use",1:nrow(thermo$protein))
          assign("default.protein",paste("PROT",(1:8),"_NEW",sep=""))
          assign("Tmin",25); assign("Tmax","")
          assign("Pmin","0"); assign("Pmax","")
          assign("line.types",1)
          assign("use.residues",TRUE)
          assign("fill.areas",TRUE)
          # species in examples 1 and 2
          assign("species.1",paste(c("GLRX1","GLRX2","GLRX3","GLRX4","GLRX5"),"YEAST",sep="_"))
          assign("states.1",rep("aq",length(species.1)))
          assign("species.2",c("sulfur","H2S","HS-","S2O3-2","HS2O3-"))
          assign("states.2",c("cr1","aq","aq","aq","aq"))
          assign("basis.1",c("CO2","H2O","NH3","e-","H2S","H+"))
          assign("basis.2",c("H2O","H2S","e-","H+"))
          assign("exampletext",c("redoxin","sulfur","iron","ATP",
            "amino","s-layer","sigma","yeast"))
          assign("whichexample",match(whichexample,exampletext))
            if(whichexample==1) {
            # glutaredoxins Eh-pH
            assign("default.basis",info(basis.1))
            assign("other.basis",info(basis.2))
            assign("default.logact",c(-5,0,-4,-7,-12,-7))
            assign("default.species",info(species.1))
            assign("default.slogact",rep(-6,length(default.species)))
            assign("other.species",info(species.2,states.2))
            assign("data.keys",c2s(as.character(c(default.basis,other.basis,
              other.species)),sep=","))
            assign("protein.use",as.numeric(rownames(protein(species.1))))
            assign("xvar","pH"); assign("yvar","Eh")
            assign("xmin",0); assign("xmax",14)
            assign("ymin",-0.6); assign("ymax",0.6)
            assign("default.coeff",c(paste("-1",protein.length("GLRX1_YEAST"),sep="/"),
              paste("1",protein.length("GLRX2_YEAST"),sep="/"),0))
            assign("Tmin",24); assign("Tmax",99)
            assign("Pmin",1); assign("Pmax","")
          } else if(whichexample==2) {
            # sulfur Eh-pH after Fig. 3 of Sato 1992 
            # (http://adsabs.harvard.edu/abs/1992GeCoA..56.3133S)
            assign("default.basis",info(basis.2))
            assign("other.basis",info(basis.1))
            assign("default.logact",c(0,99,99,-7))
            assign("default.species",info(species.2,states.2))
            assign("default.slogact",c(0,-1,-1,-1,-1))
            assign("data.keys",c2s(as.character(c(other.basis,default.basis,
              default.species)),sep=","))
            assign("protein.use",as.numeric(rownames(protein(species.1))))
            assign("xvar","pH"); assign("yvar","Eh")
            assign("xmin",0); assign("xmax",14)
            assign("ymin",-0.6); assign("ymax",0.6)
            assign("default.coeff",c(0,-1,1))
            assign("Tmin",0); assign("Tmax",100)
            assign("Pmin",0); assign("Pmax","")
          } else if(whichexample==3) {
            # Fe-S-O after Helgeson 1970
            assign("default.basis",info(c("FeS2","O2","S2")))
            assign("default.logact",c(0,-50,0))
            assign("default.species",info(c("iron","ferrous-oxide","magnetite",
              "hematite","pyrite","pyrrhotite")))
            assign("data.keys","HDN+78,WEP+82,RH95")
            assign("protein.use",0)
            assign("temperature",200); assign("pressure",1)
            assign("xvar","S2"); assign("yvar","O2")
            assign("xmin",-50); assign("xmax",0)
            assign("ymin",-90); assign("ymax",-10)
            assign("default.coeff",c(0,0,2,-3,0,0))
            assign("Tmin",0); assign("Tmax",300)
            assign("Pmin",1); assign("Pmax","")
          } else if(whichexample==4) {
            # MgATP after LaRowe and Helgeson, 2007 (and ADP too ...)
            assign("data.keys","WEP+82,SH88,SHS89,LH06a,LH06b")
            assign("default.basis",info(c("CO2","NH3","H2O","H3PO4","O2","H+","Mg+2")))
            assign("default.logact",c(99,99,0,5,99,-5,-4))
            assign("default.species",info(c("HATP-3","H2ATP-2","MgATP-2","MgHATP-")))
            assign("default.species",c(default.species,
              info(c("HADP-2","H2ADP-","MgADP-","MgHADP"))))
            assign("default.coeff",c(0,0,0,1,-1,0,0,0))
            assign("default.slogact",rep(-5,length(default.species)))
            assign("species.1",c("LYSC_CHICK","RNAS1_BOVIN","CYC_BOVIN","MYG_PHYCA"))
            assign("protein.use",as.numeric(rownames(protein(species.1))))
            assign("pressure","Psat")
            assign("xvar","T"); assign("yvar","")
            assign("xmin",0); assign("xmax",140)
            assign("ymin",0); assign("ymax",1)
            assign("Tmin",0); assign("Tmax",140)
            assign("Pmin",""); assign("Pmax","")
            assign("line.types",4)
          } else if(whichexample==5) {
            # amino acid-nucleobase interaction balanced on carbon
            assign("data.keys","WEP+82,SH88,SHS89,LH06a,DLH06")
            assign("default.basis",info(c("CO2","H2O","glutamine","H2")))
            assign("default.logact",c(-3,-1,-3,-4))
            assign("default.species",info(c('uracil','cytosine','adenine',
              'guanine','phenylalanine','proline','lysine','glycine')))
            assign("default.slogact",rep(-3,length(default.species)))
            assign("default.coeff",c(0,0,0,-0.5,0,0,0,7.5))
            assign("pressure","Psat")
            assign("xvar","H2O"); assign("yvar","H2")
            assign("xmin",-5); assign("xmax",0)
            assign("ymin",-9); assign("ymax",-4)
            assign("Tmin",""); assign("Tmax",25)
            assign("Pmin",1); assign("Pmax",2000)
            assign("fill.areas",FALSE)
          } else if(whichexample==6) {
            # surface-layer proteins H2
            assign("default.basis",info(c('CO2','H2O','NH3','H2','H2S')))
            assign("default.logact",rep(-5,length(default.basis)))
            assign("default.species",info(c("CSG_HALJP","CSG_METFE","CSG_METJA",
              "CSG_METSC","CSG_METVO","CSG_METBU","SLAP_ACEKI","SLAP_BACST",
              "SLAP_CAMFE","SLAP_LACAC")))
            assign("temperature",25); assign("pressure",1)
            assign("xvar","H2"); assign("yvar","")
            assign("xmin",-12); assign("xmax",0)
            # ylim=c(0,1) for a degree of formation (alpha) plot
            assign("ymin",0); assign("ymax",1)
            assign("line.types",2)
          } else if(whichexample==7) {
            # ecoli sigma factors T-pH
            assign("default.basis",info(c("CO2","H2O","NH3","H2S","H2","H+")))
            assign("default.logact",c(-3,0,-1,-12,-1,-7))
            assign("default.species",info(paste(c("RPOE","RPOS","RP32","FLIA",
              "RPOD","RP54"),"ECOLI",sep="_")))
            assign("default.coeff",c(paste("-1",protein.length("RPOE_ECOLI"),sep="/"),
              0,0,0,paste("1",protein.length("RPOD_ECOLI"),sep="/"),0))
            assign("use.residues",FALSE)
            assign("fill.areas",FALSE)
            assign("xvar","pH"); assign("yvar","T")
            assign("xmin",0); assign("xmax",14)
            assign("ymin",0); assign("ymax",100)
            assign("Tmin",0); assign("Tmax",100)
            assign("Pmin",0); assign("Pmax",0)
          } else if(whichexample==8) {
            # yeast subproteomes O2-H2O
            assign("data.keys","HDN+78,WEP+82,SH88,SHS89,RH95,Gur96,SSW+97,SSW01,TS01,DLH06,LH06a,LH06b")
            assign("default.basis",info(c("CO2","H2O","NH3","H2S","O2")))
            assign("default.logact",c(-3,0,-4,-7,-74))
            assign("locations",paste(c("actin","early.Golgi","ER",
              "vacuole","cell.periphery","nucleolus")))
            assign("default.protein",c(paste(locations,"YEASTGFP",sep="_"),default.protein[1:2]))
            assign("default.species",seq(nrow(thermo$obigt)+1,length.out=length(locations)))
            assign("default.coeff",c(paste("1",469.4,sep="/"),
              0,paste("-1",245.6,sep="/")))
            assign("default.slogact",rep(-3,length(default.species)))
            assign("xvar","O2"); assign("yvar","")
            assign("xmin",-83); assign("xmax",-67)
            assign("ymin",-4); assign("ymax",-1)
            assign("Tmin",25); assign("Tmax",30)
            assign("line.types",2)
          }
          sink()
          # keep only the selected proteins
          if(identical(protein.use,0)) assign("np",0) else assign("np",length(protein.use))
          cat(paste("WWW: this example provides",np,"of",nrow(thermo$protein),"proteins in the database.\n"))
        #}
      

search thermodynamic data

        if(!exists("SetupSearchNo")) assign("SetupSearchNo",1) else 
          assign("SetupSearchNo",SetupSearchNo+1)
        cat(paste("WWW: SetupSearch (",SetupSearchNo,")\n"))
        HTMLon()
        cat("search by name, formula, or species number")
        cat(paste("<br>"))
        H("input",name="mysearch",value="",rpadType="Rstring",size=20)
        cat(paste("<br>"))
        # a GO button for processing the search
        H("input",type="button",value="GO",onclick="rpad.calculateTree(dojo.byId('SearchDiv'))")
      
        if(!exists("ProcessSearchNo")) assign("ProcessSearchNo",1) else 
          assign("ProcessSearchNo",ProcessSearchNo+1)
        cat(paste("WWW: ProcessSearch (",ProcessSearchNo,")\n"))
        if(mysearch=="" | is.na(mysearch)) {
          #if(ProcessSearchNo > 1) cat("WWW: empty search string, try again\n")
        } else {
          if(can.be.numeric(mysearch)) {
            info(as.numeric(mysearch))
          } else {
            # extract a trailing aq, cr1,2,3,4, liq, gas
            assign("mystate",NULL)
            for(thisstate in c("aq","cr","cr1","cr2","cr3","cr4","liq","gas")) {
              assign("isstate",regexpr(paste(" ",thisstate,"$",sep=""),mysearch))
              if(isstate > 0) {
                assign("mysearch",substr(mysearch,1,isstate))
                assign("mystate",thisstate)
                break
              }
            }
            # show the results of a search using info()
            assign("myinfo",info(mysearch,mystate))
            # show the whole entry if available
            if(is.list(myinfo)) myinfo <- myinfo[[1]]
            if(any(!is.na(myinfo))) {
              if(length(which(!is.na(myinfo))) < 40) print(info(myinfo[!is.na(myinfo)]))
              # if it was a protein also print the length
              if(length(grep("_",mysearch) > 0)) cat(paste("WWW: the length of protein",
                mysearch,"is",protein.length(mysearch),"\n"))
            }
          }
        }
      

select thermodynamic data

        if(!exists("SelectDataNo")) assign("SelectDataNo",1) else assign("SelectDataNo",SelectDataNo+1)
        cat(paste("WWW: SelectData (",SelectDataNo,")\n"))
        cat("enter data source codes and/or species numbers\n")
        HTMLon()
        if(length(data.keys) > 1) assign("data.keys",c2s(data.keys,sep=","))
        H("input",name="data.keys",value=data.keys,rpadType="Rstring",standaloneTag=FALSE,size=42)
        cat(paste("<br>"))
        # a GO button for processing the data selection
        H("input",type="button",value="GO",onclick="rpad.calculateTree(dojo.byId('DataDiv'))")
      
        if(!exists("ProcessDataNo")) assign("ProcessDataNo",1) else assign("ProcessDataNo",ProcessDataNo+1)
        cat(paste("WWW: ProcessData (",ProcessDataNo,")\n"))
        # we use the selected sources plus entries 1 and 2 (water and electron)
        assign("data.keys",s2c(data.keys,sep=",",keep.sep=FALSE))
        assign("data.use",numeric())
        for(i in 1:length(data.keys)) {
          if(can.be.numeric(data.keys[i]))
          assign("data.use",c(data.use,as.numeric(data.keys[i]))) else assign("data.use",
            c(data.use,which(thermo$obigt$source1==data.keys[i] | thermo$obigt$source2==data.keys[i])))
        }
        assign("data.use",unique(c(1,2,data.use)))
        cat(paste("WWW:",length(data.use),"of",nrow(thermo$obigt)-length(grep("_",thermo$obigt$name)),
          "data entries activated\n"))
      
        // refresh the basis and species lists after ProcessData
        dojo.event.connect(dojo.widget.byId("ProcessData"), "onReceive",  
           dojo.widget.byId("SelectBasis"), "calculate");
        " ";
      

select basis species

        if(!exists("SelectBasisNo")) assign("SelectBasisNo",1) else assign("SelectBasisNo",SelectBasisNo+1)
        cat(paste("WWW: SelectBasis (",SelectBasisNo,")\n"))
        cat("basis species / logarithm of activity\n")
      ## the default species and logarithms of activities
        assign("n.basis",7)
        # assemble the names of species (with states)
        assign("namestuff",thermo$obigt$name)
        assign("basisnames",paste(namestuff,thermo$obigt$state))
        assign("usedbasisnames",paste(namestuff,thermo$obigt$state)[data.use])
        if(!identical(protein.use,0)) assign("usedbasisnames",c(usedbasisnames,
          paste(thermo$protein$protein[protein.use],"_",thermo$protein$organism[protein.use]," aq",sep="")))
        if(SelectBasisNo==1 | SelectBasisGo) {
          assign("default.basis",match(basisnames[default.basis],usedbasisnames))
          assign("SelectBasisGo",FALSE)
        } else {
          # keep the ones set by the user
          assign("default.basis",numeric())
          assign("default.logact",numeric())
          for(i in 1:n.basis) if(exists(paste("basis",i,sep=""))) {
            assign("mybasis",get(paste("basis",i,sep="")))
            assign("mylogact",get(paste("logact",i,sep="")))
            assign("mymatch",match(mybasis,usedbasisnames))
            if(!is.na(mymatch)) {
              assign("default.basis",c(default.basis,mymatch))
              assign("default.logact",c(default.logact,mylogact))
            }
          }
        }
        # make sure they're long enough
        assign("default.basis",c(default.basis,rep(0,n.basis)))
        assign("default.logact",c(default.logact,rep(0,n.basis)))
      ## make the form
        cat(HTMLon())
        for(i in 1:n.basis) {
          # the names (species states) pull-down menu
          assign("name",paste("basis",i,sep=""))
          cat(HTMLselect(name, c("",usedbasisnames),default=default.basis[i]+1))
          # entry field for logarithms of activities
          cat(HTMLinput(paste("logact",i,sep=""),value=default.logact[i],size=5))
          cat(paste("<br>"))
        }
        # a GO button for processing the basis species
        H("input",type="button",value="GO",onclick="rpad.calculateTree(dojo.byId('BasisDiv'))")
      
        if(!exists("ProcessBasisNo")) assign("ProcessBasisNo",1) else 
          assign("ProcessBasisNo",ProcessBasisNo+1)
        cat(paste("WWW: ProcessBasis (",ProcessBasisNo,")\n"))
        assign("mybasis",character())
        assign("mybasisname",character())
        assign("basisstates",character())
        assign("mylogact",character())
        assign("n.basis",as.numeric(n.basis))
        for(i in 1:n.basis) {
          if(exists(paste("basis",i,sep=""))) {
            assign("thisspecies",get(paste("basis",i,sep="")))
            if(!identical(thisspecies,"")) {
              assign("j",match(thisspecies,usedbasisnames))
              assign("mybasis",c(mybasis,c(thermo$obigt$formula[data.use],
                paste(thermo$protein$protein[protein.use],thermo$protein$organism[protein.use],sep="_"))[j]))
              assign("mybasisname",c(mybasisname,c(thermo$obigt$name[data.use],
                paste(thermo$protein$protein[protein.use],thermo$protein$organism[protein.use],sep="_"))[j]))
              assign("basisstates",c(basisstates,c(thermo$obigt$state[data.use],
                rep("aq",nrow(thermo$protein[protein.use,])))[j]))
              assign("mylogact",c(mylogact,get(paste("logact",i,sep=""))))
            }
          }
        }
        # now to apply the changes (and output the result)
        thermo$basis <- NULL
        trybasis <- try(basis(mybasisname,basisstates),TRUE)
        if(class(trybasis)=="try-error") {
          # cat("WWW: error in basis definition; the output was\n")
          cat(trybasis[1])
        } else basis(mybasis,mylogact)
        HTMLoff()
      
        // refresh the species list after the basis list
        dojo.event.connect(dojo.widget.byId("SelectBasis"), "onReceive",  
           dojo.widget.byId("SelectSpecies"), "calculate");
        " ";
      
        dojo.event.connect(dojo.widget.byId("ProcessBasis"), "onReceive",  
           dojo.widget.byId("SelectVariables"), "calculate");
        " ";
      

add protein sequences

        if(!exists("SelectProteinsNo")) assign("SelectProteinsNo",1) else 
          assign("SelectProteinsNo",SelectProteinsNo+1)
        cat(paste("WWW: SelectProteins (",SelectProteinsNo,")\n"))
        cat("protein sequence / name\n")
        HTMLon()
        assign("n.proteins",8)
        for(i in 1:n.proteins) {
          cat(H("input",name=paste("sequence",i,sep=""),value="",size="30"))
          cat(H("input",name=paste("protein",i,sep=""),value=default.protein[i],size="10"))
          cat("<br>")
        }
        # the GO button
        H("input",type="button",value="GO",
          onclick="javascript:rpad.calculateTree(dojo.byId('ProteinDiv'))")
      
        if(!exists("ProcessProteinsNo")) assign("ProcessProteinsNo",1) else 
          assign("ProcessProteinsNo",ProcessProteinsNo+1)
        cat(paste("WWW: ProcessProteins (",ProcessProteinsNo,")\n"))
        for(i in 1:n.proteins) {
          assign("mysequence",get(paste("sequence",i,sep="")))
          assign("myprotein",get(paste("protein",i,sep="")))
          if(mysequence=="") {
            # empty sequence: attempt to retrieve protein from
            # data files. 20080420 jmd
            if(length(grep("_",myprotein))==0) next
            assign("myprot",s2c(get(paste("protein",i,sep="")),sep="_",keep.sep=FALSE)[1])
            assign("myorg",s2c(get(paste("protein",i,sep="")),sep="_",keep.sep=FALSE)[2])
            if(myorg %in% c("ECO","SGD")) {
              # using protein compositions in thermo$ ECO or SGD
              assign("tryprotein",try(get.protein(myprot,myorg),TRUE))
              if(class(tryprotein)=="try-error") {
                cat(paste("WWW: adding protein",myprotein,"failed; the message was\n"))
                cat(tryprotein[1])
              } else {
                cat(paste("WWW: adding protein",myprotein,"\n"))
                assign("ip",add.protein(tryprotein)) 
                protein.use <- unique(c(protein.use,ip))
              }
            } else if(myorg=="YEASTGFP") {
              # make a proteolog for the given subcellular location
              if(!myprot %in% colnames(thermo$yeastgfp)[6:28]) {
                cat(paste("WWW: YEASTGFP location",myprot,"is unknown\n"))
                next
              } else {
                cat(paste("WWW: looking up proteins for YEASTGFP location",myprot,"\n"))
                assign("yprot",yeastgfp(myprot))
                assign("tp",get.protein(yprot$yORF,"SGD",yprot$abundance,myprot))
                assign("ip",add.protein(tp))
                # for the yeast example to work we manually add the protein
                # to thermo$obigt
                info(paste(myprot,"SGD",sep="_"))
                protein.use <- unique(c(protein.use,ip))
              }
            } else {
              # try looking up the protein in thermo$protein
              assign("tp",protein(myprotein,online=FALSE))
              if(nrow(tp)!=0) {
                protein.use <- unique(c(protein.use,as.numeric(rownames(tp))))
                cat(paste("WWW: added",myprotein,"\n"))
              }
              else if(myorg!="NEW") cat(paste("WWW: protein",myprotein,"is unknown\n"))
            }
          } else {
            if(length(grep("_",myprotein))==0) assign("myprotein",paste(myprotein,"NEW",sep="_"))
            assign("tryprotein",try(protein(mysequence,myprotein),TRUE))
            if(class(tryprotein)=="try-error") {
              cat(paste("WWW: adding protein",myprotein,"failed; the message was\n"))
              cat(tryprotein[1])
            } else {
              protein.use <- unique(c(protein.use,as.numeric(rownames(tryprotein))))
            }
          }
        }
      
        // refresh the species list after ProcessProteins
        dojo.event.connect(dojo.widget.byId("ProcessProteins"), "onReceive",  
           dojo.widget.byId("SelectSpecies"), "calculate");
        " ";
      

select species of interest

        if(!exists("SelectSpeciesNo")) assign("SelectSpeciesNo",1) else 
          assign("SelectSpeciesNo",SelectSpeciesNo+1)
        cat(paste("WWW: SelectSpecies (",SelectSpeciesNo,")\n"))
        cat("reaction coefficient / species / logarithm of activity\n")
      ## the default species and logarithms of activities
        assign("n.species",12)
        # assemble the names of species (with states)
        assign("speciesnames",paste(thermo$obigt$name,thermo$obigt$state))
        assign("usedspeciesnames",paste(thermo$obigt$name,thermo$obigt$state)[data.use])
        if(nrow(thermo$protein) > 0) assign("usedspeciesnames",c(usedspeciesnames,
          paste(thermo$protein$protein[protein.use],"_",thermo$protein$organism[protein.use]," aq",sep="")))
        if(SelectSpeciesNo==1 | SelectSpeciesGo) {
          assign("default.species",match(speciesnames[default.species],usedspeciesnames))
          assign("SelectSpeciesGo",FALSE)
        } else {
          # keep the ones set by the user
          assign("default.species",numeric())
          assign("default.slogact",numeric())
          for(i in 1:n.species) if(exists(paste("species",i,sep=""))) {
            assign("myspecies",get(paste("species",i,sep="")))
            assign("myslogact",get(paste("slogact",i,sep="")))
            assign("mymatch",match(myspecies,usedspeciesnames))
            if(!is.na(mymatch)) {
              assign("default.species",c(default.species,mymatch))
              assign("default.slogact",c(default.slogact,myslogact))
            }
          }
        }
        # make sure they're long enough
        assign("default.species",c(default.species,rep(0,n.species)))
        assign("default.slogact",c(default.slogact,rep(0,n.species)))
        assign("default.coeff",c(default.coeff,rep(0,n.species)))
      ## make the form
        cat(HTMLon())
        for(i in 1:n.species) {
          # entry field for reaction coefficients
          cat(HTMLinput(paste("coeff",i,sep=""),value=default.coeff[i],size=5))
          # the names (species states) pull-down menu
          assign("name",paste("species",i,sep=""))
          cat(HTMLselect(name, c("",usedspeciesnames),default=default.species[i]+1))
          # entry field for reference activities
          cat(HTMLinput(paste("slogact",i,sep=""),value=default.slogact[i],size=4))
          cat(paste("<br>"))
        }
        H("input",type="button",value="GO",onclick="rpad.calculateTree(dojo.byId('SpeciesDiv'))")
      
        if(!exists("ProcessSpeciesNo")) assign("ProcessSpeciesNo",1) else 
          assign("ProcessSpeciesNo",ProcessSpeciesNo+1)
        cat(paste("WWW: ProcessSpecies (",ProcessSpeciesNo,")\n"))
        assign("myspecies",character())
        assign("mystates",character())
        assign("myslogact",character())
        assign("mycoeff",character())
        assign("n.species",as.numeric(n.species))
        for(i in 1:n.species) {
          if(exists(paste("species",i,sep=""))) {
            assign("thisspecies",get(paste("species",i,sep="")))
            if(!identical(thisspecies,"")) {
              assign("j",match(thisspecies,usedspeciesnames))
              assign("myspecies",c(myspecies,c(thermo$obigt$formula[data.use],
                paste(thermo$protein$protein[protein.use],thermo$protein$organism[protein.use],sep="_"))[j]))
              assign("mystates",c(mystates,c(thermo$obigt$state[data.use],
                rep("aq",nrow(thermo$protein)))[j]))
              assign("myslogact",c(myslogact,get(paste("slogact",i,sep=""))))
              assign("mycoeff",c(mycoeff,get(paste("coeff",i,sep=""))))
            }
          }
        }
        # now to apply the changes (and output the result)
        thermo$species <- NULL
        tryspecies <- try(capture.output(species(myspecies,mystates)),TRUE)
        if(class(tryspecies)=="try-error") {
          # cat("WWW: error in species definition; the output was\n")
          cat(tryspecies[1])
        } else species(1:length(myspecies),as.numeric(myslogact),quiet=FALSE)
      

setup and plot diagram

        if(!exists("SelectTPNo")) assign("SelectTPNo",1) else assign("SelectTPNo",SelectTPNo+1)
        cat(paste("WWW: SelectTP (",SelectTPNo,")\n"))
        HTMLon()
        # reference temperature and pressure
        HTMLinput("temperature",value=temperature,size=5)
        cat("Temperature, degC<br>")
        HTMLinput("pressure",value=pressure,size=5)
        cat("Pressure, bar")
        # substitute pH and Eh?
        cat(paste("<br>"))
        if(!exists("use.pH")) assign("use.pH",TRUE)
        if(!exists("use.Eh")) assign("use.Eh",TRUE)
        # HTMLcheckbox("use.pH",checked=use.pH,text="use pH for H+")
        # HTMLcheckbox("use.Eh",checked=use.Eh,text="use Eh for e-")
        cat(paste("<br>"))
      
        if(!exists("SelectVariablesNo")) assign("SelectVariablesNo",1) else 
          assign("SelectVariablesNo",SelectVariablesNo+1)
        cat(paste("WWW: SelectVariables (",SelectVariablesNo,")\n"))
        HTMLon()
        if(!is.null(thermo$basis)) {
          assign("myvar",rownames(thermo$basis))
          assign("myvar",c(myvar,"T","P"))
          if("H+" %in% myvar) if(use.pH) myvar[myvar=="H+"] <- "pH"
          if("e-" %in% myvar) if(use.Eh) myvar[myvar=="e-"] <- "Eh"
          # now some pull-down boxes and range settings
          cat("X,min,max")
          assign("xdefault",match(xvar,myvar))
          if(is.na(xdefault)) assign("xdefault",1)
          cat(HTMLselect("xvar",myvar,default=xdefault,id="xvarSelect"))
          cat(HTMLinput("xmin",value=xmin,size=5,id="xminInput"))
          cat(HTMLinput("xmax",value=xmax,size=5,id="xmaxInput"))
          cat(paste("<br>"))
          cat("Y,min,max")
          assign("ydefault",match(yvar,c(myvar,"")))
          if(is.na(ydefault)) assign("ydefault",1)
          cat(HTMLselect("yvar",c(myvar,""),default=ydefault,id="yvarSelect"))
          cat(HTMLinput("ymin",value=ymin,size=5,id="yminInput"))
          cat(HTMLinput("ymax",value=ymax,size=5,id="yminInput"))
          cat(paste("<br>"))
        } else {
          cat("WWW: no basis species ... skipping variable selection")
          cat(paste("<br>"))
        }
      
        dojo.event.connect(dojo.widget.byId("SelectVariables"), "onReceive",  
           dojo.widget.byId("ProcessSpecies"), "calculate");
        " ";
      
        if(!exists("SetupDiagramNo")) assign("SetupDiagramNo",1) else 
          assign("SetupDiagramNo",SetupDiagramNo+1)
        cat(paste("WWW: SetupDiagram (",SetupDiagramNo,")\n"))
        HTMLon()
        cat("all plots: ")
        HTMLcheckbox("show.graph",checked=TRUE,text="show graph ")
        if(use.residues) HTMLcheckbox("use.residues",checked=TRUE,
          text="use residues ") else HTMLcheckbox("use.residues",text="use residues ")
        HTMLcheckbox("draw.names",checked=TRUE,text="draw names ")
        cat("<br>")
        cat("2-D plots: ")
        if(fill.areas )HTMLcheckbox("fill.areas",checked=TRUE,
          text="fill areas ") else HTMLcheckbox("fill.areas",text="fill areas ")
        HTMLselect("line.color",text=palette(),default=1)
        cat("line color")
        cat("<br>")
        cat("1-D plots: ")
        HTMLselect("line.types",text=as.character(1:12),default=line.types)
        cat("number of line types ")
        cat("<br>")
        # the CALCULATE button
        #H("input",type="button",value="Calculate",onclick="javascript:rpad.calculatePage()")
        H("input",type="button",value="Calculate",onclick="rpad.calculateTree(dojo.byId('DiagramDiv'))")
      
        if(!exists("ShowDiagramNo")) assign("ShowDiagramNo",1) else 
          assign("ShowDiagramNo",ShowDiagramNo+1)
        cat(paste("WWW: ShowDiagram (",ShowDiagramNo,")\n"))
        # arguments for affinity()
        if(yvar==xvar) assign("yvar","")
        if(xmin=="") assign("xmin",-20)
        if(xmax=="") assign("xmax",20)
        if(yvar!="") {
          if(ymin=="") assign("ymin",-20)
          if(ymax=="") assign("ymax",20)
          assign("myargs",list(c(xmin,xmax),c(ymin,ymax)))
          names(myargs) <- c(xvar,yvar)
          assign("nd",2)
        } else {
          if(ymin=="") assign("ymin",-5)
          if(ymax=="") assign("ymax",0)
          assign("myargs",list(c(xmin,xmax)))
          names(myargs) <- c(xvar)
          assign("nd",1)
        }
        if(!"T" %in% c(xvar,yvar)) assign("myargs",c(myargs,list(T=temperature)))
        if(!"P" %in% c(xvar,yvar)) assign("myargs",c(myargs,list(P=pressure)))
        if(is.null(thermo$species)) {
          cat("WWW: error: no species present\n")
        } else {
          assign("a",do.call(affinity,myargs))
          if(is.null(a)) stop("WWW: something didn't work ... try again")
          if(draw.names) assign("legend.x","topright") else assign("legend.x",NULL)
          if(draw.names) assign("names",NA) else assign("names",NULL)
          assign("addtoplot",TRUE)
          assign("trypoints",try(points(0,0,pch=""),TRUE))
          if(class(trypoints)=="try-error") assign("addtoplot",FALSE) 
          # preparing the arguments for diagram
          if(nd == 1) {
            # 1-D plot (speciation diagram)
            assign("col",NULL)
            assign("line.types",as.numeric(line.types))
            assign("color",palette())
            assign("ncols",ceiling(length(a$values)/line.types))
            if(ncols > length(color)) assign("ncols",length(color))
            assign("color",rep(color[1:ncols],each=ceiling(length(a$values)/ncols)))
            assign("lty",rep(1:line.types,length.out=length(a$values)))
            if(ymin==0 & ymax==1) {
              assign("alpha",TRUE)
              assign("ylim",c(0,1))
            } else {
              assign("alpha",FALSE)
              assign("ylim",c(ymin,ymax))
            }
            assign("cex.names",0.7)
          } else {
            # 2-D plot (predominance diagram)
            assign("col",line.color)
            if(fill.areas) assign("color","heat") else assign("color",NULL)
            assign("lty",NULL)
            assign("alpha",FALSE)
            assign("ylim",c(0,1))
            assign("cex.names",1)
          }
          if(length(grep("_",myspecies))==0) assign("use.residues",FALSE)
          assign("d",try(diagram(a,col=col,color=color,legend.x=legend.x,
            add=addtoplot,residue=use.residues,cex.names=cex.names,
            lty=lty,alpha=alpha,ylim=ylim,names=names,col.names=col)))
          if(class(d)=="try-error") {
            cat("WWW: error in diagram function; the output was\n")
            cat(d[1])
            stop()
          }
          # water lines on Eh-pH and logfO2-pH diagrams
          if(xvar=='pH' & yvar %in% c('Eh','O2')) water.lines(yaxis=yvar,
            T=convert(temperature,'K'),P=pressure)
          write(capture.output(print(d)),file="diagram.txt")
          cat("opens in new window ")
          cat(HTMLon())
          cat(HTMLlink(RpadURL("affinity.txt"),"output from affinity()",target="_blank"))
          cat(HTMLoff())
          cat("
") cat("opens in new window ") cat(HTMLon()) cat(HTMLlink(RpadURL("diagram.txt"),"output from diagram()",target="_blank")) cat(HTMLoff()) cat("
") if(!show.graph) { cat("WWW: plot device kept open ... now define an overlay\n") } else { cat(HTMLon()) cat(showgraph(link=TRUE)) } write(capture.output(print(a)),file="affinity.txt") }

calculate properties

        if(!exists("SetupPropertyNo")) assign("SetupPropertyNo",1) else 
          assign("SetupPropertyNo",SetupPropertyNo+1)
        cat(paste("WWW: SetupProperty (",SetupPropertyNo,")\n"))
        HTMLon()
        # now some range settings
        cat("Tmin,Tmax")
        cat(HTMLinput("Tmin",value=Tmin,size=5))
        cat(HTMLinput("Tmax",value=Tmax,size=5))
        cat(paste("<br>"))
        cat("Pmin,Pmax")
        cat(HTMLinput("Pmin",value=Pmin,size=5))
        cat(HTMLinput("Pmax",value=Pmax,size=5))
        cat(paste("<br>"))
        # the CALCULATE button
        H("input",type="button",value="Calculate",onclick="rpad.calculateTree(dojo.byId('PropertyDiv'))")
      
        if(!exists("CalculatePropertyNo")) assign("CalculatePropertyNo",1) else 
          assign("CalculatePropertyNo",CalculatePropertyNo+1)
        cat(paste("WWW: CalculateProperty (",CalculatePropertyNo,")\n"))
        if(length(myspecies)==0) stop("WWW: select one or more species to calculate properties")
        assign("doit",TRUE)
        if(can.be.numeric(Tmin) & can.be.numeric(Tmax)) {
          if(Tmin==Tmax) assign("Targ",as.numeric(Tmin)) else
          assign("Targ",seq(as.numeric(Tmin),as.numeric(Tmax),length.out=11))
        } else {
          if(can.be.numeric(Tmin)) assign("Targ",as.numeric(Tmin)) else {
            if(can.be.numeric(Tmax)) assign("Targ",as.numeric(Tmax)) else {
              cat("WWW: temperature must be numeric\n")
              assign("doit",FALSE)
            }
          }
        }
        if(can.be.numeric(Pmin) & can.be.numeric(Pmax)) {
          if(Pmin==Pmax) assign("Parg",as.numeric(Pmin)) else
          assign("Parg",seq(as.numeric(Pmin),as.numeric(Pmax),length.out=11))
        } else {
          if(can.be.numeric(Pmin)) assign("Parg",as.numeric(Pmin)) else {
            if(can.be.numeric(Pmax)) assign("Parg",as.numeric(Pmax)) else {
              if(Pmin=="" & Pmax=="") assign("Parg","Psat") else {
                cat("WWW: pressure must be numeric\n")
                assign("doit",FALSE)
              }
            }
          }
        }
        if(doit) {
          if(!any(Parg!=0)) assign("Parg","Psat")
          assign("grid",NULL)
          if(length(Parg) > 1 & length(Targ) > 1) assign("grid","P")
          if(any(mycoeff!=0)) {
            assign("inrxn",which(mycoeff!=0))
            assign("s.out",subcrt(species=myspecies[inrxn],state=mystates[inrxn],
              T=Targ,P=Parg,grid=grid,do.phases=FALSE,logact=as.numeric(myslogact)[inrxn],
              coeff=as.numeric(mycoeff)[inrxn])) 
          } else {
            assign("s.out",subcrt(species=myspecies,state=mystates,
              T=Targ,P=Parg,grid=grid,do.phases=FALSE,logact=as.numeric(myslogact)))
          }
          print(s.out)
        }
        HTMLoff()