Overview

This page explains how PIMMS represents a system and how a simulation works. By the end you should understand the lattice model, the Monte Carlo move set, the energy function, and how to set up and judge the convergence of a run. The exhaustive per-keyword reference lives in Keyword reference.

The lattice and polymers

PIMMS is a lattice model: space is discretised into a square (2D) or cubic (3D) grid of sites, set by the DIMENSIONS keyword (e.g. DIMENSIONS : 70 70 70 for a 70×70×70 cube). Every site is either empty (implicit solvent) or occupied by exactly one bead - PIMMS enforces hard-sphere exclusion, so two beads can never share a site.

A polymer chain is a connected string of beads. Consecutive beads along a chain must be lattice-adjacent in the Chebyshev sense: each coordinate differs by -1, 0 or +1 from the previous bead (so a bead has up to 26 valid successor sites in 3D, 8 in 2D, including diagonals). This flexible connectivity lets chains fold compactly on the lattice.

Chains are declared with the CHAIN keyword, which gives a count and a one-letter sequence and may be repeated to build a multi-component mixture:

CHAIN : 20 QQQQQQQQQQ      # 20 copies of a 10-bead poly-Q homopolymer
CHAIN : 5  EEEEKKKKEEEE     # 5 copies of a hetero-polymer
CHAIN : 100 A               # 100 single-bead "particles"

Every bead letter used in a chain must be defined in the parameter file (Energy: interactions, solvation and angles). A chain of length 1 is a single free particle; a chain where all beads are identical is a homopolymer; otherwise it is a heteropolymer. Sequences are upper-cased on read unless CASE_INSENSITIVE_CHAINS : False (which lets a and A be distinct bead types).

The simulation box is periodic by default (a chain that leaves one face re-enters the opposite face). Setting HARDWALL : True instead makes the box edges hard, reflecting walls - see Simulation setup: boundaries and box resizing.

The box need not be cubic/square: unequal axes (e.g. DIMENSIONS : 20 20 60) are fully supported in 2D and 3D, with periodic or hardwall boundaries. The only restriction is that cluster-rotation moves (MOVE_CLUSTER_ROTATE) cannot be combined with a non-cubic box under periodic boundaries, because a 90° rigid rotation is only an energy-preserving symmetry of a cube/square (or of any box under hardwall).

Monte Carlo and the move set

PIMMS samples configurations with Metropolis Monte Carlo (MC). Starting from the current configuration it repeatedly proposes a random change (a move), computes the resulting energy change ΔE, and accepts the move with probability

\[P_\text{accept} = \min\!\left(1,\; e^{-\Delta E / T}\right),\]

where T is the TEMPERATURE. Downhill moves (ΔE < 0) are always accepted; uphill moves are accepted with a temperature-dependent probability. Each move is constructed to satisfy detailed balance, which guarantees that - given enough steps - the simulation samples the correct Boltzmann distribution of configurations. Hard-sphere overlaps are rejected outright.

One outer-loop step (N_STEPS counts these) is not a single move: the core crankshaft move is a “megamove” that performs CRANKSHAFT_SUBSTEPS single-bead sub-moves in total (spread at random over all beads), and several other moves are likewise batched. The true number of accept/reject operations is therefore far larger than N_STEPS (it is reported in TOTAL_MOVES.dat).

Which moves are attempted, and how often, is set by the MOVE_* keywords - fractions that must sum to 1.0. The available moves (with their internal move codes) are:

Local / single-chain moves

  • Crankshaft (MOVE_CRANKSHAFT, 1) - the workhorse. Fast, local bead rotations applied throughout the system in optimised Cython; most of your move budget should usually go here. CRANKSHAFT_SUBSTEPS tunes how much work each crankshaft step does.

  • Chain translate / rotate / pivot (MOVE_CHAIN_TRANSLATE 2, MOVE_CHAIN_ROTATE 3, MOVE_CHAIN_PIVOT 4) - rigid-body translation or rotation of a whole chain, or a pivot of one half of a chain about a randomly chosen point.

  • Head pivot (MOVE_HEAD_PIVOT, 5) - pivots a single terminus; rarely useful (keep at 0).

  • Slither / reptation (MOVE_SLITHER, 6) - advances a chain forwards or backwards through the lattice “like a snake”. Efficiently relaxes chain conformations; SLITHER_SUBSTEPS sets how many per chain.

  • Pull (MOVE_PULL, 11) - cooperative reptation of a sub-segment: an interior bead is displaced and the following beads are “pulled” along to keep the chain connected, letting chains rearrange in dense systems where rigid moves clash. Requires chains of length ≥ 3.

  • Jump-and-relax (MOVE_JUMP_AND_RELAX, 13) - relax a chain, relocate it, and relax again.

Collective / many-chain moves

  • Cluster translate / rotate (MOVE_CLUSTER_TRANSLATE 7, MOVE_CLUSTER_ROTATE 8) - rigid-body moves of a whole connected cluster of chains. Relatively expensive; keep small (0.01-0.05).

  • Virtual-Move Monte Carlo, VMMC (MOVE_VMMC, 14) - recruits a cluster of chains by interaction-energy gradients and moves it rigidly, to escape the kinetic traps that single-chain moves hit in dense/condensed phases.

  • Temperature-switch MC, TSMMC (MOVE_CTSMMC 9, MOVE_MULTICHAIN_TSMMC 10, MOVE_SYSTEM_TSMMC 12) - take a chain, subset of chains, or the whole system on a temperature excursion to hop over energy barriers. See Advanced features.

The collective and enhanced-sampling moves (TSMMC, pull, jump-and-relax, VMMC) are powerful for assembly/condensate problems; of these only VMMC is still experimental and gated behind EXPERIMENTAL_FEATURES : True. A robust default move set for most problems is mostly crankshaft with a little translate/rotate/pivot and slither.

Energy: interactions, solvation and angles

The total potential energy is a sum of pairwise bead-bead interactions, bead-solvent (solvation) terms, and backbone-angle penalties. Pairwise interactions act over three nested length scales, defined by how far apart two beads sit on the lattice (Chebyshev distance):

  • Short range (SR) - beads in contact (Chebyshev distance 1).

  • Long range (LR) - Chebyshev distance 2.

  • Super-long range (SLR) - Chebyshev distance 3.

SR is always present; LR and SLR are optional and only act between bead types that declare them. This lets you model, e.g., a strong short-ranged “sticker” attraction plus a weak longer-ranged electrostatic-like tail.

All of this is specified in the parameter file (the PARAMETER_FILE keyword). All interaction and angle energies must be integers (floats are rejected with an error). It has a few kinds of line:

## pairwise interactions:  R1 R2  e_SR [e_LR [e_SLR]]
A  A   -8                 # A-A short-range contact energy
A  B   -3  -2             # A-B short-range AND long-range energies
B  B   -6  -3   3         # B-B short, long and super-long-range energies

## solvation (interaction with solvent, denoted 0):  R 0  e_solv
A  0   -2                 # required for EVERY bead type
B  0   -1

## backbone angle penalties:  ANGLE_PENALTY R  a1 a2 a3
ANGLE_PENALTY  A   30 10 0
ANGLE_PENALTY  B   50 20 0

## ...or temperature-normalised angle penalties (units of kT, k=1):
##   ANGLE_PENALTY_T_NORM R  a1 a2 a3
ANGLE_PENALTY_T_NORM  A   0.5 0.2 0    # multiplied by TEMPERATURE at parse time

Notes:

  • Negative energies are favourable (attractive); positive are repulsive.

  • The short-range interaction matrix must be complete and non-redundant: every pair of bead types you use (including each type with itself) needs exactly one short-range line. For types {A, B} that means A A, A B and B B - a missing or duplicated pair is an error.

  • A solvation line for every bead type is mandatory - the energy is measured relative to a fully solvated reference, so PIMMS needs to know each bead’s bead-solvent energy. Solvent is the special type 0; the solvent-solvent energy is fixed at 0. Long-range (LR/SLR) terms are for solute-solute pairs only - a solvent (0) entry in an LR/SLR line is an error. Unlike the short-range matrix, LR/SLR pairs need not be complete (any pair you omit defaults to 0).

  • Angle penalties bias the local backbone geometry (three values per residue for the distinct lattice bend angles). Use either ANGLE_PENALTY (absolute integer penalties) or ANGLE_PENALTY_T_NORM (penalties in units of \(k_BT\) with \(k_B=1\), scaled by TEMPERATURE when the file is read - handy for keeping the stiffness fixed relative to temperature). Set ANGLES_OFF : True to disable angles entirely (then no angle lines are needed).

  • Set NON_INTERACTING : True to zero all interaction energies and run a pure excluded-volume reference simulation, regardless of the parameter file.

Internally the energy decomposes into SR, LR, SLR and angle components (see Hamiltonian.evaluate_total_energy); the running total is what ENERGY.dat records, and ENERGY_CHECK periodically re-derives it from scratch as a correctness guard.

Simulation setup: boundaries and box resizing

Boundary conditions. With the default periodic boundaries a system behaves as a bulk phase. HARDWALL : True gives reflective walls instead - appropriate for a droplet in a finite container, or whenever you do not want chains wrapping across the box.

Box resizing for equilibration. Sometimes you want to condense a system at high effective concentration and then study it in a larger box. RESIZED_EQUILIBRATION runs the equilibration phase in a smaller box and then grows it (re-centring the chains) to the full DIMENSIONS for production; EQUILIBRATION_OFFSET places the small box within the large one. Resizing requires HARDWALL : True.

Centring. For single-chain runs, AUTOCENTER : True keeps the chain in the middle of the box every frame, so the saved trajectory needs no post-hoc alignment.

(Box resizing also appears when restarting from a previous configuration into a larger box - see Restart files.)

Running and converging a simulation

A simulation is driven by a keyfile: a plain-text file of KEYWORD : value lines (# starts a comment). A minimal but complete keyfile looks like this:

## --- system ---
DIMENSIONS      : 30 30 30
PARAMETER_FILE  : params.prm
CHAIN           : 50 AABBAABB     # 50 copies of an 8-bead heteropolymer
TEMPERATURE     : 60

## --- run length ---
N_STEPS         : 5000
EQUILIBRATION   : 1000

## --- moves (must sum to 1.0) ---
MOVE_CRANKSHAFT     : 0.8
CRANKSHAFT_SUBSTEPS : 20000
MOVE_CHAIN_TRANSLATE: 0.1
MOVE_SLITHER        : 0.1
SLITHER_SUBSTEPS    : 200

## --- output / analysis ---
EN_FREQ         : 10
XTC_FREQ        : 100
ANA_CLUSTER     : 100

Run it with:

PIMMS -k KEYFILE.kf

PIMMS first runs EQUILIBRATION steps (no analysis output, and trajectory frames only if SAVE_EQ : True) and then the remaining production steps, writing the requested output files (see Output files).

Judging convergence. The first thing to check is ENERGY.dat: the potential energy should fall (or rise) and then plateau with stationary fluctuations - that plateau marks equilibrium. If the energy is still drifting at the end of equilibration, increase EQUILIBRATION (and probably N_STEPS). For assembly/condensate problems, also confirm that structural observables (e.g. cluster sizes in CLUSTERS.dat, radius of gyration in RG.dat) have stabilised, and that move acceptance ratios (ACCEPTANCE.dat divided by MOVE_FREQS.dat) are reasonable - extremely low acceptance for a move means it is doing little useful work. PERFORMANCE.dat reports throughput and an estimated time-to-completion. As a correctness safeguard you can set ENERGY_CHECK so PIMMS periodically re-computes the total energy from scratch and aborts if it has drifted.