Checking and exploring umbrella sampling runs

An umbrella sampling (US) run leaves dozens or hundreds of windows on disk. WHAM turns them into a free energy profile (PMF), but WHAM does not tell you whether the windows were good: whether neighbouring windows overlap, whether each window stayed near its target, whether the run is long enough. The Umbrella Sampling Browser answers those questions and lets you load the structures of any window, a range of windows or a path on a 2D PMF as objects in EasyHybrid.

This tutorial uses a small system that runs in minutes: the gas-phase SN2 reaction Cl− + CH3Br → ClCH3 + Br− with the AM1 Hamiltonian (6 atoms, charge −1). Three runs are used:

• SN2_US_1D: reaction coordinate d(C–Br) − d(C–Cl), 31 windows, k = 600 kJ mol-1 Å-2: a well-behaved 1D run.

• SN2_US_1D_stiff: the same with k = 4000: windows that do not overlap.

• SN2_US_2D: d(C–Br) and d(C–Cl), 15 × 15 windows, k = 600: a 2D run, its PMF and the reaction path.

The starting structure is in examples/tutorials/11_umbrella_browser_sn2/; section 1 shows how to produce the three runs (a few minutes with 12 processes).

1. Running the example

Each US window starts from one frame of a relaxed PES scan, so a US run in EasyHybrid always has two steps: a scan along the reaction coordinate, then the US from the scan frames.

• The system. File > Open examples/tutorials/11_umbrella_browser_sn2/sn2_cl_ch3br.xyz (the reactant complex: Cl− 3.2 Å from the carbon, on the side opposite to Br). Select all atoms, Selection > Set as QC Atoms, and in System > Setup QC Region choose the semiempirical method AM1, charge −1, multiplicity 1.

• The scan. Simulate > Reaction Coordinate Scans. Reaction Coordinate 1: Coordinate Type multiple distance, Atom 1 = Br, Atom 2 = C, Atom 3 = Cl (pick them in the 3D view and press Import from Picking Selection), so that RC = d(C–Br) − d(C–Cl). Initial Distance −1.4 Å, Step Size 0.1 Å, Number of Steps 31, Force Constant 4000. The scan takes a few seconds; the minimum of the reactant complex is at RC ≈ −0.9 Å and the top of the profile at RC ≈ +0.3 Å.

• The umbrella sampling. Simulate > Umbrella Sampling, with the scan trajectory as Input Trajectory and the same reaction coordinate. 1D run: Force Constant 600 (and 4000 for the "stiff" run), Langevin integrator at 300 K with collision frequency 25 ps-1, time step 0.001 ps, 2000 steps of equilibration and 10 000 of data collection (ptRes every step, ptGeo every 100), 12 CPUs. 2D run: RC1 the simple distance Br–C from 1.8 Å, RC2 the simple distance C–Cl from 1.7 Å, 15 steps of 0.1 Å each, starting from a 2D scan of the two distances over the same grid; 1000 + 4000 steps. With 12 processes the 1D runs take about 1.5 min each and the 2D run about 5 min.

• Then run Analysis > WHAM on each run (100 bins for 1D, 60 × 60 for 2D, 300 K): its *_pmf.log file goes to the wham folder of the run, where the browser finds it.

Each run is a folder named after the trajectory (here SN2_US_1D) with output.log (all the parameters of the run), optmization.ptGeo/ (the starting structure of each window), equilibration/ (the equilibration frames), data_collection/ (the reaction coordinate at every step, window*.ptRes, used by WHAM), data_collection_traj/ (the data collection frames) and wham/ (the WHAM output).

2. Opening a run

Load the system first (File > Open sn2_cl_ch3br.xyz): the browser loads frames as new objects of a system with the same atoms. Then open Analysis > Umbrella Sampling Browser... (or press Browse Windows... in the WHAM window, which opens the browser on the run chosen there). Press Browse... and choose the SN2_US_1D folder. The browser shows what it found: a 1D run with 31 windows at 300 K, 10 001 restraint values and 101 frames per window, and the definition of RC1 (Br(4)–C(0) / C(0)–Cl(5) (×−1), the weighted distance d(C–Br) − d(C–Cl), with the 0-based atom indices of pDynamo).

• System: the system that receives the loaded frames (it must have the same atoms). Coordinates: data collection, equilibration or the initial structure of each window. Stride: read every n-th frame. PMF: a PMF file of pDynamo's WHAM, or a WHAM computed by the browser. Energy: kJ/mol or kcal/mol.

The left side has the views (Windows, Diagnostics, Free energy), the list of selected windows and the buttons that load or export them. The right side shows the plot of the view, a table and a Details tab with the explanation of the numbers.

3. Are the windows good?

The window map. Windows > Window map shows one cell per window (a strip for 1D, a matrix for 2D), coloured by the quantity chosen in Color by:

• Frames: coordinates saved in the window. Warning sign: a window with fewer frames than the others (a failed window).

• |mean − target| / restraint width: how far the window's mean RC is from its target, in units of the width sqrt(kT/2k) the restraint alone would give. Above about 3, the window slid away from its target.

• Drift: change of the mean RC between the first and the second half; a large drift means the window was not equilibrated.

• RC spread: standard deviation of the RC in the window.

• Lowest overlap with a neighbour: see below; a warning below 0.03.

• PMF at the target: the PMF where the window was centred.

The table under the plot has every number per window. In SN2_US_1D all windows stay within 1.2 widths of their targets (Figure 1).

Overlap of neighbouring windows. WHAM can only connect two windows if they sample some common values of the reaction coordinate. Diagnostics > RC distributions and overlap shows the distribution of the RC in every window (Figure 2), and Overlap between neighbours the overlap coefficient of each pair of neighbours: the common area of their two normalized histograms, 1 for identical distributions and 0 for disjoint ones (Figure 3). With k = 600 kJ mol-1 Å-2 each window spreads over about 0.045 Å (one standard deviation) and the targets are 0.1 Å apart: the distributions touch, and the lowest overlap of a pair is 0.22.

Mean RC against the target. Diagnostics > Mean RC against the target plots, for each window, where it actually sampled (Figure 4). The mean is pulled a few hundredths of an ångström downhill, towards the reactant and product complexes: the restraint does not cancel the slope of the PMF. That is normal; WHAM corrects for it. The Drift view shows the change of the mean between the two halves of the run. Diagnostics > RC along the run shows the RC at every step of the selected windows, with their targets as dashed lines (select windows first, section 5; with none selected it shows the middle window; Figure 5).

A run whose windows do not overlap. Open SN2_US_1D_stiff, the same run with k = 4000 kJ mol-1 Å-2. Each window now spreads over only 0.018 Å, five times less than the distance between targets (Figures 6 and 7). All 30 pairs are below 0.03. WHAM still runs and still gives a PMF (here with a barrier close to the good one), but nothing in that PMF is determined by the data between the windows: the result cannot be trusted. The map coloured by Lowest overlap with a neighbour makes the problem visible at a glance (Figure 8).

The fix is a softer force constant or more windows. A rule of thumb: the width of a window, sqrt(kT/2k), should be at least about half the distance between targets. At 300 K, kT = 2.49 kJ/mol, so for windows 0.1 Å apart k should not exceed about 500–600 kJ mol-1 Å-2.

4. The free energy

Go back to SN2_US_1D. Free energy > PMF and windows plots the PMF chosen in the PMF box (here the file of pDynamo's WHAM, wham/output_pmf.log) with the window targets on it; clicking near a window selects it (Figure 9). The reactant complex is at RC ≈ −1.1 Å, the barrier at RC ≈ +0.3 Å, about 59 kJ/mol higher, and the product complex (ClCH3···Br−) at RC ≈ +1.4 Å, about 21 kJ/mol above the reactant complex.

Convergence. Free energy > Convergence computes the PMF again from the first half and from the second half of every window (and, in 1D, from each quarter), with a WHAM implemented in the browser, and compares them. Press Compute (a few seconds in 1D, about 20 s per PMF in 2D). The two halves differ by at most 3.1 kJ/mol, mostly near the barrier: the run is reasonably converged (Figure 10). A large difference, or a systematic shift between the halves, means the windows need longer data collection or longer equilibration. The browser's own WHAM is used to compare halves and blocks with the same bins and equations; the PMF to report is the one of Analysis > WHAM (pDynamo), which can also give bootstrap errors.

5. Loading structures

Selecting windows. Click windows in the map (or in the PMF) to select them; click again to unselect. The order of the clicks is the order of the windows in a trajectory. Clear / All select nothing or every window in order; Range (1D) every window between the first and the last selected; Row and Column (2D) the windows with the same RC1 (row) or RC2 (column) index as the last selected. For the figure below, windows 4 to 28 were selected: click window 4, then window 28, then Range.

• One object per window: each selected window as an object, with all its frames (source and stride above).

• One trajectory (windows in order): all the frames of the selected windows in one object.

• Representative frames (one per window): one frame per window, the one whose RC is closest to the window's target: a "movie" of the reaction.

• Export...: the representative frames or all the frames to DCD, XYZ or PDB.

Double-clicking a window in a table loads that window alone. Representative frames with windows 4–28 gives an object of 25 frames; the mini player plays it, and the C–Cl bond forms while the C–Br bond breaks and the methyl group inverts (Figure 11).

6. A 2D run

Open SN2_US_2D. The window map is now a 15 × 15 matrix, rows ordered by the RC1 target (d(C–Br)) and columns by the RC2 target (d(C–Cl)), in the same orientation as the 2D PMF (Figure 12). Two corners are dark: the windows with both distances short (Cl and Br both bonded to carbon) or both long (a free CH3 between two ions). Those are hundreds of kJ/mol above the reaction path and a force constant of 600 cannot hold them on target; they slide towards lower energy. That does not matter for the reaction, which runs along the other diagonal, but the map shows where the grid was not sampled as planned.

Coverage. Diagnostics > Coverage shows all the samples of all windows on the RC1 × RC2 plane (log of the count per bin), with the window targets as black dots (Figure 13). Grey bins have no sample. Holes between targets are regions WHAM cannot fill; here the sampled region covers the whole reaction valley. Figure 14 shows the 2D PMF from pDynamo's WHAM with the window targets.

The lowest-barrier path. Free energy > Lowest-barrier path finds, on the 2D PMF, the path between two points that keeps the highest free energy met as low as possible: the route through the lowest saddle (Figure 15). By default it goes from the global minimum to the lowest point far from it; Set start and Set end, then a click on the PMF, move the ends. The Profile along the path view gives the free energy along it: the saddle is at d(C–Br) = 2.35 Å and d(C–Cl) = 2.10 Å, 60.8 kJ/mol above the reactant complex, in agreement with the 1D run (59 kJ/mol; Figure 16).

The table lists the windows nearest to the path, in order (25 here). Use these windows selects them; the map shows them (Figure 17), and Representative frames turns them into a movie of the reaction along the 2D surface.

Convergence in 2D. In 2D the convergence view shows the difference between the PMFs of the two halves (Figure 18). Where the PMF is within 30 kJ/mol of the minimum, they differ by 0.95 kJ/mol on average; the larger differences are at the high-energy edges, which are sampled by few points.

7. Quick reference

• Overlap below 0.03 between neighbours: force constant too large for the window spacing; lower k, or add windows in the gap.

• Windows with deviation above 3: force constant too small for the slope of the PMF, or the system escaped; higher k there, or accept it in regions off the path.

• Large drift: equilibration too short; longer equilibration.

• Halves differ near the barrier: data collection too short; longer data collection.

• Grey holes in the 2D coverage: windows missing or not sampled; add windows, or restrict the analysis to the sampled region.

• A window with fewer frames: the window failed; check the output of that window and run it again.

Every plot can be saved with Export Image... (PNG or SVG) and every table with Save Table... (CSV). In the WHAM window, Browse Windows... opens the browser on the chosen run (Figure 19).