OpenMECP

OpenMECP is a high-performance Rust implementation of the MECP (Minimum Energy Crossing Point) optimizer for locating crossing points between two potential energy surfaces (PES) in quantum chemistry calculations. OpenMECP is developed and maintained by Dr. Le Nhan Pham.

Status: Beta testing phase — OpenMECP is highly robust and ready for production MECP calculations. All supported features are actively tested and improved.

What is an MECP?

A Minimum Energy Crossing Point is the lowest-energy geometry at which two electronic states become degenerate. MECPs are essential for understanding:

• Photochemical reactions — light-driven processes crossing between electronic states
• Spin-forbidden processes — reactions changing spin multiplicity
• Intersystem crossing — transitions between singlet and triplet manifolds
• Conical intersections — degeneracy points governing ultrafast photochemistry
• Non-adiabatic dynamics — processes beyond the Born–Oppenheimer approximation

Key Algorithm

OpenMECP implements the algorithm reported by Harvey et al. in Theor. Chem. Acc. 99, 95–99 (1998).

The MECP effective gradient combines two orthogonal components:

f-vector — drives the energy difference to zero:

$$\mathbf{f}=(E_1-E_2)\hat{\mathbf{x}}$$

g-vector — minimizes the average energy perpendicular to the gradient difference:

$$\mathbf{g}={\mathbf{f}}_1-(\hat{\mathbf{x}}\cdot {\mathbf{f}}_1)\hat{\mathbf{x}}$$

where the normalized gradient difference is:

$$\hat{\mathbf{x}}=\frac{{\mathbf{f}}_1-{\mathbf{f}}_2}{|{\mathbf{f}}_1-{\mathbf{f}}_2|}$$

Feature Highlights

Quick Navigation

• New to OpenMECP? Start with Installation then Quick Start.
• Ready to run? See the Input File Format and Keywords Reference.
• Want to understand the algorithms? See MECP Algorithm.
• Looking for the full API? See API Reference.

Citation

If you use OpenMECP in your research, please cite this preprint:

Pham, Le Nhan. OpenMECP: A High-Performance Rust Implementation for the Rigorous Location of Minimum Energy Crossing Points in Chemical Dynamics. 2026, https://doi.org/10.13140/RG.2.2.21309.73443.

Installation

Prerequisites

• Rust 1.70 or later — install from rustup.rs
• One or more QM programs: Gaussian 16 (recommended) or ORCA 5+

Building from Source

# Clone the repository
git clone https://github.com/lenhanpham/OpenMECP.git
cd OpenMECP

# Build optimized release binary
cargo build --release

# The binary is placed at:
#   target/release/omecp

Installing the Binary

Copy the binary to a directory on your PATH:

# Linux / macOS
cp target/release/omecp ~/.local/bin/
chmod 700 ~/.local/bin/omecp

# Or to a system-wide location (requires sudo)
sudo cp target/release/omecp /usr/local/bin/

On HPC systems, copy the binary to a personal bin/ directory and export it in your shell profile or job script:

export PATH="$HOME/bin:$PATH"

Verifying the Installation

omecp --help

You should see the OpenMECP help output listing all available options and topics.

HPC Setup

On HPC clusters, add the following to your job submission script:

# Load your QM program module
module load gaussian/g16   # or: module load orca/5.0

# Export the directory containing omecp
export PATH="/path/to/your/bin:$PATH"

# Run the calculation
omecp input.inp > output.log

Building Documentation Locally

# Install mdBook
cargo install mdbook

# Build and serve the documentation
mdbook serve docs

Open http://localhost:3000 in your browser.

Building API Reference Locally

cargo doc --no-deps --open

Quick Start

This page walks through the two ways to set up a calculation and how to run it.

Option 1: Generate a Template (Recommended)

OpenMECP can generate a ready-to-edit input file from any existing geometry file.

# From an XYZ file
omecp ci molecule.xyz
# Creates: molecule.inp

# From a Gaussian output file
# (adds _omecp suffix to avoid overwriting your output)
omecp ci calculation.log
# Creates: calculation_omecp.inp

# From a Gaussian input file
omecp ci input.gjf
# Creates: input.inp

# Specify a custom output name
omecp ci molecule.xyz my_run.inp

Then open the generated .inp file and fill in the required keywords (program, method, nprocs, mem, charge, mult_a, mult_b).

Option 2: Write an Input File Manually

Create a file named input.inp with the following structure:

*geom
C     -4.108383   -0.269570   -0.105531
C     -3.086048   -0.004314   -0.907012
H     -4.464908    0.456695    0.621271
H     -4.622302   -1.226624   -0.151874
C     -2.326155    1.302381   -0.898892
H     -2.751025   -0.749763   -1.624593
H     -2.733108    1.996398   -0.159695
H     -1.274400    1.129883   -0.660641
H     -2.379567    1.780215   -1.879473
*

*tail_a
# Extra Gaussian route keywords for state A
*

*tail_b
# Extra Gaussian route keywords for state B
*

program  = gaussian
method   = B3LYP/6-31G*
nprocs   = 4
mem      = 4GB
charge   = 0
mult_a   = 1
mult_b   = 3

See the Input File Format chapter for the full section reference and the Keywords Reference for all available parameters.

Running OpenMECP

omecp input.inp > output.log

On an HPC cluster, put the above in a job submission script and load your QM program module before calling omecp.

Checking the Results

Built-in Help

OpenMECP includes a comprehensive built-in help system:

omecp --help               # General usage
omecp --help keywords      # All input keywords with descriptions
omecp --help methods       # Supported QM methods
omecp --help features      # Feature overview
omecp --help examples      # Worked examples
omecp ci --help            # Template generation help

Minimal Gaussian Example

*geom
Fe    0.000000    0.000000    0.000000
N     0.000000    0.000000    2.100000
N     0.000000    2.100000    0.000000
*

*tail_a
*

*tail_b
*

program  = gaussian
method   = B3LYP/def2-SVP
nprocs   = 8
mem      = 16GB
charge   = 0
mult_a   = 1
mult_b   = 3

Minimal ORCA Example

*geom
Fe    0.000000    0.000000    0.000000
N     0.000000    0.000000    2.100000
*

*tail_a
%scf
  MaxIter 500
end
*

*tail_b
*

program  = orca
method   = B3LYP def2-SVP
nprocs   = 8
mem      = 8000
charge   = 0
mult_a   = 1
mult_b   = 3

Input File Format

OpenMECP input files use a custom section-based format. Sections are delimited by a *SECTION_NAME opening marker and a bare * closing marker. Keywords are specified as key = value pairs outside sections (or they can appear after the geometry section).

File Layout

*GEOM
  <Cartesian coordinates>
*

*TAIL_A
  <extra QM keywords for state A>
*

*TAIL_B
  <extra QM keywords for state B>
*

keyword1 = value1
keyword2 = value2
...

Section names are case-insensitive (*GEOM, *geom, and *Geom are all equivalent).

Section Reference

Keywords

All parameters controlling the run are specified as key = value pairs after the geometry section. See the Keywords Reference for the full list.

Minimum required keywords:

program  = gaussian    # or orca
method   = B3LYP/6-31G*
nprocs   = 4
mem      = 4GB
charge   = 0
mult_a   = 1
mult_b   = 3

Input Sections

*GEOM Section

Provides the initial molecular geometry in Cartesian coordinates (Ångström).

Format:

*GEOM
Element  X  Y  Z
...
*

Example:

*GEOM
C     -4.108383   -0.269570   -0.105531
C     -3.086048   -0.004314   -0.907012
H     -4.464908    0.456695    0.621271
H     -4.622302   -1.226624   -0.151874
C     -2.326155    1.302381   -0.898892
H     -2.751025   -0.749763   -1.624593
*

You can also reference an external XYZ file:

*GEOM
@geometry.xyz
*

The XYZ file must use standard format (Ångström units, with or without the atom-count / comment header lines).

*TAIL_A and *TAIL_B Sections

These sections inject extra text directly into the QM program's input for state A and state B respectively. Their content is program-specific.

Gaussian — append additional route section keywords:

*TAIL_A
scf(maxcycle=500,xqc)
*

*TAIL_B
scf(maxcycle=500,xqc)
*

ORCA — append additional input blocks (e.g., %tddft, %scf):

*TAIL_A
%tddft
  nroots 5
  iroot  1
end
*

*TAIL_B
%tddft
  nroots 5
  iroot  2
end
*

If no extra keywords are needed, include empty sections:

*TAIL_A
*

*TAIL_B
*

*CONSTR Section

Specifies geometric constraints (fixed internal coordinates) or PES scan specifications. See Constraints & Scans for the full syntax reference.

Example — fix a bond and an angle:

*CONSTR
R 1 2 1.54
A 1 2 3 109.5
*

Example — 1D bond scan:

*CONSTR
s r 1 2 1.0 20 0.05
*

*LST_A and *LST_B Sections

Provide the endpoint geometry for Linear Synchronous Transit (LST) interpolation. When both sections are present, OpenMECP generates an interpolated initial guess by linearly mixing the *GEOM geometry and the LST target.

Format — same as *GEOM (Cartesian coordinates, Ångström):

*LST_A
C  0.0  0.0  0.0
H  1.0  0.0  0.0
*

*LST_B
C  0.0  0.0  0.5
H  1.2  0.0  0.5
*

Atom ordering must match between *GEOM, *LST_A, and *LST_B. OpenMECP applies the Kabsch algorithm to find the optimal rigid-body alignment before interpolating, so pre-alignment of the structures is not required.

See LST Interpolation for full details.

Keywords Reference

All keywords are case-insensitive and use key = value syntax. Boolean values accept true/false.

Required Keywords

Program-Specific Method Format

Gaussian — method goes in the route section:

method = n scf(maxcycle=500,xqc) uwb97xd/def2svpp scrf=(smd,solvent=acetonitrile)
mem    = 8GB

ORCA — method is the ! line (without the !):

method = B3LYP def2-SVP CPCM(2-octanone) VeryTightSCF
mem    = 8000

For ORCA, mem is the per-core memory in MB (%maxcore in the ORCA input).

TD-DFT and State Selection

td_state_a and td_state_b are accepted as aliases for td_a and td_b.

Run Mode and Basic Optimization

Optimizer Settings

Trust Radius (GDIIS_blend)

Convergence Thresholds

Robust DIIS (Step Validation)

These options are only active when use_robust_diis = true.

See Robust DIIS for detailed descriptions.

Optimizer Switching Control

Constraints and Fixed Atoms

Coordinate Driving

Program Commands

ONIOM (QM/MM)

Advanced / Miscellaneous

Constraints & Scans

Constraints and PES scan specifications go inside the *CONSTR section. Constraints use Lagrange multipliers via the augmented Hessian method, with analytical gradients for all constraint types.

Bond Constraint

Fixes a bond length between two atoms at a target distance (Ångström).

R atom1 atom2 target_distance

Example:

*CONSTR
R 1 2 1.54    # Fix C–C bond at 1.54 Å
*

Angle Constraint

Fixes a bond angle defined by three atoms at a target value (degrees).

A atom1 atom2 atom3 target_angle

atom2 is the central atom.

Example:

*CONSTR
A 1 2 3 109.5    # Fix angle at 109.5°
*

Dihedral Constraint

Fixes a dihedral angle defined by four atoms at a target value (degrees).

D atom1 atom2 atom3 atom4 target_dihedral

Example:

*CONSTR
D 1 2 3 4 180.0    # Fix dihedral at 180°
*

Fixed Atoms

Individual atoms can be frozen using the fixedatoms keyword (not inside *CONSTR). Ranges are supported.

fixedatoms = 1,3-5,7    # Fix atoms 1, 3, 4, 5, and 7

PES Scans

Scan specifications are placed inside *CONSTR and trigger a PES scan run instead of a single MECP optimization.

Bond Scan

s r atom1 atom2 start num_points step_size

Example — scan a C–H bond from 1.0 to 1.95 Å in 20 steps of 0.05 Å:

*CONSTR
s r 1 2 1.0 20 0.05
*

Angle Scan

s a atom1 atom2 atom3 start num_points step_size

Example — scan from 90° to 135° in 10 steps of 5°:

*CONSTR
s a 1 2 3 90 10 5
*

Dihedral Scan

s d atom1 atom2 atom3 atom4 start num_points step_size

Example — full dihedral rotation from 0° to 350° in 36 steps of 10°:

*CONSTR
s d 1 2 3 4 0 36 10
*

2D Scan

Specify two scan lines in *CONSTR for a 2D scan:

*CONSTR
s r 1 2 1.0 10 0.1    # First dimension: bond
s a 1 2 3 90 5 10     # Second dimension: angle
*

Output: Geometries saved as scan_{value1}_{value2}.xyz.

Combining Constraints

Multiple constraints of different types can be combined freely:

*CONSTR
R 1 2 1.54      # Fix bond
A 1 2 3 109.5   # Fix angle
D 1 2 3 4 60.0  # Fix dihedral
*

MECP Algorithm

OpenMECP implements the MECP optimization algorithm of Harvey et al. (Theor. Chem. Acc. 99, 95–99, 1998) with several modern enhancements.

Problem Statement

Given two electronic states with energies $E_1(\mathbf{R})$ and $E_2(\mathbf{R})$ and gradients ${\mathbf{f}}_1(\mathbf{R})$ and ${\mathbf{f}}_2(\mathbf{R})$, find the geometry ${\mathbf{R}}^{\ast }$ that simultaneously satisfies:

1. Degeneracy: $E_1({\mathbf{R}}^{\ast })=E_2({\mathbf{R}}^{\ast })$
2. Minimum energy: ${\mathbf{R}}^{\ast }$ is the lowest-energy point on the crossing seam

Effective Gradient

Harvey et al. construct an effective gradient that encodes both conditions. First, define the normalized gradient difference:

$$\hat{\mathbf{x}}=\frac{{\mathbf{f}}_1-{\mathbf{f}}_2}{|{\mathbf{f}}_1-{\mathbf{f}}_2|}$$

The effective gradient has two orthogonal components:

f-vector — drives the energy gap to zero:

$$\mathbf{f}=(E_1-E_2)\hat{\mathbf{x}}$$

g-vector — minimizes the average energy in the intersection seam:

$$\mathbf{g}={\mathbf{f}}_1-(\hat{\mathbf{x}}\cdot {\mathbf{f}}_1)\hat{\mathbf{x}}$$

The total effective gradient is:

$${\mathbf{G}}_{\text{eff}}=\mathbf{f}+\mathbf{g}$$

Convergence Criteria

All five criteria must be satisfied simultaneously:

Optimization Strategy

OpenMECP uses a hybrid strategy combining stability and speed:

Step 1–switch_step:  BFGS (builds Hessian curvature information)
Step switch_step+1:  GDIIS / GEDIIS (quadratic-convergence DIIS acceleration)

The default switch_step = 3 provides a 3-step BFGS warm-up before engaging DIIS. See Optimizers for a full description of each algorithm.

Two-State QM Calculation

At each MECP optimization step, OpenMECP:

1. Writes the QM input files for state A and state B from the current geometry.
2. Runs the QM program for state A, reads $E_1$ and ${\mathbf{f}}_1$.
3. Runs the QM program for state B, reads $E_2$ and ${\mathbf{f}}_2$.
4. Evaluates the MECP effective gradient ${\mathbf{G}}_{\text{eff}}$.
5. Updates the geometry using the selected optimizer.
6. Checks convergence; writes checkpoint file if requested.

Unit Conventions

Optimizers

OpenMECP provides several geometry optimizers. All optimizers work on the MECP effective gradient described in the MECP Algorithm chapter.

BFGS (Initial Steps)

The Broyden–Fletcher–Goldfarb–Shanno quasi-Newton method provides a stable warm-up for the first switch_step iterations (default: 3). It builds up curvature information in the approximate Hessian before DIIS takes over.

When to use BFGS exclusively: Set switch_step = 999 (or any value ≥ max_steps) to run pure BFGS throughout — most stable but slower.

Relevant keywords: hessian, bfgs_rho, max_step_size

GDIIS — Geometry DIIS (Default)

use_gediis = false (default)

The Geometry Direct Inversion in the Iterative Subspace method constructs the next geometry as an optimal linear combination of previous Newton-step corrections:

$${\mathbf{e}}_i={\mathbf{H}}^{-1}\cdot ({\mathbf{g}}_i+{\mathbf{f}}_i)$$

$${\mathbf{B}}_{jk}=⟨{\mathbf{e}}_j\mid {\mathbf{e}}_k⟩$$

$$\left(\begin{array}{cc}\mathbf{B}& 1\\ 1^{\top }& 0\end{array}\right)\left(\begin{array}{c}\mathbf{c}\\ \lambda \end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)$$

$${\mathbf{x}}^{\ast }=\sum _i{c}_i{\mathbf{x}}_i,{\mathbf{x}}_{\text{GDIIS}}={\mathbf{x}}^{\ast }-{\mathbf{H}}^{-1}\cdot P({\mathbf{x}}^{\ast })$$

GDIIS minimizes residual Newton-step forces and delivers robust quadratic convergence. Proven reliable for ~10-step convergence on production systems.

Performance: ~2–3× faster than pure BFGS.

GEDIIS — Energy-Informed DIIS

use_gediis = true, use_hybrid_gediis = false

Extends GDIIS by adding an energy-based diagonal bias to the B-matrix:

$$B_{ii}+=\alpha \cdot E_i^2,{\text{RHS}}_i=-E_i$$

where $E_i=|E_{1,i}-E_{2,i}|$ is the energy gap at step $i$. The diagonal term biases interpolation toward geometries with small energy gaps (i.e., close to the crossing seam). The constraint is $c_i>0$ — GEDIIS only interpolates, never extrapolates.

Best for: systems where the crossing seam needs to be approached from a geometry far from it.

Sequential Hybrid GEDIIS/GDIIS

use_gediis = true, use_hybrid_gediis = true

Implements the 3-phase switching strategy of Li & Frisch (JCTC 2006):

Default thresholds: gediis_switch_rms = 0.005 Ha/Å, gediis_switch_step = 0.001 Å.

Performance: ~2–4× faster than pure GDIIS on most systems.

GDIIS_blend

use_gediis = "blend"

A per-step blend of a GDIIS step and an EDIIS step, governed by a trust radius:

$${\mathbf{x}}_{\text{new}}=w\cdot {\mathbf{x}}_{\text{EDIIS}}+(1-w)\cdot {\mathbf{x}}_{\text{GDIIS}}$$

When use_hybrid_gediis = false, no EDIIS component is added (pure GDIIS with trust radius control). When use_hybrid_gediis = true, the blend weight $w$ is determined by gediis_blend_mode:

The trust radius automatically contracts when energy increases and expands when energy decreases, preventing wild steps.

Requirement: A direct Hessian update method (direct_psb, bofill, powell, or bfgs_powell_mix). inverse_bfgs is incompatible.

Optimizer Switching

The switch_step parameter controls when the warm-up BFGS steps hand off to DIIS:

switch_step = 3     # Default: 3 BFGS steps, then DIIS
switch_step = 0     # Pure DIIS from step 1 (fastest, needs good geometry)
switch_step = 10    # 10 BFGS steps, then DIIS (better for difficult cases)
switch_step = 999   # Pure BFGS (most stable, slowest)

The smart_history option retains the most useful geometry history points rather than using a simple FIFO queue, which can accelerate convergence.

Performance Comparison

Hessian Updates

The hessian keyword selects how the Hessian (or inverse Hessian) matrix is stored and updated at each optimization step.

Available Methods

direct_psb — Powell-Symmetric-Broyden (Default)

hessian = direct_psb

Stores the direct Hessian $\mathbf{H}$ (units: Ha/Ų) and solves the linear system $\mathbf{H}\cdot {\mathbf{d}}_k=-\mathbf{g}$ at each step via LU decomposition.

The PSB update formula preserves the secant condition while maintaining symmetry:

$${\mathbf{H}}_{k+1}={\mathbf{H}}_k+\frac{(\mathbf{y}-{\mathbf{H}}_k\mathbf{s}){\mathbf{s}}^{\top }+\mathbf{s}(\mathbf{y}-{\mathbf{H}}_k\mathbf{s}{)}^{\top }}{{\mathbf{s}}^{\top }\mathbf{s}}-\frac{(\mathbf{y}-{\mathbf{H}}_k\mathbf{s}{)}^{\top }\mathbf{s}}{{\left({\mathbf{s}}^{\top }\mathbf{s}\right)}^2}\mathbf{s}{\mathbf{s}}^{\top }$$

where $\mathbf{s}={\mathbf{R}}_{k+1}-{\mathbf{R}}_k$ and $\mathbf{y}={\mathbf{g}}_{k+1}-{\mathbf{g}}_k$.

Required when use_gediis = blend. Recommended for all production runs.

inverse_bfgs — Inverse Hessian BFGS

hessian = inverse_bfgs

Stores the inverse Hessian ${\mathbf{H}}^{-1}$ (units: Ų/Ha) and computes steps as ${\mathbf{d}}_k=-{\mathbf{H}}^{-1}\mathbf{g}$ via matrix-vector multiply (no linear solve required).

Uses the standard BFGS rank-2 update:

$${\mathbf{H}}_{k+1}^{-1}=\left(\mathbf{I}-{\rho }_k\mathbf{s}{\mathbf{y}}^{\top }\right){\mathbf{H}}_k^{-1}\left(\mathbf{I}-{\rho }_k\mathbf{y}{\mathbf{s}}^{\top }\right)+{\rho }_k\mathbf{s}{\mathbf{s}}^{\top },{\rho }_k=\frac{1}{{\mathbf{y}}^{\top }\mathbf{s}}$$

Note: Incompatible with use_gediis = blend. Use direct_psb instead for the blend optimizer.

bofill — Bofill Update (Experimental)

hessian = bofill

Blends the Powell (SR1) and Murtagh–Sargent updates using a Bofill-style weighting parameter $\varphi \in [0,1]$:

$${\mathbf{H}}_{k+1}=\varphi \cdot {\mathbf{H}}_{k+1}^{\text{Powell}}+(1-\varphi )\cdot {\mathbf{H}}_{k+1}^{\text{MS}}$$

The weight $\varphi$ is computed from the curvature: $\varphi$ is close to 1 when the gradient change is well-described by a symmetric rank-1 correction (Powell), and close to 0 otherwise (Murtagh–Sargent).

Recommended for: TS-like crossing points where the Hessian may have negative eigenvalues.

powell — Powell Symmetric Rank-1 (Experimental)

hessian = powell

Uses the symmetric SR1 update. Unlike BFGS, SR1 can build up negative curvature, which is useful for crossing points with Hessians that have negative eigenvalues.

bfgs_powell_mix — Adaptive Blend (Experimental)

hessian = bfgs_powell_mix

Adaptively switches between BFGS and Powell updates using Bofill-style weighting at each step. Inherits the robustness of Powell near negative curvature and the superlinear convergence of BFGS near positive-definite regions.

Choosing a Method

Robust DIIS

Standard GDIIS/GEDIIS can occasionally produce unreliable steps when the DIIS matrix is ill-conditioned — for example, near convergence or when the history contains redundant points. The Robust DIIS system adds configurable sanity checks that detect and reject bad steps before they are applied.

Enabling Robust DIIS

use_robust_diis = true

This is the master switch. All options described on this page are inactive unless use_robust_diis = true.

Cosine Check

gdiis_cosine_check = standard    # default

Computes the cosine between the proposed DIIS step and the reference gradient. A low cosine means DIIS is stepping in a direction inconsistent with the local gradient — a sign that the step is unreliable and should be rejected.

Coefficient Check

gdiis_coeff_check = regular    # default

Validates that the DIIS interpolation coefficients $c_i$ are physically reasonable. Well-conditioned DIIS produces coefficients that sum to 1 with individual values near $[0,1]$.

GEDIIS Variants

gediis_variant = auto    # default

Only active when use_robust_diis = true and use_gediis = true.

The gediis_sim_switch parameter (default: 0.0025) controls the RMS error threshold at which auto switches from energy to rfo.

Negative Eigenvalue Control

n_neg = 0    # default (minimum search)
n_neg = 1    # transition-state search (one imaginary mode)

Set n_neg = 1 for TS-like MECP searches where the crossing point is a saddle point with one negative Hessian eigenvalue.

Worked Examples

Oscillating Convergence

If the optimizer oscillates or takes wild steps near convergence:

use_robust_diis    = true
gdiis_cosine_check = standard
gdiis_coeff_check  = regular

Fragile System (Flat PES / Near-Dissociation)

use_robust_diis    = true
gdiis_cosine_check = strict

Restart from Checkpoint

use_robust_diis   = true
gdiis_coeff_check = force_recent
restart           = true

Full Robust GEDIIS

use_robust_diis    = true
use_gediis         = true
gediis_variant     = energy
gdiis_cosine_check = standard
gdiis_coeff_check  = regular

QM Program Support

OpenMECP communicates with external quantum chemistry programs through a unified interface. Each program has its own input/output format; OpenMECP handles all the file writing and parsing automatically.

Supported Programs

Selecting a Program

Set the program keyword in your input file:

program = gaussian    # or orca / custom

Common Configuration

Both programs share these keywords:

nprocs  = 8       # Number of parallel CPU cores
charge  = 0       # Molecular charge
mult_a  = 1       # Multiplicity of state A
mult_b  = 3       # Multiplicity of state B

See individual program pages for method and mem format details:

• Gaussian
• ORCA
• Custom Interface

Wavefunction Checkpoints

OpenMECP saves wavefunction files from each QM step to use as a starting guess for the next step. This dramatically speeds up SCF convergence over a multi-step optimization.

The Run Mode controls whether these checkpoints are read.

File Organization

All QM input/output files are written to a running_dir/ subdirectory, keeping your working directory clean. The final MECP geometry is written to {basename}_mecp.xyz in the current directory.

Gaussian

Setting Up a Gaussian Calculation

program = gaussian
method  = B3LYP/6-31G*
nprocs  = 8
mem     = 16GB
charge  = 0
mult_a  = 1
mult_b  = 3

method Format

The method value is inserted directly into the Gaussian route section (# line). You can include any valid Gaussian route keywords:

method = n scf(maxcycle=500,xqc) uwb97xd/def2svpp scrf=(smd,solvent=acetonitrile)

mem Format

Specify memory using Gaussian's syntax — e.g. 4GB, 8GB, 4000MB.

Executable

By default, OpenMECP calls g16. Override with:

gau_comm = g09    # or full path: /opt/g16/g16

TD-DFT and Excited States

Use td_a / td_b to inject TD-DFT keywords into the route section for each state, and state_a / state_b to select which root to use for the gradient:

td_a   = TD(nstates=5,root=1)
td_b   = TD(nstates=5,root=2)
state_a = 1
state_b = 2

MP2

Enable MP2 gradients (instead of DFT) with:

mp2 = true

CASSCF and Multi-Reference Methods

Put CASSCF keywords in the *TAIL_A / *TAIL_B sections as extra route or additional input blocks:

*TAIL_A
casscf(8,8)/cc-pVDZ
*

*TAIL_B
casscf(8,8)/cc-pVDZ
*

Extra Input via *TAIL Sections

Anything in *TAIL_A / *TAIL_B is appended to the Gaussian route section keywords for the respective state. Common uses:

*TAIL_A
scf(maxcycle=500,xqc) stable=opt
*

ONIOM (QM/MM)

Multi-layer ONIOM calculations are supported. Mark atom layers in *GEOM using Gaussian's H / M / L layer notation:

*GEOM
C  0.0  0.0  0.0  H
H  1.0  0.0  0.0  L
H  0.0  1.0  0.0  L
*

isoniom = true
chargeandmultforoniom1 = 0 1 0 1   # state A: QM and MM charges/mults
chargeandmultforoniom2 = 0 3 0 1   # state B

Run Modes and Checkpoints

OpenMECP automatically manages Gaussian checkpoint files (.chk) to pass the wavefunction between steps. The mode keyword controls this behavior — see Run Modes for details.

The inter_read mode is essential for open-shell singlet (OSS) calculations: it copies the triplet wavefunction to the singlet and adds guess=(read,mix).

Output Files

ORCA

Setting Up an ORCA Calculation

program = orca
method  = B3LYP def2-SVP
nprocs  = 8
mem     = 8000
charge  = 0
mult_a  = 1
mult_b  = 3

method Format

The method value is inserted after the ! on the ORCA simple input line. Include all keywords that would normally follow ! in your ORCA input:

method = B3LYP def2-SVP CPCM(2-octanone) VeryTightSCF

mem Format

Specify per-core memory in MB — this maps to %maxcore in the ORCA input:

mem = 8000    # 8000 MB per core

Executable

By default, OpenMECP calls orca. Override with:

orca_comm = /opt/orca5/orca    # full path recommended for ORCA

TD-DFT and Excited States

Use the *TAIL_A / *TAIL_B sections to inject %tddft blocks:

*TAIL_A
%tddft
  nroots 5
  iroot  1
end
*

*TAIL_B
%tddft
  nroots 5
  iroot  2
end
*

Combined with:

state_a = 1
state_b = 2

CASSCF

Put CASSCF blocks in the tail sections:

*TAIL_A
%casscf
  nel    8
  norb   8
  nroots 3
  iroot  0
end
*

Extra Input via *TAIL Sections

Anything in *TAIL_A / *TAIL_B is appended as an additional input block for the respective state. Common uses:

*TAIL_A
%scf
  MaxIter 500
  Convergence VeryTight
end
*

Output Files

Run Modes

ORCA reads wavefunction from .gbw files. The read and noread run modes control whether OpenMECP passes a MORead keyword to ORCA. See Run Modes for details.

Notes

• Ensure the orca executable is on your PATH or specify the full path via orca_comm.
• ORCA requires an exclusive MPI launch; do not wrap it in mpirun in your job script — ORCA handles parallel execution internally.
• The .engrad file is written only when ORCA is asked for gradients (the default for MECP runs).

Custom QM Interface

OpenMECP supports user-defined QM programs via a JSON configuration file. This allows connecting any quantum chemistry code that can produce energies and gradients to the MECP optimizer.

Enabling Custom Interface

program              = custom
custom_interface_file = my_program.json

JSON Configuration Format

Create a JSON file describing how to write inputs, run the program, and parse its outputs:

{
  "name": "my_qm_program",
  "executable": "/path/to/qm_program",
  "input_template": "template_state.inp",
  "output_extension": "out",
  "gradient_extension": "grad",
  "energy_pattern": "Total Energy:\\s+([\\-0-9\\.]+)",
  "gradient_pattern": "Gradient:\\s+([\\-0-9\\.eE+\\-]+)",
  "extra_args": []
}

Fields

Creating a Template

Generate a starting JSON template with:

omecp ci custom_interface.json

Edit the generated file to match your program's input/output format.

Notes

• The regex patterns use Rust regex syntax.
• The first capturing group in energy_pattern must match a floating-point number (hartree).
• The gradient_pattern is applied repeatedly to extract all gradient components in the order: x₁ y₁ z₁ x₂ y₂ z₂ … for N atoms.
• OpenMECP writes a separate input file for state A and state B and calls the executable twice per optimization step.

Advanced Features

This chapter covers OpenMECP's advanced capabilities beyond basic MECP optimization.

State Selection for TD-DFT

Choose specific excited states for TD-DFT calculations using state_a and state_b:

state_a = 1    # Use first excited state for state A
state_b = 2    # Use second excited state for state B
td_a = TD(nstates=5,root=1)
td_b = TD(nstates=5,root=2)

This is particularly important when:

• The ground state is contaminated by near-degeneracy
• Targeting a specific conical intersection between excited states
• Studying multi-state photochemical reactivity

Pre-point Calculations

Pre-point calculations run a single-point QM job at the starting geometry before the first optimization step. Enable with mode = stable (see Run Modes). This checks SCF stability and ensures that wavefunction checkpoints are available from step 1.

Restart from Checkpoint

restart = true

OpenMECP saves a JSON checkpoint ({basename}_checkpoint.json) after each step that records the full optimization state: geometry, DIIS history, Hessian, and configuration. Setting restart = true picks up from the last saved checkpoint.

Enable checkpoint writing:

print_checkpoint = true
restart          = true

Debug Logging

Enable file-based debug logging via omecp_config.cfg:

[logging]
file_logging = true

Log files are written to omecp_debug_{basename}.log in the working directory.

PES Scan

OpenMECP can perform 1D and 2D potential energy surface scans by running a constrained MECP optimization at each point on a grid.

1D Bond Scan

Specify a scan in the *CONSTR section:

*CONSTR
s r atom1 atom2 start num_points step_size
*

Example — scan a C–C bond from 1.0 to 1.95 Å in 20 steps of 0.05 Å:

*CONSTR
s r 1 2 1.0 20 0.05
*

Final geometry at each scan point: scan_1.00.xyz, scan_1.05.xyz, …

1D Angle Scan

*CONSTR
s a atom1 atom2 atom3 start num_points step_size
*

Example — scan a C–C–H angle from 90° to 135° in 10 steps of 5°:

*CONSTR
s a 1 2 3 90 10 5
*

1D Dihedral Scan

*CONSTR
s d atom1 atom2 atom3 atom4 start num_points step_size
*

Example — full dihedral rotation from 0° to 350° in 36 steps of 10°:

*CONSTR
s d 1 2 3 4 0 36 10
*

2D Scan

Specify two scan lines in *CONSTR:

*CONSTR
s r 1 2 1.0 10 0.1    # First dimension: C–C bond, 10 points
s a 1 2 3 90 5 10     # Second dimension: C–C–H angle, 5 points
*

This performs 10 × 5 = 50 MECP optimizations on a grid. Output geometries are named scan_{val1}_{val2}.xyz.

Output

The energy values printed to stdout can be used directly to plot the PES scan.

Tips

• For 2D scans, reduce max_steps per point to keep the total runtime manageable: max_steps = 50.
• Use print_level = 0 to suppress per-step output and only log the final energies for each scan point.
• Converged MECP geometries from one scan point are used as the starting geometry for the next, which generally speeds up convergence along the scan.

LST Interpolation

Linear Synchronous Transit (LST) interpolation generates an initial MECP guess by linearly mixing between two endpoint geometries. This is useful when the MECP lies somewhere between a known reactant geometry and a known product geometry.

Setup

Add *LST_A and *LST_B sections alongside the standard *GEOM:

*GEOM
C  -4.10   -0.27   -0.11
H  -4.46    0.46    0.62
*

*LST_A
C  0.0  0.0  0.0
H  1.0  0.0  0.0
*

*LST_B
C  0.0  0.0  0.5
H  1.2  0.0  0.5
*

When both *LST_A and *LST_B sections are present, OpenMECP:

1. Applies the Kabsch algorithm to optimally align both structures to the *GEOM reference (removes translational and rotational differences).
2. Generates a sequence of 10 interpolation points along the LST path.
3. Prints the energy profile along the interpolated path.
4. Asks for interactive confirmation before proceeding with MECP optimization.

Interpolation Methods

Note: QST currently uses a midpoint approximation for the third geometry. Full QST with an explicit user-supplied midpoint is planned for a future release.

Kabsch Alignment

The Kabsch algorithm finds the optimal rotation matrix $\mathbf{R}$ that minimizes the RMSD between corresponding atoms of the two structures. This ensures that the interpolation path follows the shortest physically meaningful route between the two geometries, avoiding spurious rotations.

No manual pre-alignment of the structures is required — OpenMECP handles this automatically.

Output

After LST interpolation:

Tips

• LST is most effective when the two endpoint geometries differ mainly in bond lengths or angles, not overall topology.
• If the energy profile shows a clear minimum-energy crossing along the LST path, the starting geometry for MECP optimization will be close to the MECP, leading to faster convergence.
• Use print_level = 2 to see the full interpolated path during the LST stage.

Reaction Path & NEB

Coordinate Driving

Coordinate driving systematically varies an internal coordinate (bond, angle, or dihedral) in a series of constrained MECP optimizations. This traces a reaction path on the crossing seam.

Keywords

drive_type       = bond        # bond | angle | dihedral
drive_atoms      = 1,2         # comma-separated atom indices
drive_start      = 1.0         # starting value (Å or degrees)
drive_end        = 2.5         # ending value
drive_steps      = 20          # number of steps

Example — C–N Bond Elongation

drive_type   = bond
drive_atoms  = 1,2
drive_start  = 1.4
drive_end    = 2.4
drive_steps  = 20

At each of the 20 steps, OpenMECP runs a constrained MECP optimization with the C–N bond fixed at the current driving value. The converged geometry is then used as the starting point for the next step.

Output: Energy profile along the driven coordinate printed to stdout; geometries saved to running_dir/.

Nudged Elastic Band (NEB)

NEB optimizes an entire reaction path simultaneously by connecting a chain of geometries ("images") with spring forces that keep them evenly spaced along the path.

Setup

drive_type   = bond
drive_atoms  = 1,2
drive_start  = 1.0
drive_end    = 2.0
drive_steps  = 10    # number of images on the path

When the run_mode is set to path_optimization, OpenMECP:

1. Creates an initial path via coordinate driving.
2. Applies simplified NEB spring forces between adjacent images.
3. Optimizes each image's position on the crossing seam.

Current Limitations: The NEB implementation uses a simplified spring model. A full NEB with perpendicular energy gradient projection is planned for a future release.

Output: Smooth energy profile along the optimized path; geometry files for each image.

Choosing Between Methods

Fixed Atoms

Selected atoms can be frozen during MECP optimization, allowing partial geometry optimization. Fixed atoms do not contribute to the MECP gradient and their coordinates are never updated by the optimizer.

Usage

Use the fixedatoms keyword (not inside *CONSTR). Specify atom indices using comma-separated values and ranges:

fixedatoms = 1,3-5,7    # Fix atoms 1, 3, 4, 5, and 7

Atom indices are 1-based and correspond to the order of atoms in the *GEOM section.

Syntax

Use Cases

Partial optimization — optimize only the active site of a large molecule:

fixedatoms = 10-50    # Freeze the molecular scaffold; optimize only atoms 1–9

QM/MM boundaries — freeze the MM region atoms while optimizing the QM core:

fixedatoms = 15-100

Surface calculations — fix the surface slab atoms while relaxing the adsorbate:

fixedatoms = 1-48    # 48-atom surface slab; adsorbate is atoms 49–54

Interaction with Constraints

Fixed atoms and geometric constraints (*CONSTR) can be combined in the same calculation. A fixed atom's coordinate is not optimized, but it can still appear as part of a constraint definition (e.g., as the anchor atom of a bond constraint).

Notes

• Fixing a large fraction of atoms significantly reduces the effective degrees of freedom, which can speed up convergence.
• The DIIS history and Hessian are built only in the subspace of free atoms.
• Fixed atoms are included in the geometry written to output files for visual inspection.

Examples Overview

This section provides complete, ready-to-run examples for locating MECPs with OpenMECP. All examples are also available in the examples/ directory of the repository.

Test System: Phenylcation (C₆H₅⁺)

The primary test case used throughout these examples is the phenylcation, a well-studied 11-atom aromatic cation. The goal is to locate the minimum energy crossing point between the triplet (S = 3) and open-shell singlet (S = 1) spin states.

This system is an excellent benchmark: it is small enough for rapid testing on any hardware, yet exercises all parts of the MECP optimizer, including gradient combination, Hessian updates, and DIIS extrapolation.

Geometry

Save the following as phenylcation.xyz (or embed it directly in the *geom section of your input):

11

C     0.000000000000      1.396613000000      0.000000000000
C     1.209503000000      0.698307000000      0.000000000000
C     1.209503000000     -0.698307000000      0.000000000000
C     0.000000000000     -1.396613000000      0.000000000000
C    -1.209503000000     -0.698307000000      0.000000000000
C    -1.209503000000      0.698307000000      0.000000000000
H     2.150882000000      1.241812000000      0.000000000000
H     2.150882000000     -1.241812000000      0.000000000000
H     0.000000000000     -2.483625000000      0.000000000000
H    -2.150882000000     -1.241812000000      0.000000000000
H    -2.150882000000      1.241812000000      0.000000000000

How to Run

omecp input.inp > output.log

OpenMECP writes a summary line to stdout at each step. Key fields:

Output Files

Available Examples

Gaussian Examples

All examples locate the triplet/singlet (S=3 / S=1) MECP of the phenylcation (C₆H₅⁺, charge = 1). See Overview for the geometry and how to run.

Example 1 — Default GDIIS

The simplest production-ready input. Uses the default GDIIS optimizer with 3 warm-up BFGS steps.

nprocs = 10
mem = 40GB
method = n b3lyp/6-31g**
charge = 1
mult_a = 3
mult_b = 1
mode = normal
program = gaussian

*geom
C     0.000000000000      1.396613000000      0.000000000000
C     1.209503000000      0.698307000000      0.000000000000
C     1.209503000000     -0.698307000000      0.000000000000
C     0.000000000000     -1.396613000000      0.000000000000
C    -1.209503000000     -0.698307000000      0.000000000000
C    -1.209503000000      0.698307000000      0.000000000000
H     2.150882000000      1.241812000000      0.000000000000
H     2.150882000000     -1.241812000000      0.000000000000
H     0.000000000000     -2.483625000000      0.000000000000
H    -2.150882000000     -1.241812000000      0.000000000000
H    -2.150882000000      1.241812000000      0.000000000000
*

*tail_a
*

*tail_b
*

Key settings: switch_step = 3 (default), hessian = direct_psb (default), use_gediis = false (default).

Example 2 — Sequential Hybrid (GEDIIS → GDIIS)

Activates the 3-phase switching strategy of Li & Frisch (JCTC 2006). Adds only two lines to Example 1 and typically converges ~2–4× faster.

nprocs = 10
mem = 40GB
method = n b3lyp/6-31g**
charge = 1
mult_a = 3
mult_b = 1
mode = normal
program = gaussian

# --- optimizer ---
use_gediis = true
use_hybrid_gediis = true

*geom
C     0.000000000000      1.396613000000      0.000000000000
C     1.209503000000      0.698307000000      0.000000000000
C     1.209503000000     -0.698307000000      0.000000000000
C     0.000000000000     -1.396613000000      0.000000000000
C    -1.209503000000     -0.698307000000      0.000000000000
C    -1.209503000000      0.698307000000      0.000000000000
H     2.150882000000      1.241812000000      0.000000000000
H     2.150882000000     -1.241812000000      0.000000000000
H     0.000000000000     -2.483625000000      0.000000000000
H    -2.150882000000     -1.241812000000      0.000000000000
H    -2.150882000000      1.241812000000      0.000000000000
*

*tail_a
*

*tail_b
*

Phase behaviour: starts with GDIIS, switches to GEDIIS when rms_grad ≤ gediis_switch_rms (0.005 Ha/Å), then back to GDIIS for the final quadratic convergence when rms_disp < gediis_switch_step (0.001 Å).

Example 3 — GDIIS_blend with Trust Radius

Uses the blend optimizer with the fixed_sequential blend mode and a live trust radius. Suitable when the default GDIIS takes excessively large or oscillating steps.

nprocs = 10
mem = 40GB
method = n b3lyp/6-31g**
charge = 1
mult_a = 3
mult_b = 1
mode = normal
program = gaussian

# --- optimizer ---
use_gediis = blend
use_hybrid_gediis = true
gediis_blend_mode = fixed_sequential

*geom
C     0.000000000000      1.396613000000      0.000000000000
C     1.209503000000      0.698307000000      0.000000000000
C     1.209503000000     -0.698307000000      0.000000000000
C     0.000000000000     -1.396613000000      0.000000000000
C    -1.209503000000     -0.698307000000      0.000000000000
C    -1.209503000000      0.698307000000      0.000000000000
H     2.150882000000      1.241812000000      0.000000000000
H     2.150882000000     -1.241812000000      0.000000000000
H     0.000000000000     -2.483625000000      0.000000000000
H    -2.150882000000     -1.241812000000      0.000000000000
H    -2.150882000000      1.241812000000      0.000000000000
*

*tail_a
*

*tail_b
*

Blend behaviour: 50 % GDIIS + 50 % EDIIS far from convergence, transitions to pure GDIIS for the final steps. Trust radius contracts when energy rises and expands when energy falls.

Note: blend mode requires a direct Hessian method (hessian = direct_psb, bofill, powell, or bfgs_powell_mix). inverse_bfgs is incompatible.

Real-World Example — Cu Complex with Solvent

For a larger, more demanding system the to-4-tep-1 examples in examples/gaussian/ show a 29-atom Cu complex optimised at the ωB97X-D/def2-SVP level with SMD solvation in benzene:

nprocs = 48
mem = 190GB
method = n scf(maxcycle=300,xqc) UWB97XD/DEF2SVPP scrf(smd,solvent=Benzene)
charge = 1
mult_a = 3
mult_b = 1
mode = normal
program = gaussian
use_gediis = blend
use_hybrid_gediis = true
gediis_blend_mode = fixed_sequential

The geometry is read from the external file to-4-tep-1-blend-gdiis.xyz (included in the repository) via *geom\n@to-4-tep-1-blend-gdiis.xyz\n*.

ORCA Examples

All examples locate the triplet/singlet (S=3 / S=1) MECP of the phenylcation (C₆H₅⁺, charge = 1). See Overview for the geometry and how to run.

ORCA memory: pass the MaxCore value in MB, not the total job memory. For example, mem = 4000 means 4 000 MB per core.

Example 1 — Default GDIIS

nprocs = 10
mem = 4000
method = B3LYP 6-31G(d,p) VeryTightSCF
charge = 1
mult_a = 3
mult_b = 1
mode = normal
program = orca
orca_comm = orca

*geom
C     0.000000000000      1.396613000000      0.000000000000
C     1.209503000000      0.698307000000      0.000000000000
C     1.209503000000     -0.698307000000      0.000000000000
C     0.000000000000     -1.396613000000      0.000000000000
C    -1.209503000000     -0.698307000000      0.000000000000
C    -1.209503000000      0.698307000000      0.000000000000
H     2.150882000000      1.241812000000      0.000000000000
H     2.150882000000     -1.241812000000      0.000000000000
H     0.000000000000     -2.483625000000      0.000000000000
H    -2.150882000000     -1.241812000000      0.000000000000
H    -2.150882000000      1.241812000000      0.000000000000
*

*tail_a
*

*tail_b
*

Key settings: default GDIIS, hessian = direct_psb, switch_step = 3.

On HPC systems, set orca_comm to the full path of the ORCA binary, e.g.:

orca_comm = /apps/orca/6.1.1/orca

Example 2 — GDIIS_blend with Gradient-Based Weight

Uses use_gediis = blend with gediis_blend_mode = gradient. The blend weight $w$ varies smoothly as $w={\text{RMS}}_g/({\text{RMS}}_g+\text{switch\_rms})$, so the optimizer is EDIIS-heavy far from convergence and GDIIS-heavy close to it.

nprocs = 10
mem = 4000
method = B3LYP 6-31G(d,p) VeryTightSCF
charge = 1
mult_a = 3
mult_b = 1
mode = normal
program = orca
orca_comm = orca

# --- optimizer ---
use_gediis = blend
use_hybrid_gediis = true
gediis_blend_mode = gradient

*geom
C     0.000000000000      1.396613000000      0.000000000000
C     1.209503000000      0.698307000000      0.000000000000
C     1.209503000000     -0.698307000000      0.000000000000
C     0.000000000000     -1.396613000000      0.000000000000
C    -1.209503000000     -0.698307000000      0.000000000000
C    -1.209503000000      0.698307000000      0.000000000000
H     2.150882000000      1.241812000000      0.000000000000
H     2.150882000000     -1.241812000000      0.000000000000
H     0.000000000000     -2.483625000000      0.000000000000
H    -2.150882000000     -1.241812000000      0.000000000000
H    -2.150882000000      1.241812000000      0.000000000000
*

*tail_a
*

*tail_b
*

TD-DFT with State Selection

For excited-state MECPs, add %tddft blocks in the tail sections and set state_a / state_b to select which excited state to track:

nprocs = 10
mem = 4000
method = B3LYP def2-SVP VeryTightSCF
charge = 0
mult_a = 1
mult_b = 1
state_a = 1
state_b = 2
mode = normal
program = orca
orca_comm = orca

*geom
...
*

*tail_a
%tddft
  nroots 5
  iroot 1
end
*

*tail_b
%tddft
  nroots 5
  iroot 2
end
*

Troubleshooting

Common Errors

"No such file or directory"

Cause: QM program not found in PATH, or input file is missing.

Solutions:

• Verify the QM program is installed and accessible in your $PATH.
• Check that the input file path is spelled correctly.

"QM calculation failed"

Cause: The QM program encountered an internal error.

Solutions:

• Open the running_dir/*.log (or *.out) files and look for program-specific error messages.
• Confirm that the method and basis set are valid for the chosen program.
• Check that nprocs and mem do not exceed available resources.
• Switch to mode = noread to force a fresh SCF each step and avoid reading a corrupted wavefunction:

mode = noread

Program-specific hints:

"Failed to parse output"

Cause: The QM output format was not recognized — often because the calculation did not finish normally.

Solutions:

• Confirm the QM calculation completed without errors before OpenMECP tries to read it.
• Check that the program version is supported.

Parsing-specific tips:

"Maximum steps exceeded"

Cause: The optimizer did not converge within max_steps (default 100).

Solutions:

• Increase the step limit:

  max_steps = 300
  

• Verify that the starting geometry is chemically reasonable.
• Use LST interpolation to generate a better initial guess.
• Check for recurring SCF convergence failures — if the QM energy/gradient is wrong the optimizer cannot converge.

SCF Convergence Issues

Symptoms: QM calculations fail repeatedly, or the energy oscillates step-to-step.

Solutions:

1. Use mode = noread — starts every SCF from scratch, avoids corrupted guess orbitals:

   mode = noread
   

2. Use mode = stable — checks wavefunction stability before proceeding:

   mode = stable
   

3. Use mode = inter_read — recommended for open-shell singlets, reads MOs from the previous step to maintain spin symmetry:

   mode = inter_read
   

4. Increase the SCF iteration limit via the *TAIL section:

   • Gaussian:

     scf=(maxcycle=200,xqc)
     

   • ORCA:

     %scf maxiter 200 end
     

See Run Modes for a full comparison of mode options.

Slow Convergence

Symptoms: Optimization requires an unusually large number of steps.

Step 1 — Increase the step limit

max_steps = 300

Step 2 — Try a different optimizer

The examples below go from fastest/least-robust to most-robust. Pick the first one that converges for your system.

Pure GDIIS (default)

No changes needed — GDIIS is the default. Explicitly:

use_gediis = false   # GDIIS only (default)
switch_step = 3      # 3 BFGS warm-up steps, then DIIS (default)

Sequential Hybrid GEDIIS/GDIIS

Three-phase switching: GDIIS → GEDIIS (near seam) → GDIIS (final). Typically 2–4× faster than pure GDIIS.

use_gediis = true
use_hybrid_gediis = true
gediis_switch_rms  = 0.005   # GDIIS → GEDIIS when RMS grad < 0.005 Ha/Å (default)
gediis_switch_step = 0.001   # GEDIIS → GDIIS when RMS disp < 0.001 Å (default)

GDIIS_blend — Pure (trust-radius GDIIS, no EDIIS)

GDIIS protected by an adaptive trust radius. Very robust; no EDIIS component.

switch_step = 3
hessian = direct_psb      # required for blend mode
use_gediis = blend
use_hybrid_gediis = false  # default for blend; can be omitted

GDIIS_blend — Fixed-Sequential Hybrid (recommended blend default)

50/50 GDIIS+EDIIS far from the minimum, transitions to pure GDIIS near convergence. Best overall blend mode for production runs.

switch_step = 3
hessian = direct_psb
use_gediis = blend
use_hybrid_gediis = true
gediis_blend_mode = fixed_sequential   # default when use_hybrid_gediis = true

GDIIS_blend — Gradient-Weighted Hybrid

Smooth sigmoid blend driven by the RMS gradient. EDIIS-heavy when forces are large, GDIIS-heavy near convergence. Useful for systems with energy plateaus.

switch_step = 3
hessian = direct_psb
use_gediis = blend
use_hybrid_gediis = true
gediis_blend_mode = gradient

GDIIS_blend — Sequential Hybrid

Per-step binary switching (GDIIS or EDIIS) based on the RMS displacement trend.

switch_step = 3
hessian = direct_psb
use_gediis = blend
use_hybrid_gediis = true
gediis_blend_mode = sequential

Optimizer Quick-Reference

Note: All blend modes require a direct Hessian update method: direct_psb (default), bofill, powell, or bfgs_powell_mix. inverse_bfgs is incompatible with blend mode.

Hessian Update Methods

The hessian keyword affects both step quality and compatibility with blend mode. Choose a method suited to your system:

For transition-state-like crossings, bofill is recommended:

hessian = bofill

Robust DIIS

For particularly difficult cases, enable the Robust DIIS sanity-check layer, which detects and rejects ill-conditioned DIIS steps:

use_robust_diis = true

Combine with any optimizer above for extra protection. See the Robust DIIS page for full options.

Run Modes

The mode keyword controls how OpenMECP manages QM wavefunction checkpoints between optimization steps. Choosing the right run mode can significantly affect SCF convergence speed and stability.

mode = normal    # default

Normal Mode (Default)

mode = normal

Standard MECP optimization. OpenMECP reads the wavefunction checkpoint from the previous step as the SCF initial guess.

• Reads .chk (Gaussian) or .gbw (ORCA) from the previous step
• Recommended for most calculations
• Fastest per-step runtime

Read Mode

mode = read

Reads existing wavefunction checkpoints from disk without running pre-point calculations. Use this when restarting a calculation where checkpoints already exist from a previous run.

• Skips pre-point calculations
• Requires existing checkpoint files
• Useful for restarting after a crash or job time limit

NoRead Mode

mode = noread

Runs a fresh SCF at every step without reading any checkpoint. This is slower but more robust for systems with difficult SCF convergence or when the wavefunction character changes significantly between steps.

• Does not read or write checkpoint files
• Slower (full SCF from scratch each step)
• Use when checkpoint reading causes SCF divergence

Stable Mode

mode = stable

Runs pre-point calculations at the initial geometry before starting the MECP optimization. This checks SCF stability (stable=opt in Gaussian) and produces reliable initial checkpoints.

• Runs stable=opt or equivalent before step 1
• Ensures stable wavefunction from the start
• Recommended for systems with unstable ground-state wavefunctions (e.g., open-shell systems)

InterRead Mode

mode = inter_read

Essential for open-shell singlet (OSS) states.

InterRead runs state B (typically the triplet) first at each step, then copies the triplet wavefunction as the initial guess for state A (the open-shell singlet) and adds guess=(read,mix) (Gaussian) or equivalent. This provides a physically correct broken-symmetry starting point for the singlet.

Example — open-shell singlet/triplet MECP:

mode   = inter_read
mult_a = 1    # Open-shell singlet
mult_b = 3    # Triplet (reference for broken-symmetry guess)

Mode Comparison

Configuration File

OpenMECP reads a configuration file omecp_config.cfg at startup to customize program-wide behavior. Settings in the config file apply to all calculations run from that directory.

Creating a Config Template

omecp ci omecp_config.cfg

This generates a fully commented template with all available options.

File Locations and Precedence

Configuration files are loaded in hierarchical order. A setting in a higher-priority file overrides the same setting in a lower-priority file.

File Format

The config file uses INI format with four sections:

[extensions]
gaussian = log
orca     = out
custom   = log

[general]
max_memory     = 4GB
default_nprocs = 4
print_level    = 0

[logging]
level        = info
file_logging = false

[cleanup]
enabled              = true
preserve_extensions  = xyz,backup
verbose              = 1
cleanup_frequency    = 5

[extensions] Section

Configures the file extensions that OpenMECP looks for when parsing QM output.

[general] Section

[logging] Section

Enabling file_logging = true writes a detailed debug log named omecp_debug_{basename}.log (where {basename} is derived from the input file name). This is useful for diagnosing convergence issues.

[cleanup] Section

OpenMECP accumulates many intermediate QM input/output files during a long run. The cleanup system automatically removes unnecessary files to prevent disk saturation.

Cleanup strategy: OpenMECP keeps the latest .engrad files and all .log / input files; it removes older intermediate files. Add custom extensions to preserve_extensions to protect specific file types.

Configuration Display at Startup

When OpenMECP starts, it prints:

• Which configuration file was loaded (or "built-in defaults" if none found)
• All input parameters and their values
• All convergence thresholds
• All configuration file settings

This ensures you always know exactly what settings are being used for a given run.

API Reference

The full Rust API documentation for the omecp crate is generated from source code comments and is available at:

View API Reference →

Module Overview

Generating API Docs Locally

cargo doc --no-deps --all-features --open

This builds and opens the full API documentation in your browser.