# 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")

# VSEPR geometry database with real molecular examples
vsepr_molecules <- data.frame(
  molecule = c("BeCl2", "BF3", "CH4", "NH3", "H2O", "SO2", "PF5", "SF6", "ClF3", "XeF4"),
  central_atom = c("Be", "B", "C", "N", "O", "S", "P", "S", "Cl", "Xe"),
  bonding_pairs = c(2, 3, 4, 3, 2, 2, 5, 6, 3, 4),
  lone_pairs = c(0, 0, 0, 1, 2, 1, 0, 0, 2, 2),
  total_pairs = c(2, 3, 4, 4, 4, 3, 5, 6, 5, 6),
  geometry = c("Linear", "Trigonal Planar", "Tetrahedral", "Trigonal Pyramidal", 
               "Bent", "Bent", "Trigonal Bipyramidal", "Octahedral", "T-shaped", "Square Planar"),
  ideal_angle = c(180, 120, 109.5, 109.5, 109.5, 120, 120, 90, 90, 90),
  actual_angle = c(180, 120, 109.5, 106.8, 104.5, 119, 120, 90, 87.5, 90),
  dipole_moment = c(0.0, 0.0, 0.0, 1.47, 1.85, 1.63, 0.0, 0.0, 0.6, 0.0),
  polarity = c("Nonpolar", "Nonpolar", "Nonpolar", "Polar", "Polar", "Polar", 
               "Nonpolar", "Nonpolar", "Polar", "Nonpolar")
)

# Calculate bond angle deviation from ideal
vsepr_molecules$angle_deviation <- vsepr_molecules$actual_angle - vsepr_molecules$ideal_angle

# Create comprehensive VSEPR analysis plot
p_vsepr <- ggplot(vsepr_molecules, aes(x = lone_pairs, y = angle_deviation)) +
  geom_point(aes(size = dipole_moment, color = geometry), alpha = 0.8) +
  geom_text(aes(label = molecule), hjust = -0.1, vjust = -0.1, size = 3) +
  geom_smooth(method = "lm", se = TRUE, alpha = 0.2, color = "purple") +
  scale_size_continuous(range = c(3, 8), name = "Dipole\nMoment (D)") +
  scale_color_brewer(type = "qual", palette = "Set3", name = "Molecular\nGeometry") +
  labs(
    title = "VSEPR Theory: Lone Pair Effects on Bond Angles",
    subtitle = "Gillespie's theory (1957) explains molecular shape through electron pair repulsion",
    x = "Number of Lone Pairs on Central Atom",
    y = "Bond Angle Deviation from Ideal (°)",
    caption = "Lone pairs cause greater repulsion than bonding pairs, compressing bond angles"
  ) +
  chemistry_theme +
  theme(legend.position = "right")

print(p_vsepr)

# Create geometry distribution analysis
geometry_summary <- vsepr_molecules %>%
  group_by(geometry, polarity) %>%
  summarise(
    count = n(),
    avg_dipole = mean(dipole_moment),
    avg_deviation = mean(angle_deviation),
    .groups = 'drop'
  ) %>%
  arrange(desc(count))

cat("\n Molecular Geometry Analysis:\n")
cat("================================\n")
print(geometry_summary)

# Correlation between lone pairs and molecular polarity
polar_correlation <- cor(vsepr_molecules$lone_pairs, vsepr_molecules$dipole_moment)
cat(sprintf("\n Lone Pairs - Dipole Moment Correlation: %.3f\n", polar_correlation))
cat(" Lone pairs create molecular asymmetry, leading to polarity!\n")