Molecular input

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 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.

Definition 3 (Spin 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:

Element   x   y   z

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

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.

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:

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

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

Example 1 (Choosing charge and 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 ()):

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:

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

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:

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:

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):

# 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. 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()).

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

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

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 1

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).

With the molecule defined, the next choice is how to expand its orbitals — the basis set. That is the next chapter.