# Setup: Load essential R packages for chemical analysis

# Avoid interactive CRAN prompts
options(repos = c(CRAN = "https://cloud.r-project.org"))

required_packages <- c("ggplot2", "dplyr", "plotly", "corrplot", "reshape2", "RColorBrewer")

# Install missing packages quietly
missing <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing)) install.packages(missing, quiet = TRUE)

# Attach packages without printing list results or masking chatter
suppressPackageStartupMessages({
  for (pkg in required_packages) {
    suppressWarnings(library(pkg, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))
  }
})
# Alternatively:
# invisible(lapply(required_packages, function(pkg)
#   suppressWarnings(library(pkg, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE))
# ))

# Theme fix: use linewidth instead of size to avoid ggplot2 deprecation warning
chemistry_theme <- ggplot2::theme_minimal() +
  ggplot2::theme(
    plot.title = ggplot2::element_text(size = 16, face = "bold", color = "#2E86AB"),
    plot.subtitle = ggplot2::element_text(size = 12, color = "#A23B72"),
    axis.title = ggplot2::element_text(size = 12, face = "bold"),
    axis.text = ggplot2::element_text(size = 10),
    legend.title = ggplot2::element_text(size = 11, face = "bold"),
    panel.grid.minor = ggplot2::element_blank(),
    panel.border = ggplot2::element_rect(color = "gray80", fill = NA, linewidth = 0.5)
  )

cat("Chemical Analysis Toolkit Loaded Successfully!\n")
cat("Ready to explore the molecular world with R\n")

# Comprehensive intermolecular forces and properties database
imf_compounds <- data.frame(
  compound = c("He", "Ne", "Ar", "Kr", "Xe", "H2", "N2", "O2", "CO2", "CH4", 
               "C2H6", "C3H8", "C4H10", "HF", "HCl", "HBr", "HI", "H2O", "NH3", "CH3OH"),
  molecular_weight = c(4.0, 20.2, 39.9, 83.8, 131.3, 2.0, 28.0, 32.0, 44.0, 16.0,
                       30.1, 44.1, 58.1, 20.0, 36.5, 80.9, 127.9, 18.0, 17.0, 32.0),
  boiling_point = c(-269, -246, -186, -153, -108, -253, -196, -183, -78, -164,
                    -89, -42, -1, 20, -85, -67, -35, 100, -33, 65), # °C
  dipole_moment = c(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                    0.0, 0.08, 0.05, 1.83, 1.08, 0.82, 0.45, 1.85, 1.47, 1.70),
  h_bonding = c("No", "No", "No", "No", "No", "No", "No", "No", "No", "No",
                "No", "No", "No", "Yes", "No", "No", "No", "Yes", "Yes", "Yes"),
  primary_imf = c("London", "London", "London", "London", "London", "London", 
                  "London", "London", "London", "London", "London", "London", 
                  "London", "H-bonding", "Dipole-Dipole", "Dipole-Dipole", 
                  "Dipole-Dipole", "H-bonding", "H-bonding", "H-bonding"),
  polarizability = c(0.2, 0.4, 1.6, 2.5, 4.0, 0.8, 1.7, 1.6, 2.6, 2.6,
                     4.5, 6.3, 8.2, 0.8, 2.6, 3.6, 5.4, 1.5, 2.3, 3.2) # ų
)

# Create comprehensive IMF analysis visualization
p_imf <- ggplot(imf_compounds, aes(x = molecular_weight, y = boiling_point)) +
  geom_point(aes(size = polarizability, color = primary_imf), alpha = 0.8) +
  geom_text(aes(label = compound), hjust = -0.1, vjust = -0.1, size = 2.8) +
  geom_smooth(method = "loess", se = TRUE, alpha = 0.2, color = "purple") +
  scale_size_continuous(range = c(2, 8), name = "Polarizability\n(ų)") +
  scale_color_manual(values = c("London" = "blue", "Dipole-Dipole" = "orange", 
                               "H-bonding" = "red"), name = "Primary\nIMF") +
  labs(
    title = "Intermolecular Forces: Structure-Property Relationships",
    subtitle = "Molecular weight, polarizability, and IMF type determine boiling points",
    x = "Molecular Weight (g/mol)",
    y = "Boiling Point (°C)",
    caption = "H-bonding creates anomalously high boiling points for small molecules"
  ) +
  chemistry_theme

print(p_imf)

# Analyze IMF effects on boiling points
imf_analysis <- imf_compounds %>%
  group_by(primary_imf) %>%
  summarise(
    count = n(),
    avg_bp = mean(boiling_point),
    avg_mw = mean(molecular_weight),
    avg_dipole = mean(dipole_moment),
    bp_range = max(boiling_point) - min(boiling_point),
    .groups = 'drop'
  ) %>%
  arrange(desc(avg_bp))

cat("\n Intermolecular Force Analysis:\n")
cat("==================================\n")
print(imf_analysis)

# Special analysis for hydrogen bonding anomalies
h_bonding_compounds <- imf_compounds[imf_compounds$h_bonding == "Yes", ]
non_h_bonding <- imf_compounds[imf_compounds$h_bonding == "No" & 
                              imf_compounds$molecular_weight < 50, ]

avg_bp_h_bonding <- mean(h_bonding_compounds$boiling_point)
avg_bp_non_h <- mean(non_h_bonding$boiling_point)

cat(sprintf("\n H-bonding Effect Analysis:\n"))
cat(sprintf("Average BP with H-bonding: %.1f°C\n", avg_bp_h_bonding))
cat(sprintf("Average BP without H-bonding (MW<50): %.1f°C\n", avg_bp_non_h))
cat(sprintf("H-bonding elevation: %.1f°C\n", avg_bp_h_bonding - avg_bp_non_h))
cat(" Hydrogen bonding dramatically increases boiling points!\n")