Force Matching

Overview

The force-matching workflow in MiMiCPy automates the parameterization of MM force field parameters by fitting to reference QM/MM data. It optimizes non-bonded interactions (charges via DRESP) and bonded interactions (bonds, angles, dihedrals) to reproduce QM forces and electrostatic properties.

Inputs

The mimicpy fm command expects:

Input Files

-top [.top]

GROMACS topology file.

-sele [.txt/.dat]

QM region selection file.

-fmdata [.json/.h5]

Force matching data file with QM/MM reference data.

-fi [.inp]

Force matching input file containing optimization parameters.

-trr [.trr]

GROMACS trajectory file for rerun force evaluation.

-ndx [.ndx]

GROMACS index file.

-coords [.gro]

Coordinate file matching the topology order.

-mdp [.mdp]

GROMACS MDP file for gmx grompp in the rerun step.

Other Options

-ff

Path to force field data directory (overrides GMXDATA/GMXLIB).

-gmx (gmx)

GROMACS executable name.

-dresp

Run only DRESP optimization (skip bonded optimization).

-opt_ff

Run only bonded optimization (skip DRESP).

-grid_search

Grid search over DRESP weights (wv, we, wh).

-compare_methods

Compare all optimization methods and choose the best one.

-n_processes

Number of processes for parallel computation (auto-detected if not set).

Force-Matching Input File

The force-matching input file lists the fitting parameters. Lines beginning with # or ! are ignored.

Trajectory sub-sampling

  • stride (mandatory) List of three integers [start, end, step]. It defines the starting snapshot, the final snapshot, and the sampling frequency of the reference data to be used for fitting.

Chemical equivalence definitions

  • eq_atoms (mandatory) Accepted values:

    • use_atomtypes to automatically determine equivalence from atom types.

    • eq_atoms_global followed by a list of tuples with equivalent GROMACS indices, e.g. [(1,2),(3,4)].

    • eq_atoms_local followed by a dictionary of molecule names and local equivalent atoms, e.g. MOL: [(1,2),(3,4)].

D-RESP fitting parameters

  • wv (mandatory) Float or list of floats. Weight for electrostatic potential fitting.

  • we (mandatory) Float or list of floats. Weight for electric field fitting.

  • wh (mandatory) Float or list of floats. Weight for Hirshfeld charge restraints.

  • wq (optional, default: 100000) Weight for total charge constraint.

  • reference_charges (optional) Use ff_charges to start from the existing force field charges, or provide an explicit list of floats that follows the MiMiCPy ordering of QM atoms. The QM atom indices can be printed with mimicpy qminfo.

  • qm_total_charge (mandatory) Total charge constraint for the QM region.

  • num_bonds_away (optional, default: 2) Number of bonds away from the boundary to include in charge redistribution.

  • charge_group_constraints (optional, default: True) Enforce charge group constraints.

  • fixed_charge_indices (optional, default: none) Atom indices with fixed charges.

  • weights_to_fix_charges (optional, default: 100000.0) Weight to fix charges to original values.

Solvent handling

  • skip_solvent_optimization (optional, default: True) Skip optimization of solvent molecules.

  • solvent_resnames (optional) Residue names to treat as solvent.

  • solvent_molecules (optional) Molecule names to treat as solvent.

Bonded parameter optimization controls

  • optimize_bond_length (optional, default: True) Optimize equilibrium bond lengths.

  • optimize_bond_force (optional, default: True) Optimize bond force constants.

  • optimize_angle_value (optional, default: True) Optimize equilibrium angle values.

  • optimize_angle_force (optional, default: True) Optimize angle force constants.

  • optimize_dihedral_force (optional, default: True) Optimize dihedral force constants.

Excluded interactions

  • exclude_bonds (optional) Atom pairs [(atom1, atom2), ...] to exclude from fitting.

  • exclude_angles (optional) Atom triplets [(atom1, atom2, atom3), ...] to exclude from fitting.

  • exclude_dihedrals (optional) Atom quartets [(atom1, atom2, atom3, atom4), ...] to exclude from fitting.

  • exclude_hydrogen_bonds (optional, default: False) Exclude all bonds involving hydrogen.

  • exclude_hydrogen_angles (optional, default: False) Exclude all angles involving hydrogen.

  • exclude_hydrogen_dihedrals (optional, default: False) Exclude all dihedrals involving hydrogen.

Regularization settings

  • regularization (optional, default: True) Enable L2 regularization.

  • regularization_alpha (optional, default: 0.1) L2 regularization strength.

Optimization method selection

  • optimization_method (optional, default: hierarchical) Accepted values: hierarchical or simultaneous.

Minimal example of an input file

stride: [50, 1103, 1]
wv: 1
we: 0.1
wh: 100
wq: 100000
eq_atoms: use_atomtypes
optimize_dihedral_force: False
reference_charges: ff_charges
optimization_method: hierarchical

Force-Matching Data Format

The data file can be JSON or HDF5. JSON is converted to HDF5 automatically on first use. Each configuration contains an atoms array with:

  • id (atom ID)

  • region (1 for QM, 2 for MM)

  • coordinate (3-vector)

  • electric_potential (scalar)

  • electric_field (3-vector)

  • hirshfeld_charge (scalar)

  • force (3-vector)

Use mimicpy json2h5 to convert large JSON files:

$ mimicpy json2h5 -i fm_data.json -o fm_data.h5

Inspect HDF5 contents with:

$ mimicpy h5info -i fm_data.h5

Command Reference

Force matching

$ mimicpy fm -top system.top -sele qm_selection.txt -fmdata fm_data.json \
             -fi fm_input.dat -trr trajectory.trr -ndx index.ndx \
             -coords system.gro -mdp run.mdp

Outputs

The force matching workflow writes the following files:

  • fm.log with detailed progress and diagnostics.

  • DRESP results: resp_*.top (when running -dresp only).

  • Non-bonded rerun files: non_bonded_*.itp, rerun_*.top, rerun_nonb.tpr, rerun_nonb.trr.

  • Bonded optimization results: optimization_results.txt, opt_*.top.

  • If -grid_search is used: dresp_grid_search_results.csv.