r"""
Debye's permanent dipole moments: polar and nonpolar molecules from geometry
============================================================================

Peter Debye (1912) explained why some gases have a dielectric constant
that falls with temperature: their molecules carry a permanent electric
dipole moment :math:`\boldsymbol\mu = \sum_i q_i \mathbf r_i`, which thermal
motion partly disorders. The measured moment tests a proposed shape.
Water's 1.85 D rules out a linear H-O-H, while carbon dioxide's zero
moment requires a linear O=C=O.

This example builds five molecules, adds up their bond dipoles with
:func:`~chemistrykit.structure.systems.dipole.bond_dipole_sum`, and
checks which ones cancel. It then uses
:func:`~chemistrykit.structure.systems.dipole.dipole_moment` to show how
water's moment would change with its bond angle. The unit, the debye
(:math:`10^{-18}` esu cm, about 0.208 e Å), is named after him.
"""

# %%
import matplotlib.pyplot as plt
import numpy as np

from chemistrykit.structure.systems.dipole import E_ANGSTROM_IN_DEBYE, bond_dipole_sum, dipole_moment
from chemistrykit.structure.systems.vsepr import build_vsepr_molecule

# Approximate bond-dipole magnitudes in debye. Every bond in a molecule is
# the same kind, so each is taken to point along its bond from the central
# atom; the direction convention does not change whether they cancel.
molecules = {
    "H2O (bent)": (build_vsepr_molecule(4, 2, bond_length=0.958, central_symbol="O", ligand_symbol="H"), 1.51),
    "NH3 (pyramidal)": (build_vsepr_molecule(4, 1, bond_length=1.012, central_symbol="N", ligand_symbol="H"), 1.31),
    "CO2 (linear)": (build_vsepr_molecule(2, 0, bond_length=1.16, central_symbol="C", ligand_symbol="O"), 2.3),
    "BF3 (trigonal planar)": (build_vsepr_molecule(3, 0, bond_length=1.31, central_symbol="B", ligand_symbol="F"), 1.6),
    "CH4 (tetrahedral)": (build_vsepr_molecule(4, 0, bond_length=1.09, central_symbol="C", ligand_symbol="H"), 0.4),
}
moments = {}
for name, (mol, mu_bond) in molecules.items():
    vectors = mol.coordinates[1:] - mol.coordinates[0]
    moments[name] = float(np.linalg.norm(bond_dipole_sum(vectors, [mu_bond] * len(vectors))))
    print(f"{name:22s}: {len(vectors)} bond dipoles of {mu_bond} D -> molecular dipole {moments[name]:.2f} D")
assert moments["CO2 (linear)"] < 1e-9 and moments["BF3 (trigonal planar)"] < 1e-9 and moments["CH4 (tetrahedral)"] < 1e-9

# %%
# VSEPR places water's hydrogens at the ideal 109.47 degrees, which gives
# 1.74 D. At the measured 104.5 degrees the vector sum rises to 1.85 D, the
# experimental value. Water's moment as a function of angle, from partial
# charges of -2q on O and +q on each H chosen to give 1.51 D per O-H bond:

r_oh = 0.958
q = 1.51 / (E_ANGSTROM_IN_DEBYE * r_oh)
angles = np.linspace(90.0, 180.0, 91)
mu = []
for theta in np.radians(angles):
    coords = [[0, 0, 0], [r_oh * np.sin(theta / 2), 0, r_oh * np.cos(theta / 2)], [-r_oh * np.sin(theta / 2), 0, r_oh * np.cos(theta / 2)]]
    mu.append(np.linalg.norm(dipole_moment(coords, [-2 * q, q, q])))
mu_at_1045 = np.interp(104.5, angles, mu)
print(f"\nwater at 104.5 deg: {mu_at_1045:.2f} D (measured 1.85 D); linear water would have {mu[-1]:.2f} D")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
ax1.barh(list(moments), list(moments.values()), color=["C3" if v > 1e-6 else "C7" for v in moments.values()])
ax1.set_xlabel("dipole moment (D)")
ax1.set_title("Bond dipoles cancel in symmetric shapes")
ax2.plot(angles, mu, color="C0")
ax2.axvline(104.5, color="gray", linestyle="--", linewidth=0.8)
ax2.axhline(1.85, color="C3", linestyle=":", linewidth=1, label="measured 1.85 D")
ax2.set_xlabel("H-O-H angle (degrees)")
ax2.set_ylabel("dipole moment (D)")
ax2.set_title("Water's dipole depends on its shape")
ax2.legend()
fig.tight_layout()
plt.show()