# PICMI - Particle-In-Cell Modeling Interface

The Particle-In-Cell Modeling Interface (PICMI) standard establishes conventions for the naming and structuring of input files for Particle-In-Cell (PIC) simulations. PICMI provides a unified Python interface that allows researchers to write simulation scripts that are portable across different PIC codes such as WarpX, FBPIC, and Warp. The standard defines a common language for specifying grids, particles, fields, lasers, diagnostics, and other simulation components.

PICMI is designed to be code-agnostic while allowing code-specific extensions through prefixed keyword arguments. The base classes provided in this package serve as the foundation that implementing codes inherit from to provide their specific functionality. This approach enables users to switch between PIC codes with minimal changes to their simulation scripts, promoting reproducibility and collaboration in plasma physics research.

## Installation

```bash
pip install picmistandard
```

## Simulation Setup

The `PICMI_Simulation` class is the main container that brings together all simulation components including the field solver, particle species, lasers, applied fields, and diagnostics.

```python
from picmistandard import picmi

# Create a 3D Cartesian grid
grid = picmi.Cartesian3DGrid(
    number_of_cells=[64, 64, 480],
    lower_bound=[-30e-6, -30e-6, -38e-6],
    upper_bound=[30e-6, 30e-6, 10e-6],
    lower_boundary_conditions=['periodic', 'periodic', 'open'],
    upper_boundary_conditions=['periodic', 'periodic', 'open'],
    moving_window_velocity=[0., 0., picmi.constants.c]
)

# Create electromagnetic solver with Cole-Karkkainen-Cowan stencil
smoother = picmi.BinomialSmoother(n_pass=[1, 1, 1], compensation=[True, True, True])
solver = picmi.ElectromagneticSolver(
    grid=grid,
    cfl=1.0,
    method='CKC',  # Options: 'Yee', 'CKC', 'Lehe', 'PSTD', 'PSATD', 'DS', 'ECT'
    source_smoother=smoother
)

# Initialize the simulation
sim = picmi.Simulation(
    solver=solver,
    time_step_size=1e-15,  # seconds
    max_steps=1000,
    verbose=1,
    particle_shape='linear'  # Options: 'NGP', 'linear', 'quadratic', 'cubic'
)

# Run the simulation
sim.step(1000)  # Run for 1000 time steps
```

## Grid Classes

PICMI provides multiple grid types for different simulation geometries including 1D, 2D, 3D Cartesian, and cylindrical (RZ) grids.

```python
from picmistandard import picmi

# 3D Cartesian Grid - using individual parameters
grid_3d = picmi.Cartesian3DGrid(
    nx=64, ny=64, nz=480,
    xmin=-30e-6, xmax=30e-6,
    ymin=-30e-6, ymax=30e-6,
    zmin=-38e-6, zmax=10e-6,
    bc_xmin='periodic', bc_xmax='periodic',
    bc_ymin='periodic', bc_ymax='periodic',
    bc_zmin='open', bc_zmax='open',
    moving_window_velocity=[0., 0., 3e8],
    guard_cells=[8, 8, 8],
    pml_cells=[10, 10, 10]
)

# Cylindrical (RZ) Grid for axisymmetric simulations
grid_rz = picmi.CylindricalGrid(
    nr=32, nz=480,
    rmin=0., rmax=30e-6,
    zmin=-38e-6, zmax=10e-6,
    bc_rmin=None,  # Axis boundary
    bc_rmax='open',
    bc_zmin='open', bc_zmax='open',
    n_azimuthal_modes=2,
    moving_window_velocity=[0., 3e8]
)

# 2D Cartesian Grid
grid_2d = picmi.Cartesian2DGrid(
    number_of_cells=[64, 480],
    lower_bound=[-30e-6, -38e-6],
    upper_bound=[30e-6, 10e-6],
    lower_boundary_conditions=['periodic', 'open'],
    upper_boundary_conditions=['periodic', 'open']
)

# Add mesh refinement region
grid_3d.add_refined_region(
    level=1,
    lo=[-10e-6, -10e-6, -20e-6],
    hi=[10e-6, 10e-6, 5e-6],
    refinement_factor=[2, 2, 2]
)
```

## Particle Species

PICMI supports defining particle species with various particle types, distributions, and physical properties.

```python
from picmistandard import picmi

# Define an electron species with analytic density distribution
plasma_distribution = picmi.AnalyticDistribution(
    density_expression="1.e23*(1+tanh((z - 20.e-6)/10.e-6))/2.",
    lower_bound=[-20e-6, -20e-6, 0.0],
    upper_bound=[20e-6, 20e-6, 1e-3],
    rms_velocity=[0., 0., 0.],
    directed_velocity=[0., 0., 0.],
    fill_in=True  # Fill in particles as moving window advances
)

electrons = picmi.Species(
    particle_type='electron',
    name='electrons',
    initial_distribution=plasma_distribution,
    particle_shape='linear',
    method='Boris'  # Options: 'Boris', 'Vay', 'Higuera-Cary', 'Li', 'free-streaming'
)

# Define an ion species with explicit charge and mass
protons = picmi.Species(
    particle_type='proton',
    name='protons',
    charge_state=1,
    initial_distribution=plasma_distribution
)

# Multi-species definition (e.g., for plasma with multiple ion types)
multi_species = picmi.MultiSpecies(
    particle_types=['He', 'Ar', 'electron'],
    names=['He+', 'Argon', 'e-'],
    charge_states=[1, 5, None],
    proportions=[0.2, 0.8, 0.2 + 5*0.8],
    initial_distribution=plasma_distribution
)

# Access individual species from multi-species
helium = multi_species['He+']       # By name
argon = multi_species[1]            # By index
```

## Particle Distributions

PICMI provides various particle distribution classes for initializing particle positions and velocities.

```python
from picmistandard import picmi
import numpy as np

# Uniform density distribution
uniform_dist = picmi.UniformDistribution(
    density=1e23,  # particles per m^3
    lower_bound=[-10e-6, -10e-6, 0.],
    upper_bound=[10e-6, 10e-6, 100e-6],
    rms_velocity=[1e5, 1e5, 1e5],
    directed_velocity=[0., 0., 1e7],
    fill_in=True
)

# Gaussian bunch distribution (for beams)
gaussian_bunch = picmi.GaussianBunchDistribution(
    n_physical_particles=1e8,
    rms_bunch_size=[1e-6, 1e-6, 1e-6],
    rms_velocity=[0., 0., 10.*3e8],
    centroid_position=[0., 0., -22e-6],
    centroid_velocity=[0., 0., 1000.*3e8],
    velocity_divergence=[0., 0., 0.]
)

# Analytic distribution with custom expression
analytic_dist = picmi.AnalyticDistribution(
    density_expression='((x**2+y**2)<rmax**2)*n0',
    rmax=5e-6,  # Parameter used in expression
    n0=1e20,    # Parameter used in expression
    momentum_expressions=[None, None, 'vz0'],
    vz0=1e6,    # Parameter for z-momentum
    lower_bound=[-10e-6, -10e-6, 0.],
    upper_bound=[10e-6, 10e-6, 50e-6]
)

# Foil distribution with exponential pre/post plasma
foil_dist = picmi.FoilDistribution(
    density=1e28,
    front=10e-6,
    thickness=1e-6,
    exponential_pre_plasma_length=2e-6,
    exponential_pre_plasma_cutoff=10e-6,
    exponential_post_plasma_length=1e-6,
    lower_bound=[-20e-6, -20e-6, None],
    upper_bound=[20e-6, 20e-6, None]
)

# Particle list distribution (explicit particle data)
particle_list = picmi.ParticleListDistribution(
    x=[0., 1e-6, 2e-6],
    y=[0., 0., 0.],
    z=[0., 0., 0.],
    ux=[0., 0., 0.],
    uy=[0., 0., 0.],
    uz=[1e8, 1e8, 1e8],
    weight=1e6
)

# Flux distribution for particle injection from a plane
flux_dist = picmi.AnalyticFluxDistribution(
    flux='1e20*exp(-t/tau)',
    tau=1e-12,  # Parameter used in expression
    flux_normal_axis='z',
    surface_flux_position=0.,
    flux_direction=1,
    lower_bound=[-10e-6, -10e-6, None],
    upper_bound=[10e-6, 10e-6, None],
    rms_velocity=[1e5, 1e5, 1e6],
    flux_tmin=0.,
    flux_tmax=1e-12
)
```

## Particle Layouts

Particle layouts define how particles are distributed spatially in the simulation domain.

```python
from picmistandard import picmi

# Gridded layout: uniform particles per cell
gridded_layout = picmi.GriddedLayout(
    n_macroparticle_per_cell=[2, 2, 2],  # 2x2x2 = 8 particles per cell
    grid=grid  # Optional: use specific grid, defaults to simulation grid
)

# Pseudo-random layout with total particle count
random_layout = picmi.PseudoRandomLayout(
    n_macroparticles=100000,
    seed=42
)

# Pseudo-random layout with particles per cell
random_layout_ppc = picmi.PseudoRandomLayout(
    n_macroparticles_per_cell=8,
    seed=42,
    grid=grid
)

# Add species to simulation with layout
sim.add_species(species=electrons, layout=gridded_layout)
sim.add_species(
    species=beam,
    layout=random_layout,
    initialize_self_field=True  # Calculate initial space-charge fields
)
```

## Field Solvers

PICMI supports electromagnetic, electrostatic, and magnetostatic field solvers with various numerical methods.

```python
from picmistandard import picmi

# Electromagnetic solver (Maxwell's equations)
em_solver = picmi.ElectromagneticSolver(
    grid=grid,
    method='PSATD',  # Pseudo-Spectral Analytical Time Domain
    stencil_order=[-1, -1, -1],  # -1 = infinite order (spectral)
    cfl=0.99,
    source_smoother=picmi.BinomialSmoother(n_pass=[1, 1, 1]),
    galilean_velocity=[0., 0., 3e8],  # For Galilean PSATD
    divE_cleaning=True,
    divB_cleaning=True
)

# Electrostatic solver (Poisson equation)
es_solver = picmi.ElectrostaticSolver(
    grid=grid,
    method='FFT',  # Options: 'FFT', 'Multigrid'
    required_precision=1e-10,
    maximum_iterations=1000
)

# Magnetostatic solver
ms_solver = picmi.MagnetostaticSolver(
    grid=grid,
    method='FFT'
)

# Binomial smoother for current/charge density
smoother = picmi.BinomialSmoother(
    n_pass=[1, 1, 1],
    compensation=[True, True, True],
    stride=[1, 1, 1],
    alpha=[0.5, 0.5, 0.5]
)
```

## Laser Pulses

PICMI provides classes for defining and injecting laser pulses into simulations.

```python
from picmistandard import picmi
import math

# Gaussian laser pulse
gaussian_laser = picmi.GaussianLaser(
    wavelength=800e-9,           # 800 nm wavelength
    waist=5e-6,                  # 5 micron waist at focus
    duration=30e-15,             # 30 fs duration
    propagation_direction=[0, 0, 1],
    polarization_direction=[1, 0, 0],
    focal_position=[0., 0., 100e-6],
    centroid_position=[0., 0., 0.],
    a0=4.0,                      # Normalized vector potential
    phi0=0.,                     # Carrier-envelope phase
    zeta=0.,                     # Spatial chirp
    beta=0.,                     # Angular dispersion
    phi2=0.,                     # Temporal chirp
    fill_in=True,
    name='main_laser'
)

# Alternative: specify E0 instead of a0
laser_with_E0 = picmi.GaussianLaser(
    wavelength=800e-9,
    waist=5e-6,
    duration=30e-15,
    propagation_direction=[0, 0, 1],
    polarization_direction=[math.cos(math.pi/4), math.sin(math.pi/4), 0],
    focal_position=[0., 0., 50e-6],
    centroid_position=[0., 0., -50e-6],
    E0=1e12  # V/m
)

# Analytic laser with custom field expression
analytic_laser = picmi.AnalyticLaser(
    field_expression='Emax*exp(-(X**2+Y**2)/w0**2)*cos(2*pi*c*t/wavelength)',
    Emax=1e12,
    w0=5e-6,
    wavelength=800e-9,
    c=3e8,
    propagation_direction=[0, 0, 1],
    polarization_direction=[1, 0, 0],
    amax=2.0
)

# Laser antenna for injection
antenna = picmi.LaserAntenna(
    position=[0., 0., 9e-6],
    normal_vector=[0., 0., 1.]
)

# Add laser to simulation
sim.add_laser(laser=gaussian_laser, injection_method=antenna)
```

## Applied Fields

PICMI supports constant, analytic, and file-loaded applied electromagnetic fields.

```python
from picmistandard import picmi

# Constant applied field
constant_field = picmi.ConstantAppliedField(
    Ex=0., Ey=0., Ez=1e6,  # V/m
    Bx=0., By=0.5, Bz=0.,  # Tesla
    lower_bound=[-10e-6, -10e-6, 0.],
    upper_bound=[10e-6, 10e-6, 100e-6]
)

# Analytic applied field with expressions
analytic_field = picmi.AnalyticAppliedField(
    Ex_expression='E0*sin(k*z - omega*t)',
    Ey_expression=None,
    Ez_expression=None,
    Bx_expression=None,
    By_expression='E0/c*sin(k*z - omega*t)',
    Bz_expression=None,
    E0=1e9,      # Parameter used in expressions
    k=7.85e6,    # Parameter used in expressions
    omega=2.36e15,  # Parameter used in expressions
    c=3e8,       # Parameter used in expressions
    lower_bound=[None, None, 0.],
    upper_bound=[None, None, 1e-3]
)

# Mirror (perfectly reflecting surface)
mirror = picmi.Mirror(
    z_front_location=50e-6,
    depth=1e-6,
    number_of_cells=4
)

# Load field from OpenPMD file (applied directly to particles)
loaded_applied_field = picmi.LoadAppliedField(
    read_fields_from_path='external_fields.h5',
    load_B=True,
    load_E=True
)

# Load field from file to initialize grid
loaded_gridded_field = picmi.LoadGriddedField(
    read_fields_from_path='initial_fields.h5',
    load_B=True,
    load_E=True
)

# Add applied fields to simulation
sim.add_applied_field(constant_field)
sim.add_applied_field(analytic_field)
```

## Diagnostics

PICMI provides comprehensive diagnostic classes for outputting field and particle data.

```python
from picmistandard import picmi

# Field diagnostic in simulation frame
field_diag = picmi.FieldDiagnostic(
    grid=grid,
    period=100,  # Output every 100 time steps
    data_list=['E', 'B', 'J', 'rho'],  # Fields to output
    write_dir='./diags',
    step_min=0,
    step_max=10000,
    number_of_cells=[32, 32, 240],  # Downsampled output
    lower_bound=[-15e-6, -15e-6, -20e-6],
    upper_bound=[15e-6, 15e-6, 5e-6],
    parallelio=True,
    name='field_output'
)

# Electrostatic field diagnostic
es_field_diag = picmi.ElectrostaticFieldDiagnostic(
    grid=grid,
    period=50,
    data_list=['rho', 'E', 'phi'],
    write_dir='./diags_es'
)

# Particle diagnostic
particle_diag = picmi.ParticleDiagnostic(
    period=100,
    species=[electrons, beam],
    data_list=['position', 'momentum', 'weighting'],
    write_dir='./diags',
    step_min=0,
    step_max=10000,
    parallelio=True,
    name='particle_output'
)

# Boundary scraping diagnostic (particles absorbed at boundaries)
scraping_diag = picmi.ParticleBoundaryScrapingDiagnostic(
    period=100,
    species=[electrons],
    data_list=['position', 'momentum', 'weighting'],
    write_dir='./boundary_particles',
    parallelio=True,
    name='boundary_scraping'
)

# Lab frame field diagnostic (for boosted frame simulations)
lab_field_diag = picmi.LabFrameFieldDiagnostic(
    grid=grid,
    num_snapshots=100,
    dt_snapshots=1e-13,
    data_list=['E', 'B', 'rho'],
    z_subsampling=2,
    time_start=0.,
    write_dir='./lab_diags',
    name='lab_field_output'
)

# Lab frame particle diagnostic
lab_particle_diag = picmi.LabFrameParticleDiagnostic(
    grid=grid,
    num_snapshots=100,
    dt_snapshots=1e-13,
    species=[electrons],
    data_list=['position', 'momentum'],
    time_start=0.,
    write_dir='./lab_diags',
    name='lab_particle_output'
)

# Add diagnostics to simulation
sim.add_diagnostic(field_diag)
sim.add_diagnostic(particle_diag)
```

## Interactions

PICMI supports physics interactions such as field ionization.

```python
from picmistandard import picmi

# Define ion and electron species
nitrogen = picmi.Species(
    particle_type='N',
    name='nitrogen',
    charge_state=0,
    initial_distribution=uniform_dist
)

ionization_electrons = picmi.Species(
    particle_type='electron',
    name='ionization_electrons',
    initial_distribution=None
)

# Field ionization using ADK model
ionization = picmi.FieldIonization(
    model='ADK',  # Ammosov-Delone-Krainov model
    ionized_species=nitrogen,
    product_species=ionization_electrons
)

# Add species and interaction to simulation
sim.add_species(species=nitrogen, layout=gridded_layout)
sim.add_species(species=ionization_electrons, layout=None)
sim.add_interaction(ionization)
```

## Complete Laser-Plasma Acceleration Example

A comprehensive example demonstrating laser wakefield acceleration simulation setup.

```python
#!/usr/bin/env python3
from plasmacode import picmi  # Replace with your PIC code's PICMI import
import math

# Physical constants
c = picmi.constants.c

# Laser parameters
laser_wavelength = 800e-9      # 800 nm
laser_waist = 5e-6             # 5 micron waist
laser_duration = 30e-15        # 30 fs
laser_a0 = 4.0                 # Normalized amplitude
laser_focal_distance = 100e-6  # Focus position

# Plasma parameters
plasma_density = "1.e23*(1+tanh((z - 20.e-6)/10.e-6))/2."
plasma_min = [-20e-6, -20e-6, 0.]
plasma_max = [20e-6, 20e-6, 1e-3]

# Grid setup
grid = picmi.Cartesian3DGrid(
    number_of_cells=[64, 64, 480],
    lower_bound=[-30e-6, -30e-6, -38e-6],
    upper_bound=[30e-6, 30e-6, 10e-6],
    lower_boundary_conditions=['periodic', 'periodic', 'open'],
    upper_boundary_conditions=['periodic', 'periodic', 'open'],
    moving_window_velocity=[0., 0., c]
)

# Solver setup
smoother = picmi.BinomialSmoother(n_pass=[1, 1, 1], compensation=[True, True, True])
solver = picmi.ElectromagneticSolver(grid=grid, cfl=1., method='CKC', source_smoother=smoother)

# Laser setup
laser = picmi.GaussianLaser(
    wavelength=laser_wavelength,
    waist=laser_waist,
    duration=laser_duration,
    focal_position=[0., 0., laser_focal_distance + 9e-6],
    centroid_position=[0., 0., 9e-6 - c*30e-15],
    polarization_direction=[0., 1., 0.],
    propagation_direction=[0., 0., 1.],
    a0=laser_a0
)

antenna = picmi.LaserAntenna(position=[0., 0., 9e-6], normal_vector=[0., 0., 1.])

# Plasma species
plasma_dist = picmi.AnalyticDistribution(
    density_expression=plasma_density,
    lower_bound=plasma_min,
    upper_bound=plasma_max,
    fill_in=True
)

electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=plasma_dist)
ions = picmi.Species(particle_type='H', name='protons', charge_state=1, initial_distribution=plasma_dist)

# Electron beam (witness bunch)
beam_dist = picmi.GaussianBunchDistribution(
    n_physical_particles=1e8,
    rms_bunch_size=[1e-6, 1e-6, 1e-6],
    rms_velocity=[0., 0., 10.*c],
    centroid_position=[0., 0., -22e-6],
    centroid_velocity=[0., 0., 1000.*c]
)
beam = picmi.Species(particle_type='electron', name='beam', initial_distribution=beam_dist)

# Particle layouts
plasma_layout = picmi.GriddedLayout(n_macroparticle_per_cell=[2, 2, 2], grid=grid)
beam_layout = picmi.PseudoRandomLayout(n_macroparticles=100000, seed=0)

# Diagnostics
field_diag = picmi.FieldDiagnostic(grid=grid, period=100, name='fields')
particle_diag = picmi.ParticleDiagnostic(period=100, species=[electrons, beam], name='particles')

# Assemble simulation
sim = picmi.Simulation(solver=solver, verbose=1)
sim.add_laser(laser, injection_method=antenna)
sim.add_species(species=electrons, layout=plasma_layout)
sim.add_species(species=ions, layout=plasma_layout)
sim.add_species(species=beam, layout=beam_layout, initialize_self_field=True)
sim.add_diagnostic(field_diag)
sim.add_diagnostic(particle_diag)

# Run simulation
sim.step(1000)
```

## Summary

PICMI serves as a powerful abstraction layer for Particle-In-Cell simulations, enabling researchers to write portable simulation scripts that work across multiple PIC codes. The standard covers all essential aspects of PIC simulations including grid definition (1D, 2D, 3D Cartesian, and cylindrical geometries), particle species with various distribution types (uniform, Gaussian, analytic, foil), field solvers (electromagnetic, electrostatic, magnetostatic with multiple numerical methods), laser pulse injection (Gaussian and analytic profiles), applied fields (constant, analytic, file-loaded), and comprehensive diagnostics for both simulation and lab frame outputs.

The key integration pattern involves creating a `Simulation` object with a field solver, then progressively adding species with their layouts, lasers with injection methods, applied fields, interactions like ionization, and diagnostics. Code-specific extensions can be passed as keyword arguments prefixed with the code name (e.g., `warpx_max_grid_size`). Implementing codes inherit from the PICMI base classes and override the `init` method to translate the standard parameters into code-specific configurations. This design allows users to develop simulation scripts once and run them on different PIC codes with minimal modifications, greatly facilitating code benchmarking, verification, and collaboration in the plasma physics community.