Hillclimber

https://github.com/zincware/hillclimber
Admin
Hillclimber is a Python framework for enhanced sampling with PLUMED that provides Pythonic...
Tokens:15,990
Snippets:120
Trust Score:7.4
License:LGPL-3.0
Update:10 months ago
Show doc for...
Tokens:
Raw
### Python: Common Real-World Examples

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Presents common use-case examples for HillClimber, including calculating distances between specific atoms, molecule COMs, and coordination numbers. Shows the corresponding PLUMED output.

```python
# Distance between two specific atoms
hc.DistanceCV(x1=sel1[0][0], x2=sel2[0][0])

# Distance between molecule COMs
hc.DistanceCV(x1=hc.VirtualAtom(sel1[0],"com"), x2=hc.VirtualAtom(sel2[0],"com"))

# Distance from one COM to all COMs
hc.DistanceCV(x1=va1[0], x2=va2)

# Coordination between all COMs
hc.CoordinationCV(x1=va, x2=va, r_0=5.0)

# COM of COMs (e.g., cluster center)
hc.VirtualAtom(hc.VirtualAtom(sel, "com"), "com")
```

--------------------------------

### Python: Common VirtualAtom Patterns

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Provides common patterns for using VirtualAtom, including calculating the COM of individual groups, nested COMs, and combining multiple VirtualAtom sets.

```python
# COM for each water
water_coms = hc.VirtualAtom(water_sel, "com")

# COM for ethanol[0]
ethanol_0_com = hc.VirtualAtom(ethanol_sel[0], "com")

# Combine both
all_coms = ethanol_0_com + water_coms
```

--------------------------------

### Install Hillclimber using uv

Source: https://github.com/zincware/hillclimber/blob/main/README.md

Installs the hillclimber package using the `uv` package manager. This is a straightforward command-line operation.

```bash
uv add hillclimber
```

--------------------------------

### Python: Indexing Atom Groups with Selectors

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Demonstrates how to access specific groups or atoms within groups selected by a Selector object. Supports slicing and advanced indexing to retrieve subsets of selected atoms.

```python
# First oxygen of each water
oxygens = water_sel[:][0]

# First two waters only
first_two = water_sel[0:2]
```

--------------------------------

### Python: CoordinationCV Calculation

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Shows how to use CoordinationCV to measure the degree of solvation or clustering between groups of atoms or virtual sites. Requires defining a cutoff radius (r_0).

```python
cn = hc.CoordinationCV(
    x1=hc.VirtualAtom(water_sel, "com"),
    x2=hc.VirtualAtom(water_sel, "com"),
    r_0=5.0,
    prefix="cn_water"
)
# → One CV: coordination among all water COMs
```

--------------------------------

### Python: DistanceCV Calculation

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Demonstrates the use of DistanceCV to calculate the distance between two groups or virtual sites. Supports single-to-single, one-to-many, and many-to-many distance calculations.

```python
# Distance between ethanol[0] COM and water[0] COM
dist = hc.DistanceCV(
    x1=hc.VirtualAtom(ethanol_sel[0], "com"),
    x2=hc.VirtualAtom(water_sel[0], "com"),
    prefix="d_com"
)

# Distance from one ethanol to each water
dist = hc.DistanceCV(
    x1=hc.VirtualAtom(ethanol_sel[0], "com"),
    x2=hc.VirtualAtom(water_sel, "com"),
)
# → Creates 3 CVs: d_0, d_1, d_2
```

--------------------------------

### Python: Create Virtual Atoms (COM/COG)

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Illustrates the creation of VirtualAtom objects, which represent virtual sites like the center of mass (COM) or center of geometry (COG) for selected groups. Each virtual site corresponds to a group in the input selector.

```python
water_coms = hc.VirtualAtom(water_sel, "com")  # 3 COMs
ethanol_coms = hc.VirtualAtom(ethanol_sel, "com")  # 2 COMs
```

--------------------------------

### Python: Combine Selectors

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Shows how to combine multiple Selector objects using the '+' operator. This creates a new selector that includes all atoms or groups matched by the individual selectors.

```python
all_mols = water_sel + ethanol_sel
```

--------------------------------

### PLUMED Equivalent for OPES Metadynamics

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This example shows the equivalent PLUMED input for setting up OPES Metadynamics. It defines collective variables (phi and psi) using the TORSION command and then applies OPES_METAD to these variables with specified PACE and BARRIER parameters.

```plumed
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
opes: OPES_METAD ARG=phi,psi PACE=500 BARRIER=40
```

--------------------------------

### Hierarchical Center of Mass (Python)

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Builds a hierarchical center of mass by first creating a virtual atom for a selection and then computing the center of mass of that virtual atom. This allows for multi-level structural analysis.

```python
VirtualAtom(VirtualAtom(sel, "com"), "com")
```

--------------------------------

### Install Hillclimber using pip

Source: https://github.com/zincware/hillclimber/blob/main/README.md

Installs the hillclimber package using the standard `pip` package manager. This command ensures the library is available in your Python environment.

```bash
pip install hillclimber
```

--------------------------------

### Combining Selectors (Python)

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Merges two existing selectors or virtual atom sets into a single, larger set. This is useful for building complex selections from simpler components.

```python
sel1 + sel2
```

--------------------------------

### Virtual Atom Creation (Python)

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Creates a virtual atom representing the center of mass (COM) or center of geometry (COG) for a given selection of atoms. This is typically used for COM/COG calculations and not for indexing.

```python
hc.VirtualAtom(water_sel, "com")
```

--------------------------------

### Python: Select Atoms with SMARTSSelector

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Uses SMARTSSelector to identify groups of atoms based on SMARTS patterns. Returns lists of atom indices for each matched group. Useful for selecting specific molecule types or functional groups.

```python
ethanol_sel = hc.SMARTSSelector("CCO")  # Matches 2 ethanols
water_sel = hc.SMARTSSelector("O")      # Matches 3 waters
```

--------------------------------

### Coordination Number Calculation (Python)

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Computes the coordination number between two sets of atoms or virtual atoms within a specified cutoff radius (r_0). This CV is useful for analyzing local chemical environments.

```python
hc.CoordinationCV(x1=va, x2=va, r_0=5.0)
```

--------------------------------

### Coordination Number Collective Variable with Python

Source: https://context7.com/zincware/hillclimber/llms.txt

This example demonstrates the creation of a Coordination Number Collective Variable (CV) using Hillclimber. It allows for calculating self-coordination within a group of molecules or solvation of one molecule type around another. Options include using virtual atoms or directly flattened atom selectors, with `to_plumed` generating the necessary PLUMED commands.

```python
import hillclimber as hc

water_sel = hc.SMARTSSelector("O")
ethanol_sel = hc.SMARTSSelector("CCO")

# Water self-coordination (clustering)
cn = hc.CoordinationNumberCV(
    x1=hc.VirtualAtom(water_sel, "com"),
    x2=hc.VirtualAtom(water_sel, "com"),
    r_0=5.0,  # Reference distance in Angstroms
    nn=6,
    mm=0,
    d_0=0.0,
    prefix="cn_water"
)

# Solvation: water molecules around ethanol
cn = hc.CoordinationNumberCV(
    x1=hc.VirtualAtom(ethanol_sel[0], "com"),
    x2=hc.VirtualAtom(water_sel, "com"),
    r_0=4.0,
    prefix="solvation"
)

# Using atom selectors directly (flattened)
cn = hc.CoordinationNumberCV(
    x1=water_sel,
    x2=ethanol_sel,
    r_0=3.5,
    flatten=True,
    prefix="cn"
)

labels, commands = cn.to_plumed(atoms)
```

--------------------------------

### Python: Flattening Option in CVs

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Explains the `flatten` option in CVs, which controls how multi-atom groups are represented in the output. `flatten=True` merges atoms into a single list, while `flatten=False` creates separate groups for each input group.

```python
dist = hc.DistanceCV(
    x1=water_sel[0],
    x2=ethanol_sel[0],
    flatten=False
)
```

--------------------------------

### Distance Calculation (Python)

Source: https://github.com/zincware/hillclimber/blob/main/Getting_started.md

Calculates the distance between two virtual atoms or selections. This is a common collective variable (CV) used to measure spatial relationships between molecular entities.

```python
hc.DistanceCV(x1=va1, x2=va2)
```

--------------------------------

### Complete Alanine Dipeptide Metadynamics Workflow

Source: https://context7.com/zincware/hillclimber/llms.txt

Demonstrates a full molecular dynamics simulation workflow using hillclimber, ipsuite, and zntrack. It covers system preparation (SMILES to atoms, solvation, optimization), Metadynamics configuration (CVs, bias, actions), simulation execution, and analysis of free energy surfaces.

```python
import hillclimber as hc
import ipsuite as ips
import zntrack

TEMPERATURE = 300.0  # Kelvin
TIMESTEP = 0.5       # fs
N_STEPS = 2_000_000

project = zntrack.Project()

with project:
    # Create system
    alanine = ips.Smiles2Atoms(smiles="CC(=O)NC(C)C(=O)NC")
    water = ips.Smiles2Conformers(smiles="O", numConfs=100)

    box = ips.Packmol(
        data=[alanine.frames, water.frames],
        count=[1, 64],
        density=1000
    )

    geoopt = ips.ASEGeoOpt(
        data=box.frames,
        model=force_field, # Assuming force_field is defined elsewhere
        optimizer="LBFGS",
        run_kwargs={"fmax": 0.05}
    )

    # Define CVs
    phi_selector = hc.SMARTSSelector("[C:1](=O)[N:2][C:3][C:4](=O)")
    psi_selector = hc.SMARTSSelector("C(=O)[N:1][C:2][C:3](=O)[N:4]")

    phi_cv = hc.TorsionCV(atoms=phi_selector, prefix="phi")
    psi_cv = hc.TorsionCV(atoms=psi_selector, prefix="psi")

    # Configure metadynamics
    phi_bias = hc.MetadBias(cv=phi_cv, sigma=0.35, grid_min="-pi", grid_max="pi")
    psi_bias = hc.MetadBias(cv=psi_cv, sigma=0.35, grid_min="-pi", grid_max="pi")

    metad_config = hc.MetaDynamicsConfig(
        height=0.1,
        pace=500,
        biasfactor=6,
        temp=TEMPERATURE
    )

    print_action = hc.PrintAction(cvs=[phi_cv, psi_cv], stride=100, file="COLVAR")

    metad_model = hc.MetaDynamicsModel(
        config=metad_config,
        bias_cvs=[phi_bias, psi_bias],
        actions=[print_action],
        data=geoopt.frames,
        model=force_field, # Assuming force_field is defined elsewhere
        timestep=TIMESTEP
    )

    # Run simulation
    md = ips.ASEMD(
        data=geoopt.frames,
        model=metad_model,
        thermostat=ips.LangevinThermostat(
            temperature=TEMPERATURE,
            friction=0.01,
            time_step=TIMESTEP
        ),
        steps=N_STEPS,
        sampling_rate=1000
    )

project.repro()

# Analysis
import metadynminer
fig, axes = hc.plot_cv_time_series("nodes/ASEMD/model/-1/COLVAR")
hillsfile = metadynminer.Hills("nodes/ASEMD/model/-1/HILLS")
fes = metadynminer.Fes(hillsfile)
fes.plot()
```

--------------------------------

### Configure Histogram and Reweighting Actions

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This Python code outlines the configuration for analysis actions. It shows how to set up a HistogramAction by specifying collective variables, grid boundaries, bin numbers, kernel type, stride, and output file. It also demonstrates initializing a ReweightBiasAction with a temperature parameter for reweighting biased simulation data.

```python
histogram = hc.HistogramAction(
    cvs=[cv1, cv2],
    grid_min=[0.0, -3.14],
    grid_max=[5.0, 3.14],
    grid_bins=[100, 100],
    kernel="discrete",
    stride=1000,
    file="histogram.dat"
)

reweight = hc.ReweightBiasAction(temp=300.0)
```

--------------------------------

### Running Tests with pytest

Source: https://github.com/zincware/hillclimber/blob/main/AGENTS.md

To execute the unit tests for the Hillclimber package, use the `uv run pytest tests/` command. Ensure that all new functionality includes corresponding unit tests that verify the complete output, not just partial matches.

```bash
uv run pytest tests/
```

--------------------------------

### Configure and Model OPES (On-the-fly Probability Enhanced Sampling)

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This snippet demonstrates how to configure and instantiate an OPES model in Python. OPES is a modern alternative to traditional metadynamics offering better convergence and automatic variance adaptation. It requires defining an OPESConfig object with parameters like barrier, pace, and temperature, and an OPESModel object that links the configuration to collective variables, data, and a specific model type.

```python
opes_config = hc.OPESConfig(
    barrier=40.0,  # kJ/mol
    pace=500,
    temp=300.0
)

opes_model = hc.OPESModel(
    config=opes_config,
    cvs=[cv1, cv2],
    data=data.frames,
    model=ips.MACEMPModel(),
    timestep=0.5
)
```

--------------------------------

### Atom Selection with SMARTS Patterns using Python

Source: https://context7.com/zincware/hillclimber/llms.txt

This snippet demonstrates how to select atoms in a molecular system using SMARTS patterns with the Hillclimber framework. It supports various selection criteria, including atom mapping for specific interactions, and allows for indexing and combining selectors. The `select` method returns groups of atom indices.

```python
import hillclimber as hc
from ase import Atoms

# Create a molecular system
atoms = Atoms(...)  # Your molecular structure

# Select all water molecules
water_sel = hc.SMARTSSelector(pattern="O")

# Select ethanol molecules ethanol_sel = hc.SMARTSSelector(pattern="CCO")

# Select carboxylic acid groups with hydrogens
acid_sel = hc.SMARTSSelector(pattern="C(=O)O", hydrogens="include")

# Select specific atoms using atom mapping for torsion angles
phi_sel = hc.SMARTSSelector(pattern="[C:1](=O)[N:2][C:3][C:4](=O)")

# Get groups of atom indices
groups = water_sel.select(atoms)  # [[0, 1, 2], [3, 4, 5], ...]

# Index selectors: first two waters
first_two = water_sel[0:2]

# Combine selectors
all_molecules = water_sel + ethanol_sel
```

--------------------------------

### Configure OPES Model with ZnTrack

Source: https://context7.com/zincware/hillclimber/llms.txt

Sets up an OPES (One-Particle Effective sampling) model configuration, including barrier height, pace, temperature, and explore mode. It also defines actions like printing CVs to a file and integrates with the ZnTrack workflow for reproducibility.

```python
import zntrack
import hillclimber as hc
import ipsuite as ips

# Assuming phi_cv, psi_cv, phi_bias, psi_bias, atoms_list, and force_field_model are defined

opes_config = hc.OPESConfig(
    barrier=40.0,       # eV (highest barrier to overcome)
    pace=500,           # Kernel deposition frequency
    temp=300.0,         # Kelvin
    explore_mode=False, # False=OPES_METAD, True=OPES_METAD_EXPLORE
    compression_threshold=1.0,
    file="KERNELS"
)

print_action = hc.PrintAction(cvs=[phi_cv, psi_cv], stride=100, file="COLVAR")

opes_model = hc.OPESModel(
    config=opes_config,
    bias_cvs=[phi_bias, psi_bias],
    actions=[print_action],
    data=atoms_list,
    data_idx=-1,
    model=force_field_model,
    timestep=0.5
)

with zntrack.Project():
    md = ips.ASEMD(
        model=opes_model,
        data=atoms_list,
 E        steps=1_000_000
    )

project.repro()
```

--------------------------------

### Import Libraries for Hillclimber Project

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Imports necessary libraries for the hillclimber project, including shutil for file operations, pathlib for path manipulation, ASE for atomic simulations, matplotlib for plotting, numpy for numerical operations, and specific modules from ASE and the hillclimber package.

```python
import shutil
from pathlib import Path

import ase
import ase.units
import matplotlib.pyplot as plt
import numpy as np
from ase.calculators.lj import LennardJones
from ase.md.langevin import Langevin
from tqdm import tqdm

import hillclimber as hc
```

--------------------------------

### Configure Metadynamics and Wall Biases

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Sets up metadynamics and upper wall biases using the hillclimber package. It defines collective variables (distance between atoms), metadynamics parameters (height, pace, temperature, bias factor), and wall parameters (position and stiffness).

```python
atom1 = hc.IndexSelector(indices=[[0]])
atom2 = hc.IndexSelector(indices=[[1]])

distance_cv = hc.DistanceCV(x1=atom1, x2=atom2, prefix="dist", flatten=True)

metad_bias = hc.MetadBias(
    cv=distance_cv,
    sigma=0.2,  # Width in Angstrom
    grid_min=0.0,
    grid_max=5.0,
)

wall = hc.UpperWallBias(
    cv=distance_cv,
    at=4,
    kappa=100.0,
)

metad_config = hc.MetaDynamicsConfig(
    height=0.005,  # Small height in eV
    pace=10,  # Deposit Gaussian every 10 steps
    temp=120.0,  # Temperature in Kelvin (Argon melts at ~84 K)
    biasfactor=10.0,  # Well-tempered metadynamics
    file="HILLS",
    flush=10,
)

class LJModel:
    """Minimal model for testing that wraps LJ calculator."""

    def get_calculator(self, directory=None, **kwargs):
        return lj_calc


metad_model = hc.MetaDynamicsModel(
    config=metad_config,
    data=[system.copy()],
    data_idx=0,
    bias_cvs=[metad_bias],
    actions=[hc.PrintAction(cvs=[distance_cv], stride=5), wall],
    timestep=2.0,  # 2 fs timestep
    model=LJModel(),
)

```

--------------------------------

### PLUMED 2 Syntax for Restraints and Walls

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This snippet shows the equivalent PLUMED 2 syntax for defining harmonic restraints and wall potentials. It mirrors the functionality provided by the hillclimber library, specifying the collective variable (ARG), force constant (KAPPA), target value (AT), and exponent (EXP) for walls.

```plumed
# Restraint
RESTRAINT ARG=cv KAPPA=200.0 AT=2.5

# Walls
UPPER_WALLS ARG=cv AT=3.0 KAPPA=100.0 EXP=2
LOWER_WALLS ARG=cv AT=1.0 KAPPA=100.0 EXP=2
```

--------------------------------

### Run OPES Enhanced Sampling Simulation

Source: https://context7.com/zincware/hillclimber/llms.txt

Sets up an On-the-fly Probability Enhanced Sampling (OPES) simulation using the OPESBias. This involves defining collective variables (e.g., TorsionCV) and configuring OPES biases with adaptive sigma. Requires hillclimber library.

```python
import hillclimber as hc

# Define collective variables
phi_cv = hc.TorsionCV(atoms=phi_selector, prefix="phi")
psi_cv = hc.TorsionCV(atoms=psi_selector, prefix="psi")

# Configure OPES biases with adaptive sigma
phi_bias = hc.OPESBias(cv=phi_cv, sigma="ADAPTIVE")
psi_bias = hc.OPESBias(cv=psi_cv, sigma="ADAPTIVE")

```

--------------------------------

### Run Molecular Dynamics Simulation with Langevin Integrator

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Executes a molecular dynamics simulation using the Langevin integrator. It sets the timestep, temperature, friction coefficient, and trajectory file. The simulation runs for a specified number of steps, depositing metadynamics Gaussians.

```python
dyn = Langevin(
    atoms=system,
    timestep=2.0 * ase.units.fs,
    temperature_K=120.0,
    friction=0.01,  # 1/fs
    trajectory=str(work_dir / "md.traj"),
)

n_steps = 100_000
for _ in tqdm(dyn.irun(n_steps), total=n_steps):
    pass

```

--------------------------------

### Run Metadynamics Simulation with Multiple Collective Variables

Source: https://context7.com/zincware/hillclimber/llms.txt

Configures and runs a well-tempered metadynamics simulation using multiple collective variables (CVs). It involves defining CVs (e.g., TorsionCV), configuring biases (MetadBias) with deposition parameters and grid settings, setting up the metadynamics configuration (height, pace, biasfactor, temp), and defining actions like printing CV values (PrintAction). The simulation is executed using ASE MD within a zntrack project context. Requires hillclimber, ipsuite, and zntrack libraries.

```python
import hillclimber as hc
import ipsuite as ips
import zntrack

# Define collective variables
phi_selector = hc.SMARTSSelector(pattern="[C:1](=O)[N:2][C:3][C:4](=O)")
psi_selector = hc.SMARTSSelector(pattern="C(=O)[N:1][C:2][C:3](=O)[N:4]")

phi_cv = hc.TorsionCV(atoms=phi_selector, prefix="phi")
psi_cv = hc.TorsionCV(atoms=psi_selector, prefix="psi")

# Configure biases for each CV
phi_bias = hc.MetadBias(
    cv=phi_cv,
    sigma=0.35,      # radians
    grid_min="-pi",
    grid_max="pi",
    grid_bin=100
)

psi_bias = hc.MetadBias(
    cv=psi_cv,
    sigma=0.35,
    grid_min="-pi",
    grid_max="pi",
    grid_bin=100
)

# Metadynamics configuration
metad_config = hc.MetaDynamicsConfig(
    height=0.1,        # eV (Gaussian height)
    pace=500,          # Deposition frequency (MD steps)
    biasfactor=6,      # Well-tempered factor
    temp=300.0,        # Kelvin
    file="HILLS",
    adaptive="GEOM"    # Adaptive sigma (optional)
)

# Print action to save CV values
print_action = hc.PrintAction(
    cvs=[phi_cv, psi_cv],
    stride=100,
    file="COLVAR"
)

# Create metadynamics model
metad_model = hc.MetaDynamicsModel(
    config=metad_config,
    bias_cvs=[phi_bias, psi_bias],
    actions=[print_action],
    data=atoms_list,
    data_idx=-1,
    model=force_field_model,
    timestep=0.5  # fs
)

# Run MD with metadynamics
with zntrack.Project():
    md = ips.ASEMD(
        model=metad_model,
        data=atoms_list,
        steps=1_000_000,
        sampling_rate=1000
    )

project.repro()
```

--------------------------------

### Generate PLUMED Calculator and Assign to System

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Generates a PLUMED-wrapped calculator from the metadynamics model and assigns it to the Argon system. This prepares the system for simulation with the defined biases and calculator.

```python
work_dir = Path("argon_metad")
work_dir.mkdir(exist_ok=True)

calc = metad_model.get_calculator(directory=work_dir)
system.calc = calc

```

--------------------------------

### Analyze Free Energy Surfaces with Hillclimber

Source: https://context7.com/zincware/hillclimber/llms.txt

Provides Python functions to reconstruct and visualize free energy surfaces from simulation data. It supports both 1D and 2D surfaces, reading simulation output, and plotting time series of collective variables.

```python
import hillclimber as hc

# Reconstruct 1D free energy surface
hc.sum_hills(
    hills_file="HILLS",
    bin=1000,
    outfile="fes.dat"
)

# 2D free energy surface (e.g., Ramachandran plot)
hc.sum_hills(
    hills_file="HILLS",
    bin=[250, 250],
    min_bounds=[-3.14, -3.14],
    max_bounds=[3.14, 3.14],
    idw=["phi", "psi"],
    outfile="ramachandran.dat"
)

# Read COLVAR file
data = hc.read_colvar("COLVAR")
print(data.keys())  # dict_keys(['time', 'phi', 'psi'])
time = data['time']
phi_values = data['phi']

# Plot CV time series with KDE
fig, (ax, ax_kde) = hc.plot_cv_time_series(
    colvar_file="COLVAR",
    cv_names=["phi", "psi"],
    figsize=(10, 6),
    colors=["blue", "red"],
    alpha=0.7
)

# Auto-detect CVs (excludes time, sigma_*, height, biasf)
fig, axes = hc.plot_cv_time_series("COLVAR")
```

--------------------------------

### Define Argon Dimer System with ASE

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Defines a simple Argon dimer system using ASE (Atomic Simulation Environment). It initializes two Argon atoms with specified positions, forming the basis for simulation.

```python
system = ase.Atoms(
    "Ar2",
    positions=[[0.0, 0.0, 0.0], [3.8, 0.0, 0.0]],
)

```

--------------------------------

### Set up Lennard-Jones Calculator for Argon

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Configures a Lennard-Jones (LJ) potential calculator for Argon atoms. It specifies the epsilon (well depth), sigma (collision diameter), and cutoff radius, along with enabling smoothing for the potential.

```python
epsilon_ar = 0.0104  # eV
sigma_ar = 3.4  # Angstrom
cutoff = 10.0  # Angstrom

lj_calc = LennardJones(sigma=sigma_ar, epsilon=epsilon_ar, rc=cutoff, smooth=True)

```

--------------------------------

### Creating Virtual Atoms (COM/COG) with Python

Source: https://context7.com/zincware/hillclimber/llms.txt

This code shows how to create virtual atoms representing the center of mass (COM) or center of geometry (COG) for groups of atoms using Hillclimber. It supports creating COMs for individual molecules, multiple molecules, or nested COMs. The `to_plumed` method generates corresponding PLUMED commands.

```python
import hillclimber as hc

# Create COM for each water molecule
water_sel = hc.SMARTSSelector("O")
water_coms = hc.VirtualAtom(water_sel, reduction="com")

# Single COM from first ethanol
ethanol_sel = hc.SMARTSSelector("CCO")
ethanol_0_com = hc.VirtualAtom(ethanol_sel[0], reduction="com")

# Center of geometry instead of mass
cog = hc.VirtualAtom(water_sel[0], reduction="cog")

# Flatten all groups into one virtual site
all_waters_com = hc.VirtualAtom(water_sel, reduction="flatten")

# Nested: COM of COMs (e.g., cluster center)
water_coms = hc.VirtualAtom(water_sel, "com")
cluster_center = hc.VirtualAtom(water_coms, reduction="com")

# Generate PLUMED commands
labels, commands = water_coms.to_plumed(atoms)
# labels: ['vsite_123_0', 'vsite_123_1', 'vsite_123_2']
# commands: ['vsite_123_0: COM ATOMS=1,2,3', ...]
```

--------------------------------

### Clean Up Simulation Directory

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Removes the working directory created for the simulation, including all generated files, to clean up the project space.

```python
shutil.rmtree(work_dir)

```

--------------------------------

### Load and Plot Free Energy Surface

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Loads the computed free energy surface data from a file and plots it using matplotlib. The plot shows the free energy as a function of the collective variable (distance).

```python
data = np.loadtxt("argon_metad/fes.dat", comments="#")

cv = data[:, 0]  # Collective variable (distance)
fes = data[:, 1]  # Free energy

_, (ax1) = plt.subplots(1, 1)

ax1.plot(cv, fes, "b-", linewidth=2)
ax1.set_ylabel("Free Energy (kJ/mol)", fontsize=12)
ax1.grid(True, alpha=0.3)

```

--------------------------------

### Implement Restraints and Walls using hillclimber Python API

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This snippet demonstrates how to define and use harmonic restraints, upper walls, and lower walls within the hillclimber Python library. These are crucial for enhanced sampling methods like umbrella sampling and steered MD. The API allows defining collective variables (CVs) using selectors and then configuring bias potentials with parameters like kappa (force constant) and at (target value). The defined biases are then passed to the MetaDynamicsModel.

```python
import hillclimber as hc

# Define a distance CV
distance_cv = hc.DistanceCV(
    x1=hc.SMARTSSelector(pattern="[OH]"),
    x2=hc.SMARTSSelector(pattern="[O]"),
    prefix="d"
)

# Harmonic restraint
restraint = hc.RestraintBias(
    cv=distance_cv,
    kappa=200.0,  # force constant (kJ/mol)
    at=2.5        # target value
)

# Upper wall
upper_wall = hc.UpperWallBias(
    cv=distance_cv,
    at=3.0,
    kappa=100.0,
    exp=2  # exponent (2=harmonic, higher=steeper)
)

# Lower wall
lower_wall = hc.LowerWallBias(
    cv=distance_cv,
    at=1.0,
    kappa=100.0,
    exp=2
)

# Use in MetaDynamicsModel via actions parameter
# Assuming config, metad_bias, and data are defined elsewhere
# model = hc.MetaDynamicsModel(
#     config=config,
#     bias_cvs=[metad_bias],
#     actions=[restraint, upper_wall, lower_wall],  # Add biases here!
#     data=data.frames,
#     model=ips.MACEMPModel()
# )
```

--------------------------------

### PLUMED 2 Syntax for Combined and Custom CVs

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This snippet presents the PLUMED 2 syntax for defining combined (linear combination) and custom (mathematical function) collective variables. It shows how to specify the input CVs (ARG), coefficients (COEFFICIENTS), functions (FUNC), and periodicity (PERIODIC).

```plumed
# Combine
diff: COMBINE ARG=d1,d2 COEFFICIENTS=1,-1 PERIODIC=NO

# Custom
custom: CUSTOM ARG=d1,d2 FUNC=x^2+y^2 PERIODIC=NO
```

--------------------------------

### PLUMED 2 Syntax for RMSD Collective Variables

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This snippet shows the PLUMED 2 syntax for defining RMSD-related collective variables. It corresponds to the proposed hillclimber API, demonstrating how to specify reference structures, types of RMSD calculation, and parameters like cutoffs and residue selections.

```plumed
rmsd: RMSD REFERENCE=native.pdb TYPE=OPTIMAL
drmsd: DRMSD REFERENCE=native.pdb LOWER_CUTOFF=0.1 UPPER_CUTOFF=0.8
alpha: ALPHARMSD RESIDUES=all
```

--------------------------------

### Expected PLUMED Output Verification

Source: https://github.com/zincware/hillclimber/blob/main/AGENTS.md

Unit tests should validate the full expected output from PLUMED commands. This involves comparing the generated PLUMED string against a list of exact command strings, rather than checking for the presence of substrings.

```python
expected = [
        "d12_g1_0_com: COM ATOMS=1,2,3,4,5,6,7,8,9",
        "d12_g2_0_com: COM ATOMS=19,20,21",
        "d12: DISTANCE ATOMS=d12_g1_0_com,d12_g2_0_com",
    ]

assert plumed_str == expected
```

--------------------------------

### Plot CV Time Series and Sum Hills for Free Energy Surface

Source: https://github.com/zincware/hillclimber/blob/main/examples/LJ.ipynb

Generates plots of the collective variable (CV) time series and sums the metadynamics hills to compute the free energy surface. The results are saved to specified output files.

```python
_ = hc.plot_cv_time_series("argon_metad/COLVAR")
hc.sum_hills(
    hills_file="argon_metad/HILLS",
    bin=500,
    outfile="argon_metad/fes.dat",
)

```

--------------------------------

### Define Torsion Collective Variable using SMARTSSelector

Source: https://context7.com/zincware/hillclimber/llms.txt

Defines a Torsion Collective Variable (CV) by selecting atoms based on SMARTS patterns. The 'all' strategy processes all matching groups. Requires a SMARTSSelector and a TorsionCV object.

```python
psi_selector = hc.SMARTSSelector(pattern="C(=O)[N:1][C:2][C:3](=O)[N:4]")

psi_cv = hc.TorsionCV(
    atoms=psi_selector,
    prefix="psi"
)

multi_torsion = hc.TorsionCV(
    atoms=phi_selector,
    prefix="torsion",
    strategy="all"  # Creates N CVs for N groups
)

labels, commands = phi_cv.to_plumed(atoms)
```

--------------------------------

### Define Angle Collective Variable using SMARTSSelector and VirtualAtoms

Source: https://context7.com/zincware/hillclimber/llms.txt

Measures angles between three atoms or virtual sites, using SMARTSSelectors to define atom positions. It supports direct atom selection or virtual atoms created from selectors. The 'first' strategy selects the first match, while 'flatten=True' combines all atoms into a single CV. Requires SMARTSSelector, VirtualAtom, and AngleCV objects.

```python
import hillclimber as hc

sel1 = hc.SMARTSSelector("C")
sel2 = hc.SMARTSSelector("O")
sel3 = hc.SMARTSSelector("N")

# Angle with x2 as vertex
angle = hc.AngleCV(
    x1=hc.VirtualAtom(sel1[0], "com"),
    x2=hc.VirtualAtom(sel2[0], "com"),  # vertex
    x3=hc.VirtualAtom(sel3[0], "com"),
    prefix="angle",
    strategy="first"
)

# Using atom selectors directly
angle = hc.AngleCV(
    x1=sel1[0],
    x2=sel2[0],  # vertex
    x3=sel3[0],
    prefix="angle",
    flatten=True
)

labels, commands = angle.to_plumed(atoms)
```

--------------------------------

### Define Radius of Gyration Collective Variable

Source: https://context7.com/zincware/hillclimber/llms.txt

Measures molecular size and compactness using SMARTSSelectors. Supports calculating Radius of Gyration ('RADIUS') or Asphericity ('ASPHERICITY'). 'flatten=True' combines all atoms into one CV, while 'strategy="all"' processes multiple groups separately. Requires SMARTSSelector and RadiusOfGyrationCV objects.

```python
import hillclimber as hc

# Select protein backbone
protein_sel = hc.SMARTSSelector("C(=O)N")

# Radius of gyration
rg = hc.RadiusOfGyrationCV(
    atoms=protein_sel,
    prefix="rg",
    flatten=True,  # Combine all atoms into one
    type="RADIUS"
)

# Asphericity (shape parameter)
asph = hc.RadiusOfGyrationCV(
    atoms=protein_sel,
    prefix="asph",
    flatten=True,
    type="ASPHERICITY"
)

# Process multiple groups separately
rg_multi = hc.RadiusOfGyrationCV(
    atoms=protein_sel,
    prefix="rg",
    flatten=False,
    strategy="all"
)

labels, commands = rg.to_plumed(atoms)
```

--------------------------------

### Define Combined and Custom CVs using hillclimber Python API

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This snippet illustrates the proposed Python API for creating combined and custom Collective Variables (CVs) in hillclimber. The `CombineCV` allows for linear combinations of existing CVs with specified coefficients and powers, while `CustomCV` enables defining CVs based on arbitrary mathematical functions of other CVs. These are essential for constructing complex reaction coordinates.

```python
# Assuming cv1, cv2 are defined Collective Variables
# Linear combination
# combined_cv = hc.CombineCV(
#     cvs=[cv1, cv2],
#     coefficients=[1.0, -1.0],
#     powers=[1, 1],  # optional
#     periodic=False,
#     prefix="combined"
# )

# Custom mathematical function
# custom_cv = hc.CustomCV(
#     cvs=[cv1, cv2],
#     function="x**2 + y**2",  # or lambda x,y: x**2 + y**2
#     periodic=False,
#     prefix="custom"
# )
```

--------------------------------

### Define Path Collective Variable for Reaction Path Sampling

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This Python code defines a Path Collective Variable (PathCV), which is crucial for transition path sampling and committor analysis. It requires a reference PDB file containing waypoints and allows specifying the path type (optimal or euclidean) and a lambda parameter. The CV outputs progress along the path (s) and deviation from the path (z).

```python
path_cv = hc.PathCV(
    reference="path.pdb",  # multi-frame PDB with waypoints
    lambda_param=15100.0,
    path_type="optimal",  # or "euclidean"
    prefix="path"
)
```

--------------------------------

### Define and Use Virtual Atoms in CVs

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This Python snippet illustrates the definition and usage of virtual atoms. Virtual atoms, such as center of mass (COM) or center of geometry (COG), can be explicitly defined using the VirtualAtom class. These defined virtual atoms can then be used as coordinates in other collective variable definitions, enabling hierarchical CV construction and reuse.

```python
center1 = hc.VirtualAtom(
    atoms=selector1,
    type="com",  # or "cog"
    label="c1"
)

center2 = hc.VirtualAtom(
    atoms=selector2,
    type="com",
    label="c2"
)

distance_cv = hc.DistanceCV(
    x1=center1,
    x2=center2,
    prefix="d_centers"
)
```

--------------------------------

### PLUMED Equivalent for Virtual Atoms and Distance CV

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This PLUMED input demonstrates defining virtual atoms for centers of mass (COM) and then using them in a DISTANCE collective variable. It shows how PLUMED supports defining intermediate sites (like c1 and c2) and referencing them in subsequent definitions.

```plumed
c1: COM ATOMS=1-100
c2: COM ATOMS=101-200
d: DISTANCE ATOMS=c1,c2
```

--------------------------------

### Torsion (Dihedral) Angle Collective Variable using Python

Source: https://context7.com/zincware/hillclimber/llms.txt

This snippet shows how to define a Torsion (Dihedral) Angle Collective Variable (CV) with Hillclimber. It is particularly useful for measuring backbone or sidechain torsion angles by specifying the atoms involved using selectors. The `strategy` parameter allows control over which torsion angle is used if multiple are found.

```python
import hillclimber as hc

# Phi angle in alanine dipeptide (requires 4 atoms)
phi_selector = hc.SMARTSSelector(pattern="[C:1](=O)[N:2][C:3][C:4](=O)")

phi_cv = hc.TorsionCV(
    atoms=phi_selector,
    prefix="phi",
    strategy="first"  # Use first matching group
)
```

--------------------------------

### Apply Upper and Lower Wall Biases to Collective Variable

Source: https://context7.com/zincware/hillclimber/llms.txt

Prevents a collective variable (CV) from exceeding upper or lower bounds using UpperWallBias and LowerWallBias. These biases require a CV, a boundary value (at), a force constant (kappa), and an exponent (exp) to define the potential shape. The `to_plumed` method generates the PLUMED commands. Requires DistanceCV, UpperWallBias, and LowerWallBias.

```python
import hillclimber as hc

distance_cv = hc.DistanceCV(...)

# Prevent distance from exceeding 3.0 Å
upper_wall = hc.UpperWallBias(
    cv=distance_cv,
    at=3.0,      # Å
    kappa=100.0,  # eV
    exp=2,       # Harmonic (exp=2), steeper for higher values
    label="d_uwall"
)

# Prevent distance from going below 1.0 Å
lower_wall = hc.LowerWallBias(
    cv=distance_cv,
    at=1.0,
    kappa=100.0,
    exp=2
)

# Generate PLUMED commands
commands = upper_wall.to_plumed(atoms)
```

--------------------------------

### PLUMED Equivalent for Path Collective Variable

Source: https://github.com/zincware/hillclimber/blob/main/TODO.md

This PLUMED input defines a PATH collective variable using a reference PDB file and specifies the type and lambda parameter. It indicates that the output will be stored in path.spath and path.zpath files.

```plumed
path: PATH REFERENCE=path.pdb TYPE=OPTIMAL LAMBDA=15100.
```

--------------------------------

### Distance Collective Variable Calculation with Python

Source: https://context7.com/zincware/hillclimber/llms.txt

This snippet illustrates how to define a Distance Collective Variable (CV) using Hillclimber to measure distances between atoms or virtual sites. It supports various configurations, including distances between specific atoms, COMs of molecules, and one-to-many or diagonal pairings. The `to_plumed` method generates PLUMED commands for distance calculations.

```python
import hillclimber as hc

ethanol_sel = hc.SMARTSSelector("CCO")
water_sel = hc.SMARTSSelector("O")

# Distance between two specific atoms
dist = hc.DistanceCV(
    x1=ethanol_sel[0][0],  # First atom of first ethanol
    x2=water_sel[0][0],    # First atom of first water
    prefix="d_atoms"
)

# Distance between molecule COMs
dist = hc.DistanceCV(
    x1=hc.VirtualAtom(ethanol_sel[0], "com"),
    x2=hc.VirtualAtom(water_sel[0], "com"),
    prefix="d_com"
)

# One-to-many: First ethanol COM to all water COMs (creates 3 CVs)
dist = hc.DistanceCV(
    x1=hc.VirtualAtom(ethanol_sel[0], "com"),
    x2=hc.VirtualAtom(water_sel, "com"),
    prefix="d",
    pairwise="all"
)

# Diagonal pairing (3 waters × 2 ethanols → only 2 CVs)
dist = hc.DistanceCV(
    x1=hc.VirtualAtom(water_sel, "com"),
    x2=hc.VirtualAtom(ethanol_sel, "com"),
    prefix="d",
    pairwise="diagonal"
)

# Generate PLUMED commands
labels, commands = dist.to_plumed(atoms)
```

--------------------------------

### Apply Harmonic Restraint Bias to Collective Variable

Source: https://context7.com/zincware/hillclimber/llms.txt

Applies harmonic restraints to a collective variable (CV) using the RestraintBias. It requires a defined CV, a force constant (kappa), and a target value (at). The `to_plumed` method generates both the CV and restraint definitions for PLUMED. Requires DistanceCV and RestraintBias.

```python
import hillclimber as hc

# Create a CV
distance_cv = hc.DistanceCV(
    x1=hc.VirtualAtom(water_sel[0], "com"),
    x2=hc.VirtualAtom(ethanol_sel[0], "com"),
    prefix="d"
)

# Restrain distance around 2.5 Å with force constant 200 eV/Å²
restraint = hc.RestraintBias(
    cv=distance_cv,
    kappa=200.0,  # eV/Å² (ASE units)
    at=2.5,       # Å
    label="d_restraint"
)

# Generate PLUMED commands (includes CV definition + restraint)
commands = restraint.to_plumed(atoms)
```
Tokens:
Raw