Geometry optimisation

Relax a molecular structure to its energy minimum. vibe-qc uses ASE’s BFGS-line-search optimiser under the hood, driven by analytic HF/DFT forces from the C++ core. run_job(..., optimize=True) is the one-liner; the standard .out / .molden / .traj files all get written as usual, plus the trajectory animation.

Water at HF/6-31G*

Starting from an approximate geometry:

from pathlib import Path
from vibeqc import Atom, Molecule, run_job

mol = Molecule([
    Atom(8, [0.0,   0.0,   0.0  ]),
    Atom(1, [0.0,   1.60, -1.10]),   # slightly off-equilibrium
    Atom(1, [0.0,  -1.60, -1.10]),
])

run_job(
    mol,
    basis="6-31g*",
    method="rhf",
    output="water_opt",
    optimize=True,
    fmax=0.05,            # eV/Å — converged when |F|_max drops below this
    max_opt_steps=50,
)

Files produced:

  • water_opt.out — initial geometry, optimiser parameters, final optimised geometry, then the full SCF trace at the minimum

  • water_opt.molden — MOs at the minimum, for Avogadro / Jmol

  • water_opt.traj — one frame per BFGS step (ASE binary); view with ase gui water_opt.traj

Reading the output

After the banner and atom table, water_opt.out contains:

  Geometry optimisation (ASE/BFGS) -> water_opt.traj
  fmax = 0.05 eV/A, max_steps = 50

  Optimised geometry
  Atoms (bohr)
  ------------------------------------------------------
     1  Z=  8       0.00000000      0.00000000      0.07390216
     2  Z=  1       0.00000000      1.44014211     -1.00364326
     3  Z=  1       0.00000000     -1.44014211     -1.00364326
  charge=0  multiplicity=1  n_electrons=10

  iter     energy (Ha)  ...
  (SCF trace at the optimised geometry)
  ...
  converged in 7 iterations; E = -76.0108432849 Ha

The OH bond length and HOH angle at this HF/6-31G* minimum: r(OH) = 0.947 Å, HOH = 105.5° (experimental: 0.957 Å, 104.5°). HF slightly under-shoots the bond and opens the angle.

Parameters

Argument

Default

Meaning

fmax

0.05

Force-convergence tolerance (eV/Å)

max_opt_steps

200

Hard iteration limit for the optimiser

rhf_options

RHFOptions()

Applied to every intermediate SCF; bump max_iter here if BFGS visits hard-to-converge geometries

rhf_options / rks_options forward through to every frame of the optimisation, not just the final single point. This matters for H-bonded systems and floppy molecules where one BFGS step sometimes lands on a particularly stubborn geometry — bump max_iter = 200 or add damping = 0.3 to the options struct and the whole trajectory stays stable.

A harder case: the water dimer

H-bonded clusters sit on a very flat potential-energy surface — plain BFGS Hessian extrapolation can push the waters apart before the quadratic model is trustworthy. run_job uses BFGSLineSearch which rejects uphill steps by construction, and pair this with a realistic fmax:

from vibeqc import Atom, InitialGuess, Molecule, RHFOptions, run_job

mol = Molecule.from_xyz("h2o_dimer.xyz")  # see examples/h2o_dimer.xyz

rhf_opts = RHFOptions()
rhf_opts.max_iter = 200                   # H-bonded intermediates need it
rhf_opts.initial_guess = InitialGuess.SAD

run_job(
    mol, basis="6-31g*", method="rhf",
    output="water_dimer_opt",
    optimize=True,
    fmax=0.2,                             # force precision floor is ~0.1 eV/Å
    max_opt_steps=25,
    rhf_options=rhf_opts,
)

You’ll get R(OO) ≈ 2.98 Å (experimental: 2.976 Å) — the well-known HF/6-31G* result. For a real binding-energy number add a correlation method and dispersion; for a chemistry demonstration this is enough.

Visualising the trajectory

The .traj file holds positions, energies, and forces for every frame. ASE’s built-in GUI animates it:

ase gui water_opt.traj

Or inspect programmatically:

import numpy as np
from ase.io import read

frames = read("water_opt.traj", index=":")
energies = np.array([f.get_potential_energy() for f in frames])
print(f"Energy dropped {energies[0] - energies[-1]:.4f} eV over "
      f"{len(frames) - 1} steps")

DFT geometries

Same call, different method:

run_job(
    mol, basis="def2-tzvp", method="rks", functional="PBE",
    output="water_pbe_opt",
    optimize=True,
)

PBE lengthens O-H bonds by ~0.02 Å relative to HF and pulls the HOH angle back toward 104°, generally matching experiment more closely (with the well-known caveat that GGA over-binds H-bonds somewhat).

Next