From 58579e69f1edb731f6afd82e592206d32a6a1b67 Mon Sep 17 00:00:00 2001 From: Tony Kan Date: Mon, 20 Jul 2026 07:45:00 +0800 Subject: [PATCH 1/4] feat(forge): Add MPSim compatibility adapter for `H___A` typing and neutral termini --- src/bin/dforge/config/forge.rs | 2 + src/forge/charge/hybrid.rs | 79 ++++-- src/forge/charge/mod.rs | 20 +- src/forge/config/mod.rs | 10 + src/forge/config/mpsim.rs | 73 +++++ src/forge/mod.rs | 69 ++++- src/forge/mpsim.rs | 499 +++++++++++++++++++++++++++++++++ src/lib.rs | 2 +- 8 files changed, 729 insertions(+), 25 deletions(-) create mode 100644 src/forge/config/mpsim.rs create mode 100644 src/forge/mpsim.rs diff --git a/src/bin/dforge/config/forge.rs b/src/bin/dforge/config/forge.rs index b84a320..d236b4a 100644 --- a/src/bin/dforge/config/forge.rs +++ b/src/bin/dforge/config/forge.rs @@ -24,6 +24,7 @@ pub fn build_bio_forge_config( bond_potential: potential.bond_potential.into(), angle_potential: potential.angle_potential.into(), vdw_potential: potential.vdw_potential.into(), + mpsim: None, }) } @@ -42,6 +43,7 @@ pub fn build_chem_forge_config( bond_potential: potential.bond_potential.into(), angle_potential: potential.angle_potential.into(), vdw_potential: potential.vdw_potential.into(), + mpsim: None, }) } diff --git a/src/forge/charge/hybrid.rs b/src/forge/charge/hybrid.rs index ab0b3b7..bf75266 100644 --- a/src/forge/charge/hybrid.rs +++ b/src/forge/charge/hybrid.rs @@ -26,6 +26,8 @@ const C_TERMINAL_PKA: f64 = 3.1; /// /// * `system` — Mutable reference to the intermediate system /// * `config` — Hybrid charge configuration +/// * `neutral_termini` — When `true` (MPSim mode), protein N/C termini use the +/// neutral charge sets (`NH₂` / `COOH`) regardless of pH /// /// # Errors /// @@ -35,6 +37,7 @@ const C_TERMINAL_PKA: f64 = 3.1; pub fn assign_hybrid_charges( system: &mut IntermediateSystem, config: &HybridConfig, + neutral_termini: bool, ) -> Result<(), Error> { if !system.has_bio_metadata() { return Err(Error::MissingBioMetadata); @@ -44,7 +47,14 @@ pub fn assign_hybrid_charges( let metadata = system.bio_metadata.as_ref().unwrap().clone(); let classification = classify_atoms(&metadata); - assign_fixed_charges(system, &metadata, config, ph, &classification)?; + assign_fixed_charges( + system, + &metadata, + config, + ph, + &classification, + neutral_termini, + )?; let ligand_groups = identify_ligand_groups(&metadata, &classification); if !ligand_groups.is_empty() { @@ -130,10 +140,11 @@ fn assign_fixed_charges( config: &HybridConfig, ph: f64, classification: &[AtomClass], + neutral_termini: bool, ) -> Result<(), Error> { for (idx, (&class, info)) in classification.iter().zip(&metadata.atom_info).enumerate() { let charge = match class { - AtomClass::Protein => lookup_protein_charge(config, info, ph)?, + AtomClass::Protein => lookup_protein_charge(config, info, ph, neutral_termini)?, AtomClass::NucleicAcid => lookup_nucleic_charge(config, info)?, AtomClass::Water => lookup_water_charge(config, info)?, AtomClass::Ion => lookup_ion_charge(info)?, @@ -145,20 +156,25 @@ fn assign_fixed_charges( } /// Maps our ResiduePosition to ffcharge::Position for proteins. -fn map_residue_position(info: &AtomResidueInfo, ph: f64) -> FfPosition { +/// +/// When `neutral_termini` is `true` (MPSim mode), protein N/C termini map to +/// their neutral charge sets (`NH₂` / `COOH`) regardless of pH. This must be +/// paired with the matching neutral terminal topology (see +/// [`crate::forge::mpsim`]) so every terminal residue stays net-neutral. +fn map_residue_position(info: &AtomResidueInfo, ph: f64, neutral_termini: bool) -> FfPosition { use crate::model::metadata::ResiduePosition; match info.position { ResiduePosition::NTerminal => { - if ph < N_TERMINAL_PKA { - FfPosition::NTerminal // Protonated NH3+ - } else { + if neutral_termini || ph >= N_TERMINAL_PKA { FfPosition::NTerminalDeprotonated // Neutral NH2 + } else { + FfPosition::NTerminal // Protonated NH3+ } } ResiduePosition::CTerminal => { - if ph < C_TERMINAL_PKA { - FfPosition::CTerminalProtonated // Protonated COOH + if neutral_termini || ph < C_TERMINAL_PKA { + FfPosition::CTerminalProtonated // Protonated COOH (neutral) } else { FfPosition::CTerminal // Deprotonated COO- } @@ -174,8 +190,9 @@ fn lookup_protein_charge( config: &HybridConfig, info: &AtomResidueInfo, ph: f64, + neutral_termini: bool, ) -> Result { - let position = map_residue_position(info, ph); + let position = map_residue_position(info, ph, neutral_termini); config .protein_scheme @@ -196,7 +213,8 @@ fn lookup_protein_charge( /// Looks up nucleic acid charge from ffcharge. fn lookup_nucleic_charge(config: &HybridConfig, info: &AtomResidueInfo) -> Result { - let position = map_residue_position(info, 7.0); // Assuming pH 7.0 for nucleic acids (no effect) + // Nucleic acids only use 5'/3' positions, which neutral_termini does not affect. + let position = map_residue_position(info, 7.0, false); config .nucleic_scheme @@ -479,7 +497,7 @@ mod tests { let info = AtomResidueInfo::builder("N", "ALA", 1, "A") .position(ResiduePosition::NTerminal) .build(); - let pos = map_residue_position(&info, 7.0); + let pos = map_residue_position(&info, 7.0, false); assert_eq!(pos, FfPosition::NTerminal); } @@ -488,7 +506,7 @@ mod tests { let info = AtomResidueInfo::builder("N", "ALA", 1, "A") .position(ResiduePosition::NTerminal) .build(); - let pos = map_residue_position(&info, 9.0); + let pos = map_residue_position(&info, 9.0, false); assert_eq!(pos, FfPosition::NTerminalDeprotonated); } @@ -497,7 +515,7 @@ mod tests { let info = AtomResidueInfo::builder("C", "ALA", 1, "A") .position(ResiduePosition::CTerminal) .build(); - let pos = map_residue_position(&info, 7.0); + let pos = map_residue_position(&info, 7.0, false); assert_eq!(pos, FfPosition::CTerminal); } @@ -506,16 +524,45 @@ mod tests { let info = AtomResidueInfo::builder("C", "ALA", 1, "A") .position(ResiduePosition::CTerminal) .build(); - let pos = map_residue_position(&info, 2.0); + let pos = map_residue_position(&info, 2.0, false); assert_eq!(pos, FfPosition::CTerminalProtonated); } + #[test] + fn map_protein_termini_neutral_regardless_of_ph_in_mpsim_mode() { + let n_info = AtomResidueInfo::builder("N", "ALA", 1, "A") + .position(ResiduePosition::NTerminal) + .build(); + let c_info = AtomResidueInfo::builder("C", "ALA", 1, "A") + .position(ResiduePosition::CTerminal) + .build(); + + // At physiological pH, neutral_termini forces the neutral charge sets. + assert_eq!( + map_residue_position(&n_info, 7.0, true), + FfPosition::NTerminalDeprotonated + ); + assert_eq!( + map_residue_position(&c_info, 7.0, true), + FfPosition::CTerminalProtonated + ); + // Even at extreme pH values the neutral state is kept. + assert_eq!( + map_residue_position(&n_info, 1.0, true), + FfPosition::NTerminalDeprotonated + ); + assert_eq!( + map_residue_position(&c_info, 13.0, true), + FfPosition::CTerminalProtonated + ); + } + #[test] fn map_nucleic_position_five_prime() { let info = AtomResidueInfo::builder("P", "DA", 1, "B") .position(ResiduePosition::FivePrime) .build(); - let pos = map_residue_position(&info, 7.0); + let pos = map_residue_position(&info, 7.0, false); assert_eq!(pos, FfPosition::FivePrime); } @@ -524,7 +571,7 @@ mod tests { let info = AtomResidueInfo::builder("O3'", "DA", 1, "B") .position(ResiduePosition::ThreePrime) .build(); - let pos = map_residue_position(&info, 7.0); + let pos = map_residue_position(&info, 7.0, false); assert_eq!(pos, FfPosition::ThreePrime); } diff --git a/src/forge/charge/mod.rs b/src/forge/charge/mod.rs index 5b21a32..e61d17c 100644 --- a/src/forge/charge/mod.rs +++ b/src/forge/charge/mod.rs @@ -22,17 +22,27 @@ use super::intermediate::IntermediateSystem; /// * `system` — Mutable reference to the intermediate system /// * `method` — Charge calculation method to use /// +/// The `neutral_termini` flag forces MPSim-compatible neutral N/C protein +/// termini in the hybrid method (see [`hybrid::assign_hybrid_charges`]); it has +/// no effect on the `None` or `Qeq` methods. +/// /// # Errors /// /// Returns [`Error`] if: /// - QEq solver fails to converge ([`Error::ChargeCalculation`]) /// - Hybrid method is used without biological metadata ([`Error::MissingBioMetadata`]) /// - Classical charge lookup fails ([`Error::HybridChargeAssignment`]) -pub fn assign_charges(system: &mut IntermediateSystem, method: &ChargeMethod) -> Result<(), Error> { +pub fn assign_charges( + system: &mut IntermediateSystem, + method: &ChargeMethod, + neutral_termini: bool, +) -> Result<(), Error> { match method { ChargeMethod::None => Ok(()), ChargeMethod::Qeq(config) => qeq::assign_qeq_charges(system, config), - ChargeMethod::Hybrid(config) => hybrid::assign_hybrid_charges(system, config), + ChargeMethod::Hybrid(config) => { + hybrid::assign_hybrid_charges(system, config, neutral_termini) + } } } @@ -66,7 +76,7 @@ mod tests { let water = make_water(); let mut int = IntermediateSystem::from_system(&water).unwrap(); - assign_charges(&mut int, &ChargeMethod::None).unwrap(); + assign_charges(&mut int, &ChargeMethod::None, false).unwrap(); for atom in &int.atoms { assert_eq!(atom.charge, 0.0); @@ -79,7 +89,7 @@ mod tests { let mut int = IntermediateSystem::from_system(&water).unwrap(); let qeq_config = QeqConfig::default(); - assign_charges(&mut int, &ChargeMethod::Qeq(qeq_config)).unwrap(); + assign_charges(&mut int, &ChargeMethod::Qeq(qeq_config), false).unwrap(); assert!(int.atoms[0].charge < 0.0); assert!(int.atoms[1].charge > 0.0); @@ -98,7 +108,7 @@ mod tests { total_charge: -1.0, ..Default::default() }; - assign_charges(&mut int, &ChargeMethod::Qeq(qeq_config)).unwrap(); + assign_charges(&mut int, &ChargeMethod::Qeq(qeq_config), false).unwrap(); let total: f64 = int.atoms.iter().map(|a| a.charge).sum(); assert!((total + 1.0).abs() < 1e-9); diff --git a/src/forge/config/mod.rs b/src/forge/config/mod.rs index 77e232f..3ee3243 100644 --- a/src/forge/config/mod.rs +++ b/src/forge/config/mod.rs @@ -16,6 +16,7 @@ //! - [`VdwPotentialType`] — Van der Waals potential selection mod charge; +mod mpsim; mod potential; pub use charge::{ @@ -23,6 +24,7 @@ pub use charge::{ LigandQeqMethod, NucleicScheme, ProteinScheme, QeqConfig, ResidueSelector, SolverOptions, WaterScheme, }; +pub use mpsim::MpsimConfig; pub use potential::{AnglePotentialType, BondPotentialType, VdwPotentialType}; /// Main configuration for DREIDING force field parameterization. @@ -70,6 +72,13 @@ pub struct ForgeConfig { /// Type of van der Waals non-bonded potential to generate. pub vdw_potential: VdwPotentialType, + + /// Optional MPSim (legacy EM/MM engine) compatibility adapter. + /// + /// When `Some`, the pipeline applies MPSim's DREIDING conventions: + /// renaming `H_HB` hydrogens to `H___A` and forcing neutral N/C protein + /// termini. See [`MpsimConfig`]. Defaults to `None` (disabled). + pub mpsim: Option, } impl Default for ForgeConfig { @@ -81,6 +90,7 @@ impl Default for ForgeConfig { bond_potential: BondPotentialType::Harmonic, angle_potential: AnglePotentialType::Cosine, vdw_potential: VdwPotentialType::LennardJones, + mpsim: None, } } } diff --git a/src/forge/config/mpsim.rs b/src/forge/config/mpsim.rs new file mode 100644 index 0000000..2ec9b7b --- /dev/null +++ b/src/forge/config/mpsim.rs @@ -0,0 +1,73 @@ +//! MPSim compatibility configuration. +//! +//! MPSim is a legacy EM/MM engine used within the group. Its DREIDING +//! implementation differs from the modern convention in two ways that this +//! adapter reconciles: +//! +//! 1. **Hydrogen-bond hydrogen naming** — MPSim expects the force-field type +//! `H___A` where modern DREIDING assigns `H_HB` to polar (donor) hydrogens. +//! 2. **Protein chain termini** — MPSim inputs use the neutral, uncharged +//! protonation state for both the N-terminus (–NH₂) and the C-terminus +//! (–COOH), regardless of the modeled pH. +//! +//! Enabling [`MpsimConfig`] on [`ForgeConfig`](super::ForgeConfig) applies +//! these conventions as a self-contained post-/pre-processing layer without +//! altering the underlying DREIDING parameterization. + +/// Compatibility adapter for the legacy MPSim EM/MM engine. +/// +/// When attached to [`ForgeConfig`](super::ForgeConfig) via +/// [`mpsim`](super::ForgeConfig::mpsim), the [`forge`](crate::forge) pipeline +/// applies MPSim's DREIDING conventions. Both behaviors are enabled by +/// default and can be toggled independently. +/// +/// # Examples +/// +/// ``` +/// use dreid_forge::{ForgeConfig, MpsimConfig}; +/// +/// // Enable the full MPSim adapter with default behavior. +/// let config = ForgeConfig { +/// mpsim: Some(MpsimConfig::default()), +/// ..Default::default() +/// }; +/// ``` +#[derive(Debug, Clone)] +pub struct MpsimConfig { + /// Rename hydrogen-bond hydrogens from DREIDING `H_HB` to MPSim `H___A`. + /// + /// The rename is applied only to the emitted force-field type names; it + /// does not affect parameter lookup, hydrogen-bond term generation, or + /// atom typing internally. + pub rename_hb_hydrogen: bool, + + /// Force neutral, uncharged N/C protein termini. + /// + /// When `true`, the C-terminus is normalized to the protonated carboxyl + /// (`–COOH`, adding `HOXT`) and the N-terminus to the neutral amine + /// (`–NH₂`, removing the extra `H3`), independent of pH. The hybrid + /// charge assignment likewise uses the neutral terminal charge sets so + /// each terminal residue carries no net formal charge. + pub neutral_termini: bool, +} + +impl Default for MpsimConfig { + fn default() -> Self { + Self { + rename_hb_hydrogen: true, + neutral_termini: true, + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn default_enables_all_behaviors() { + let config = MpsimConfig::default(); + assert!(config.rename_hb_hydrogen); + assert!(config.neutral_termini); + } +} diff --git a/src/forge/mod.rs b/src/forge/mod.rs index 6ebf0da..03b79d5 100644 --- a/src/forge/mod.rs +++ b/src/forge/mod.rs @@ -9,13 +9,14 @@ mod charge; mod config; mod error; mod intermediate; +mod mpsim; mod paramgen; mod params; mod typer; pub use config::{ AnglePotentialType, BasisType, BondPotentialType, ChargeMethod, DampingStrategy, - EmbeddedQeqConfig, ForgeConfig, HybridConfig, LigandChargeConfig, LigandQeqMethod, + EmbeddedQeqConfig, ForgeConfig, HybridConfig, LigandChargeConfig, LigandQeqMethod, MpsimConfig, NucleicScheme, ProteinScheme, QeqConfig, ResidueSelector, SolverOptions, VdwPotentialType, WaterScheme, }; @@ -69,13 +70,32 @@ use crate::model::topology::ForgedSystem; pub fn forge(system: &System, config: &ForgeConfig) -> Result { let ff_params = params::load_parameters(config.params.as_deref())?; + let neutral_termini = config.mpsim.as_ref().is_some_and(|m| m.neutral_termini); + + // MPSim neutral-termini requires the terminal topology (NH2 / COOH) to match + // the neutral charge sets, so normalize a working copy up front. The + // normalized structure is also what the ForgedSystem carries into output. + let owned_system; + let system: &System = if neutral_termini { + let mut normalized = system.clone(); + mpsim::normalize_terminal_hydrogens(&mut normalized); + owned_system = normalized; + &owned_system + } else { + system + }; + let mut intermediate = intermediate::IntermediateSystem::from_system(system)?; typer::assign_atom_types(&mut intermediate, config.rules.as_deref())?; - charge::assign_charges(&mut intermediate, &config.charge_method)?; + charge::assign_charges(&mut intermediate, &config.charge_method, neutral_termini)?; + + let mut forged = paramgen::generate_parameters(system, &intermediate, &ff_params, config)?; - let forged = paramgen::generate_parameters(system, &intermediate, &ff_params, config)?; + if config.mpsim.as_ref().is_some_and(|m| m.rename_hb_hydrogen) { + mpsim::rename_hb_hydrogens(&mut forged.atom_types); + } Ok(forged) } @@ -259,6 +279,49 @@ mod tests { assert!(matches!(result, Err(Error::EmptySystem))); } + #[test] + fn mpsim_renames_hb_hydrogen_to_h_a() { + use crate::forge::config::MpsimConfig; + let water = make_water(); + + // Baseline: default DREIDING emits H_HB for the polar hydrogens. + let baseline = forge(&water, &ForgeConfig::default()).unwrap(); + assert!(baseline.atom_types.contains(&"H_HB".to_string())); + assert!(!baseline.atom_types.contains(&"H___A".to_string())); + + // MPSim mode renames H_HB -> H___A in the emitted type table. + let config = ForgeConfig { + mpsim: Some(MpsimConfig::default()), + ..Default::default() + }; + let forged = forge(&water, &config).unwrap(); + assert!(forged.atom_types.contains(&"H___A".to_string())); + assert!(!forged.atom_types.contains(&"H_HB".to_string())); + + // Type indexing stays consistent and H-bond terms still generate. + assert_eq!(forged.atom_types.len(), baseline.atom_types.len()); + assert_eq!( + forged.potentials.h_bonds.len(), + baseline.potentials.h_bonds.len() + ); + } + + #[test] + fn mpsim_rename_can_be_disabled() { + use crate::forge::config::MpsimConfig; + let water = make_water(); + let config = ForgeConfig { + mpsim: Some(MpsimConfig { + rename_hb_hydrogen: false, + neutral_termini: true, + }), + ..Default::default() + }; + let forged = forge(&water, &config).unwrap(); + assert!(forged.atom_types.contains(&"H_HB".to_string())); + assert!(!forged.atom_types.contains(&"H___A".to_string())); + } + #[test] fn errors_on_invalid_custom_params() { let water = make_water(); diff --git a/src/forge/mpsim.rs b/src/forge/mpsim.rs new file mode 100644 index 0000000..77c823c --- /dev/null +++ b/src/forge/mpsim.rs @@ -0,0 +1,499 @@ +//! MPSim compatibility adapter. +//! +//! Implements the two structural/naming conventions required by the legacy +//! MPSim EM/MM engine, driven by [`MpsimConfig`](super::config::MpsimConfig): +//! +//! - [`rename_hb_hydrogens`] rewrites emitted `H_HB` force-field type names to +//! MPSim's `H___A`. +//! - [`normalize_terminal_hydrogens`] rebuilds protein chain termini into their +//! neutral protonation state — stripping the extra N-terminal hydrogen +//! (`NH₃⁺ → NH₂`) and building the C-terminal carboxyl proton +//! (`COO⁻ → COOH`) — so that the neutral terminal charge sets applied during +//! hybrid charge assignment stay self-consistent. +//! +//! The neutral-terminal *charge* selection itself lives in the hybrid charge +//! module; this module owns the matching topology transform. Both are gated on +//! the same [`MpsimConfig::neutral_termini`] flag. + +use crate::model::atom::Atom; +use crate::model::metadata::{AtomResidueInfo, ResidueCategory, ResiduePosition, StandardResidue}; +use crate::model::system::{Bond, System}; +use crate::model::types::{BondOrder, Element}; +use std::collections::{HashMap, HashSet}; + +/// DREIDING force-field type for polar (hydrogen-bond) hydrogens. +pub const HB_HYDROGEN_TYPE: &str = "H_HB"; +/// MPSim force-field type replacing [`HB_HYDROGEN_TYPE`]. +pub const MPSIM_HB_HYDROGEN_TYPE: &str = "H___A"; + +/// PDB atom name of the extra N-terminal ammonium hydrogen (`NH₃⁺`). +const N_TERM_EXTRA_HYDROGEN: &str = "H3"; +/// PDB atom name of the C-terminal carboxyl hydrogen (`COOH`). +const C_TERM_HYDROGEN: &str = "HOXT"; + +/// O–H bond length for the carboxyl proton (Å); mirrors bio-forge. +const COOH_BOND_LENGTH: f64 = 0.97; +/// sp³ tetrahedral bond angle for hydroxyl placement (degrees). +const SP3_ANGLE_DEG: f64 = 109.5; +/// Dihedral offset used when placing the carboxyl hydrogen (degrees). +const HOXT_DIHEDRAL_DEG: f64 = 60.0; + +/// Rewrites `H_HB` force-field type names to MPSim's `H___A`. +/// +/// Operates only on the emitted type-name table. Because potentials reference +/// atom types by index and the rename is one-to-one, all downstream references +/// (per-atom type index, hydrogen-bond terms, BGF output) remain valid. +pub fn rename_hb_hydrogens(atom_types: &mut [String]) { + for atom_type in atom_types.iter_mut() { + if atom_type == HB_HYDROGEN_TYPE { + *atom_type = MPSIM_HB_HYDROGEN_TYPE.to_string(); + } + } +} + +/// Normalizes protein chain termini to their neutral protonation state. +/// +/// For every standard protein residue at a terminus: +/// +/// - **N-terminal** — removes the extra `H3` hydrogen, leaving the neutral +/// amine (`–NH₂`, or `–NH` for proline). +/// - **C-terminal** — adds the `HOXT` carboxyl hydrogen (bonded to `OXT`), +/// forming the neutral acid (`–COOH`), unless it is already present. +/// +/// Atom indices, bonds, and biological metadata are kept consistent. Systems +/// without biological metadata (e.g. small molecules) are left untouched, as +/// are termini that already carry the neutral hydrogen set. +pub fn normalize_terminal_hydrogens(system: &mut System) { + if system.bio_metadata.is_none() { + return; + } + + let removals = plan_nterm_removals(system); + let additions = plan_cterm_additions(system); + + if removals.is_empty() && additions.is_empty() { + return; + } + + let old_len = system.atoms.len(); + + // Rebuild the atom / metadata arrays, dropping removed atoms and recording + // the old-index -> new-index mapping (None for removed atoms). + let mut remap = vec![None; old_len]; + { + let mut new_atoms = Vec::with_capacity(old_len); + let mut new_info = Vec::with_capacity(old_len); + let info = &system.bio_metadata.as_ref().unwrap().atom_info; + for i in 0..old_len { + if removals.contains(&i) { + continue; + } + remap[i] = Some(new_atoms.len()); + new_atoms.push(system.atoms[i].clone()); + new_info.push(info[i].clone()); + } + system.atoms = new_atoms; + system.bio_metadata.as_mut().unwrap().atom_info = new_info; + } + + // Remap surviving bonds, dropping any that touched a removed atom. + let mut new_bonds = Vec::with_capacity(system.bonds.len()); + for bond in &system.bonds { + if let (Some(i), Some(j)) = (remap[bond.i], remap[bond.j]) { + new_bonds.push(Bond::new(i, j, bond.order)); + } + } + system.bonds = new_bonds; + + // Append the newly built C-terminal carboxyl hydrogens. + for add in additions { + let Some(partner) = remap[add.partner] else { + continue; + }; + let new_idx = system.atoms.len(); + system.atoms.push(Atom::new(Element::H, add.position)); + system + .bio_metadata + .as_mut() + .unwrap() + .atom_info + .push(add.info); + system + .bonds + .push(Bond::new(partner, new_idx, BondOrder::Single)); + } +} + +/// Collects atom indices of extra N-terminal hydrogens to remove. +fn plan_nterm_removals(system: &System) -> HashSet { + let meta = system.bio_metadata.as_ref().unwrap(); + let mut removals = HashSet::new(); + + for (idx, info) in meta.atom_info.iter().enumerate() { + if info.position == ResiduePosition::NTerminal + && info.category == ResidueCategory::Standard + && is_protein_residue(info.standard_name) + && info.atom_name == N_TERM_EXTRA_HYDROGEN + { + removals.insert(idx); + } + } + + removals +} + +/// A carboxyl hydrogen to be appended, referencing its `OXT` partner by the +/// atom index in the *original* (pre-removal) system. +struct PendingHydrogen { + position: [f64; 3], + info: AtomResidueInfo, + partner: usize, +} + +/// Plans C-terminal `HOXT` additions for all neutral-capped protein termini. +fn plan_cterm_additions(system: &System) -> Vec { + let meta = system.bio_metadata.as_ref().unwrap(); + + // Group C-terminal protein residue atom indices by residue identity. + let mut groups: HashMap<(&str, i32, Option), Vec> = HashMap::new(); + for (idx, info) in meta.atom_info.iter().enumerate() { + if info.position == ResiduePosition::CTerminal + && info.category == ResidueCategory::Standard + && is_protein_residue(info.standard_name) + { + let key = (info.chain_id.as_str(), info.residue_id, info.insertion_code); + groups.entry(key).or_default().push(idx); + } + } + + let mut additions = Vec::new(); + for atoms in groups.values() { + let find = |name: &str| { + atoms + .iter() + .copied() + .find(|&i| meta.atom_info[i].atom_name == name) + }; + + // Skip if already protonated. + if find(C_TERM_HYDROGEN).is_some() { + continue; + } + + let (Some(c_idx), Some(oxt_idx)) = (find("C"), find("OXT")) else { + continue; + }; + let Some(ref_idx) = find("CA").or_else(|| find("O")) else { + continue; + }; + + let position = place_hydroxyl_hydrogen( + system.atoms[oxt_idx].position, + system.atoms[c_idx].position, + system.atoms[ref_idx].position, + ); + + let template = &meta.atom_info[oxt_idx]; + let info = AtomResidueInfo::builder( + C_TERM_HYDROGEN, + template.residue_name.clone(), + template.residue_id, + template.chain_id.clone(), + ) + .insertion_code_opt(template.insertion_code) + .standard_name(template.standard_name) + .category(template.category) + .position(template.position) + .build(); + + additions.push(PendingHydrogen { + position, + info, + partner: oxt_idx, + }); + } + + // Deterministic append order regardless of HashMap iteration order. + additions.sort_by_key(|a| a.partner); + additions +} + +/// Returns `true` for the 20 standard amino-acid residues. +fn is_protein_residue(name: Option) -> bool { + use StandardResidue::*; + matches!( + name, + Some( + ALA | ARG + | ASN + | ASP + | CYS + | GLN + | GLU + | GLY + | HIS + | ILE + | LEU + | LYS + | MET + | PHE + | PRO + | SER + | THR + | TRP + | TYR + | VAL + ) + ) +} + +/// Places a hydroxyl hydrogen on `oxt` using sp³ tetrahedral geometry, mirroring +/// bio-forge's carboxyl-proton construction (bond length, angle, dihedral). +fn place_hydroxyl_hydrogen(oxt: [f64; 3], attached: [f64; 3], reference: [f64; 3]) -> [f64; 3] { + let (x, y, z) = build_sp3_frame(oxt, attached, reference); + + let theta = SP3_ANGLE_DEG.to_radians(); + let phi = HOXT_DIHEDRAL_DEG.to_radians(); + let (sin_theta, cos_theta) = (theta.sin(), theta.cos()); + + let local = [sin_theta * phi.cos(), sin_theta * phi.sin(), -cos_theta]; + let dir = add( + add(scale(x, local[0]), scale(y, local[1])), + scale(z, local[2]), + ); + + add(oxt, scale(dir, COOH_BOND_LENGTH)) +} + +/// Builds an orthonormal sp³ frame at `center` relative to `attached`, oriented +/// by `reference`. +fn build_sp3_frame( + center: [f64; 3], + attached: [f64; 3], + reference: [f64; 3], +) -> ([f64; 3], [f64; 3], [f64; 3]) { + let z = normalize(sub(center, attached)); + let ref_vec = normalize(sub(reference, attached)); + let x = normalize(sub(ref_vec, scale(z, dot(z, ref_vec)))); + let y = cross(z, x); + (x, y, z) +} + +#[inline] +fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} + +#[inline] +fn add(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [a[0] + b[0], a[1] + b[1], a[2] + b[2]] +} + +#[inline] +fn scale(a: [f64; 3], s: f64) -> [f64; 3] { + [a[0] * s, a[1] * s, a[2] * s] +} + +#[inline] +fn dot(a: [f64; 3], b: [f64; 3]) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} + +#[inline] +fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] { + [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ] +} + +#[inline] +fn normalize(a: [f64; 3]) -> [f64; 3] { + let n = dot(a, a).sqrt(); + if n < 1e-12 { a } else { scale(a, 1.0 / n) } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::model::metadata::BioMetadata; + + #[test] + fn renames_only_hb_hydrogens() { + let mut types = vec![ + "C_3".to_string(), + "H_HB".to_string(), + "H_".to_string(), + "O_3".to_string(), + ]; + rename_hb_hydrogens(&mut types); + assert_eq!(types, vec!["C_3", "H___A", "H_", "O_3"]); + } + + fn atom( + name: &str, + res: StandardResidue, + pos: ResiduePosition, + xyz: [f64; 3], + ) -> (Atom, AtomResidueInfo) { + let element = if name.starts_with('H') { + Element::H + } else if name.starts_with('O') { + Element::O + } else if name.starts_with('N') { + Element::N + } else { + Element::C + }; + let info = AtomResidueInfo::builder(name, format!("{:?}", res), 1, "A") + .standard_name(Some(res)) + .category(ResidueCategory::Standard) + .position(pos) + .build(); + (Atom::new(element, xyz), info) + } + + fn build_system(entries: Vec<(Atom, AtomResidueInfo)>, bonds: Vec) -> System { + let mut atoms = Vec::new(); + let mut atom_info = Vec::new(); + for (a, i) in entries { + atoms.push(a); + atom_info.push(i); + } + System { + atoms, + bonds, + box_vectors: None, + bio_metadata: Some(BioMetadata { + atom_info, + target_ph: None, + }), + } + } + + #[test] + fn strips_nterminal_h3() { + use ResiduePosition::NTerminal; + use StandardResidue::ALA; + let entries = vec![ + atom("N", ALA, NTerminal, [0.0, 0.0, 0.0]), + atom("H1", ALA, NTerminal, [-0.5, 0.8, 0.0]), + atom("H2", ALA, NTerminal, [-0.5, -0.8, 0.0]), + atom("H3", ALA, NTerminal, [-0.5, 0.0, 0.8]), + atom("CA", ALA, NTerminal, [1.4, 0.0, 0.0]), + ]; + let bonds = vec![ + Bond::new(0, 1, BondOrder::Single), + Bond::new(0, 2, BondOrder::Single), + Bond::new(0, 3, BondOrder::Single), + Bond::new(0, 4, BondOrder::Single), + ]; + let mut system = build_system(entries, bonds); + normalize_terminal_hydrogens(&mut system); + + let names: Vec<&str> = system + .bio_metadata + .as_ref() + .unwrap() + .atom_info + .iter() + .map(|i| i.atom_name.as_str()) + .collect(); + assert_eq!(names, vec!["N", "H1", "H2", "CA"]); + assert_eq!(system.atoms.len(), 4); + // The N–H3 bond is gone; three bonds remain, all valid indices. + assert_eq!(system.bonds.len(), 3); + for b in &system.bonds { + assert!(b.i < 4 && b.j < 4); + } + } + + #[test] + fn builds_cterminal_hoxt() { + use ResiduePosition::CTerminal; + use StandardResidue::ALA; + let entries = vec![ + atom("CA", ALA, CTerminal, [0.0, 0.0, 0.0]), + atom("C", ALA, CTerminal, [1.5, 0.0, 0.0]), + atom("O", ALA, CTerminal, [2.1, 1.0, 0.0]), + atom("OXT", ALA, CTerminal, [2.1, -1.0, 0.0]), + ]; + let bonds = vec![ + Bond::new(0, 1, BondOrder::Single), + Bond::new(1, 2, BondOrder::Double), + Bond::new(1, 3, BondOrder::Single), + ]; + let mut system = build_system(entries, bonds); + normalize_terminal_hydrogens(&mut system); + + let meta = system.bio_metadata.as_ref().unwrap(); + assert_eq!(system.atoms.len(), 5); + let hoxt_idx = meta + .atom_info + .iter() + .position(|i| i.atom_name == "HOXT") + .expect("HOXT added"); + assert_eq!(system.atoms[hoxt_idx].element, Element::H); + + // Bond OXT–HOXT exists. + let oxt_idx = meta + .atom_info + .iter() + .position(|i| i.atom_name == "OXT") + .unwrap(); + assert!( + system + .bonds + .iter() + .any(|b| (b.i == oxt_idx && b.j == hoxt_idx) || (b.j == oxt_idx && b.i == hoxt_idx)) + ); + + // Placed near OXT at ~0.97 Å. + let d = dot( + sub( + system.atoms[hoxt_idx].position, + system.atoms[oxt_idx].position, + ), + sub( + system.atoms[hoxt_idx].position, + system.atoms[oxt_idx].position, + ), + ) + .sqrt(); + assert!((d - COOH_BOND_LENGTH).abs() < 1e-6, "HOXT bond length {d}"); + } + + #[test] + fn idempotent_when_already_neutral() { + use ResiduePosition::CTerminal; + use StandardResidue::ALA; + let entries = vec![ + atom("CA", ALA, CTerminal, [0.0, 0.0, 0.0]), + atom("C", ALA, CTerminal, [1.5, 0.0, 0.0]), + atom("O", ALA, CTerminal, [2.1, 1.0, 0.0]), + atom("OXT", ALA, CTerminal, [2.1, -1.0, 0.0]), + atom("HOXT", ALA, CTerminal, [3.0, -1.0, 0.0]), + ]; + let bonds = vec![ + Bond::new(1, 3, BondOrder::Single), + Bond::new(3, 4, BondOrder::Single), + ]; + let mut system = build_system(entries, bonds); + let before = system.atoms.len(); + normalize_terminal_hydrogens(&mut system); + assert_eq!(system.atoms.len(), before, "no duplicate HOXT"); + } + + #[test] + fn no_metadata_is_noop() { + let mut system = System { + atoms: vec![Atom::new(Element::O, [0.0, 0.0, 0.0])], + bonds: vec![], + box_vectors: None, + bio_metadata: None, + }; + normalize_terminal_hydrogens(&mut system); + assert_eq!(system.atoms.len(), 1); + } +} diff --git a/src/lib.rs b/src/lib.rs index 76e58c3..b3a2627 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -143,7 +143,7 @@ pub use model::metadata::{ pub use forge::{ AnglePotentialType, BasisType, BondPotentialType, ChargeMethod, DampingStrategy, - EmbeddedQeqConfig, ForgeConfig, HybridConfig, LigandChargeConfig, LigandQeqMethod, + EmbeddedQeqConfig, ForgeConfig, HybridConfig, LigandChargeConfig, LigandQeqMethod, MpsimConfig, NucleicScheme, ProteinScheme, QeqConfig, ResidueSelector, SolverOptions, VdwPotentialType, WaterScheme, forge, }; From f043abcc3646acaffdc6b8e79f29bcadc5ad6518 Mon Sep 17 00:00:00 2001 From: Tony Kan Date: Mon, 20 Jul 2026 07:45:24 +0800 Subject: [PATCH 2/4] feat(cli): Add `--mpsim` flag for MPSim-compatible output --- src/bin/dforge/cli.rs | 20 ++++++++++++++++++++ src/bin/dforge/commands/bio.rs | 13 +++++++++++-- src/bin/dforge/commands/chem.rs | 7 ++++++- src/bin/dforge/config/forge.rs | 20 +++++++++++++++----- src/bin/dforge/util/convert.rs | 12 +++++++++++- 5 files changed, 63 insertions(+), 9 deletions(-) diff --git a/src/bin/dforge/cli.rs b/src/bin/dforge/cli.rs index 8c2de1a..1741145 100644 --- a/src/bin/dforge/cli.rs +++ b/src/bin/dforge/cli.rs @@ -138,6 +138,20 @@ pub struct QeqSolverOptions { pub damping: DampingStrategy, } +/// MPSim compatibility options shared by bio and chem commands. +#[derive(Args)] +#[command(next_help_heading = "MPSim Compatibility")] +pub struct MpsimOptions { + /// Emit input for the legacy MPSim EM/MM engine. + /// + /// Renames the DREIDING `H_HB` hydrogen type to MPSim's `H___A` and, for + /// protein systems, forces neutral N/C termini (–NH₂ / –COOH) in both + /// topology and charge, independent of pH. For `bio`, this selects the + /// hybrid charge method unless `--charge` is set explicitly. + #[arg(long = "mpsim")] + pub enabled: bool, +} + /// Potential function options shared by bio and chem commands. #[derive(Args)] #[command(next_help_heading = "Potential Functions")] @@ -199,6 +213,9 @@ pub struct BioArgs { #[command(flatten)] pub potential: PotentialOptions, + + #[command(flatten)] + pub mpsim: MpsimOptions, } #[derive(Args)] @@ -322,6 +339,9 @@ pub struct ChemArgs { #[command(flatten)] pub potential: PotentialOptions, + + #[command(flatten)] + pub mpsim: MpsimOptions, } #[derive(Clone, Copy, ValueEnum)] diff --git a/src/bin/dforge/commands/bio.rs b/src/bin/dforge/commands/bio.rs index 8e5eb83..4c48710 100644 --- a/src/bin/dforge/commands/bio.rs +++ b/src/bin/dforge/commands/bio.rs @@ -44,8 +44,13 @@ pub fn run_bio(args: BioArgs, ctx: DisplayContext) -> Result<()> { } progress.step("Running DREIDING parameterization"); - let forge_config = - build_bio_forge_config(&args.charge, &args.hybrid, &args.qeq, &args.potential)?; + let forge_config = build_bio_forge_config( + &args.charge, + &args.hybrid, + &args.qeq, + &args.potential, + &args.mpsim, + )?; let forged = forge(&system, &forge_config).context("Parameterization failed")?; let param_substeps = build_param_substeps(&args); @@ -226,6 +231,10 @@ fn build_param_substeps(args: &BioArgs) -> Vec { steps.push("Compute H-bond pairs".to_string()); + if args.mpsim.enabled { + steps.push("Apply MPSim adapter (H___A, neutral N/C termini)".to_string()); + } + steps } diff --git a/src/bin/dforge/commands/chem.rs b/src/bin/dforge/commands/chem.rs index 301e4f2..1cacb70 100644 --- a/src/bin/dforge/commands/chem.rs +++ b/src/bin/dforge/commands/chem.rs @@ -39,7 +39,8 @@ pub fn run_chem(args: ChemArgs, ctx: DisplayContext) -> Result<()> { } progress.step("Running DREIDING parameterization"); - let forge_config = build_chem_forge_config(&args.charge, &args.qeq, &args.potential)?; + let forge_config = + build_chem_forge_config(&args.charge, &args.qeq, &args.potential, &args.mpsim)?; let forged = forge(&system, &forge_config).context("Parameterization failed")?; let param_substeps = build_param_substeps(&args); @@ -122,6 +123,10 @@ fn build_param_substeps(args: &ChemArgs) -> Vec { steps.push("Apply custom parameters".to_string()); } + if args.mpsim.enabled { + steps.push("Apply MPSim adapter (H___A type naming)".to_string()); + } + steps } diff --git a/src/bin/dforge/config/forge.rs b/src/bin/dforge/config/forge.rs index d236b4a..5161c64 100644 --- a/src/bin/dforge/config/forge.rs +++ b/src/bin/dforge/config/forge.rs @@ -3,9 +3,11 @@ use std::path::Path; use anyhow::{Context, Result}; -use dreid_forge::ForgeConfig; +use dreid_forge::{ForgeConfig, MpsimConfig}; -use crate::cli::{ChargeOptions, HybridChargeOptions, PotentialOptions, QeqSolverOptions}; +use crate::cli::{ + ChargeOptions, HybridChargeOptions, MpsimOptions, PotentialOptions, QeqSolverOptions, +}; use crate::util::convert::{build_bio_charge_method, build_chem_charge_method}; pub fn build_bio_forge_config( @@ -13,6 +15,7 @@ pub fn build_bio_forge_config( hybrid: &HybridChargeOptions, qeq: &QeqSolverOptions, potential: &PotentialOptions, + mpsim: &MpsimOptions, ) -> Result { let rules = load_optional_file(&potential.rules, "typing rules")?; let params = load_optional_file(&potential.params, "force field parameters")?; @@ -20,11 +23,13 @@ pub fn build_bio_forge_config( Ok(ForgeConfig { rules, params, - charge_method: build_bio_charge_method(charge, hybrid, qeq), + // MPSim's neutral termini need force-field charges to be meaningful, so + // default to the hybrid method when the user left --charge unset. + charge_method: build_bio_charge_method(charge, hybrid, qeq, mpsim.enabled), bond_potential: potential.bond_potential.into(), angle_potential: potential.angle_potential.into(), vdw_potential: potential.vdw_potential.into(), - mpsim: None, + mpsim: build_mpsim_config(mpsim), }) } @@ -32,6 +37,7 @@ pub fn build_chem_forge_config( charge: &ChargeOptions, qeq: &QeqSolverOptions, potential: &PotentialOptions, + mpsim: &MpsimOptions, ) -> Result { let rules = load_optional_file(&potential.rules, "typing rules")?; let params = load_optional_file(&potential.params, "force field parameters")?; @@ -43,10 +49,14 @@ pub fn build_chem_forge_config( bond_potential: potential.bond_potential.into(), angle_potential: potential.angle_potential.into(), vdw_potential: potential.vdw_potential.into(), - mpsim: None, + mpsim: build_mpsim_config(mpsim), }) } +fn build_mpsim_config(mpsim: &MpsimOptions) -> Option { + mpsim.enabled.then(MpsimConfig::default) +} + fn load_optional_file( path: &Option, description: &str, diff --git a/src/bin/dforge/util/convert.rs b/src/bin/dforge/util/convert.rs index e106924..a803a14 100644 --- a/src/bin/dforge/util/convert.rs +++ b/src/bin/dforge/util/convert.rs @@ -31,8 +31,18 @@ pub fn build_bio_charge_method( charge: &cli::ChargeOptions, hybrid: &cli::HybridChargeOptions, qeq: &cli::QeqSolverOptions, + mpsim: bool, ) -> LibChargeMethod { - match charge.method { + // MPSim's neutral-termini charges only exist in the hybrid (force-field) + // path. When --mpsim is set and no charge method was chosen, default to + // hybrid so the requested neutral termini are actually assigned. + let method = if mpsim && matches!(charge.method, cli::ChargeMethod::None) { + cli::ChargeMethod::Hybrid + } else { + charge.method + }; + + match method { cli::ChargeMethod::None => LibChargeMethod::None, cli::ChargeMethod::Qeq => LibChargeMethod::Qeq(QeqConfig { total_charge: charge.total_charge, From 0ceebdfcd49ac7c68ec6e36d83bb7fee522cb885 Mon Sep 17 00:00:00 2001 From: Tony Kan Date: Mon, 20 Jul 2026 07:45:33 +0800 Subject: [PATCH 3/4] docs(cli): Document `--mpsim` compatibility flag in the manual --- MANUAL.md | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/MANUAL.md b/MANUAL.md index 2b82b2c..139a96f 100644 --- a/MANUAL.md +++ b/MANUAL.md @@ -287,6 +287,36 @@ These affect the DREIDING parameterization stage. | `--rules ` | path | none | Custom typing rules (TOML file). | | `--params ` | path | none | Custom force field parameters (TOML file). | +### MPSim Compatibility + +| Flag | Type | Default | Meaning | +| --------- | ---- | ------- | --------------------------------------------- | +| `--mpsim` | flag | off | Emit input for the legacy MPSim EM/MM engine. | + +`--mpsim` adapts the output to MPSim's DREIDING conventions: + +- **Atom types** — the hydrogen-bond hydrogen type `H_HB` is renamed to `H___A`. + This is a pure naming change; parameter lookup and H-bond terms are unaffected. +- **Protein termini** — both chain termini are forced into their neutral, + uncharged protonation state, independent of pH: the N-terminus becomes `–NH₂` + (the extra `H3` is removed) and the C-terminus becomes `–COOH` (the `HOXT` + proton is built onto `OXT`). The matching neutral terminal charge sets are + applied, so each terminal residue's backbone carries no net formal charge. + Nucleic-acid 5′/3′ termini are untouched. + +Because the neutral terminal charges live in the force-field (hybrid) charge +path, `--mpsim` selects `--charge hybrid` automatically when `--charge` is left +unset. An explicit `--charge` value is respected — the terminal _topology_ is +still normalized, but the charges come from the method you chose. + +```bash +# Protein → MPSim-ready BGF in one step (hybrid charges, neutral termini, H___A) +dforge bio -i protein.pdb -o protein.bgf --mpsim +``` + +For `dforge chem`, `--mpsim` only applies the `H_HB → H___A` rename (small +molecules have no protein termini). + --- ## Command: `dforge chem` From 834e2a1295f115940f93a974c7125b31ceb196e2 Mon Sep 17 00:00:00 2001 From: Tony Kan Date: Mon, 20 Jul 2026 17:34:38 +0800 Subject: [PATCH 4/4] refactor(cli): Use `sort_by_key` for descending distribution sort --- src/bin/dforge/display/tables.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/bin/dforge/display/tables.rs b/src/bin/dforge/display/tables.rs index a2259fe..28d836a 100644 --- a/src/bin/dforge/display/tables.rs +++ b/src/bin/dforge/display/tables.rs @@ -64,7 +64,7 @@ fn print_type_distribution(out: &mut impl Write, forged: &ForgedSystem) { (name, count) }) .collect(); - sorted.sort_by(|a, b| b.1.cmp(&a.1)); + sorted.sort_by_key(|(_, count)| std::cmp::Reverse(*count)); print_distribution_table(out, "Atom Type Distribution", &sorted, total); } @@ -80,7 +80,7 @@ fn print_element_distribution(out: &mut impl Write, system: &System) { .into_iter() .map(|(e, c)| (format!("{:?}", e), c)) .collect(); - sorted.sort_by(|a, b| b.1.cmp(&a.1)); + sorted.sort_by_key(|(_, count)| std::cmp::Reverse(*count)); print_distribution_table(out, "Element Distribution", &sorted, total); }