porereax.simulate¶

Simulation setup and management for ReaxFF molecular dynamics simulations.

This module provides functionality to convert GROMACS structure files (.gro) to LAMMPS data files and generate complete simulation workflows for ReaxFF force field calculations. It handles atom type mapping, charge assignment, and creates all necessary input files for running molecular dynamics simulations on HPC systems.

The module automatically generates: - LAMMPS data files from GROMACS structures - LAMMPS input scripts for equilibration and production runs - Job submission scripts for HPC clusters - Analysis scripts for post-processing simulation results

Example

>>> from porereax.simulate import Simulate
>>> gro_lib = {"Si": "Si", "O": "O", "OM": "O", "HW": "H", "MW": ""}
>>> gro_charges = {'Si': 2.4, 'O': -1.2, 'H': 0.6, 'OM': -0.8}
>>> atom_masses = {'Si': 28.085, 'O': 15.999, 'H': 1.008}
>>> sim = Simulate(gro_lib, gro_charges, atom_masses, 'system.gro')
>>> sim.set_force_field('reax.ffield')
>>> sim.add_sim('nvt', nsteps=100000, temp=300)
>>> sim.generate()

Classes¶

Simulate

Main class for setting up and managing ReaxFF molecular dynamics simulations.

Module Contents¶

class porereax.simulate.Simulate(gro_lib: dict, gro_charges: dict, atom_masses: dict, structure_file: str = None)[source]¶

Main class for setting up and managing ReaxFF molecular dynamics simulations.

This class handles the conversion of GROMACS structure files to LAMMPS format, manages atom type mappings and charges, and generates all necessary files for running ReaxFF simulations on HPC systems.

Example

>>> gro_lib = {"Si": "Si", "O": "O", "OM": "O", "HW": "H", "MW": ""}
>>> gro_charges = {'Si': 2.4, 'O': -1.2, 'H': 0.6, 'OM': -0.8}
>>> atom_masses = {'Si': 28.085, 'O': 15.999, 'H': 1.008}
>>> sim = Simulate(gro_lib, gro_charges, atom_masses, 'system.gro')

set_job_file(job_file: str = None, submit_command: str = None, lammps_command: str = None)[source]¶

Specify a custom job submission template file and command.

This method allows you to provide a custom job script template for HPC job submission. The template should use Jinja2 syntax and include placeholders for job parameters like node count, wall time, and LAMMPS commands.

Parameters:

job_file (str) – Path to the job submission template file. Must be a valid file that exists. The template should be compatible with the HPC scheduler (SLURM, PBS, etc.).
submit_command (str) – Command to submit jobs to the scheduler (e.g., “sbatch” for SLURM, “qsub” for PBS/Torque). Must be a non-empty string.
lammps_command (str, optional) – Custom command to run LAMMPS. If None, a default MPI command will be used. Use placeholders {input_file} and {log_file} for input and log file names.

Raises:

FileNotFoundError – If the specified job file does not exist.
ValueError – If submit_command is not a non-empty string.

Notes

If this method is not called, a default SLURM template will be used with ‘sbatch’ as the submission command. If lammps_command is not provided, a default MPI command will be used.

Example

>>> sim.set_job_file('/path/to/custom.job', 'sbatch', 'mpirun lmp -in {input_file} -log {log_file}')

set_force_field(force_field: str = None)[source]¶

Specify a custom ReaxFF force field parameter file.

Parameters:: force_field (str) – Path to the ReaxFF force field parameter file (ffield). This file contains all the reactive force field parameters for the atom types in the system.
Raises:: FileNotFoundError – If the specified force field file does not exist.

Notes

If this method is not called, a default Si/O/H force field from https://doi.org/10.1063/1.3407433 will be used.

Example

>>> sim.set_force_field('/path/to/reax.ffield')

add_image_dump(plane='xy', dump_freq=None, zoom=1.5, image_width=1200, image_height=1200, atom_colors=None, atom_sizes=None, map_by_charge=None, kwargs=None)[source]¶

Add an image rendering during LAMMPS simulations. It allows for the generation of image snapshots of the simulation at specified intervals, with customizable viewing planes, zoom levels, and atom visualizations. Multiple image dumps can be added to the simulation workflow with different settings.

Parameters:

plane (str or None, optional) – Viewing plane for the snapshot. Supported values are “xy”, “xz”, “yz”. Default is “xy”.
dump_freq (int, optional) – Frequency (in steps) for writing image snapshots. If None, uses the same frequency as the trajectory dump.
zoom (float, optional) – Camera zoom factor for rendered images.
image_width (int, optional) – Output image width in pixels.
image_height (int, optional) – Output image height in pixels.
atom_colors (dict or None, optional) – Optional mapping of atom names to LAMMPS color names (example: {“Si”: “yellow”, “O”: “red”}).
atom_sizes (dict or None, optional) – Optional mapping of atom names to relative sphere sizes for rendering (example: {“Si”: 1.0, “O”: 0.8}).
map_by_charge (str or None, optional) – If not None, must be a string representing an amap string to map atom colors by their partial charges using a color gradient. If provided, this overrides atom_colors. Example: “-1 2 ca 0.0 3 min royalblue 0 green max orangered”
kwargs (str or None, optional) – Additional keyword arguments for the LAMMPS dump command. This can be used to pass extra options to the dump command.

add_sim(type: str, nsteps: int, temp: float, pressure: float = 1.0, dt: float = 0.5, nodes: int = 1, tasks_per_node: int = 64, wall_time: str = '20:00:00', dump_freq: int = 100, thermo_freq: int = 100) → None[source]¶

Add a simulation step to the workflow.

This method adds a molecular dynamics simulation step with specified parameters. Multiple simulation steps can be added sequentially to create a multi-stage workflow (e.g., equilibration followed by production runs).

Parameters:

type (str) – Ensemble type for the simulation. Common values: - ‘nvt’: Constant number of particles, volume, and temperature - ‘npt’: Constant number of particles, pressure, and temperature - ‘nve’: Constant number of particles, volume, and energy
nsteps (int) – Number of MD steps to run in this simulation stage.
temp (float) – Temperature in Kelvin for the simulation.
pressure (float, optional) – Pressure in atmospheres for NPT simulations. Default is 1 atm.
dt (float, optional) – Time step in femtoseconds. Default is 0.5 fs.
nodes (int, optional) – Number of compute nodes to request for this job. Default is 1.
tasks_per_node (int, optional) – Number of MPI tasks per node. Default is 64.
wall_time (str, optional) – Maximum wall time for the job in HH:MM:SS format. Default is “20:00:00”.
dump_freq (int, optional) – Frequency (in steps) for writing trajectory snapshots. Default is 100.
thermo_freq (int, optional) – Frequency (in steps) for writing thermodynamic output. Default is 100.

Notes

Simulation steps are executed in the order they are added. Each step will automatically submit the next step upon completion if multiple steps exist.

Example

>>> # Add equilibration run
>>> sim.add_sim('npt', nsteps=500000, temp=300, pressure=1, dt=0.5)
>>> # Add production run
>>> sim.add_sim('nvt', nsteps=50000, temp=300, dt=0.5)

freeze_region(region: collections.abc.Callable[[float, float, float], bool])[source]¶

Freeze atoms in a specified region of the simulation box.

Parameters:: region (callable) – A function that takes three arguments (x, y, z) representing the coordinates of an atom in Angstroms and returns True if the atom should be frozen, or False otherwise.
Raises:: ValueError – If the provided region is not callable.

Notes

This method modifies the atom types to create frozen versions of the specified atom types. Frozen atoms will not move during the simulation, allowing for the study of surface interactions or confinement effects.

auto_freeze(reactive_radius: float = 10.0)[source]¶

Automatically freeze some atoms of the pore structure and leave a reactive layer between frozen and solvent.

Parameters:: reactive_radius (float, optional) – The thickness of the reactive layer in Angstroms. Atoms within this distance from the pore wall will remain unfrozen. Default is 10.0 Angstroms.
Raises:: NotImplementedError – If the pore type specified in the system YAML file is not “cylinder”. Currently, only cylindrical pores are supported for automatic freezing.

generate()[source]¶

Generate all simulation files and scripts.

This method is the main entry point for generating a complete simulation workflow. It creates all necessary files for running ReaxFF simulations:

Converts GROMACS structure to LAMMPS data file (system.data)
Copies or uses default ReaxFF force field file (reax.ffield)
Generates LAMMPS input scripts for initial equilibration and all simulation steps
Creates job submission scripts for each simulation stage
Generates analysis script (ana.py) for post-processing

The method automatically chains job submissions so that each simulation step submits the next one upon completion.

Raises:: ValueError – If atom names in the structure file are not found in gro_lib or if charges are missing from gro_charges.

Notes

If no force field is specified, uses default Si/O/H parameters
If no job template is specified, uses default SLURM template
Virtual sites and excluded atoms (mapped to “” in gro_lib) are filtered out
Prints detailed information about system composition and total charge

Example

>>> sim = Simulate(gro_lib, gro_charges, atom_masses, 'system.gro')
>>> sim.set_force_field('reax.ffield')
>>> sim.add_sim('nvt', nsteps=100000, temp=300)
>>> sim.generate()
Using force field from https://doi.org/10.1063/1.3407433 for Si/O/H systems.
Box dimensions (Angstroms): [50.0, 50.0, 50.0]
Total charge in system: 0.0000e
Atom counts by type:
  Type Si: 100 atoms
  Type O: 200 atoms
LAMMPS data file written to /path/to/system.data