# Load default settings for CHNOSZ
reset()

# 20220206
info <- "Database is set up correctly"
# Activate the AD model for aqueous nonelectrolytes
iAD <- add.OBIGT("AD")
# Get the indices of the modified aqueous species
iaq <- iAD[info(iAD)$state == "aq"]
# Each aqueous species should be associated with the AD model
expect_equal(unique(info(iaq)$model), "AD", info = info)
# Each aqueous species should have a gaseous counterpart
igas <- info(info(iAD)$name, "gas")
expect_true(!any(is.na(igas)), info = info)

# Test added on 20260106
# NOTE: This and the analogous code in demo/AD.R appears to be
# the only code in CHNOSZ that uses LVSeqn, LVSsat, etc. in H2O92D.f
info <- "We can get close to the critical point without NAs"
expect_silent(sres <- subcrt("CO2", "aq", T = 370:373, P = "Psat"), info = info)
expect_false(any(is.na(sres$out$CO2$G)), info = info)
expect_equal(round(sres$out$CO2$G[4], 1), -459810.9, info = info)

# First version of tests used a circular reference (values previously calculated in CHNOSZ) 20190220
# Tests now use values calculated with the AD_full program provided by N. Akinfiev and E. Bastrakov 20210206
info <- "AD produces correct values for CO2 along saturation curve"
P <- "Psat"
T <- c(50, 150, 250, 350)
# J mol-1
G_ref1 <- c(-389128.8, -405167.1, -425423.6, -450572.9)
# J K-1 mol-1
S_ref1 <- c(135.082, 183.295, 226.32, 366.626)
Cp_ref1 <- c(189.791, 178.452, 273.662, 7231.723)
V_ref1 <- c(32.57, 38.125, 58.269, 298.659)
# Calculate values using AD model in CHNOSZ
sout1 <- subcrt("CO2", T = T, P = P)$out[[1]]
expect_equal(sout1$G, G_ref1, tolerance = 13, scale = 1, info = info)
expect_equal(sout1$S, S_ref1, tolerance = 0.8, scale = 1, info = info)
expect_equal(sout1$V[1:3], V_ref1[1:3], tolerance = 0.11, scale = 1, info = info)
# Cp and V get much larger, and so do the differences, near the critical point
expect_equal(sout1$V[4], V_ref1[4], tolerance = 24, scale = 1, info = info)

info <- "AD produces correct values for CO2 at 1000 bar"
P <- 1000
T <- c(300, 400, 500, 600)
G_ref2 <- c(-431672.1, -455065.1, -480692, -507005.3)
S_ref2 <- c(221.255, 246.139, 263.759, 259.812)
Cp_ref2 <- c(150.613, 157.014, 54.489, -76.384)
V_ref2 <- c(45.831, 67.408, 107.656, 129.264)
# Calculate values using AD model in CHNOSZ
sout2 <- subcrt("CO2", T = T, P = P)$out[[1]]
expect_equal(sout2$G, G_ref2, tolerance = 11, scale = 1, info = info)
expect_equal(sout2$S, S_ref2, tolerance = 0.1, scale = 1, info = info)
expect_equal(sout2$V, V_ref2, tolerance = 0.4, scale = 1, info = info)

info <- "AD gives consistent values of G, H, and S"
T_in_Kelvin <- convert(T, "K")
S_of_elements_in_Joules <- entropy("CO2")
expect_equal(sout2$H - T_in_Kelvin * sout2$S + 298.15 * S_of_elements_in_Joules, sout2$G, info = info)

# 20220206
info <- "Fugacity, density, and density derivatives of H2O are close to values in Akinfiev and Diamond (2003)"
# This is the value of fugacity given in the paper
expect_equal(log(CHNOSZ:::.f1(298.15, 1, TRUE)), -3.4773, tolerance = 0.026, scale = 1, info = info)
# This is the value of fugacity output by AD_full
expect_equal(log(CHNOSZ:::.f1(298.15, 1, TRUE)), -3.4524, tolerance = 0.0007, scale = 1, info = info)
# The following values are from the paper
# g / cm3
expect_equal(CHNOSZ:::.rho1(298.15, 1),        0.9970,      tolerance = 0.0001, scale = 1, info = info)
# g / cm3 / K
expect_equal(CHNOSZ:::.drho1_dT(298.15, 1, FALSE),   -0.0002571,   tolerance = 0.000002, scale = 1, info = info)
# g / cm3 / bar
expect_equal(CHNOSZ:::.drho1_dP(298.15, 1, FALSE),    0.00004511,  tolerance = 0.0000001, scale = 1, info = info)

# 20190220 Compare Gibbs energies at 25 degrees calculated with AD model to default OBIGT database
info <- "Gibbs energies at 25 degree C are comparable between AD and default OBIGT database"
# Remove some hydroxides that aren't in the default database
iaq <- iaq[!info(iaq)$name %in% c("Si(OH)4", "As(OH)3")]
# This would produce an error if any calculations failed
# (e.g. because gases corresponding to any aqueous species were unavailable)
sAD <- subcrt(iaq, T = 25)
GAD <- do.call(rbind, sAD$out)$G
# Now calculate the parameters using default OBIGT database
OBIGT()
sOBIGT <- subcrt(iaq, T = 25)
GOBIGT <- do.call(rbind, sOBIGT$out)$G
# Mean absolute differences are less than 280 J / mol
expect_equal(GAD, GOBIGT, tolerance = 280, scale = 1, info = info)
# The largest differences are for HCl, ethane, and B(OH)3
expect_equal(sort(info(iaq[abs(GAD - GOBIGT) > 900])$name), sort(c("HCl", "ethane", "B(OH)3")))

## The remaining tests fail on some platforms (including CRAN machines)
exit_file("Skipping tests that fail on some platforms")

# Results from R-Universe MacOS arm64  20260616
#  ----- FAILED[attr]: test-AD.R<98--98>
#   call| expect_equal(sout1$Cp[4], Cp_ref1[4], tolerance = 800, scale = 1,
#   call| -->    info = "AD produces correct values for CO2 along saturation curve")
#   diff| Expected '7231.723', got '6072.37863187773'
#   info| AD produces correct values for CO2 along saturation curve
#  ----- FAILED[data]: test-AD.R<104--105>
#   call| expect_equal(CHNOSZ:::.d2rho1_dT2(298.15, 1, FALSE), -9.503e-06,
#   call| -->    tolerance = 7e-06, scale = 1, info = "Fugacity, density, and density derivatives of H2O are close to values in Akinfiev and Diamond (2003)")
#   diff| Expected '-9.503e-06', got '0.000176838824242332'
#   info| Fugacity, density, and density derivatives of H2O are close to values in Akinfiev and Diamond (2003)

# Fails on ATLAS and M1mac on CRAN 20220210
expect_equal(sout1$Cp[4], Cp_ref1[4], tolerance = 800, scale = 1, info = "AD produces correct values for CO2 along saturation curve")
# Fails on ATLAS (also on GitHub actions 20251111)
expect_equal(sout2$Cp, Cp_ref2, tolerance = 14, scale = 1, info = "AD produces correct values for CO2 at 1000 bar")
# Fails on M1mac (units: g / cm3 / K^2)
expect_equal(CHNOSZ:::.d2rho1_dT2(298.15, 1, FALSE), -0.000009503, tolerance = 0.000007, scale = 1,
             info = "Fugacity, density, and density derivatives of H2O are close to values in Akinfiev and Diamond (2003)")
# Fails on Windows (especially on R-Forge and winbuilder) (tolerance = 9 works) 20220208
expect_equal(sout1$Cp[1:3], Cp_ref1[1:3], tolerance = 3, scale = 1, info = "AD produces correct values for CO2 along saturation curve")

