main.tex

% Computational Chemistry LaTeX Template for CoCalc
% Optimized for molecular dynamics, quantum chemistry, and drug discovery
% Keywords: computational chemistry latex template, molecular dynamics latex, quantum chemistry template
% DFT calculations latex, HOMO LUMO latex template, pymol figures latex, gaussian output latex

\documentclass[11pt,a4paper]{article}

% Essential packages for computational chemistry manuscripts
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{graphicx}
\usepackage{xcolor}
\usepackage{booktabs}
\usepackage{siunitx} % For proper scientific units
\usepackage{chemfig} % Chemical structure drawing
\usepackage[version=4]{mhchem} % Chemical equations and formulas
\usepackage{listings} % Code listings for computational methods
\usepackage{subcaption} % Multiple subfigures
\usepackage{float} % Better figure placement
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{cleveref}
\usepackage{natbib} % Bibliography management

% PythonTeX for embedded calculations and molecular visualization
\usepackage{pythontex}

% Geometry setup
\geometry{margin=2.5cm}

% Color definitions for molecular representations
\definecolor{carboncolor}{RGB}{64,64,64}
\definecolor{oxygencolor}{RGB}{255,13,13}
\definecolor{nitrogencolor}{RGB}{48,80,248}
\definecolor{hydrogencolor}{RGB}{255,255,255}
\definecolor{sulfurcolor}{RGB}{255,200,50}

% siunitx setup for chemistry units
\sisetup{
    per-mode=symbol,
    detect-all,
    group-separator={,},
    group-minimum-digits=4
}

% Hyperref setup
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    filecolor=magenta,
    urlcolor=cyan,
    citecolor=green,
    pdftitle={Computational Chemistry LaTeX Template},
    pdfauthor={CoCalc Scientific Templates},
    pdfsubject={Molecular Dynamics, Quantum Chemistry, Drug Discovery},
    pdfkeywords={computational chemistry, molecular dynamics, quantum chemistry, DFT, HOMO LUMO, drug discovery, gaussian, amber, gromacs, pymol, rdkit}
}

% Custom commands for common chemistry notations
\newcommand{\homo}{\text{HOMO}}
\newcommand{\lumo}{\text{LUMO}}
\newcommand{\dft}{\text{DFT}}
\newcommand{\mptwo}{\text{MP2}}
\newcommand{\ccsd}{\text{CCSD(T)}}
\newcommand{\hartree}{\text{E}_\text{h}}
\newcommand{\kcalmol}{\si{\kilo\cal\per\mol}}
\newcommand{\kjmol}{\si{\kJ\per\mol}}
\DeclareSIUnit\angstrom{\text{\AA}}

\title{Computational Chemistry Research Template:\\
Molecular Dynamics, Quantum Chemistry \& Drug Discovery}

\author{Your Name\thanks{Department of Chemistry, Institution Name} \\
\small \texttt{[e-mail]}}

\date{\today}

\begin{document}

\maketitle

\begin{abstract}
This computational chemistry template demonstrates advanced molecular modeling techniques including density functional theory (DFT) calculations, molecular dynamics (MD) simulations, and drug discovery applications. The template showcases quantum chemical property predictions, protein-ligand interactions, orbital visualizations, and energy landscape analysis using embedded Python calculations optimized for reproducible research in CoCalc.

\textbf{Keywords:} computational chemistry, molecular dynamics, quantum chemistry, DFT calculations, HOMO-LUMO gap, drug discovery, gaussian calculations, amber simulations, pymol visualization
\end{abstract}

\section{Introduction}
% SEO: computational chemistry latex template, molecular dynamics latex

Computational chemistry has revolutionized our understanding of molecular systems, enabling researchers to predict chemical properties, design new drugs, and understand complex biological processes at the atomic level. This template provides a comprehensive framework for presenting computational chemistry research, including:

\begin{itemize}
\item Quantum chemical calculations (\dft, \mptwo, \ccsd)
\item Molecular dynamics simulations (AMBER, GROMACS, CHARMM)
\item Drug discovery and virtual screening
\item Protein-ligand interaction analysis
\item Molecular visualization and property prediction
\end{itemize}

\section{Computational Methods}
% SEO: gaussian output latex template, dft calculations latex

\subsection{Quantum Chemical Calculations}

All quantum chemical calculations were performed using Gaussian 16~\cite{gaussian16} at the B3LYP/6-31G(d,p) level of theory. Geometry optimizations were followed by frequency calculations to confirm stationary points.

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import seaborn as sns

# Set style for publication-quality figures
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("husl")

# Example: HOMO-LUMO energy gap calculation
molecules = ['benzene', 'naphthalene', 'anthracene', 'phenanthrene']
homo_energies = np.array([-5.78, -5.42, -5.15, -5.25])  # eV
lumo_energies = np.array([-1.23, -2.15, -2.85, -2.45])  # eV
gap_energies = homo_energies - lumo_energies

# Create HOMO-LUMO gap plot
fig, ax = plt.subplots(figsize=(10, 6))
x_pos = np.arange(len(molecules))

# Plot HOMO and LUMO levels
homo_bars = ax.bar(x_pos - 0.2, homo_energies, 0.4, label='HOMO',
                   color='#FF6B6B', alpha=0.8)
lumo_bars = ax.bar(x_pos + 0.2, lumo_energies, 0.4, label='LUMO',
                   color='#4ECDC4', alpha=0.8)

# Add gap annotations
for i, gap in enumerate(gap_energies):
    ax.annotate(f'{gap:.2f} eV', xy=(i, (homo_energies[i] + lumo_energies[i])/2),
                ha='center', va='center', fontweight='bold', fontsize=10,
                bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0.8))

ax.set_xlabel('Polycyclic Aromatic Hydrocarbons', fontsize=12, fontweight='bold')
ax.set_ylabel('Energy (eV)', fontsize=12, fontweight='bold')
ax.set_title('HOMO-LUMO Energy Gaps in PAH Series', fontsize=14, fontweight='bold')
ax.set_xticks(x_pos)
ax.set_xticklabels(molecules, rotation=45, ha='right')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/homo_lumo_gaps.pdf', dpi=300, bbox_inches='tight')
plt.show()

print(f"Average HOMO-LUMO gap: {np.mean(gap_energies):.2f} ± {np.std(gap_energies):.2f} eV")
\end{pycode}

The \homo-\lumo gap is a crucial parameter for understanding electronic properties and reactivity. As shown in \cref{fig:homo_lumo}, the gap systematically decreases with increasing conjugation length in polycyclic aromatic hydrocarbons.

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{assets/homo_lumo_gaps.pdf}
\caption{HOMO-LUMO energy gaps calculated at the B3LYP/6-31G(d,p) level of theory for polycyclic aromatic hydrocarbons. The systematic decrease in gap energy with increasing molecular size demonstrates the effect of extended conjugation.}
\label{fig:homo_lumo}
\end{figure}

\subsection{Molecular Dynamics Simulations}
% SEO: amber molecular dynamics latex, protein structure latex

Molecular dynamics simulations were performed using AMBER 20~\cite{amber20} with the ff19SB force field for proteins and TIP3P water model. The system was equilibrated for \SI{10}{ns} followed by \SI{100}{ns} production runs.

\begin{pycode}
# Generate sample MD analysis data
import numpy as np
import matplotlib.pyplot as plt

# Simulate RMSD trajectory data
time = np.linspace(0, 100, 1000)  # 100 ns simulation
np.random.seed(42)  # Reproducible results

# RMSD with equilibration and fluctuations
rmsd_protein = 1.5 + 0.8 * (1 - np.exp(-time/10)) + 0.3 * np.random.normal(0, 1, len(time)) * np.exp(-time/50)
rmsd_ligand = 2.0 + 1.2 * (1 - np.exp(-time/8)) + 0.5 * np.random.normal(0, 1, len(time)) * np.exp(-time/30)

# Apply smoothing
from scipy.ndimage import gaussian_filter1d
rmsd_protein = gaussian_filter1d(rmsd_protein, sigma=2)
rmsd_ligand = gaussian_filter1d(rmsd_ligand, sigma=2)

# Create MD analysis plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

# RMSD plot
ax1.plot(time, rmsd_protein, linewidth=2, label='Protein Backbone', color='#2E86AB')
ax1.plot(time, rmsd_ligand, linewidth=2, label='Ligand Heavy Atoms', color='#A23B72')
ax1.set_ylabel('RMSD (Å)', fontsize=12, fontweight='bold')
ax1.set_title('Molecular Dynamics Trajectory Analysis', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(0, 4)

# Energy components
potential_energy = -45000 + 2000 * np.sin(time/10) + 500 * np.random.normal(0, 1, len(time))
kinetic_energy = 15000 + 1000 * np.sin(time/12 + np.pi/4) + 300 * np.random.normal(0, 1, len(time))

# Apply smoothing to energies
potential_energy = gaussian_filter1d(potential_energy, sigma=3)
kinetic_energy = gaussian_filter1d(kinetic_energy, sigma=3)

ax2.plot(time, potential_energy/1000, linewidth=2, label='Potential Energy', color='#F18F01')
ax2.plot(time, kinetic_energy/1000, linewidth=2, label='Kinetic Energy', color='#C73E1D')
ax2.set_xlabel('Time (ns)', fontsize=12, fontweight='bold')
ax2.set_ylabel('Energy (kcal/mol × 10³)', fontsize=12, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('assets/md_analysis.pdf', dpi=300, bbox_inches='tight')
plt.show()

# Calculate and print statistics
print(f"Final RMSD (protein): {rmsd_protein[-1]:.2f} ± {np.std(rmsd_protein[-100:]):.2f} Å")
print(f"Final RMSD (ligand): {rmsd_ligand[-1]:.2f} ± {np.std(rmsd_ligand[-100:]):.2f} Å")
print(f"Average potential energy: {np.mean(potential_energy)/1000:.1f} kcal/mol")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.9\textwidth]{assets/md_analysis.pdf}
\caption{Molecular dynamics trajectory analysis showing (top) root-mean-square deviation (RMSD) of protein backbone and ligand heavy atoms, and (bottom) potential and kinetic energy components over a \SI{100}{ns} simulation period.}
\label{fig:md_analysis}
\end{figure}

\section{Drug Discovery Applications}
% SEO: drug discovery latex template, virtual screening results

\subsection{Virtual Screening and Docking Analysis}

High-throughput virtual screening was performed against a library of \num{10000} drug-like compounds using AutoDock Vina~\cite{vina}. The binding affinities and poses were analyzed for structure-activity relationships.

\begin{pycode}
# Generate virtual screening results
import matplotlib.pyplot as plt
import numpy as np

np.random.seed(123)
n_compounds = 1000

# Generate docking scores with realistic distribution
binding_affinities = np.random.gamma(2, 2, n_compounds) * -1 - 4  # Gamma distribution for realistic scores
binding_affinities = np.clip(binding_affinities, -12, -4)  # Realistic range

# Create a subset of high-affinity compounds
high_affinity_indices = np.where(binding_affinities < -8)[0]
n_hits = len(high_affinity_indices)

# Generate ADMET properties for top hits
lipophilicity = np.random.normal(3.5, 1.2, n_hits)  # LogP
molecular_weight = np.random.normal(350, 80, n_hits)  # Da
polar_surface_area = np.random.normal(70, 25, n_hits)  # Ų

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))

# Docking score distribution
ax1.hist(binding_affinities, bins=50, alpha=0.7, color='#6A994E', edgecolor='black', linewidth=0.5)
ax1.axvline(x=-8, color='red', linestyle='--', linewidth=2, label='Hit Threshold')
ax1.set_xlabel('Binding Affinity (kcal/mol)', fontweight='bold')
ax1.set_ylabel('Number of Compounds', fontweight='bold')
ax1.set_title('Virtual Screening Results Distribution', fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# ADMET analysis - Lipophilicity vs Molecular Weight
scatter = ax2.scatter(molecular_weight, lipophilicity, c=binding_affinities[high_affinity_indices],
                     cmap='viridis_r', alpha=0.7, s=60, edgecolors='black', linewidth=0.5)
ax2.set_xlabel('Molecular Weight (Da)', fontweight='bold')
ax2.set_ylabel('Lipophilicity (LogP)', fontweight='bold')
ax2.set_title('ADMET Properties of High-Affinity Hits', fontweight='bold')
ax2.grid(True, alpha=0.3)
cbar = plt.colorbar(scatter, ax=ax2)
cbar.set_label('Binding Affinity (kcal/mol)', fontweight='bold')

# Rule of Five compliance
ro5_compliant = ((molecular_weight <= 500) &
                (lipophilicity <= 5) &
                (polar_surface_area <= 140))

labels = ['Compliant', 'Non-Compliant']
sizes = [np.sum(ro5_compliant), np.sum(~ro5_compliant)]
colors = ['#A8DADC', '#E63946']

ax3.pie(sizes, labels=labels, colors=colors, autopct='%1.1f%%', startangle=90)
ax3.set_title('Lipinski Rule of Five Compliance', fontweight='bold')

# SAR analysis - Structure-Activity Relationship
pharmacophore_features = np.random.poisson(3, n_hits)  # Number of pharmacophore features
ax4.scatter(pharmacophore_features, binding_affinities[high_affinity_indices],
           alpha=0.7, s=60, color='#F77F00', edgecolors='black', linewidth=0.5)
ax4.set_xlabel('Number of Pharmacophore Features', fontweight='bold')
ax4.set_ylabel('Binding Affinity (kcal/mol)', fontweight='bold')
ax4.set_title('Structure-Activity Relationship', fontweight='bold')
ax4.grid(True, alpha=0.3)

# Add trendline
z = np.polyfit(pharmacophore_features, binding_affinities[high_affinity_indices], 1)
p = np.poly1d(z)
ax4.plot(sorted(pharmacophore_features), p(sorted(pharmacophore_features)),
         "r--", alpha=0.8, linewidth=2)

plt.tight_layout()
plt.savefig('assets/virtual_screening.pdf', dpi=300, bbox_inches='tight')
plt.show()

print(f"Total compounds screened: {len(binding_affinities):,}")
print(f"Hits identified (≤ -8 kcal/mol): {n_hits}")
print(f"Hit rate: {(n_hits/len(binding_affinities)*100):.1f}%")
print(f"Rule of Five compliance: {np.sum(ro5_compliant)}/{n_hits} ({np.sum(ro5_compliant)/n_hits*100:.1f}%)")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.95\textwidth]{assets/virtual_screening.pdf}
\caption{Virtual screening analysis: (a) Distribution of binding affinities from molecular docking, (b) ADMET properties correlation for high-affinity hits, (c) Lipinski Rule of Five compliance, and (d) Structure-activity relationship analysis.}
\label{fig:virtual_screening}
\end{figure}

\subsection{Protein-Ligand Interaction Analysis}
% SEO: protein structure latex, solvation energy latex

The protein-ligand interactions were analyzed using hydrogen bonding, hydrophobic contacts, and electrostatic interactions. Key binding residues were identified through interaction frequency analysis.

\begin{pycode}
# Analyze protein-ligand interactions
import matplotlib.pyplot as plt
import numpy as np

# Simulation data for interaction analysis
residue_names = ['GLU123', 'ARG156', 'TYR189', 'PHE234', 'ASP267', 'LYS301', 'TRP345', 'HIS378']
hydrogen_bonds = np.array([85, 72, 45, 12, 78, 56, 23, 89])  # % occupancy
hydrophobic_contacts = np.array([23, 45, 89, 92, 34, 67, 95, 56])  # % occupancy
electrostatic = np.array([78, 89, 12, 5, 82, 76, 15, 67])  # % occupancy

# Create interaction heatmap
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Interaction frequency heatmap
interaction_data = np.array([hydrogen_bonds, hydrophobic_contacts, electrostatic])
interaction_types = ['H-bonds', 'Hydrophobic', 'Electrostatic']

im = ax1.imshow(interaction_data, cmap='YlOrRd', aspect='auto', vmin=0, vmax=100)
ax1.set_xticks(range(len(residue_names)))
ax1.set_xticklabels(residue_names, rotation=45, ha='right')
ax1.set_yticks(range(len(interaction_types)))
ax1.set_yticklabels(interaction_types)
ax1.set_title('Protein-Ligand Interaction Frequencies (%)', fontweight='bold', fontsize=12)

# Add percentage labels
for i in range(len(interaction_types)):
    for j in range(len(residue_names)):
        text = ax1.text(j, i, f'{interaction_data[i, j]:.0f}%',
                       ha="center", va="center", color="black" if interaction_data[i, j] < 50 else "white",
                       fontweight='bold', fontsize=9)

cbar1 = plt.colorbar(im, ax=ax1, shrink=0.8)
cbar1.set_label('Interaction Frequency (%)', fontweight='bold')

# Binding energy decomposition
residue_contributions = -(hydrogen_bonds * 0.05 + hydrophobic_contacts * 0.03 + electrostatic * 0.04)
total_binding_energy = np.sum(residue_contributions)

bars = ax2.bar(range(len(residue_names)), residue_contributions,
               color=['#E74C3C', '#3498DB', '#2ECC71', '#F39C12', '#9B59B6', '#1ABC9C', '#E67E22', '#34495E'])
ax2.set_xlabel('Binding Site Residues', fontweight='bold')
ax2.set_ylabel('Energy Contribution (kcal/mol)', fontweight='bold')
ax2.set_title('Per-Residue Binding Energy Decomposition', fontweight='bold')
ax2.set_xticks(range(len(residue_names)))
ax2.set_xticklabels(residue_names, rotation=45, ha='right')
ax2.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for i, v in enumerate(residue_contributions):
    ax2.text(i, v - 0.1, f'{v:.1f}', ha='center', va='top', fontweight='bold', fontsize=9)

plt.tight_layout()
plt.savefig('assets/protein_ligand_interactions.pdf', dpi=300, bbox_inches='tight')
plt.show()

print(f"Total calculated binding energy: {total_binding_energy:.1f} kcal/mol")
print(f"Most important residue: {residue_names[np.argmin(residue_contributions)]} ({residue_contributions[np.argmin(residue_contributions)]:.1f} kcal/mol)")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.95\textwidth]{assets/protein_ligand_interactions.pdf}
\caption{Protein-ligand interaction analysis: (left) Interaction frequency heatmap showing hydrogen bonding, hydrophobic contacts, and electrostatic interactions for key binding site residues, and (right) per-residue binding energy decomposition calculated from interaction frequencies.}
\label{fig:protein_interactions}
\end{figure}

\section{Chemical Reactions and Mechanisms}
% SEO: reaction pathway latex, transition state searches

\subsection{Reaction Pathway Analysis}

The reaction mechanism was investigated using transition state theory and intrinsic reaction coordinate (IRC) calculations. The energy profile reveals the rate-determining step and intermediate stability.

\begin{pycode}
# Generate reaction pathway energy profile
import matplotlib.pyplot as plt
import numpy as np

# Reaction coordinate and energies (in kcal/mol)
reaction_coord = np.array([0, 1, 2, 3, 4, 5, 6])
energies = np.array([0, 12.5, 25.8, 18.3, 22.1, 8.5, -15.2])  # Reactants to products

# Labels for reaction steps
step_labels = ['Reactants', 'TS1', 'Int1', 'TS2', 'Int2', 'TS3', 'Products']

fig, ax = plt.subplots(figsize=(12, 8))

# Plot energy profile
ax.plot(reaction_coord, energies, 'o-', linewidth=3, markersize=8, color='#2E86AB')

# Fill areas to highlight energy barriers
for i in range(len(reaction_coord)-1):
    if 'TS' in step_labels[i] or 'TS' in step_labels[i+1]:
        ax.fill_between(reaction_coord[i:i+2], energies[i:i+2], alpha=0.3, color='red')
    else:
        ax.fill_between(reaction_coord[i:i+2], energies[i:i+2], alpha=0.3, color='lightblue')

# Add energy labels
for i, (x, y) in enumerate(zip(reaction_coord, energies)):
    ax.annotate(f'{y:.1f}', xy=(x, y), xytext=(0, 15),
               textcoords='offset points', ha='center', fontweight='bold',
               bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0.8))

# Customize plot
ax.set_xlabel('Reaction Coordinate', fontsize=14, fontweight='bold')
ax.set_ylabel('Relative Energy (kcal/mol)', fontsize=14, fontweight='bold')
ax.set_title('Reaction Energy Profile: Diels-Alder Cycloaddition', fontsize=16, fontweight='bold')
ax.set_xticks(reaction_coord)
ax.set_xticklabels(step_labels, rotation=45, ha='right')
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color='black', linestyle='--', alpha=0.5)

# Add reaction scheme using chemical equations
reaction_text = r'$\mathrm{C_4H_6 + C_2H_4 \rightarrow C_6H_{10}}$'
ax.text(0.02, 0.98, reaction_text, transform=ax.transAxes, fontsize=14,
        verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

plt.tight_layout()
plt.savefig('assets/reaction_pathway.pdf', dpi=300, bbox_inches='tight')
plt.show()

# Calculate thermodynamic and kinetic parameters
activation_energy = max(energies) - energies[0]
reaction_energy = energies[-1] - energies[0]

print(f"Activation energy (forward): {activation_energy:.1f} kcal/mol")
print(f"Reaction energy: {reaction_energy:.1f} kcal/mol")
print(f"Rate-determining step: {step_labels[np.argmax(energies)]} ({max(energies):.1f} kcal/mol)")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.85\textwidth]{assets/reaction_pathway.pdf}
\caption{Reaction energy profile for the Diels-Alder cycloaddition calculated at the B3LYP/6-31G(d,p) level. Transition states (TS) and intermediates (Int) are labeled with their relative energies. The reaction is thermodynamically favorable with $\Delta E = \SI{-15.2}{\kcalmol}$.}
\label{fig:reaction_pathway}
\end{figure}

\subsection{Molecular Orbital Analysis}

The frontier molecular orbitals (\homo\ and \lumo) control the reactivity and electronic properties. The orbital energies and symmetries determine the feasibility of chemical reactions.

\begin{pycode}
# Create molecular orbital energy diagram
import matplotlib.pyplot as plt
import numpy as np

fig, ax = plt.subplots(figsize=(10, 8))

# Molecular orbital energies (eV) for reactants and product
reactant1_orbitals = [-8.2, -7.1, -6.8, -5.9, -5.2, -1.8, -0.9, 0.3, 1.2]  # Diene
reactant2_orbitals = [-9.1, -8.3, -6.2, -2.1, 0.8, 1.5, 2.3]  # Dienophile
product_orbitals = [-9.8, -8.9, -7.8, -7.2, -6.5, -5.8, -5.1, -2.5, -1.9, -0.8, 0.5, 1.4, 2.1]

# Plot energy levels
x_positions = [1, 3, 5]
x_labels = ['Diene\n(C₄H₆)', 'Dienophile\n(C₂H₄)', 'Product\n(C₆H₁₀)']

# Draw energy levels
for i, orbitals in enumerate([reactant1_orbitals, reactant2_orbitals, product_orbitals]):
    x = x_positions[i]
    for j, energy in enumerate(orbitals):
        line_color = 'red' if energy > 0 else 'blue'  # LUMO vs HOMO
        line_width = 3 if abs(energy) < 2.5 else 2  # Highlight frontier orbitals
        ax.hlines(energy, x-0.3, x+0.3, colors=line_color, linewidth=line_width)

        # Label HOMO and LUMO
        if i < 2:  # Only for reactants
            if energy == max([e for e in orbitals if e < 0]):  # HOMO
                ax.text(x+0.4, energy, 'HOMO', fontsize=10, fontweight='bold', va='center')
            elif energy == min([e for e in orbitals if e > 0]):  # LUMO
                ax.text(x+0.4, energy, 'LUMO', fontsize=10, fontweight='bold', va='center')

# Add arrows showing orbital interactions
ax.annotate('', xy=(2.7, -1.8), xytext=(1.3, -5.2),  # HOMO-LUMO interaction
            arrowprops=dict(arrowstyle='<->', color='green', lw=2))
ax.text(2, -3.5, 'HOMO-LUMO\nInteraction', ha='center', fontsize=10,
        color='green', fontweight='bold')

# Customize plot
ax.set_xlim(0, 6)
ax.set_ylim(-12, 3)
ax.set_ylabel('Energy (eV)', fontsize=14, fontweight='bold')
ax.set_title('Molecular Orbital Energy Diagram: Diels-Alder Reaction', fontsize=16, fontweight='bold')
ax.set_xticks(x_positions)
ax.set_xticklabels(x_labels, fontsize=12, fontweight='bold')

# Add energy scale
ax.axhline(y=0, color='black', linestyle='--', alpha=0.5, linewidth=1)
ax.text(5.5, 0.2, 'Vacuum Level', fontsize=10, ha='right')

# Add legend
ax.plot([], [], 'r-', linewidth=3, label='LUMO (unoccupied)')
ax.plot([], [], 'b-', linewidth=3, label='HOMO (occupied)')
ax.legend(loc='upper right', fontsize=11)

plt.tight_layout()
plt.savefig('assets/molecular_orbitals.pdf', dpi=300, bbox_inches='tight')
plt.show()

# Calculate orbital gap
diene_gap = min([e for e in reactant1_orbitals if e > 0]) - max([e for e in reactant1_orbitals if e < 0])
dienophile_gap = min([e for e in reactant2_orbitals if e > 0]) - max([e for e in reactant2_orbitals if e < 0])

print(f"Diene HOMO-LUMO gap: {diene_gap:.1f} eV")
print(f"Dienophile HOMO-LUMO gap: {dienophile_gap:.1f} eV")
\end{pycode}

\begin{figure}[H]
\centering
\includegraphics[width=0.85\textwidth]{assets/molecular_orbitals.pdf}
\caption{Molecular orbital energy diagram for the Diels-Alder cycloaddition showing the frontier orbital interactions between diene and dienophile. The HOMO-LUMO gap determines the activation energy and reaction feasibility.}
\label{fig:molecular_orbitals}
\end{figure}

\section{Results and Discussion}

\subsection{Quantum Chemical Properties}

The calculated properties demonstrate excellent agreement with experimental values where available. The \dft\ calculations predict:

\begin{itemize}
\item \homo-\lumo\ gaps ranging from \SIrange{3.5}{4.8}{\eV} for the aromatic series
\item Binding energies of \SIrange{-8.5}{-12.3}{\kcalmol} for high-affinity ligands
\item Activation barriers of \SI{25.8}{\kcalmol} for the cycloaddition reaction
\end{itemize}

\subsection{Molecular Dynamics Insights}

The MD simulations reveal stable protein-ligand complexes with average RMSD values below \SI{2.5}{\angstrom}. Key findings include:

\begin{enumerate}
\item Hydrogen bonding dominates binding affinity (GLU123, HIS378)
\item Hydrophobic interactions provide selectivity (PHE234, TRP345)
\item Electrostatic complementarity stabilizes the complex
\end{enumerate}

\subsection{Drug Design Implications}

The virtual screening identified promising lead compounds with:
\begin{itemize}
\item High binding affinity ($\leq$ \SI{-8}{\kcalmol})
\item Favorable ADMET properties
\item Good Lipinski Rule of Five compliance (\SI{78}{\percent})
\end{itemize}

\section{Computational Details}

\begin{table}[H]
\centering
\caption{Summary of computational methods and software used.}
\label{tab:computational_methods}
\begin{tabular}{@{}lll@{}}
\toprule
\textbf{Method} & \textbf{Software} & \textbf{Key Parameters} \\
\midrule
DFT Calculations & Gaussian 16 & B3LYP/6-31G(d,p), Freq \\
MD Simulations & AMBER 20 & ff19SB, TIP3P, 100 ns \\
Molecular Docking & AutoDock Vina & Exhaustiveness = 8 \\
Virtual Screening & Schrödinger Suite & HTVS, SP, XP protocols \\
Visualization & PyMOL, VMD & Ray tracing, publication quality \\
Analysis & Python/NumPy & matplotlib, seaborn, scipy \\
\bottomrule
\end{tabular}
\end{table}

\section{Conclusion}

This computational chemistry template demonstrates the integration of quantum chemical calculations, molecular dynamics simulations, and drug discovery applications. The reproducible workflows and embedded calculations enable transparent and verifiable research outcomes.

Future extensions could include:
\begin{itemize}
\item Machine learning property prediction
\item Free energy perturbation calculations
\item QM/MM hybrid methods
\item Enhanced sampling techniques
\end{itemize}

\section*{Data and Code Availability}

All calculation inputs, outputs, and analysis scripts are available in the project repository. The embedded Python code ensures full reproducibility of results and figures.

\bibliographystyle{unsrt}
\begin{thebibliography}{99}

\bibitem{gaussian16}
M. J. Frisch, G. W. Trucks, H. B. Schlegel, et al.,
\textit{Gaussian 16, Revision C.01},
Gaussian, Inc., Wallingford CT, 2016.

\bibitem{amber20}
D.A. Case, H.M. Aktulga, K. Belfon, et al.,
\textit{AMBER 2020},
University of California, San Francisco, 2020.

\bibitem{vina}
O. Trott and A. J. Olson,
\textit{AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading},
J. Comput. Chem. \textbf{31}, 455-461 (2010).

\end{thebibliography}

\end{document}