# Molecular input

Welcome to **Part III — the User guide**. Parts I and II got you installed and gave you the theory; from
here on we work feature by feature, the way you will in practice. Every calculation begins in the same
place: you **describe a molecule**. This chapter is about that first step — how you tell qc-rs *what* you
are computing, and how to check that it understood you.

## Theory: what defines a molecule for a quantum-chemistry calculation

In [Part II](../10-foundations/many-electron-and-bo.md) we met the **Born–Oppenheimer approximation**: the
nuclei are so much heavier than the electrons that we **clamp them in place** and solve the electronic
problem in the field they create. That single idea dictates *exactly* what a calculation needs as input.
Three things, and only three:

1. **The nuclear geometry** — which atoms, and where. Because the nuclei are clamped, their positions
   $\{\mathbf R_A\}$ are fixed numbers you supply. They set the external potential
   $v_{\text{ext}}(\mathbf r) = -\sum_A Z_A / |\mathbf r - \mathbf R_A|$ that the electrons move in, and the
   constant **nuclear repulsion** $E_{\text{nuc}} = \sum_{A<B} Z_A Z_B / R_{AB}$ that is added to the
   electronic energy.
2. **The total charge** — how many electrons there are. A neutral molecule has $N_e = \sum_A Z_A$ electrons;
   a charge $q$ removes ($q>0$) or adds ($q<0$) that many. This fixes $N_e$, the number the SCF fills.
3. **The spin state** — how those electrons' spins are arranged, given as the **multiplicity** $2S+1$.

Everything else — the basis set, the method, convergence settings — is a *choice about how to solve* the
problem. Geometry, charge, and spin are the problem itself. In qc-rs you provide all three to
`qc.chk.new(...)`, which builds the **checkpoint** from [Core concepts](../00-intro/concepts.md).

:::{prf:definition} Spin multiplicity
:label: def-multiplicity

For a state with total spin quantum number $S$, the **multiplicity** is $2S+1$ — the number of $M_S$
projections ($-S, \dots, +S$). Each **unpaired electron** contributes spin $\tfrac12$, so with $n$ unpaired
electrons $S = n/2$ and the multiplicity is $n+1$:

| unpaired electrons $n$ | $S$ | multiplicity $2S+1$ | name |
|---|---|---|---|
| 0 | 0 | **1** | singlet |
| 1 | ½ | **2** | doublet |
| 2 | 1 | **3** | triplet |

This is the **Gaussian convention**, and it is what qc-rs's `spin=` expects — pass `spin=1` for a
closed-shell singlet, `spin=2` for a radical, `spin=3` for a triplet.
:::

## Usage: the `atom=` string

You describe the geometry with the **`atom=`** argument — a string with **one atom per line**, each line an
element label followed by three Cartesian coordinates:

```text
Element   x   y   z
```

Atoms are separated by a **newline** *or* a **semicolon** `;`, so both of these are the same water molecule:

```python
import qc

# semicolon-separated (compact, good for one-liners)
w1 = qc.chk.new(atom="O 0 0 0; H 0 0.757 0.587; H 0 -0.757 0.587", ao="sto-3g")

# newline-separated (readable, good for real molecules) — a triple-quoted string
w2 = qc.chk.new(atom="""
    O   0.000   0.000   0.117
    H   0.000   0.757  -0.469
    H   0.000  -0.757  -0.469
""", ao="sto-3g")
```

The **element label** can be a chemical symbol (`H`, `C`, `Fe` — case-insensitive), an **atomic number**
(`6` means carbon), or a symbol with a **label suffix** to distinguish otherwise-identical atoms
(`C1`, `C2`, `Fe3`). The label matters when you later want to give *one specific atom* its own basis set —
see [Basis sets & AO representation](basis-and-ao.md).

:::{tip} A Python note: triple-quoted strings
The `"""… """` form is a **multi-line string**. It is the natural way to write a real geometry — one atom
per line, aligned columns, just as you would in an input file. Leading/trailing blank lines and indentation
are fine; qc-rs ignores them.
:::

### Units: ångström or bohr

Coordinates are in **ångström by default** (`unit="angstrom"`). Set `unit="bohr"` to give them in atomic
units instead:

```python
# these are the *same* H2 molecule (bond length ≈ 0.741 Å ≈ 1.400 bohr)
h2_ang  = qc.chk.new(atom="H 0 0 0; H 0 0 0.741", ao="sto-3g", unit="angstrom")
h2_bohr = qc.chk.new(atom="H 0 0 0; H 0 0 1.400", ao="sto-3g", unit="bohr")
```

Internally qc-rs **always stores coordinates in bohr** (atomic units), converting on input with the single
constant $\text{1 Å} = 1/0.52917721092 \approx 1.8897$ bohr. That is why the `coordinates()` accessor below
always reports bohr, whatever unit you typed.

:::{important} ångström is the default — do not assume bohr
A geometry copied from a PDB file, a paper, or most GUIs is in **ångström**, so the default is usually what
you want. But if you paste atomic-unit coordinates and forget `unit="bohr"`, every bond will be ~1.9× too
long and the calculation will be nonsense (or fail to converge). When in doubt, check with `coordinates()`.
:::

### Charge and spin

`charge=` and `spin=` default to a **neutral singlet** (`charge=0, spin=1`). Change them for ions and
open-shell species. For example, the neutral **hydroxyl radical** OH has 9 electrons — an odd number, so it
*cannot* be a singlet; it is a doublet (`spin=2`):

```python
oh = qc.chk.new(atom="O 0 0 0; H 0 0 0.97", ao="sto-3g", unit="angstrom",
                charge=0, spin=2)          # neutral doublet radical
print(oh.nelectron(), oh.spin)             # 9 2
```

:::{prf:example} Choosing charge and spin
:label: ex-charge-spin

- **Water**, H₂O — neutral, all electrons paired → `charge=0, spin=1`.
- **Hydroxide**, OH⁻ — one extra electron → `charge=-1, spin=1` (10 electrons, closed shell).
- **Methyl radical**, CH₃ — neutral with one unpaired electron → `charge=0, spin=2`.
- **Dioxygen**, O₂ — its ground state is a **triplet** (two unpaired electrons) → `charge=0, spin=3`.

A quick sanity check: the electron count and the multiplicity must have consistent parity. An even electron
count cannot be a doublet; an odd count cannot be a singlet. If they disagree, `run()` will reject the state.
:::

## Inspecting what qc-rs parsed

Before running anything heavy, **look at what the checkpoint understood**. This catches unit mistakes,
miscounted electrons, and typos immediately. The relevant accessors (some are properties, some are methods
— note the `()`):

```python
w = qc.chk.new(atom="O 0 0 0; H 0 0.757 0.587; H 0 -0.757 0.587", ao="sto-3g", unit="angstrom")

w.natom               # 3          — number of (real + ghost) atoms
w.symbols             # ['O','H','H']
w.charge, w.spin      # (0, 1)
w.nelectron()         # 10         — electrons the SCF will fill
w.coordinates()       # (3, 3) array, always in BOHR
w.nuclear_energy()    # 9.188258   — E_nuc (hartree), the constant added to the electronic energy
```

That is the whole molecular specification, echoed back to you. If `nelectron()` or `nuclear_energy()` looks
wrong, fix the input *now* — every downstream number depends on it.

## Special atoms

Beyond ordinary atoms, the `atom=` grammar supports three special kinds. You will not need them for routine
work, but they are essential for a few important techniques.

### Ghost atoms (`Element-Bq`) — for counterpoise / BSSE

A **ghost atom**, written `Element-Bq` (e.g. `O-Bq`, `H-Bq`), sits at a point carrying that element's
**basis functions but no nucleus and no electrons** — nuclear charge 0, zero electrons, zero nuclear
repulsion. Its purpose is the **counterpoise correction** for **basis-set superposition error (BSSE)**: when
two molecules approach, each borrows the *other's* basis functions and its energy drops artificially. To
measure that artefact, you recompute one monomer with the *partner's* basis present as ghosts:

```python
# a monomer computed in the FULL dimer basis: the partner's atoms are ghosts
mono = qc.chk.new(atom="O 0 0 0; H 0 0 1.8; H-Bq 0 0 3.5", ao="sto-3g", unit="bohr")
print(mono.natom, mono.nelectron())     # 3 9   — ghost counts as an atom, adds NO electrons
```

Ghosts **are** counted by `natom` (they carry AO functions) but contribute **nothing** to `nelectron()` or
`nuclear_energy()`. The basis for a ghost is looked up under its label first, then the base element, so
`H-Bq` automatically inherits the `H` basis (details in [Basis sets & AO representation](basis-and-ao.md)).

### Dummy atoms (`X`, `Xx`) — geometric reference points

A **dummy atom** (`X`, `Xx`, or a labelled `X1`, `X-ref`) is a pure **geometric marker**: no nuclear charge,
no electrons, **and no basis functions**. It is a coordinate you can reference (e.g. to define a symmetry
axis or a measurement point) without adding anything physical:

```python
d = qc.chk.new(atom="C 0 0 0; H 0 0 1.089; X 0 0 5.0", ao="sto-3g", unit="angstrom")
print(d.natom, d.nelectron())     # 2 7   — the dummy is NOT counted, adds nothing
```

Unlike ghosts, dummies are **not** counted by `natom`; they live in `dummy_atoms()`.

### Translation vectors (`TV`) — periodic systems

A line beginning `TV` supplies a **lattice translation vector** for a periodic system (one `TV` for a
polymer, two for a sheet, three for a crystal). They are stored in `translation_vectors()` and, like
dummies, are not counted as atoms:

```python
poly = qc.chk.new(atom="""
    C  -0.574 -0.143  0.376
    C   0.579  0.022 -0.301
    TV  4.848  0.171  0.511
""", ao="sto-3g", unit="angstrom")
```

## Advanced: per-atom nuclear parameters

You can attach **nuclear parameters** to an atom in parentheses after its label — most usefully the
**isotope** (`Iso=`) and a **fragment index** (`Fragment=`, used to group atoms for counterpoise jobs):

```python
# a ¹³C isotopologue, atoms tagged into two fragments
m = qc.chk.new(atom="C(Iso=13,Fragment=1) 0 0 0; H(Fragment=1) 0 0 1.089",
               ao="sto-3g", unit="angstrom")
```

The full set (`Iso`, `Spin`, `ZNuc`, `ZEff`, `QMom`, `NMagM`, `RadNuclear`, `Fragment`, …), the MM
force-field fields, and the optimizer **freeze code** are documented in the reference chapter,
[Molecule specification](../40-reference/mol-spec.md). Most of this metadata is parsed and stored but does
not yet affect the Hamiltonian or the optimizer — the reference notes which fields are currently active.

:::{note} What is *not* accepted
qc-rs takes the molecular geometry as a **string** in the grammar above. A Python **list of tuples** (as
some other packages accept) is **not** a valid `atom=` value — use the string form. And there is no
automatic geometry lookup by name: you always supply coordinates (or optimize them later with
[`.opt()`](gradients-geomopt.md)).
:::

## Worked example: build, inspect, and hand off to the SCF

Putting it together — describe water, verify it, and it is ready for a calculation:

```python
import qc

water = qc.chk.new(
    atom="""
        O   0.000000   0.000000   0.117300
        H   0.000000   0.757200  -0.469200
        H   0.000000  -0.757200  -0.469200
    """,
    ao="cc-pvdz",
    unit="angstrom",
    charge=0,
    spin=1,
)

# --- sanity check before doing any work ---
print("atoms     :", water.natom, water.symbols)   # 3 ['O', 'H', 'H']
print("electrons :", water.nelectron())            # 10
print("E_nuc     :", round(water.nuclear_energy(), 6))

# ready to compute — the checkpoint now flows into the rest of the workflow
water = water.scf(ref="r").run()
print("E(RHF)    :", round(water.scf.energy, 6))   # -76.026772
```

The molecule is now fully specified and the checkpoint carries it into every later step — the basis set
(next chapter), the SCF, properties, and gradients. Get the input right and verified here, and everything
downstream rests on solid ground.

:::{exercise}
:label: ex-mol-input

Write the `qc.chk.new(...)` call for each species, then check `nelectron()` and `spin`:

1. The **ammonium cation** NH₄⁺ (neutral NH₄ would have 11 electrons).
2. The **oxygen molecule** O₂ in its ground state (bond length 1.208 Å), *including the correct spin*.
3. A helium atom computed in the basis of a **helium dimer** (the second He is a ghost 3 Å away).
:::

:::{solution} ex-mol-input
:class: dropdown

```python
# 1. NH4+ : remove one electron from neutral NH4 (11 e-) -> 10 e-, closed-shell singlet
nh4 = qc.chk.new(atom="N 0 0 0; H 0.63 0.63 0.63; H -0.63 -0.63 0.63; "
                      "H -0.63 0.63 -0.63; H 0.63 -0.63 -0.63",
                 ao="sto-3g", unit="angstrom", charge=1, spin=1)
# nh4.nelectron() -> 10 ,  nh4.spin -> 1

# 2. O2 ground state is a TRIPLET (two unpaired electrons) -> spin=3
o2 = qc.chk.new(atom="O 0 0 0; O 0 0 1.208", ao="sto-3g", unit="angstrom",
                charge=0, spin=3)
# o2.nelectron() -> 16 ,  o2.spin -> 3

# 3. He in the He-dimer basis: the partner is a ghost (Bq)
he = qc.chk.new(atom="He 0 0 0; He-Bq 0 0 3.0", ao="cc-pvdz", unit="angstrom")
# he.nelectron() -> 2   (the ghost adds basis functions but no electrons)
```
:::

With the molecule defined, the next choice is **how to expand its orbitals** — the basis set. That is the
[next chapter](basis-and-ao.md).