chemistrykit.structure#

chemistrykit.structure: molecular structure and bonding.

A lightweight Molecule container (atoms + 3D coordinates + bond list, no external file-format parsing); VSEPR geometry prediction from steric number, generating real 3D coordinates for the five idealized electron-domain polyhedra; point- group determination from 3D coordinates (genuine symmetry-element detection via distance/angle checks, not a formula lookup) plus character tables for the common point groups; bond order from the Pauling bond-length correlation and from Huckel-theory MO coefficients (the Coulson bond order); and formal-charge/oxidation-state assignment from a Lewis structure, using electronegativities from chemistrykit.periodic_table; Kekule-structure enumeration, Baeyer ring angle strain, and dipole moments from point charges or bond dipoles.

class chemistrykit.structure.LewisStructure(symbols, bond_orders, lone_pairs=<factory>)[source]#

Bases: object

A Lewis structure: element symbols, bond orders, and lone-pair counts.

Parameters:
  • symbols (list) – Element symbol of each atom (main-group elements only – see chemistrykit.periodic_table.valence_electrons()), 0-indexed.

  • bond_orders (dict) – Bond order for each bonded atom pair, e.g. {(0, 1): 1, (0, 2): 2} for one single and one double bond from atom 0. Each pair should appear once (order doesn’t matter, (i, j) and (j, i) are treated identically).

  • lone_pairs (dict) – Number of lone pairs (not electrons) on each atom, defaulting to 0 for any atom not given explicitly.

Examples

Water: O has 2 single bonds (to atoms 1, 2) and 2 lone pairs.

>>> water = LewisStructure(symbols=["O", "H", "H"], bond_orders={(0, 1): 1, (0, 2): 1}, lone_pairs={0: 2})
>>> water.formal_charges()
{0: 0.0, 1: 0.0, 2: 0.0}
>>> water.oxidation_states()
{0: -2.0, 1: 1.0, 2: 1.0}
bond_orders: dict#
formal_charges()[source]#

Compute the formal charge \(FC=V-N-B/2\) for every atom.

Return type:

dict

Returns:

dict of int to float

Examples

Ammonium, \(NH_4^+\): nitrogen has 4 single bonds and no lone pairs, giving the textbook +1 formal charge:

>>> nh4_plus = LewisStructure(symbols=["N", "H", "H", "H", "H"], bond_orders={(0, k): 1 for k in (1, 2, 3, 4)})
>>> nh4_plus.formal_charges()[0]
1.0
lone_pairs: dict#
oxidation_states()[source]#

Compute the oxidation state of every atom via the electronegativity-based algorithm.

For each bond, both bonding electrons (per unit of bond order) are assigned entirely to the more electronegative atom (split evenly for a homonuclear bond); nonbonding electrons always stay with their own atom. See the module docstring for the full formula.

Return type:

dict

Returns:

dict of int to float

Raises:

KeyError – If any atom’s element has no tabulated Pauling electronegativity (see chemistrykit.periodic_table.electronegativity()).

Examples

Hydrogen peroxide, \(H_2O_2\) (H-O-O-H): the O-O bond is homonuclear (no net electron shift), each O keeps its lone pairs and wins its one O-H bond, giving the well-known -1 oxidation state (intermediate between water’s -2 and O2’s 0):

>>> h2o2 = LewisStructure(symbols=["H", "O", "O", "H"], bond_orders={(0, 1): 1, (1, 2): 1, (2, 3): 1}, lone_pairs={1: 2, 2: 2})
>>> h2o2.oxidation_states()
{0: 1.0, 1: -1.0, 2: -1.0, 3: 1.0}
symbols: list#
total_formal_charge()[source]#

Sum of all atoms’ formal charges – must equal the molecule’s net charge.

A useful self-consistency check on a proposed Lewis structure: any valid structure’s formal charges must sum to the actual net charge of the species (0 for a neutral molecule).

Return type:

float

Returns:

float

class chemistrykit.structure.Molecule(symbols, coordinates, bonds=<factory>)[source]#

Bases: object

A minimal molecular-geometry container: element symbols, 3D coordinates, and a bond list.

No file-format parsing (PDB/XYZ/SMILES/…) is implemented or depended on anywhere in chemistrykit, per the package’s “no cheminformatics dependencies” scope decision – a Molecule is built directly from plain Python/numpy data, e.g. by chemistrykit.structure.systems.vsepr.build_vsepr_molecule() or by hand for a worked example.

Parameters:
  • symbols (list) – Element symbol of each atom, length n_atoms.

  • coordinates (ndarray) – Cartesian coordinates, in angstrom (Å) by convention throughout this subpackage (chosen because bond lengths and van der Waals radii are most naturally tabulated in Å; nothing here depends on that unit beyond consistent internal use).

  • bonds (list) – 0-indexed atom-index pairs describing the bond connectivity (an adjacency/bond list, not a bond-order-aware structure – see chemistrykit.structure.systems.lewis.LewisStructure for a representation that also carries bond orders and lone pairs).

Examples

A bent water molecule (experimental geometry: r(O-H) = 0.958 Å, angle(H-O-H) = 104.5 degrees):

>>> import numpy as np
>>> angle = np.radians(104.5)
>>> r = 0.958
>>> water = Molecule(
...     symbols=["O", "H", "H"],
...     coordinates=[
...         [0.0, 0.0, 0.0],
...         [r * np.sin(angle / 2), 0.0, r * np.cos(angle / 2)],
...         [-r * np.sin(angle / 2), 0.0, r * np.cos(angle / 2)],
...     ],
...     bonds=[(0, 1), (0, 2)],
... )
>>> round(water.bond_length(0, 1), 3)
0.958
>>> round(water.bond_angle(1, 0, 2), 1)
104.5
adjacency()[source]#

Build a neighbor list from bonds.

Return type:

dict

Returns:

dict of int to list of int – Maps each atom index with at least one bond to the sorted list of its bonded-neighbor indices.

bond_angle(i, j, k)[source]#

Angle \(\angle ijk\) at vertex atom j, in degrees.

\[\theta = \arccos\left(\frac{\vec{u}_{ji}\cdot\vec{u}_{jk}}{|\vec{u}_{ji}||\vec{u}_{jk}|}\right)\]
Parameters:
  • i (int) – 0-indexed atom indices; j is the vertex (central atom).

  • j (int) – 0-indexed atom indices; j is the vertex (central atom).

  • k (int) – 0-indexed atom indices; j is the vertex (central atom).

Return type:

float

Returns:

float – Angle, in degrees, in \([0, 180]\).

bond_length(i, j)[source]#

Distance between atoms i and j, in the same units as coordinates.

Parameters:
  • i (int) – 0-indexed atom indices (need not be an entry in bonds).

  • j (int) – 0-indexed atom indices (need not be an entry in bonds).

Return type:

float

Returns:

float

bonds: list#
center_of_mass()[source]#

Mass-weighted centroid of the atoms, using standard atomic weights.

Return type:

ndarray

Returns:

ndarray, shape (3,)

centered_coordinates(weighted=False)[source]#

Coordinates translated so the centroid (or center of mass) sits at the origin.

Parameters:

weighted (bool) – If True, center on center_of_mass(); otherwise on the unweighted centroid(). For a molecule with genuine point- group symmetry, every symmetry element passes through both (equivalent atoms share a mass, so the two centers coincide; see chemistrykit.structure.systems.point_group), so either choice is valid for symmetry-element detection.

Return type:

ndarray

Returns:

ndarray, shape (n_atoms, 3)

centroid()[source]#

Unweighted (geometric) centroid of the atoms.

Return type:

ndarray

Returns:

ndarray, shape (3,)

coordinates: ndarray#
is_linear(tol=0.0001)[source]#

Whether all atoms lie on a single straight line.

Parameters:

tol (float) – Tolerance on the normalized cross product used to detect collinearity.

Return type:

bool

Returns:

bool

Examples

>>> co2 = Molecule(symbols=["O", "C", "O"], coordinates=[[0, 0, -1.16], [0, 0, 0], [0, 0, 1.16]])
>>> co2.is_linear()
True
property n_atoms: int#

Number of atoms.

Type:

int

symbols: list#
class chemistrykit.structure.PointGroupCharacterTable(name, operations, irreps, characters)[source]#

Bases: object

A finite point group’s character table: irreducible representations, operation classes, and characters.

This is reference data – transcribed from the standard tables (Cotton, Chemical Applications of Group Theory, 3rd ed., Appendix A), in the same spirit as chemistrykit.periodic_table being plain data rather than a dependency – not something determine_point_group() derives computationally.

Parameters:
character(irrep, operation)[source]#

Look up a single character by irrep and operation-class label.

Parameters:
Return type:

float

Returns:

float

Examples

>>> table = get_character_table("C2v")
>>> table.character("A1", "E")
1.0
>>> table.character("B1", "C2")
-1.0
characters: list#

characters[i][j] is the character of irrep i under operation class j.

Type:

list of list of float

property class_sizes: list#

Number of operations in each class, read from the leading count of each label.

"8C3" is a class of 8 operations, "E" or "sigma_v(xz)" a class of 1. Raises ValueError for the infinite groups (C_inf_v, D_inf_h), whose classes are continuous.

Type:

list of int

irreps: list#

Irreducible-representation (Mulliken) labels, e.g. ["A1", "A2", "B1", "B2"].

Type:

list of str

name: str#

Point-group symbol (matches a determine_point_group() output where finite).

Type:

str

operations: list#

Symmetry-operation class labels, e.g. ["E", "C2", "sigma_v(xz)", "sigma_v'(yz)"].

Type:

list of str

property order: int#

The group order h (total number of symmetry operations).

Type:

int

reduce(reducible_characters)[source]#

Decompose a reducible representation into irreducible ones.

Uses the reduction formula (the “great orthogonality theorem” applied to characters; Cotton, Chemical Applications of Group Theory, 3rd ed., Ch. 4.3):

\[a_i = \frac{1}{h}\sum_{R} g_R\,\chi(R)\,\chi_i(R)\]

where the sum runs over operation classes of size \(g_R\) and \(h\) is the group order.

Parameters:

reducible_characters (sequence of float) – Characters of the reducible representation, one per operation class, in the order of operations.

Return type:

dict

Returns:

dict of str to int – Multiplicity of each irrep that occurs (irreps with zero multiplicity are left out).

Raises:

ValueError – If the input length does not match the number of classes, or the multiplicities are not integers (the characters do not form a representation of this group).

Examples

The five d orbitals in an octahedral field (Bethe, 1929) split into a doubly degenerate \(e_g\) and a triply degenerate \(t_{2g}\) set:

>>> oh = get_character_table("Oh")
>>> oh.reduce([5, -1, 1, -1, 1, 5, -1, -1, 1, 1])
{'Eg': 1, 'T2g': 1}
class chemistrykit.structure.PointGroupResult(group_name, is_linear, principal_axis_order, n_perpendicular_c2, has_sigma_h, has_sigma_v, has_inversion_center, n_c3_axes, n_mirror_planes, extra=<factory>)[source]#

Bases: object

The detected symmetry elements of a molecule, and the resulting point-group assignment.

Returned by determine_point_group(); every boolean/count field reflects an element that was actually found by direct geometric testing against the input coordinates (see the module docstring), not inferred from the final group_name.

Parameters:
  • group_name (str)

  • is_linear (bool)

  • principal_axis_order (int)

  • n_perpendicular_c2 (int)

  • has_sigma_h (bool)

  • has_sigma_v (bool)

  • has_inversion_center (bool)

  • n_c3_axes (int)

  • n_mirror_planes (int)

  • extra (dict)

extra: dict#

Free-form slot for additional diagnostics.

Type:

dict

group_name: str#

The assigned point-group symbol (e.g. "C2v", "Td", "D_inf_h" – the infinity symbol is spelled out ASCII-style since it is not a valid Python/LaTeX-free identifier component).

Type:

str

has_inversion_center: bool#

Whether an inversion center was found.

Type:

bool

has_sigma_h: bool#

Whether a mirror plane perpendicular to the principal axis was found.

Type:

bool

has_sigma_v: bool#

Whether a mirror plane containing the principal axis was found.

Type:

bool

is_linear: bool#

Whether all atoms are collinear.

Type:

bool

n_c3_axes: int#

Number of distinct C3 axes found (>= 4 triggers the cubic T/Td/Th/O/Oh branch of the decision tree).

Type:

int

n_mirror_planes: int#

Total number of distinct mirror planes found.

Type:

int

n_perpendicular_c2: int#

Number of C2 axes found perpendicular to the principal axis (the condition for the D point-group families).

Type:

int

principal_axis_order: int#

The highest-order proper rotation axis found, 1 if none (beyond the trivial identity) was detected. Not meaningful when is_linear is True (the true axis order is infinite).

Type:

int

class chemistrykit.structure.VSEPRGeometry(steric_number, lone_pairs)[source]#

Bases: object

The predicted shape for a given steric number and lone-pair count.

Parameters:
  • steric_number (int) – Number of electron domains (sigma bonds + lone pairs).

  • lone_pairs (int) – Number of lone pairs on the central atom, <= steric_number.

Examples

>>> VSEPRGeometry(steric_number=4, lone_pairs=2).shape_name
'bent'
>>> VSEPRGeometry(steric_number=4, lone_pairs=0).shape_name
'tetrahedral'
property bonding_positions: ndarray#

Unit-vector directions of just the bonded ligands.

Lone pairs preferentially occupy the least sterically crowded domain positions (see domain_positions()), so these are the remaining positions after that assignment – e.g. for steric_number=5, lone_pairs=1 (seesaw), one equatorial position is given up to the lone pair and the four bonding positions are the other two equatorial plus both axial sites.

Type:

ndarray, shape (n_bonding_domains, 3)

property electron_domain_positions: ndarray#

All domain unit-vector directions, bonds and lone pairs alike.

Type:

ndarray, shape (steric_number, 3)

lone_pairs: int#
property n_bonding_domains: int#

Number of sigma-bond (ligand) positions, steric_number - lone_pairs.

Type:

int

property shape_name: str#

The AXE-method molecular-shape name (Gillespie’s nomenclature).

Type:

str

steric_number: int#
chemistrykit.structure.angle_between(v1, v2)[source]#

Angle between two vectors, in degrees.

Parameters:
Return type:

float

Returns:

float

Examples

>>> round(angle_between([1, 0, 0], [0, 1, 0]), 6)
90.0
chemistrykit.structure.baeyer_angle_strain(ring_size)[source]#

Baeyer’s angle strain per bond for a planar ring, in degrees.

\[\delta(n) = \tfrac{1}{2}\left[109.47^\circ - \frac{180^\circ\,(n-2)}{n}\right]\]
Parameters:

ring_size (int) – Number of ring atoms, at least 3.

Return type:

float

Returns:

float – Positive when the ring angle is squeezed below tetrahedral, negative when a planar ring would force it wider.

Examples

Baeyer’s own values: 24 deg 44 min for cyclopropane, 9 deg 44 min for cyclobutane and only 0 deg 44 min for cyclopentane:

>>> [round(baeyer_angle_strain(n), 2) for n in (3, 4, 5, 6)]
[24.74, 9.74, 0.74, -5.26]
chemistrykit.structure.bond_dipole_sum(bond_vectors, bond_moments)[source]#

Molecular dipole as the vector sum of bond dipoles.

\[\boldsymbol\mu = \sum_b \mu_b\,\hat{\mathbf u}_b\]
Parameters:
  • bond_vectors (array-like, shape (n_bonds, 3)) – Direction of each bond dipole (any length; normalized here).

  • bond_moments (array-like, shape (n_bonds,)) – Magnitude of each bond dipole, in debye.

Return type:

ndarray

Returns:

ndarray, shape (3,) – Molecular dipole vector, in debye.

Examples

Two O-H bond dipoles of 1.51 D at water’s 104.5 degree angle give \(2\mu_{OH}\cos(\theta/2) \approx 1.85\) D, the measured value:

>>> half = np.radians(104.5 / 2)
>>> u = [[np.sin(half), 0, np.cos(half)], [-np.sin(half), 0, np.cos(half)]]
>>> round(float(np.linalg.norm(bond_dipole_sum(u, [1.51, 1.51]))), 2)
1.85
chemistrykit.structure.bond_length_from_order(single_bond_length, bond_order, c=0.71)[source]#

Invert bond_order_from_length(): predict a bond length from a bond order.

\[D(n) = D(1) - c\log_{10} n\]
Parameters:
  • single_bond_length (float) – Reference single-bond length \(D(1)\).

  • bond_order (float) – Bond order n (> 0).

  • c (float)

Return type:

float

Returns:

float – Predicted bond length, same units as single_bond_length.

Raises:

ValueError – If bond_order is not positive.

Examples

A round trip through bond_order_from_length() recovers the original length exactly (the two functions are exact inverses):

>>> D1, Dn = 1.54, 1.34
>>> n = bond_order_from_length(D1, Dn)
>>> round(bond_length_from_order(D1, n), 9) == round(Dn, 9)
True

Higher bond order predicts a shorter bond – triple shorter than double shorter than single:

>>> lengths = [bond_length_from_order(1.54, n) for n in (1, 2, 3)]
>>> bool(lengths[0] > lengths[1] > lengths[2])
True
chemistrykit.structure.bond_order_from_length(single_bond_length, observed_length, c=0.71)[source]#

Estimate a bond order from an observed bond length, via Pauling’s empirical correlation.

\[D(n) = D(1) - c\log_{10} n \quad\Longrightarrow\quad n = 10^{(D(1)-D(n))/c}\]

where \(D(1)\) is the reference single-bond length and c is an empirical, bond-type-specific constant (L. Pauling, J. Am. Chem. Soc. 69, 542 (1947)). Shorter bonds correspond to a higher (order > 1) bond order; this is an empirical correlation, not a first- principles bonding calculation – see the module docstring.

Parameters:
  • single_bond_length (float) – Reference single-bond (n=1) length \(D(1)\), in any consistent length unit (e.g. angstrom).

  • observed_length (float) – The bond length whose order is being estimated, \(D(n)\), same units as single_bond_length.

  • c (float) – The empirical correlation constant; must be re-fit for bond types other than carbon-carbon.

Return type:

float

Returns:

float – Estimated (generally non-integer) bond order.

Examples

Benzene’s C-C bond (1.397 Å) is intermediate between a single (1.54 Å) and double (1.34 Å) bond, giving an estimated order between 1 and 2 (the classic evidence for delocalized aromatic bonding):

>>> n = bond_order_from_length(single_bond_length=1.54, observed_length=1.397)
>>> bool(1.0 < n < 2.0)
True

A bond exactly at the reference single-bond length has order 1:

>>> round(bond_order_from_length(1.54, 1.54), 6)
1.0

With Pauling’s C-C constant (0.71 Å per decade of bond order), ethylene’s 1.33 Å double bond comes out very close to order 2:

>>> round(bond_order_from_length(1.54, 1.33), 2)
1.98
chemistrykit.structure.build_vsepr_molecule(steric_number, lone_pairs, bond_length=1.0, central_symbol='C', ligand_symbol='H')[source]#

Build a Molecule from a VSEPR prediction.

Parameters:
  • steric_number (int) – Number of electron domains (sigma bonds + lone pairs).

  • lone_pairs (int) – Number of lone pairs on the central atom.

  • bond_length (float) – Central-atom-to-ligand distance, in angstrom.

  • central_symbol (str) – Element symbol placed at the central atom.

  • ligand_symbol (str) – Element symbol placed at every ligand position.

Return type:

Molecule

Returns:

Molecule – Central atom at the origin plus n_bonding_domains identical ligand atoms; lone pairs are not represented as atoms (they carry no nuclear position), only via their effect on which directions the ligands occupy.

Examples

Methane (AX4E0): a perfect tetrahedron, every H-C-H angle exactly 109.47 degrees:

>>> methane = build_vsepr_molecule(steric_number=4, lone_pairs=0, bond_length=1.09, central_symbol="C", ligand_symbol="H")
>>> round(methane.bond_angle(1, 0, 2), 4)
109.4712

Water (AX2E2): the idealized VSEPR angle is the parent tetrahedron’s 109.47 degrees (see the module docstring for why this is an idealization – the real H-O-H angle is compressed to 104.5 degrees by extra lone-pair repulsion):

>>> water = build_vsepr_molecule(steric_number=4, lone_pairs=2, bond_length=0.96, central_symbol="O", ligand_symbol="H")
>>> round(water.bond_angle(1, 0, 2), 4)
109.4712
chemistrykit.structure.chair_cyclohexane_coordinates(bond_length=1.54)[source]#

Carbon coordinates of an ideal chair cyclohexane with exactly tetrahedral C-C-C angles.

The six carbons sit alternately at heights \(\pm z\) on a circle of radius \(\rho\), 60 degrees apart. Requiring bond length d and bond angle \(\theta_T\) fixes \(\rho^2 = 2d^2(1-\cos\theta_T)/3\) (the 1-3 distance is \(\rho\sqrt3\)) and then \(4z^2 = d^2 - \rho^2\).

Parameters:

bond_length (float) – C-C bond length, in angstrom.

Return type:

ndarray

Returns:

ndarray, shape (6, 3) – Ring atoms in order around the ring.

Examples

Every C-C-C angle of the puckered chair is tetrahedral, so Baeyer’s predicted strain for a six-membered ring vanishes once the ring is not forced to be planar:

>>> from chemistrykit.structure.core.base_system import angle_between
>>> x = chair_cyclohexane_coordinates()
>>> round(angle_between(x[0] - x[1], x[2] - x[1]), 4)
109.4712
chemistrykit.structure.coulson_pi_bond_order(coefficients, occupations, i, j)[source]#

The Coulson pi bond order between atoms i and j from Huckel molecular-orbital coefficients.

\[p_{ij} = \sum_k n_k c_{ik} c_{jk}\]

summed over molecular orbitals k with occupation number \(n_k\in\{0,1,2\}\) (C. A. Coulson, Proc. R. Soc. Lond. A 169, 413 (1939); Streitwieser, Molecular Orbital Theory for Organic Chemists, Ch. 2). This measures the pi-bonding contribution only (Huckel theory does not model the sigma framework at all, per chemistrykit.quantum.systems.huckel’s module docstring), so a formally “single” sigma-bonded pair with a fractional \(p_{ij}\) (e.g. benzene’s 2/3) has a total bond order of \(1+p_{ij}\).

Parameters:
  • coefficients (ndarray) – MO coefficient matrix, columns are individual MOs – e.g. chemistrykit.quantum.core.base_system.EigenstateResult.coefficients from a chemistrykit.quantum.systems.huckel.HuckelSystem.solve() call.

  • occupations (array-like of float, length n_mo) – Occupation number of each MO (0, 1, or 2 electrons), same column order as coefficients.

  • i (int) – 0-indexed atom (basis-function) indices.

  • j (int) – 0-indexed atom (basis-function) indices.

Return type:

float

Returns:

float

Examples

Ethene: the single pi bond is fully formed, \(p_{01}=1\):

>>> from chemistrykit.quantum.systems.huckel import HuckelSystem
>>> ethene = HuckelSystem(n_atoms=2, bonds=[(0, 1)])
>>> result = ethene.solve()
>>> order = np.argsort(result.energies)
>>> occ = np.zeros(2)
>>> occ[order[0]] = 2.0
>>> round(coulson_pi_bond_order(result.coefficients, occ, 0, 1), 6)
1.0

Benzene: the textbook Coulson bond order of 2/3 for every adjacent carbon pair (delocalization spreads the pi bonding evenly around the ring, weaker than a localized double bond):

>>> benzene = HuckelSystem.cyclic_polyene(n_atoms=6)
>>> result = benzene.solve()
>>> order = np.argsort(result.energies)
>>> occ = np.zeros(6)
>>> occ[order[:3]] = 2.0
>>> round(coulson_pi_bond_order(result.coefficients, occ, 0, 1), 4)
0.6667
chemistrykit.structure.count_kekule_structures(n_atoms, bonds)[source]#

Number of Kekulé structures \(K\) of a conjugated skeleton.

Parameters:
Return type:

int

Returns:

int

Examples

Naphthalene (two fused six-membered rings, 10 carbons) has three:

>>> naphthalene = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 9), (9, 0),
...                (4, 5), (5, 6), (6, 7), (7, 8), (8, 9)]
>>> count_kekule_structures(10, naphthalene)
3
chemistrykit.structure.determine_point_group(molecule, tol=0.001)[source]#

Determine a molecule’s point group from its 3D coordinates.

Implements the standard point-group decision tree (Cotton, Chemical Applications of Group Theory, 3rd ed., Ch. 4, flowchart Fig. 4.1) on top of symmetry elements found by direct geometric testing (module docstring):

  1. Linear molecules (all atoms collinear) are \(D_{\infty h}\) if they also have an inversion center, else \(C_{\infty v}\) (a finite rotation-order search cannot find the true \(C_\infty\) axis, so linearity is checked first and handled as a special case).

  2. Cubic groups: four or more distinct \(C_3\) axes (the hallmark of a tetrahedral/octahedral arrangement) route to \(O_h\) (if an inversion center is also present) or \(T_d\) otherwise.

  3. Otherwise, the principal axis is the highest-order proper rotation axis found (if none, the group is \(C_i\), \(C_s\), or \(C_1\) depending on whether an inversion center or a single mirror plane was found). Given a principal \(C_n\) axis: finding n (or more) \(C_2\) axes perpendicular to it selects the \(D_n\) family (further split into \(D_{nh}\)/\(D_{nd}\)/\(D_n\) by \(\sigma_h\)/\(\sigma_v\)); otherwise the \(C_n\) family (split into \(C_{nh}\)/\(C_{nv}\)/\(C_n\)).

Parameters:
  • molecule (Molecule)

  • tol (float) – Absolute distance tolerance (in the molecule’s coordinate units) for treating two atomic positions as coincident when testing a candidate symmetry operation.

Return type:

PointGroupResult

Returns:

PointGroupResult

Examples

Water is \(C_{2v}\), ammonia is \(C_{3v}\), and methane is \(T_d\) – the three textbook reference cases this classifier is validated against:

>>> import numpy as np
>>> angle = np.radians(104.5)
>>> r = 0.958
>>> water = Molecule(
...     symbols=["O", "H", "H"],
...     coordinates=[[0, 0, 0], [r * np.sin(angle / 2), 0, r * np.cos(angle / 2)], [-r * np.sin(angle / 2), 0, r * np.cos(angle / 2)]],
... )
>>> determine_point_group(water).group_name
'C2v'
>>> verts = np.array([[1, 1, 1], [1, -1, -1], [-1, 1, -1], [-1, -1, 1]], dtype=float)
>>> verts = verts / np.linalg.norm(verts[0]) * 1.09
>>> methane = Molecule(symbols=["C", "H", "H", "H", "H"], coordinates=np.vstack([[0, 0, 0], verts]))
>>> determine_point_group(methane).group_name
'Td'

Carbon dioxide is linear and centrosymmetric, \(D_{\infty h}\):

>>> co2 = Molecule(symbols=["O", "C", "O"], coordinates=[[0, 0, -1.16], [0, 0, 0], [0, 0, 1.16]])
>>> determine_point_group(co2).group_name
'D_inf_h'
chemistrykit.structure.dipole_moment(coordinates, charges)[source]#

Dipole-moment vector of a set of point charges, in debye.

\[\boldsymbol\mu = \sum_i q_i\,\mathbf r_i\]
Parameters:
  • coordinates (array-like, shape (n, 3)) – Positions in angstrom (e.g. Molecule.coordinates).

  • charges (array-like, shape (n,)) – Partial charges in units of the elementary charge. If they do not sum to zero the result depends on the origin.

Return type:

ndarray

Returns:

ndarray, shape (3,) – Dipole vector in debye (pointing from negative to positive charge, the physics convention).

Examples

Charges of +e and -e one angstrom apart make a dipole of 4.803 D:

>>> mu = dipole_moment([[0, 0, 0], [0, 0, 1.0]], [-1.0, 1.0])
>>> round(float(mu[2]), 3)
4.803

Linear CO2 with partial charges has no net dipole:

>>> mu = dipole_moment([[0, 0, -1.16], [0, 0, 0], [0, 0, 1.16]], [-0.4, 0.8, -0.4])
>>> float(np.linalg.norm(mu))
0.0
chemistrykit.structure.domain_positions(steric_number)[source]#

Return unit-vector electron-domain positions for the ideal N-domain polyhedron.

Real 3D coordinate generation (not a lookup table): each polyhedron’s vertices are placed from their defining symmetry, then normalized to the unit sphere.

Parameters:

steric_number (int) – Number of electron domains (sigma bonds + lone pairs), 2 through 6.

Return type:

ndarray

Returns:

ndarray, shape (steric_number, 3) – Unit vectors from the central atom, ordered so that – for a steric number with structurally inequivalent sites – the sites least favorable to a lone pair (see build_vsepr_molecule()) come last: equatorial-before-axial for steric_number=5 (equatorial positions have only two 90-degree neighbors vs. an axial position’s three, so a lone pair placed equatorially incurs less strong lone-pair/bonding-pair repulsion – Gillespie & Nyholm, Q. Rev. Chem. Soc. 11, 339 (1957)).

Examples

Every pair of tetrahedral vertices subtends exactly the tetrahedral angle, 109.47 degrees:

>>> import numpy as np
>>> from chemistrykit.structure.core.base_system import angle_between
>>> verts = domain_positions(4)
>>> angles = [angle_between(verts[i], verts[j]) for i in range(4) for j in range(i + 1, 4)]
>>> bool(np.allclose(angles, 109.4712206))
True
chemistrykit.structure.get_character_table(name)[source]#

Look up a built-in character table by point-group symbol.

Parameters:

name (str) – E.g. "C2v", "Td", "D_inf_h".

Return type:

PointGroupCharacterTable

Returns:

PointGroupCharacterTable

Raises:

KeyError – If name is not one of the point groups tabulated in CHARACTER_TABLES.

Examples

Every irrep’s character under the identity operation equals its dimension, and (for a real, unitary representation) the sum over all operations, weighted by class size, of the identity character times itself equals the group order – here just checking the totally symmetric irrep is trivially present:

>>> table = get_character_table("Td")
>>> table.character("A1", "E")
1.0
>>> table.character("T2", "E")
3.0
chemistrykit.structure.kekule_structures(n_atoms, bonds)[source]#

Enumerate every Kekulé structure (perfect matching) of a conjugated skeleton.

Parameters:
  • n_atoms (int) – Number of conjugated atoms (e.g. the carbons of a benzenoid hydrocarbon), labelled 0 .. n_atoms - 1.

  • bonds (sequence of tuple(int, int)) – The sigma-bond connectivity of the skeleton.

Return type:

list

Returns:

list of list of tuple(int, int) – Each entry is one Kekulé structure: the sorted list of bonds that are double bonds in it. Every atom appears in exactly one of them. An odd-membered or otherwise unmatchable skeleton returns [].

Examples

Benzene has exactly Kekulé’s two structures:

>>> ring = [(i, (i + 1) % 6) for i in range(6)]
>>> for s in kekule_structures(6, ring):
...     print(s)
[(0, 1), (2, 3), (4, 5)]
[(0, 5), (1, 2), (3, 4)]
chemistrykit.structure.pauling_electronegativity_difference(d_ab, d_aa, d_bb)[source]#

Pauling’s electronegativity difference from bond dissociation energies.

Pauling (1932) noticed that a heteronuclear bond A-B is stronger than the arithmetic mean of the homonuclear bonds A-A and B-B. He assigned this “extra ionic energy”

\[\Delta = D(\mathrm{A{-}B}) - \tfrac{1}{2}\left[D(\mathrm{A{-}A}) + D(\mathrm{B{-}B})\right]\]

to the bond’s partial ionic character and defined the electronegativity difference through \(|\chi_A - \chi_B| = \sqrt{\Delta/\mathrm{eV}}\) (L. Pauling, J. Am. Chem. Soc. 54, 3570 (1932)). With energies in kJ/mol this becomes \(|\chi_A-\chi_B| = 0.102\sqrt{\Delta}\), since 1 eV is 96.485 kJ/mol.

Parameters:
  • d_ab (float) – Bond dissociation energies of A-B, A-A and B-B, in kJ/mol.

  • d_aa (float) – Bond dissociation energies of A-B, A-A and B-B, in kJ/mol.

  • d_bb (float) – Bond dissociation energies of A-B, A-A and B-B, in kJ/mol.

Return type:

float

Returns:

float – \(|\chi_A - \chi_B|\) on the Pauling scale (0 if the extra ionic energy is negative).

Examples

H-Cl (432 kJ/mol) against H-H (436) and Cl-Cl (242): the extra ionic energy of 93 kJ/mol gives a difference of about 0.98, close to the tabulated Pauling values 3.16 - 2.20 = 0.96:

>>> round(pauling_electronegativity_difference(432.0, 436.0, 242.0), 2)
0.98

Equal bond energies mean no ionic contribution:

>>> pauling_electronegativity_difference(200.0, 200.0, 200.0)
0.0
chemistrykit.structure.planar_ring_angle(ring_size)[source]#

Interior angle of a regular planar polygon with ring_size vertices, in degrees.

Parameters:

ring_size (int) – Number of ring atoms, at least 3.

Return type:

float

Returns:

float

Examples

>>> planar_ring_angle(3), planar_ring_angle(6)
(60.0, 120.0)
chemistrykit.structure.unit_vector(v)[source]#

Normalize a vector to unit length.

Parameters:

v (ndarray)

Return type:

ndarray

Returns:

ndarray