Running Simulations¶

This guide provides progressively advanced examples of running cavity-coupled molecular dynamics simulations with Cavity HOOMD.

Basic Cavity Simulations ¶

Using the 05_advanced_run.py Script ¶

The primary interface for running simulations is the examples/05_advanced_run.py script, which provides a comprehensive command-line interface for cavity simulations.

Basic Syntax:

python examples/05_advanced_run.py [OPTIONS]

Essential Parameters:

--coupling VALUE: Cavity-molecule coupling strength (atomic units)
--temperature VALUE: System temperature in Kelvin
--frequency VALUE: Cavity frequency in cm⁻¹
--runtime VALUE: Simulation time in picoseconds

Choosing Coupling Strength ¶

The coupling strength determines the interaction between molecules and the cavity mode.

Weak Coupling (g < 10⁻⁴):

python examples/05_advanced_run.py --coupling 1e-5 --runtime 1000

Perturbative regime
Minimal modification of molecular dynamics
Useful for baseline comparisons

Moderate Coupling (10⁻⁴ ≤ g ≤ 10⁻³):

python examples/05_advanced_run.py --coupling 5e-4 --runtime 1000

Observable cavity-molecule interactions
Energy exchange between modes
Typical for initial studies

Strong Coupling (g > 10⁻³):

python examples/05_advanced_run.py --coupling 1e-3 --runtime 1000

Polariton formation regime
Significant modification of dynamics
Collective strong coupling effects

Physical Interpretation:

The effective collective coupling scales as:

\[g_{\text{eff}} = g \sqrt{N}\]

where N is the number of molecules. For N=512 molecules with g=1×10⁻³:

\[g_{\text{eff}} = 10^{-3} \times \sqrt{512} \approx 0.0226\]

This collective enhancement enables strong coupling even with modest single-molecule coupling.

Temperature and Frequency Selection ¶

Temperature Guidelines:

# Low temperature (minimize thermal fluctuations)
python examples/05_advanced_run.py --coupling 1e-3 --temperature 50 --runtime 1000

# Room temperature
python examples/05_advanced_run.py --coupling 1e-3 --temperature 300 --runtime 1000

# High temperature (study thermal effects)
python examples/05_advanced_run.py --coupling 1e-3 --temperature 500 --runtime 1000

Cavity Frequency Selection:

The cavity frequency should be chosen relative to molecular vibrational frequencies:

# On-resonance with O-O stretch (~1550 cm⁻¹)
python examples/05_advanced_run.py --coupling 1e-3 --frequency 1550 --runtime 1000

# Red-detuned
python examples/05_advanced_run.py --coupling 1e-3 --frequency 1200 --runtime 1000

# Blue-detuned
python examples/05_advanced_run.py --coupling 1e-3 --frequency 2000 --runtime 1000

Resonance Condition:

Maximum coupling effects occur when:

\[\omega_{\text{cavity}} \approx \omega_{\text{molecular}}\]

For O₂ molecules, the main vibrational mode is around 1550 cm⁻¹.

Thermostat Options ¶

Cavity HOOMD supports multiple thermostats for both molecular and cavity degrees of freedom.

Molecular Thermostats ¶

Bussi Thermostat (Default, Recommended):

python examples/05_advanced_run.py --molecular-bath bussi --coupling 1e-3 --runtime 1000

Advantages:

Correct canonical ensemble sampling
Fast equilibration
Minimal perturbation to dynamics
Suitable for both equilibrium and dynamical studies

Theory: Bussi-Donadio-Parrinello stochastic velocity rescaling thermostat samples the canonical (NVT) ensemble rigorously.

Langevin Thermostat:

python examples/05_advanced_run.py --molecular-bath langevin --coupling 1e-3 --runtime 1000

Advantages:

Simple and robust
Well-established in MD community
Good for strongly coupled systems

Considerations: Introduces friction that may affect fast dynamics.

No Thermostat (NVE):

python examples/05_advanced_run.py --molecular-bath none --coupling 1e-3 --runtime 1000

Use Cases:

Microcanonical ensemble studies
Energy conservation testing
Short-time dynamics analysis

Warning: Temperature will drift if system is not at equilibrium initially.

Cavity Thermostats ¶

Langevin (Default, Recommended):

python examples/05_advanced_run.py --cavity-bath langevin --coupling 1e-3 --runtime 1000

Advantages:

Represents cavity dissipation physically
Prevents energy accumulation in cavity mode
Mimics realistic cavity loss

Friction Coefficient: Default γ=0.1 corresponds to cavity lifetime τ ≈ 10 ps.

Bussi Thermostat:

python examples/05_advanced_run.py --cavity-bath bussi --coupling 1e-3 --runtime 1000

Use Cases:

Canonical ensemble for cavity mode
Comparison with molecular Bussi thermostat
Studying cavity thermal fluctuations

No Thermostat (Conservative):

python examples/05_advanced_run.py --cavity-bath none --coupling 1e-3 --runtime 1000

Use Cases:

Undamped cavity dynamics
Energy conservation testing
Short-time coherent dynamics

Warning: Cavity mode may accumulate large energy without dissipation.

Recommended Thermostat Combinations ¶

Standard Setup (Most Common):

python examples/05_advanced_run.py \
    --molecular-bath bussi \
    --cavity-bath langevin \
    --coupling 1e-3 \
    --runtime 1000

Molecules: Canonical ensemble (Bussi)
Cavity: Dissipative Langevin (realistic loss)

Pure Canonical Ensemble:

python examples/05_advanced_run.py \
    --molecular-bath bussi \
    --cavity-bath bussi \
    --coupling 1e-3 \
    --runtime 1000

Both molecules and cavity in canonical ensemble
Useful for equilibrium thermodynamics

Conservative Dynamics:

python examples/05_advanced_run.py \
    --molecular-bath none \
    --cavity-bath none \
    --coupling 1e-3 \
    --runtime 100

Microcanonical (NVE) ensemble
Total energy conserved
Best for short-time dynamics and energy conservation tests

Cavity Modes ¶

Cavity HOOMD supports two types of cavity modes.

q=0 Mode (Uniform, Default)¶

Description: Uniform electric field throughout the simulation box.

python examples/05_advanced_run.py --coupling 1e-3 --runtime 1000
# No additional flags needed (default mode)

Properties:

Cavity mode has no spatial variation
k-vector: \(\vec{k} = (0, 0, 0)\)
All molecules couple with same phase
Maximum collective enhancement
Two cavity particles (x and y polarizations)

When to Use:

Strong coupling studies
Maximum collective effects
Initial explorations

Cavity Hamiltonian:

\[H_{\text{cavity}} = \frac{1}{2}K(q_x^2 + q_y^2) + g(q_x D_x + q_y D_y) + \frac{g^2}{2K}(D_x^2 + D_y^2)\]

where \(D_x, D_y\) are total system dipole moments.

Finite-q Mode ¶

Description: Standing wave with spatial modulation.

python examples/05_advanced_run.py --coupling 1e-3 --finite-q --runtime 1000

Properties:

k-vector: \(\vec{k} = (0, 0, k_z)\)
Spatial phase variation: molecules at different positions couple with different phases
Reduces collective enhancement
Cavity particles positioned at finite coordinates

When to Use:

Realistic cavity geometries
Wave vector dependence studies
Reduced collective coupling

Phase Factor:

Each molecule at position \(\vec{r}_n\) couples with phase:

\[\phi_n = \vec{k} \cdot \vec{r}_n = k_z z_n\]

Coupling Strength:

\[H_{\text{coupling}} = g q \sum_{n=1}^{N} d_n \cos(\phi_n)\]

The cosine factor reduces collective enhancement compared to q=0 mode.

Cavity Particle Positioning:

In finite-q mode, cavity particles are positioned in the simulation box such that their displacement represents the mode amplitude and maintains energy conservation.

Multiple Replicas ¶

Running Independent Replicas ¶

For statistical analysis, run multiple independent simulations with different random seeds:

Sequential Replica Numbering:

python examples/05_advanced_run.py --coupling 1e-3 --replicas 1-5 --runtime 1000

This runs 5 replicas (replica IDs: 1, 2, 3, 4, 5) sequentially.

Custom Replica IDs:

python examples/05_advanced_run.py --coupling 1e-3 --replicas 1,5,10 --runtime 1000

This runs 3 replicas with IDs 1, 5, and 10.

Large Replica Sets:

python examples/05_advanced_run.py --coupling 1e-3 --replicas 1-50 --runtime 1000

For comprehensive statistics (50 replicas).

Output Organization ¶

Each replica creates a separate output file:

cavity_coupling_1e-03/
├── prod-1.gsd
├── prod-1_energy_tracker.txt
├── prod-1_cavity_mode.txt
├── prod-2.gsd
├── prod-2_energy_tracker.txt
├── prod-2_cavity_mode.txt
├── ...
├── prod-5.gsd
├── prod-5_energy_tracker.txt
└── prod-5_cavity_mode.txt

Replica Differences:

Random number seed (determined by replica ID)
Initial velocities (Boltzmann distribution)
Stochastic thermostat evolution

Replica Similarities:

Same initial positions
Same parameters (coupling, temperature, frequency)
Same forces and integration scheme

Statistical Analysis ¶

Analyze replica statistics using Python:

import numpy as np
import glob

# Load energy from all replicas
energy_files = glob.glob('cavity_coupling_1e-03/prod-*_energy_tracker.txt')

energies = []
for file in sorted(energy_files):
    data = np.loadtxt(file, skiprows=1, delimiter='\t')
    energies.append(data[:, -1])  # Total energy column

energies = np.array(energies)

# Calculate statistics
mean_energy = np.mean(energies, axis=0)
std_energy = np.std(energies, axis=0)
sem_energy = std_energy / np.sqrt(len(energies))

print(f"Final energy: {mean_energy[-1]:.2f} ± {sem_energy[-1]:.2f}")

Standard Error of the Mean:

For N independent replicas:

\[\text{SEM} = \frac{\sigma}{\sqrt{N}}\]

where σ is the standard deviation across replicas.

Device Selection ¶

CPU Execution ¶

python examples/05_advanced_run.py --device CPU --coupling 1e-3 --runtime 500

When to Use:

No GPU available
Small systems (<1000 particles)
Debugging and testing

Performance: ~100-1000 steps/second for typical systems.

GPU Execution ¶

python examples/05_advanced_run.py --device GPU --coupling 1e-3 --runtime 1000

Requirements:

NVIDIA GPU with CUDA 11.0+
GPU-enabled HOOMD-blue build
GPU-enabled Cavity HOOMD build

When to Use:

Large systems (>1000 particles)
Long simulations
Production runs

Performance: ~10,000-100,000 steps/second depending on GPU.

GPU Selection:

For multi-GPU systems:

CUDA_VISIBLE_DEVICES=0 python examples/05_advanced_run.py --device GPU --coupling 1e-3
CUDA_VISIBLE_DEVICES=1 python examples/05_advanced_run.py --device GPU --coupling 1e-3