Version: 1.0 Status: Proposed Date: 2025-11-15
This document specifies the Chemistry and Chemical Engineering domain for Morphogen. It defines eight sub-domains covering molecular simulation, reaction kinetics, quantum chemistry, transport phenomena, multiphase systems, thermodynamics, catalysis, and electrochemistry.
The chemistry domain leverages Morphogen's existing infrastructure:
- Operator graph paradigm to wrap external solvers
- Field operators for PDEs (diffusion, advection, reaction)
- Agent operators for particle-based methods
- Optimization domain for design discovery
- ML domain for surrogate models and inverse design
- Visualization domain for molecules, orbitals, and fields
- Molecular Domain
- Reaction Kinetics Domain
- Quantum Chemistry Domain
- Transport Phenomena Domain
- Multiphase Domain
- Thermodynamics Domain
- Catalysis Domain
- Electrochemistry Domain
Represent, manipulate, and analyze molecular structures. Perform classical molecular mechanics calculations, molecular dynamics simulations, and conformer generation.
struct Molecule {
atoms: Vec<Atom>,
bonds: Vec<Bond>,
positions: Field<Vec3<f32 [Angstrom]>>,
velocities: Field<Vec3<f32 [Angstrom/fs]>>,
forces: Field<Vec3<f32 [kcal/(mol·Angstrom)]>>,
masses: Field<f32 [amu]>,
charges: Field<f32 [e]>
}
struct Atom {
element: String,
atomic_number: u32,
mass: f32 [amu],
charge: f32 [e]
}
struct Bond {
atom1: u32,
atom2: u32,
order: f32 # 1.0=single, 2.0=double, 1.5=aromatic
}
struct Trajectory {
frames: Vec<Molecule>,
times: Vec<f32 [fs]>,
energies: Vec<f32 [kcal/mol]>
}
struct ForceField {
name: String, # "amber", "charmm", "uff", "gaff"
parameters: Map<String, f32>
}
# Load from various formats
molecule = molecular.load_smiles("CCO") # SMILES string
molecule = molecular.load_pdb("protein.pdb") # PDB file
molecule = molecular.load_xyz("structure.xyz") # XYZ file
molecule = molecular.load_mol2("ligand.mol2") # MOL2 file
molecule = molecular.load_sdf("database.sdf") # SDF file
# Convert between formats
smiles = molecular.to_smiles(molecule)
pdb = molecular.to_pdb(molecule)
xyz = molecular.to_xyz(molecule)
# Generate 3D coordinates from 2D
molecule = molecular.generate_3d(molecule, force_field="uff")
# Basic properties
mw = molecular.molecular_weight(molecule) -> f32 [g/mol]
formula = molecular.molecular_formula(molecule) -> String
charge = molecular.total_charge(molecule) -> f32 [e]
multiplicity = molecular.spin_multiplicity(molecule) -> u32
# Geometric properties
com = molecular.center_of_mass(molecule) -> Vec3<f32 [Angstrom]>
moi = molecular.moment_of_inertia(molecule) -> Mat3<f32 [amu·Angstrom²]>
rg = molecular.radius_of_gyration(molecule) -> f32 [Angstrom]
# Electronic properties
dipole = molecular.dipole_moment(molecule) -> Vec3<f32 [Debye]>
polarizability = molecular.polarizability(molecule) -> f32 [Angstrom³]
# Topological properties
rings = molecular.find_rings(molecule) -> Vec<Vec<u32>>
aromaticity = molecular.is_aromatic(molecule, ring_id) -> bool
hbond_donors = molecular.count_hbond_donors(molecule) -> u32
hbond_acceptors = molecular.count_hbond_acceptors(molecule) -> u32
# Energy calculation
energy = molecular.compute_energy(
molecule,
force_field="amber", # "amber", "charmm", "uff", "gaff", "oplsaa"
include_terms=["bond", "angle", "dihedral", "vdw", "electrostatic"]
) -> f32 [kcal/mol]
# Force calculation
forces = molecular.compute_forces(
molecule,
force_field="amber"
) -> Field<Vec3<f32 [kcal/(mol·Angstrom)]>>
# Energy components
bond_energy = molecular.bond_energy(molecule, force_field)
angle_energy = molecular.angle_energy(molecule, force_field)
dihedral_energy = molecular.dihedral_energy(molecule, force_field)
vdw_energy = molecular.vdw_energy(molecule, force_field)
electrostatic = molecular.electrostatic_energy(molecule, force_field)
# Minimize energy
molecule_opt = molecular.optimize_geometry(
molecule,
force_field="uff",
method="bfgs", # "steepest_descent", "bfgs", "conjugate_gradient"
max_iterations=1000,
convergence=1e-6 [kcal/(mol·Angstrom)]
)
# Constrained optimization
molecule_opt = molecular.optimize_constrained(
molecule,
force_field="amber",
constraints=[
molecular.fix_atom(0), # Fix atom 0
molecular.fix_distance(1, 2, 1.5 [Angstrom]) # Fix bond length
]
)
# Generate conformers
conformers = molecular.generate_conformers(
molecule,
n=100,
method="rdkit", # "rdkit", "omega", "random"
energy_window=10.0 [kcal/mol],
rms_threshold=0.5 [Angstrom]
) -> Vec<Molecule>
# Cluster conformers
clusters = molecular.cluster_conformers(
conformers,
method="rmsd",
threshold=1.0 [Angstrom]
)
# Neighbor list (for efficient force calculation)
neighbors = molecular.compute_neighbor_list(
molecule,
cutoff=10.0 [Angstrom],
method="cell_list" # "cell_list", "verlet_list", "naive"
)
# MD integrators
positions_new = molecular.velocity_verlet(
positions,
velocities,
forces,
masses,
dt=1.0 [fs]
)
positions_new, velocities_new = molecular.langevin_integrator(
positions,
velocities,
forces,
masses,
temp=300.0 [K],
friction=1.0 [1/ps],
dt=1.0 [fs]
)
# Thermostats
velocities = molecular.berendsen_thermostat(velocities, masses, temp_target=300.0 [K], tau=0.1 [ps])
velocities = molecular.nose_hoover_thermostat(velocities, masses, temp_target=300.0 [K], Q=1.0)
# Barostats
positions, box = molecular.berendsen_barostat(positions, box, pressure_target=1.0 [atm], tau=1.0 [ps])
# Run MD simulation
trajectory = molecular.md_simulate(
molecule,
force_field="amber",
temp=300.0 [K],
pressure=1.0 [atm],
time=10.0 [ns],
dt=2.0 [fs],
ensemble="npt", # "nve", "nvt", "npt"
thermostat="nose_hoover",
barostat="berendsen"
)
# RMSD
rmsd = molecular.rmsd(molecule1, molecule2, align=true) -> f32 [Angstrom]
rmsd_traj = molecular.rmsd_trajectory(trajectory, reference=trajectory.frames[0])
# RMSF (root-mean-square fluctuation)
rmsf = molecular.rmsf(trajectory) -> Field<f32 [Angstrom]>
# Diffusion coefficient
D = molecular.diffusion_coefficient(trajectory) -> f32 [cm²/s]
# Radial distribution function
g_r = molecular.rdf(trajectory, atom_type_1="O", atom_type_2="H", r_max=10.0 [Angstrom])
# Hydrogen bonds
hbonds = molecular.hydrogen_bonds(trajectory, donor_acceptor_distance=3.5 [Angstrom], angle_cutoff=30.0 [deg])
Model chemical reaction rates, reactor behavior, and reaction networks.
struct Reaction {
reactants: Map<String, f32>, # species -> stoichiometry
products: Map<String, f32>,
rate_law: RateLaw,
reversible: bool
}
enum RateLaw {
Arrhenius(A: f32, Ea: f32 [J/mol]),
ModifiedArrhenius(A: f32, n: f32, Ea: f32 [J/mol]),
Custom(fn(conc: Map<String, f32>, temp: f32) -> f32)
}
struct Reactor {
type: ReactorType,
volume: f32 [m³],
temp: f32 [K],
pressure: f32 [Pa]
}
enum ReactorType {
Batch,
CSTR,
PFR,
PBR # Packed bed reactor
}
# Arrhenius kinetics
k = kinetics.arrhenius(
temp=350.0 [K],
A=1e6 [1/s],
Ea=50000.0 [J/mol]
) -> f32 [1/s]
# Modified Arrhenius
k = kinetics.modified_arrhenius(
temp=350.0 [K],
A=1e6,
n=0.5,
Ea=50000.0 [J/mol]
) -> f32 [1/s]
# Temperature-dependent equilibrium constant
K_eq = kinetics.vant_hoff(
temp=350.0 [K],
delta_H=-50000.0 [J/mol],
delta_S=-100.0 [J/(mol·K)]
)
# Define reaction network
reactions = [
kinetics.reaction(
reactants={"A": 1.0},
products={"B": 1.0},
rate_law=kinetics.arrhenius(A=1e6, Ea=50000.0)
),
kinetics.reaction(
reactants={"B": 1.0, "C": 1.0},
products={"D": 1.0},
rate_law=kinetics.arrhenius(A=1e8, Ea=60000.0)
)
]
# Compute reaction rates
rates = kinetics.reaction_rates(
conc={"A": 1.0 [mol/L], "B": 0.5 [mol/L], "C": 0.3 [mol/L]},
temp=350.0 [K],
reactions=reactions
) -> Map<String, f32 [mol/(L·s)]>
# Time evolution (ODE integration)
conc_new = kinetics.integrate_ode(
conc_initial={"A": 1.0, "B": 0.0, "C": 1.0, "D": 0.0},
reactions=reactions,
temp=350.0 [K],
time=3600.0 [s],
method="bdf" # "euler", "rk4", "bdf", "lsoda"
)
# Batch reactor
result = kinetics.batch_reactor(
initial_conc={"A": 1.0 [mol/L]},
reactions=reactions,
temp=350.0 [K],
time=3600.0 [s]
) -> Map<String, f32 [mol/L]>
# CSTR (continuous stirred-tank reactor)
result = kinetics.cstr(
feed_conc={"A": 1.0 [mol/L]},
feed_flow=0.1 [L/s],
volume=10.0 [L],
reactions=reactions,
temp=350.0 [K]
) -> Map<String, f32 [mol/L]>
# PFR (plug flow reactor)
conc_profile = kinetics.pfr(
feed_conc={"A": 1.0 [mol/L]},
reactions=reactions,
length=10.0 [m],
area=0.01 [m²],
flow_velocity=1.0 [m/s],
temp=350.0 [K]
) -> Field1D<Map<String, f32 [mol/L]>>
# Mass-transfer-limited reaction
k_eff = kinetics.mass_transfer_limited(
k_intrinsic=1e6 [1/s],
k_mass_transfer=1e3 [1/s]
) -> f32 [1/s]
# Dispersion in PFR
conc_profile = kinetics.pfr_with_dispersion(
feed_conc={"A": 1.0 [mol/L]},
reactions=reactions,
length=10.0 [m],
velocity=1.0 [m/s],
dispersion_coeff=0.01 [m²/s]
)
Interface with quantum chemistry codes (DFT, ab initio) and ML surrogate models.
# Single-point energy
energy = qchem.dft_energy(
molecule,
method="B3LYP", # "B3LYP", "PBE", "M06-2X", "wB97X-D"
basis="6-31G*", # "sto-3g", "6-31G*", "def2-TZVP", "cc-pVTZ"
code="orca" # "orca", "psi4", "gaussian", "qchem"
) -> f32 [Hartree]
# Geometry optimization
molecule_opt, energy = qchem.dft_optimize(
molecule,
method="B3LYP",
basis="6-31G*",
code="orca"
)
# Forces (gradient)
forces = qchem.dft_forces(
molecule,
method="B3LYP",
basis="6-31G*"
) -> Field<Vec3<f32 [Hartree/Bohr]>>
# Frequency calculation
frequencies = qchem.dft_frequencies(
molecule,
method="B3LYP",
basis="6-31G*"
) -> Vec<f32 [cm⁻¹]>
# Transition state search
ts_molecule, energy = qchem.find_transition_state(
reactant,
product,
method="B3LYP",
basis="6-31G*"
)
# Faster, less accurate methods
energy = qchem.semi_empirical(
molecule,
method="PM7" # "PM3", "PM6", "PM7", "AM1"
) -> f32 [kcal/mol]
# Predict energy with neural network
energy = qchem.ml_pes(
molecule,
model="SchNet" # "SchNet", "MPNN", "DimeNet", "PaiNN"
) -> f32 [kcal/mol]
# Predict forces
forces = qchem.ml_forces(
molecule,
model="SchNet"
) -> Field<Vec3<f32 [kcal/(mol·Angstrom)]>>
# Train ML PES from data
model = qchem.train_pes(
training_data=[
(molecule1, energy1),
(molecule2, energy2),
...
],
architecture="SchNet",
epochs=1000,
learning_rate=1e-4
)
Heat transfer, mass diffusion, convection, and porous media transport.
# Conduction (Fourier's law)
q = transport.conduction(
temp_gradient=10.0 [K/m],
thermal_conductivity=50.0 [W/(m·K)],
area=0.1 [m²]
) -> f32 [W]
# Convection (Newton's law of cooling)
q = transport.convection(
temp_surface=350.0 [K],
temp_bulk=300.0 [K],
h=100.0 [W/(m²·K)], # Heat transfer coefficient
area=1.0 [m²]
) -> f32 [W]
# Radiation (Stefan-Boltzmann)
q = transport.radiation(
temp_surface=500.0 [K],
temp_ambient=300.0 [K],
emissivity=0.8,
area=1.0 [m²]
) -> f32 [W]
# Heat transfer coefficient correlations
h = transport.nusselt_correlation(
Re=10000.0, # Reynolds number
Pr=0.7, # Prandtl number
geometry="pipe"
) -> f32 [W/(m²·K)]
# Fickian diffusion
flux = transport.fickian_diffusion(
conc_gradient=0.1 [mol/(m³·m)],
diffusivity=1e-9 [m²/s],
area=0.01 [m²]
) -> f32 [mol/s]
# Knudsen diffusion (in pores)
D_knudsen = transport.knudsen_diffusion(
pore_diameter=10.0 [nm],
temp=300.0 [K],
molecular_weight=28.0 [g/mol]
) -> f32 [m²/s]
# Convective mass transfer
flux = transport.convective_mass_transfer(
conc_bulk=1.0 [mol/m³],
conc_surface=0.5 [mol/m³],
k_mass=1e-3 [m/s],
area=0.1 [m²]
) -> f32 [mol/s]
# Mass transfer coefficient
k_L = transport.sherwood_correlation(
Re=10000.0,
Sc=1000.0, # Schmidt number
geometry="sphere"
) -> f32 [m/s]
# Effective diffusivity
D_eff = transport.effective_diffusivity(
D_molecular=1e-9 [m²/s],
porosity=0.4,
tortuosity=2.0
) -> f32 [m²/s]
# Darcy flow
velocity = transport.darcy_flow(
pressure_gradient=1000.0 [Pa/m],
permeability=1e-12 [m²],
viscosity=1e-3 [Pa·s]
) -> f32 [m/s]
Vapor-liquid equilibrium, gas-liquid reactions, multiphase flow.
# Flash calculation
y_vapor, x_liquid = multiphase.vle_flash(
feed_composition=[0.5, 0.5],
temp=350.0 [K],
pressure=1e5 [Pa],
thermo_model="peng_robinson" # "ideal", "peng_robinson", "srk"
) -> (Vec<f32>, Vec<f32>)
# Bubble point
temp_bubble = multiphase.bubble_point(
liquid_composition=[0.5, 0.5],
pressure=1e5 [Pa]
) -> f32 [K]
# Dew point
temp_dew = multiphase.dew_point(
vapor_composition=[0.5, 0.5],
pressure=1e5 [Pa]
) -> f32 [K]
# Volumetric mass transfer coefficient
k_L_a = multiphase.volumetric_mass_transfer(
bubble_diameter=1.0 [mm],
gas_holdup=0.1,
diffusivity=1e-9 [m²/s]
) -> f32 [1/s]
# Gas-liquid reaction rate
rate = multiphase.gas_liquid_reaction(
conc_liquid=0.5 [mol/L],
pressure_gas=1e5 [Pa],
k_L_a=10.0 [1/s],
k_rxn=1e3 [1/s],
henry_constant=1e-3 [mol/(L·Pa)]
) -> f32 [mol/(L·s)]
Equations of state, activity coefficients, thermodynamic properties.
# Peng-Robinson EOS
Z = thermo.peng_robinson(
temp=350.0 [K],
pressure=10e5 [Pa],
critical_temp=647.0 [K],
critical_pressure=220.6e5 [Pa],
acentric_factor=0.344
) -> f32 # Compressibility factor
# Activity coefficient (UNIFAC, NRTL, Wilson)
gamma = thermo.activity_coefficient(
composition=[0.5, 0.5],
temp=350.0 [K],
model="nrtl",
parameters=nrtl_params
) -> Vec<f32>
# Heat capacity
Cp = thermo.heat_capacity(
species="water",
temp=350.0 [K]
) -> f32 [J/(mol·K)]
# Enthalpy of reaction
delta_H = thermo.enthalpy_of_reaction(
reactants={"A": 1.0, "B": 1.0},
products={"C": 1.0},
temp=298.15 [K]
) -> f32 [J/mol]
Heterogeneous catalysis, surface reactions, catalyst characterization.
# Langmuir-Hinshelwood mechanism
rate = catalysis.langmuir_hinshelwood(
coverage_A=0.5,
coverage_B=0.3,
k_surface=1e6 [1/s]
) -> f32 [mol/(m²·s)]
# Eley-Rideal mechanism
rate = catalysis.eley_rideal(
coverage_A=0.5,
pressure_B=1e5 [Pa],
k_surface=1e6 [mol/(m²·s·Pa)]
) -> f32 [mol/(m²·s)]
# Surface coverage evolution
coverage_new = catalysis.surface_coverage_step(
coverage,
r_adsorption=0.1 [1/s],
r_desorption=0.05 [1/s],
r_reaction=0.02 [1/s],
dt=0.1 [s]
)
# BET surface area
surface_area = catalysis.bet_surface_area(
adsorption_isotherm=[...],
adsorbate="N2"
) -> f32 [m²/g]
# Pore size distribution
psd = catalysis.pore_size_distribution(
adsorption_isotherm=[...],
method="bjh" # "bjh", "dft"
) -> Field1D<f32>
Batteries, fuel cells, electrolysis, corrosion.
# Butler-Volmer kinetics
i = electrochem.butler_volmer(
overpotential=0.1 [V],
i0=1e-3 [A/m²], # Exchange current density
alpha=0.5, # Transfer coefficient
n=2 # Electrons transferred
) -> f32 [A/m²]
# Nernst equation
E = electrochem.nernst(
E_standard=0.34 [V],
conc_oxidized=0.1 [mol/L],
conc_reduced=0.01 [mol/L],
n=2,
temp=298.15 [K]
) -> f32 [V]
# Battery simulation
voltage = electrochem.battery_discharge(
capacity=3000.0 [mAh],
current=1.0 [A],
time=3600.0 [s],
model="equivalent_circuit"
) -> f32 [V]
Chemistry leverages existing field operators for spatially-resolved simulations:
use field, kinetics
# Reaction-diffusion in 3D reactor
@state conc : Field3D<Vec<f32 [mol/m³]>> = initial_conc()
@state temp : Field3D<f32 [K]> = uniform(350.0)
flow(dt=0.1 [s]) {
# Diffusion
conc = diffuse(conc, rate=1e-9 [m²/s], dt)
# Reaction
reaction_rate = kinetics.arrhenius_field(temp, conc)
conc = conc - reaction_rate * dt
# Convection (if velocity field exists)
conc = advect(conc, velocity_field, dt)
}
Chemistry uses existing optimization operators for design:
use optimization, molecular
# Optimize catalyst structure
best_molecule = optimize.ga(
objective=catalyst_activity,
bounds=molecular_parameter_space,
population=100,
generations=50
)
# Multi-objective: yield vs. cost vs. toxicity
pareto_front = optimize.nsga2(
objectives=[
maximize(yield),
minimize(cost),
minimize(toxicity)
],
bounds=reaction_conditions
)
Chemistry adds ML operators for surrogate models and generative design:
# Train surrogate model
surrogate = ml.train_surrogate(
expensive_dft_calculation,
input_bounds=molecule_space,
n_samples=1000
)
# Generative design
generator = ml.train_generative(
molecules=training_set,
architecture="vae"
)
new_molecules = generator.sample(n=100)
- ✅ Molecule loading (SMILES, PDB, XYZ)
- ✅ Basic properties (MW, formula)
- ✅ Force fields (AMBER, UFF)
- ✅ MD integrators (Verlet, Langevin)
- Dependency: RDKit or OpenBabel
- ✅ Arrhenius kinetics
- ✅ ODE integration
- ✅ Ideal reactors (batch, CSTR, PFR)
- Dependency: SciPy for ODE solvers
- ✅ Heat transfer operators
- ✅ Mass diffusion
- ✅ Correlations (Nusselt, Sherwood)
- Integration with existing field ops
- ✅ DFT wrappers (ORCA, Psi4)
- ✅ ML PES (SchNet, DimeNet)
- Dependency: External QM codes
- ✅ VLE flash calculations
- ✅ Thermodynamic models (PR, SRK)
- Dependency: Thermo library (CoolProp or custom)
- ✅ Surface kinetics
- ✅ Butler-Volmer
- Advanced applications
- RDKit — Molecular informatics, SMILES parsing
- OpenBabel — Format conversion
- Cantera — Reaction kinetics (can wrap or reimplement)
- CoolProp — Thermodynamic properties
- ORCA/Psi4/Q-Chem — Quantum chemistry (external executables)
- SchNet/TorchMD — ML potential energy surfaces
- Field operators (diffusion, advection, Laplacian)
- Agent operators (particle systems)
- Optimization domain (GA, CMA-ES, Bayesian)
- ML domain (surrogate models, training)
- Visualization domain (molecules, fields, trajectories)
use molecular, qchem, ml, optimization
# Define objective: catalyst activity
fn catalyst_activity(molecule: Molecule) -> f32 {
# Compute binding energy (expensive DFT)
binding_energy = qchem.dft_energy(molecule, method="B3LYP", basis="6-31G*")
# Lower binding energy = better catalyst
return -binding_energy
}
# Train surrogate model
surrogate = ml.train_surrogate(
catalyst_activity,
input_space=molecular_descriptors,
n_samples=500
)
# Optimize with genetic algorithm
best_catalyst = optimize.ga(
objective=surrogate,
population=100,
generations=50
)
# Validate with real DFT
final_activity = catalyst_activity(best_catalyst)
use field, kinetics, transport, visual
@state conc : Field3D<f32 [mol/m³]> = uniform(1.0)
@state temp : Field3D<f32 [K]> = uniform(300.0)
@state velocity : Field3D<Vec3<f32 [m/s]>> = inlet_profile()
flow(dt=0.01 [s], steps=10000) {
# Fluid flow
velocity = navier_stokes_step(velocity, temp, dt)
# Heat transfer
reaction_heat = kinetics.reaction_heat(conc, temp)
temp = heat_equation_step(temp, velocity, source=reaction_heat, dt)
# Mass transfer + reaction
conc = diffuse(conc, rate=1e-9 [m²/s], dt)
conc = advect(conc, velocity, dt)
reaction_rate = kinetics.arrhenius_field(temp, conc)
conc = conc - reaction_rate * dt
# Visualize
output visual.volume_render(temp, palette="fire")
output visual.isosurface(conc, level=0.5)
}
- Specification: COMPLETE
- Implementation: NOT STARTED
- Dependencies: RDKit, SciPy, optional external QM codes
- Target: Morphogen v0.9.0+
- ADR-006: Chemistry Domain
- use-cases/chemistry-unified-framework.md
- architecture/domain-architecture.md
Last Updated: 2025-11-15 Version: 1.0