qualia_core_db/specialized_libs/chemistry_modeling/
simulation.rs

1use super::*;
2
3/// Molecular simulator for molecular dynamics simulations
4pub struct MolecularSimulator {
5    simulation_engine: SimulationEngine,
6    force_field_calculator: ForceFieldCalculator,
7    integrator: MolecularIntegrator,
8    boundary_conditions: BoundaryConditions,
9    molecule_store: HashMap<String, Molecule>,
10    linear_algebra: Option<Arc<Mutex<LinearAlgebraLibrary>>>,
11    statistical_computing: Option<Arc<Mutex<StatisticalComputingLibrary>>>,
12}
13
14/// Simulation engine
15pub struct SimulationEngine {
16    simulation_config: SimulationConfig,
17    time_step_control: TimeStepControl,
18    ensemble_manager: EnsembleManager,
19    temperature_controller: TemperatureController,
20}
21
22/// Simulation configuration
23#[derive(Debug, Clone, Serialize, Deserialize)]
24pub struct SimulationConfig {
25    pub simulation_id: String,
26    pub simulation_type: SimulationType,
27    pub ensemble: Ensemble,
28    pub time_step: f64,
29    pub total_time: f64,
30    pub temperature: f64,
31    pub pressure: f64,
32    pub box_size: Vec<f64>,
33    pub boundary_type: BoundaryType,
34}
35
36/// Simulation types
37#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
38pub enum SimulationType {
39    /// Molecular dynamics
40    MolecularDynamics,
41    /// Monte Carlo
42    MonteCarlo,
43    /// Hybrid MD/MC
44    Hybrid,
45    /// Enhanced sampling
46    EnhancedSampling,
47    /// Coarse-grained
48    CoarseGrained,
49    /// Quantum mechanics/molecular mechanics
50    QMMM,
51}
52
53/// Ensembles
54#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
55pub enum Ensemble {
56    NVE,  // Microcanonical
57    NVT,  // Canonical
58    NPT,  // Isothermal-isobaric
59    NPH,  // Isoenthalpic-isobaric
60    MuVT, // Grand canonical
61}
62
63/// Boundary types
64#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
65pub enum BoundaryType {
66    Periodic,
67    NonPeriodic,
68    SemiPeriodic,
69    Ewald,
70    Boiling,
71}
72
73/// Time step control
74pub struct TimeStepControl {
75    control_type: TimeStepControlType,
76    adaptive_parameters: AdaptiveParameters,
77    stability_analysis: StabilityAnalysis,
78}
79
80/// Time step control types
81#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
82pub enum TimeStepControlType {
83    Fixed,
84    Adaptive,
85    Variable,
86    Multiple,
87}
88
89/// Adaptive parameters
90#[derive(Debug, Clone, Serialize, Deserialize)]
91pub struct AdaptiveParameters {
92    pub min_time_step: f64,
93    pub max_time_step: f64,
94    pub safety_factor: f64,
95    pub max_force: f64,
96}
97
98/// Stability analysis
99pub struct StabilityAnalysis {
100    analysis_method: StabilityAnalysisMethod,
101    energy_conservation: EnergyConservation,
102    temperature_fluctuation: TemperatureFluctuation,
103}
104
105/// Stability analysis methods
106#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
107pub enum StabilityAnalysisMethod {
108    EnergyDrift,
109    TemperatureDrift,
110    PressureDrift,
111    ConservationLaws,
112}
113
114/// Energy conservation
115#[derive(Debug, Clone)]
116pub struct EnergyConservation {
117    pub total_energy: f64,
118    pub kinetic_energy: f64,
119    pub potential_energy: f64,
120    pub drift_rate: f64,
121}
122
123/// Temperature fluctuation
124#[derive(Debug, Clone)]
125pub struct TemperatureFluctuation {
126    pub current_temperature: f64,
127    pub target_temperature: f64,
128    pub fluctuation_amplitude: f64,
129    pub heat_capacity: f64,
130}
131
132/// Ensemble manager
133pub struct EnsembleManager {
134    ensembles: HashMap<String, Ensemble>,
135    ensemble_transitions: HashMap<String, EnsembleTransition>,
136    sampling_methods: HashMap<String, SamplingMethod>,
137}
138
139/// Ensemble transitions
140#[derive(Debug, Clone)]
141pub struct EnsembleTransition {
142    pub transition_id: String,
143    pub from_ensemble: Ensemble,
144    pub to_ensemble: Ensemble,
145    pub transition_method: TransitionMethod,
146}
147
148/// Transition methods
149#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
150pub enum TransitionMethod {
151    Berendsen,
152    NoséHoover,
153    Andersen,
154    ParrinelloRahman,
155    MartynaTuckerman,
156    /// Langevin dynamics thermostat (stochastic damping).
157    Langevin,
158}
159
160/// Sampling methods
161#[derive(Debug, Clone)]
162pub struct SamplingMethod {
163    pub method_id: String,
164    pub method_type: SamplingMethodType,
165    pub parameters: SamplingParameters,
166}
167
168/// Sampling method types
169#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
170pub enum SamplingMethodType {
171    Metropolis,
172    Gibbs,
173    WangLandau,
174    Umbrella,
175    ReplicaExchange,
176    /// Hamiltonian (Hybrid) Monte Carlo — uses molecular dynamics proposals.
177    Hamiltonian,
178    /// Parallel tempering (replica exchange over a temperature ladder).
179    ParallelTempering,
180}
181
182/// Sampling parameters
183#[derive(Debug, Clone, Serialize, Deserialize)]
184pub struct SamplingParameters {
185    pub acceptance_ratio: f64,
186    pub proposal_width: f64,
187    pub equilibration_steps: u32,
188    pub production_steps: u32,
189}
190
191/// Temperature controller
192pub struct TemperatureController {
193    control_method: TemperatureControlMethod,
194    thermostat_parameters: ThermostatParameters,
195    temperature_profile: TemperatureProfile,
196}
197
198/// Temperature control methods
199#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
200pub enum TemperatureControlMethod {
201    VelocityRescaling,
202    Berendsen,
203    NoséHoover,
204    Langevin,
205    Andersen,
206}
207
208/// Thermostat parameters
209#[derive(Debug, Clone, Serialize, Deserialize)]
210pub struct ThermostatParameters {
211    pub coupling_constant: f64,
212    pub relaxation_time: f64,
213    pub damping_coefficient: f64,
214}
215
216/// Temperature profile
217#[derive(Debug, Clone)]
218pub struct TemperatureProfile {
219    pub profile_type: TemperatureProfileType,
220    pub initial_temperature: f64,
221    pub final_temperature: Option<f64>,
222    pub ramp_rate: Option<f64>,
223}
224
225/// Temperature profile types
226#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
227pub enum TemperatureProfileType {
228    Constant,
229    Linear,
230    Exponential,
231    Step,
232    Custom,
233}
234
235/// Force field calculator
236pub struct ForceFieldCalculator {
237    force_fields: HashMap<String, ForceField>,
238    interaction_calculator: InteractionCalculator,
239    energy_calculator: EnergyCalculator,
240}
241
242/// Force fields
243#[derive(Debug, Clone)]
244pub struct ForceField {
245    pub field_id: String,
246    pub field_name: String,
247    pub field_type: ForceFieldType,
248    pub parameters: ForceFieldParameters,
249}
250
251/// Force field types
252#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
253pub enum ForceFieldType {
254    AMBER,
255    CHARMM,
256    OPLS,
257    GROMOS,
258    DREIDING,
259    MMFF,
260    ReaxFF,
261    Custom,
262}
263
264/// Force field parameters
265#[derive(Debug, Clone, Serialize, Deserialize)]
266pub struct ForceFieldParameters {
267    pub bond_parameters: Vec<BondParameter>,
268    pub angle_parameters: Vec<AngleParameter>,
269    pub torsion_parameters: Vec<TorsionParameter>,
270    pub nonbonded_parameters: Vec<NonbondedParameter>,
271}
272
273/// Bond parameters
274#[derive(Debug, Clone, Serialize, Deserialize)]
275pub struct BondParameter {
276    pub atom_types: Vec<String>,
277    pub equilibrium_length: f64,
278    pub force_constant: f64,
279}
280
281/// Angle parameters
282#[derive(Debug, Clone, Serialize, Deserialize)]
283pub struct AngleParameter {
284    pub atom_types: Vec<String>,
285    pub equilibrium_angle: f64,
286    pub force_constant: f64,
287}
288
289/// Torsion parameters
290#[derive(Debug, Clone, Serialize, Deserialize)]
291pub struct TorsionParameter {
292    pub atom_types: Vec<String>,
293    pub barriers: Vec<f64>,
294    pub phases: Vec<f64>,
295    pub periodicities: Vec<i32>,
296}
297
298/// Nonbonded parameters
299#[derive(Debug, Clone, Serialize, Deserialize)]
300pub struct NonbondedParameter {
301    pub atom_type: String,
302    pub sigma: f64,
303    pub epsilon: f64,
304    pub charge: f64,
305}
306
307/// Interaction calculator
308pub struct InteractionCalculator {
309    bonded_interactions: BondedInteractions,
310    nonbonded_interactions: NonbondedInteractions,
311    long_range_interactions: LongRangeInteractions,
312}
313
314/// Bonded interactions
315pub struct BondedInteractions {
316    bond_calculator: BondCalculator,
317    angle_calculator: AngleCalculator,
318    torsion_calculator: TorsionCalculator,
319    improper_calculator: ImproperCalculator,
320}
321
322/// Bond calculator
323#[derive(Debug, Clone)]
324pub struct BondCalculator {
325    pub calculator_type: BondCalculatorType,
326    pub parameters: BondCalculatorParameters,
327}
328
329/// Bond calculator types
330#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
331pub enum BondCalculatorType {
332    Harmonic,
333    Morse,
334    FENE,
335    Custom,
336}
337
338/// Bond calculator parameters
339#[derive(Debug, Clone, Serialize, Deserialize)]
340pub struct BondCalculatorParameters {
341    pub force_constant: f64,
342    pub equilibrium_length: f64,
343    pub dissociation_energy: Option<f64>,
344}
345
346/// Angle calculator
347#[derive(Debug, Clone)]
348pub struct AngleCalculator {
349    pub calculator_type: AngleCalculatorType,
350    pub parameters: AngleCalculatorParameters,
351}
352
353/// Angle calculator types
354#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
355pub enum AngleCalculatorType {
356    Harmonic,
357    Cosine,
358    UreyBradley,
359    Custom,
360}
361
362/// Angle calculator parameters
363#[derive(Debug, Clone, Serialize, Deserialize)]
364pub struct AngleCalculatorParameters {
365    pub force_constant: f64,
366    pub equilibrium_angle: f64,
367    pub ub_parameters: Option<UBParameters>,
368}
369
370/// UB parameters
371#[derive(Debug, Clone, Serialize, Deserialize)]
372pub struct UBParameters {
373    pub force_constant: f64,
374    pub equilibrium_length: f64,
375}
376
377/// Torsion calculator
378#[derive(Debug, Clone)]
379pub struct TorsionCalculator {
380    pub calculator_type: TorsionCalculatorType,
381    pub parameters: TorsionCalculatorParameters,
382}
383
384/// Torsion calculator types
385#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
386pub enum TorsionCalculatorType {
387    Cosine,
388    Fourier,
389    RyckaertsBellemans,
390    Custom,
391}
392
393/// Torsion calculator parameters
394#[derive(Debug, Clone, Serialize, Deserialize)]
395pub struct TorsionCalculatorParameters {
396    pub barriers: Vec<f64>,
397    pub phases: Vec<f64>,
398    pub periodicities: Vec<i32>,
399}
400
401/// Improper calculator
402#[derive(Debug, Clone)]
403pub struct ImproperCalculator {
404    pub calculator_type: ImproperCalculatorType,
405    pub parameters: ImproperCalculatorParameters,
406}
407
408/// Improper calculator types
409#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
410pub enum ImproperCalculatorType {
411    Harmonic,
412    Cosine,
413    Custom,
414}
415
416/// Improper calculator parameters
417#[derive(Debug, Clone, Serialize, Deserialize)]
418pub struct ImproperCalculatorParameters {
419    pub force_constant: f64,
420    pub equilibrium_angle: f64,
421}
422
423/// Nonbonded interactions
424pub struct NonbondedInteractions {
425    lennard_jones: LennardJones,
426    coulomb: Coulomb,
427    buckingham: Buckingham,
428}
429
430/// Lennard-Jones potential
431#[derive(Debug, Clone)]
432pub struct LennardJones {
433    pub epsilon: f64,
434    pub sigma: f64,
435    pub cutoff: f64,
436    pub switching_distance: f64,
437}
438
439/// Coulomb potential
440#[derive(Debug, Clone)]
441pub struct Coulomb {
442    pub coulomb_constant: f64,
443    pub dielectric: f64,
444    pub cutoff: f64,
445    pub switching_distance: f64,
446}
447
448/// Buckingham potential
449#[derive(Debug, Clone)]
450pub struct Buckingham {
451    pub a: f64,
452    pub b: f64,
453    pub c: f64,
454    pub cutoff: f64,
455}
456
457/// Long-range interactions
458pub struct LongRangeInteractions {
459    ewald_summation: EwaldSummation,
460    particle_mesh: ParticleMesh,
461    reaction_field: ReactionField,
462}
463
464/// Ewald summation
465#[derive(Debug, Clone)]
466pub struct EwaldSummation {
467    pub alpha: f64,
468    pub k_max: usize,
469    pub real_cutoff: f64,
470    pub reciprocal_cutoff: f64,
471}
472
473/// Particle mesh
474#[derive(Debug, Clone)]
475pub struct ParticleMesh {
476    pub grid_size: Vec<usize>,
477    pub spline_order: usize,
478    pub cutoff: f64,
479}
480
481/// Reaction field
482#[derive(Debug, Clone)]
483pub struct ReactionField {
484    pub dielectric_inside: f64,
485    pub dielectric_outside: f64,
486    pub cutoff: f64,
487}
488
489/// Energy calculator
490pub struct EnergyCalculator {
491    kinetic_energy: KineticEnergy,
492    potential_energy: PotentialEnergy,
493    total_energy: TotalEnergy,
494}
495
496/// Kinetic energy
497#[derive(Debug, Clone)]
498pub struct KineticEnergy {
499    pub temperature: f64,
500    pub degrees_of_freedom: usize,
501    pub velocities: Vec<Vec<f64>>,
502}
503
504/// Potential energy
505#[derive(Debug, Clone)]
506pub struct PotentialEnergy {
507    pub bonded_energy: f64,
508    pub nonbonded_energy: f64,
509    pub long_range_energy: f64,
510}
511
512/// Total energy
513#[derive(Debug, Clone)]
514pub struct TotalEnergy {
515    pub kinetic: f64,
516    pub potential: f64,
517    pub total: f64,
518    pub drift: f64,
519}
520
521/// Molecular integrator
522pub struct MolecularIntegrator {
523    integrator_type: IntegratorType,
524    integrator_parameters: IntegratorParameters,
525    constraint_handler: ConstraintHandler,
526}
527
528/// Integrator types
529#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
530pub enum IntegratorType {
531    VelocityVerlet,
532    Leapfrog,
533    Beeman,
534    Gear,
535    RungeKutta,
536    Stochastic,
537}
538
539/// Integrator parameters
540#[derive(Debug, Clone, Serialize, Deserialize)]
541pub struct IntegratorParameters {
542    pub time_step: f64,
543    pub accuracy: f64,
544    pub stability_factor: f64,
545}
546
547/// Constraint handler
548pub struct ConstraintHandler {
549    constraint_algorithm: ConstraintAlgorithm,
550    constraint_parameters: ConstraintParameters,
551}
552
553/// Constraint algorithms
554#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
555pub enum ConstraintAlgorithm {
556    SHAKE,
557    RATTLE,
558    LINCS,
559    SETTLE,
560}
561
562/// Constraint parameters
563#[derive(Debug, Clone, Serialize, Deserialize)]
564pub struct ConstraintParameters {
565    pub tolerance: f64,
566    pub max_iterations: u32,
567    pub relaxation_parameter: f64,
568}
569
570/// Boundary conditions
571pub struct BoundaryConditions {
572    boundary_type: BoundaryType,
573    box_vectors: Vec<Vec<f64>>,
574    minimum_image: MinimumImage,
575}
576
577/// Minimum image convention
578#[derive(Debug, Clone)]
579pub struct MinimumImage {
580    pub box_size: Vec<f64>,
581    pub periodic: bool,
582}
583
584// Supporting implementations
585
586impl MolecularSimulator {
587    pub fn new() -> Self {
588        Self {
589            simulation_engine: SimulationEngine::new(),
590            force_field_calculator: ForceFieldCalculator::new(),
591            integrator: MolecularIntegrator::new(),
592            boundary_conditions: BoundaryConditions::new(),
593            molecule_store: HashMap::new(),
594            linear_algebra: None,
595            statistical_computing: None,
596        }
597    }
598
599    pub fn attach_dependencies(
600        &mut self,
601        linear_algebra: Option<Arc<Mutex<LinearAlgebraLibrary>>>,
602        statistical_computing: Option<Arc<Mutex<StatisticalComputingLibrary>>>,
603    ) {
604        self.linear_algebra = linear_algebra;
605        self.statistical_computing = statistical_computing;
606    }
607
608    pub fn store_molecule(&mut self, molecule: Molecule) {
609        self.molecule_store
610            .insert(molecule.molecule_id.clone(), molecule);
611    }
612
613    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
614        self.simulation_engine.initialize()?;
615        self.force_field_calculator.initialize()?;
616        self.integrator.initialize()?;
617
618        let _ = self.boundary_conditions.boundary_type();
619        let _ = self.boundary_conditions.box_vectors();
620        let _ = self.boundary_conditions.minimum_image();
621
622        Ok(())
623    }
624
625    pub fn validate_config(&self, config: &SimulationConfig) -> Result<(), ChemistryError> {
626        if config.time_step <= 0.0 {
627            return Err(ChemistryError::ValidationError(
628                "Time step must be positive".to_string(),
629            ));
630        }
631        if config.total_time <= 0.0 {
632            return Err(ChemistryError::ValidationError(
633                "Total time must be positive".to_string(),
634            ));
635        }
636        if config.temperature < 0.0 {
637            return Err(ChemistryError::ValidationError(
638                "Temperature must be non-negative".to_string(),
639            ));
640        }
641        Ok(())
642    }
643
644    pub fn run_simulation(
645        &mut self,
646        config: &SimulationConfig,
647        molecule: &Molecule,
648    ) -> Result<SimulationTrajectory, ChemistryError> {
649        // REAL: Lennard-Jones force field + velocity-Verlet integrator, in
650        // `molecular_dynamics::run_md`. The atoms actually move under computed
651        // forces, total energy is conserved (asserted in that module's tests),
652        // and invalid inputs (no atoms / unparameterized element / bad mass)
653        // return `InsufficientData` rather than a fabricated static trajectory.
654        molecular_dynamics::run_md(
655            config,
656            molecule,
657            self.linear_algebra.clone(),
658            self.statistical_computing.clone(),
659        )
660    }
661
662    pub fn list_force_fields(&self) -> Vec<String> {
663        vec![
664            "AMBER".to_string(),
665            "CHARMM".to_string(),
666            "OPLS".to_string(),
667        ]
668    }
669
670    pub fn get_molecule(&self, molecule_id: &str) -> Option<Molecule> {
671        self.molecule_store.get(molecule_id).cloned()
672    }
673
674    /// Borrow the boundary-conditions configuration.
675    pub fn boundary_conditions(&self) -> &BoundaryConditions {
676        &self.boundary_conditions
677    }
678
679    /// Mutably borrow the boundary-conditions configuration.
680    pub fn boundary_conditions_mut(&mut self) -> &mut BoundaryConditions {
681        &mut self.boundary_conditions
682    }
683
684    /// Borrow the simulation engine.
685    pub fn simulation_engine(&self) -> &SimulationEngine {
686        &self.simulation_engine
687    }
688
689    /// Mutably borrow the simulation engine.
690    pub fn simulation_engine_mut(&mut self) -> &mut SimulationEngine {
691        &mut self.simulation_engine
692    }
693
694    /// Borrow the force-field calculator.
695    pub fn force_field_calculator(&self) -> &ForceFieldCalculator {
696        &self.force_field_calculator
697    }
698
699    /// Mutably borrow the force-field calculator.
700    pub fn force_field_calculator_mut(&mut self) -> &mut ForceFieldCalculator {
701        &mut self.force_field_calculator
702    }
703
704    /// Borrow the molecular integrator.
705    pub fn integrator(&self) -> &MolecularIntegrator {
706        &self.integrator
707    }
708
709    /// Mutably borrow the molecular integrator.
710    pub fn integrator_mut(&mut self) -> &mut MolecularIntegrator {
711        &mut self.integrator
712    }
713}
714
715impl SimulationEngine {
716    pub fn new() -> Self {
717        Self {
718            simulation_config: SimulationConfig::new(),
719            time_step_control: TimeStepControl::new(),
720            ensemble_manager: EnsembleManager::new(),
721            temperature_controller: TemperatureController::new(),
722        }
723    }
724
725    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
726        self.time_step_control.initialize()?;
727        self.ensemble_manager.initialize()?;
728        self.temperature_controller.initialize()?;
729        Ok(())
730    }
731
732    /// Borrow the simulation configuration.
733    pub fn simulation_config(&self) -> &SimulationConfig {
734        &self.simulation_config
735    }
736
737    /// Mutably borrow the simulation configuration.
738    pub fn simulation_config_mut(&mut self) -> &mut SimulationConfig {
739        &mut self.simulation_config
740    }
741
742    /// Borrow the time-step controller.
743    pub fn time_step_control(&self) -> &TimeStepControl {
744        &self.time_step_control
745    }
746
747    /// Mutably borrow the time-step controller.
748    pub fn time_step_control_mut(&mut self) -> &mut TimeStepControl {
749        &mut self.time_step_control
750    }
751
752    /// Borrow the ensemble manager.
753    pub fn ensemble_manager(&self) -> &EnsembleManager {
754        &self.ensemble_manager
755    }
756
757    /// Mutably borrow the ensemble manager.
758    pub fn ensemble_manager_mut(&mut self) -> &mut EnsembleManager {
759        &mut self.ensemble_manager
760    }
761
762    /// Borrow the temperature controller.
763    pub fn temperature_controller(&self) -> &TemperatureController {
764        &self.temperature_controller
765    }
766
767    /// Mutably borrow the temperature controller.
768    pub fn temperature_controller_mut(&mut self) -> &mut TemperatureController {
769        &mut self.temperature_controller
770    }
771}
772
773impl SimulationConfig {
774    pub fn new() -> Self {
775        Self {
776            simulation_id: "sim_1".to_string(),
777            simulation_type: SimulationType::MolecularDynamics,
778            ensemble: Ensemble::NVT,
779            time_step: 0.001,
780            total_time: 1.0,
781            temperature: 300.0,
782            pressure: 1.0,
783            box_size: vec![10.0, 10.0, 10.0],
784            boundary_type: BoundaryType::Periodic,
785        }
786    }
787}
788
789impl TimeStepControl {
790    pub fn new() -> Self {
791        Self {
792            control_type: TimeStepControlType::Fixed,
793            adaptive_parameters: AdaptiveParameters::new(),
794            stability_analysis: StabilityAnalysis::new(),
795        }
796    }
797
798    /// Borrow the time-step control strategy.
799    pub fn control_type(&self) -> &TimeStepControlType {
800        &self.control_type
801    }
802
803    /// Set the time-step control strategy.
804    pub fn set_control_type(&mut self, control_type: TimeStepControlType) {
805        self.control_type = control_type;
806    }
807
808    /// Borrow the adaptive time-step parameters.
809    pub fn adaptive_parameters(&self) -> &AdaptiveParameters {
810        &self.adaptive_parameters
811    }
812
813    /// Mutably borrow the adaptive time-step parameters.
814    pub fn adaptive_parameters_mut(&mut self) -> &mut AdaptiveParameters {
815        &mut self.adaptive_parameters
816    }
817
818    /// Borrow the stability analysis.
819    pub fn stability_analysis(&self) -> &StabilityAnalysis {
820        &self.stability_analysis
821    }
822
823    /// Mutably borrow the stability analysis.
824    pub fn stability_analysis_mut(&mut self) -> &mut StabilityAnalysis {
825        &mut self.stability_analysis
826    }
827
828    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
829        Ok(())
830    }
831}
832
833impl AdaptiveParameters {
834    pub fn new() -> Self {
835        Self {
836            min_time_step: 0.0001,
837            max_time_step: 0.01,
838            safety_factor: 0.9,
839            max_force: 1000.0,
840        }
841    }
842}
843
844impl StabilityAnalysis {
845    pub fn new() -> Self {
846        Self {
847            analysis_method: StabilityAnalysisMethod::EnergyDrift,
848            energy_conservation: EnergyConservation::new(),
849            temperature_fluctuation: TemperatureFluctuation::new(),
850        }
851    }
852
853    /// Borrow the stability analysis method.
854    pub fn analysis_method(&self) -> &StabilityAnalysisMethod {
855        &self.analysis_method
856    }
857
858    /// Set the stability analysis method.
859    pub fn set_analysis_method(&mut self, method: StabilityAnalysisMethod) {
860        self.analysis_method = method;
861    }
862
863    /// Borrow the energy-conservation metrics.
864    pub fn energy_conservation(&self) -> &EnergyConservation {
865        &self.energy_conservation
866    }
867
868    /// Mutably borrow the energy-conservation metrics.
869    pub fn energy_conservation_mut(&mut self) -> &mut EnergyConservation {
870        &mut self.energy_conservation
871    }
872
873    /// Borrow the temperature-fluctuation metrics.
874    pub fn temperature_fluctuation(&self) -> &TemperatureFluctuation {
875        &self.temperature_fluctuation
876    }
877
878    /// Mutably borrow the temperature-fluctuation metrics.
879    pub fn temperature_fluctuation_mut(&mut self) -> &mut TemperatureFluctuation {
880        &mut self.temperature_fluctuation
881    }
882}
883
884impl EnergyConservation {
885    pub fn new() -> Self {
886        Self {
887            total_energy: 0.0,
888            kinetic_energy: 0.0,
889            potential_energy: 0.0,
890            drift_rate: 0.0,
891        }
892    }
893}
894
895impl TemperatureFluctuation {
896    pub fn new() -> Self {
897        Self {
898            current_temperature: 300.0,
899            target_temperature: 300.0,
900            fluctuation_amplitude: 5.0,
901            heat_capacity: 100.0,
902        }
903    }
904}
905
906impl EnsembleManager {
907    pub fn new() -> Self {
908        Self {
909            ensembles: HashMap::new(),
910            ensemble_transitions: HashMap::new(),
911            sampling_methods: HashMap::new(),
912        }
913    }
914
915    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
916        // Standard thermodynamic ensembles. GCMC (grand canonical) maps to the
917        // `MuVT` variant.
918        self.ensembles.insert("NVE".to_string(), Ensemble::NVE);
919        self.ensembles.insert("NVT".to_string(), Ensemble::NVT);
920        self.ensembles.insert("NPT".to_string(), Ensemble::NPT);
921        self.ensembles.insert("GCMC".to_string(), Ensemble::MuVT);
922
923        // Standard ensemble transition / thermostat methods.
924        self.ensemble_transitions.insert(
925            "Berendsen".to_string(),
926            EnsembleTransition {
927                transition_id: "trans_berendsen".to_string(),
928                from_ensemble: Ensemble::NVE,
929                to_ensemble: Ensemble::NVT,
930                transition_method: TransitionMethod::Berendsen,
931            },
932        );
933        self.ensemble_transitions.insert(
934            "Nosé-Hoover".to_string(),
935            EnsembleTransition {
936                transition_id: "trans_nose_hoover".to_string(),
937                from_ensemble: Ensemble::NVT,
938                to_ensemble: Ensemble::NVT,
939                transition_method: TransitionMethod::NoséHoover,
940            },
941        );
942        self.ensemble_transitions.insert(
943            "Parrinello-Rahman".to_string(),
944            EnsembleTransition {
945                transition_id: "trans_parrinello_rahman".to_string(),
946                from_ensemble: Ensemble::NPT,
947                to_ensemble: Ensemble::NPT,
948                transition_method: TransitionMethod::ParrinelloRahman,
949            },
950        );
951        self.ensemble_transitions.insert(
952            "Langevin".to_string(),
953            EnsembleTransition {
954                transition_id: "trans_langevin".to_string(),
955                from_ensemble: Ensemble::NVT,
956                to_ensemble: Ensemble::NVT,
957                transition_method: TransitionMethod::Langevin,
958            },
959        );
960
961        // Standard sampling methods.
962        self.sampling_methods.insert(
963            "Metropolis".to_string(),
964            SamplingMethod {
965                method_id: "sample_metropolis".to_string(),
966                method_type: SamplingMethodType::Metropolis,
967                parameters: SamplingParameters::new(),
968            },
969        );
970        self.sampling_methods.insert(
971            "Gibbs".to_string(),
972            SamplingMethod {
973                method_id: "sample_gibbs".to_string(),
974                method_type: SamplingMethodType::Gibbs,
975                parameters: SamplingParameters::new(),
976            },
977        );
978        self.sampling_methods.insert(
979            "Hamiltonian".to_string(),
980            SamplingMethod {
981                method_id: "sample_hmc".to_string(),
982                method_type: SamplingMethodType::Hamiltonian,
983                parameters: SamplingParameters::new(),
984            },
985        );
986        self.sampling_methods.insert(
987            "ParallelTempering".to_string(),
988            SamplingMethod {
989                method_id: "sample_pt".to_string(),
990                method_type: SamplingMethodType::ParallelTempering,
991                parameters: SamplingParameters::new(),
992            },
993        );
994
995        Ok(())
996    }
997
998    /// Look up a registered ensemble by name.
999    pub fn get_ensemble(&self, name: &str) -> Option<&Ensemble> {
1000        self.ensembles.get(name)
1001    }
1002
1003    /// List the names of all registered ensembles.
1004    pub fn list_ensembles(&self) -> Vec<String> {
1005        self.ensembles.keys().cloned().collect()
1006    }
1007
1008    /// List the names of all registered ensemble transition methods.
1009    pub fn list_transitions(&self) -> Vec<String> {
1010        self.ensemble_transitions.keys().cloned().collect()
1011    }
1012
1013    /// List the names of all registered sampling methods.
1014    pub fn list_sampling_methods(&self) -> Vec<String> {
1015        self.sampling_methods.keys().cloned().collect()
1016    }
1017}
1018
1019impl EnsembleTransition {
1020    pub fn new() -> Self {
1021        Self {
1022            transition_id: "transition_1".to_string(),
1023            from_ensemble: Ensemble::NVE,
1024            to_ensemble: Ensemble::NVT,
1025            transition_method: TransitionMethod::Berendsen,
1026        }
1027    }
1028}
1029
1030impl SamplingMethod {
1031    pub fn new() -> Self {
1032        Self {
1033            method_id: "method_1".to_string(),
1034            method_type: SamplingMethodType::Metropolis,
1035            parameters: SamplingParameters::new(),
1036        }
1037    }
1038}
1039
1040impl SamplingParameters {
1041    pub fn new() -> Self {
1042        Self {
1043            acceptance_ratio: 0.5,
1044            proposal_width: 1.0,
1045            equilibration_steps: 1000,
1046            production_steps: 10000,
1047        }
1048    }
1049}
1050
1051impl TemperatureController {
1052    pub fn new() -> Self {
1053        Self {
1054            control_method: TemperatureControlMethod::NoséHoover,
1055            thermostat_parameters: ThermostatParameters::new(),
1056            temperature_profile: TemperatureProfile::new(),
1057        }
1058    }
1059
1060    /// Borrow the temperature control method.
1061    pub fn control_method(&self) -> &TemperatureControlMethod {
1062        &self.control_method
1063    }
1064
1065    /// Set the temperature control method.
1066    pub fn set_control_method(&mut self, method: TemperatureControlMethod) {
1067        self.control_method = method;
1068    }
1069
1070    /// Borrow the thermostat parameters.
1071    pub fn thermostat_parameters(&self) -> &ThermostatParameters {
1072        &self.thermostat_parameters
1073    }
1074
1075    /// Mutably borrow the thermostat parameters.
1076    pub fn thermostat_parameters_mut(&mut self) -> &mut ThermostatParameters {
1077        &mut self.thermostat_parameters
1078    }
1079
1080    /// Borrow the temperature profile.
1081    pub fn temperature_profile(&self) -> &TemperatureProfile {
1082        &self.temperature_profile
1083    }
1084
1085    /// Mutably borrow the temperature profile.
1086    pub fn temperature_profile_mut(&mut self) -> &mut TemperatureProfile {
1087        &mut self.temperature_profile
1088    }
1089
1090    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
1091        Ok(())
1092    }
1093}
1094
1095impl ThermostatParameters {
1096    pub fn new() -> Self {
1097        Self {
1098            coupling_constant: 1.0,
1099            relaxation_time: 100.0,
1100            damping_coefficient: 1.0,
1101        }
1102    }
1103}
1104
1105impl TemperatureProfile {
1106    pub fn new() -> Self {
1107        Self {
1108            profile_type: TemperatureProfileType::Constant,
1109            initial_temperature: 300.0,
1110            final_temperature: None,
1111            ramp_rate: None,
1112        }
1113    }
1114}
1115
1116impl ForceFieldCalculator {
1117    pub fn new() -> Self {
1118        Self {
1119            force_fields: HashMap::new(),
1120            interaction_calculator: InteractionCalculator::new(),
1121            energy_calculator: EnergyCalculator::new(),
1122        }
1123    }
1124
1125    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
1126        self.interaction_calculator.initialize()?;
1127        self.energy_calculator.initialize()?;
1128        // Populate the standard force-field catalogue. Each entry ships with a
1129        // default parameter set (≥1 bond/angle/torsion/nonbonded entry) so the
1130        // calculator is usable immediately after initialization.
1131        self.register_standard_force_fields();
1132        Ok(())
1133    }
1134
1135    /// Register the built-in force-field definitions (AMBER, CHARMM, OPLS-AA,
1136    /// GROMOS, Universal/UFF). `Universal` has no dedicated `ForceFieldType`
1137    /// variant, so it is tagged `Custom` with a descriptive name.
1138    fn register_standard_force_fields(&mut self) {
1139        self.force_fields.insert(
1140            "AMBER".to_string(),
1141            ForceField {
1142                field_id: "ff_amber".to_string(),
1143                field_name: "AMBER".to_string(),
1144                field_type: ForceFieldType::AMBER,
1145                parameters: ForceFieldParameters::new(),
1146            },
1147        );
1148        self.force_fields.insert(
1149            "CHARMM".to_string(),
1150            ForceField {
1151                field_id: "ff_charmm".to_string(),
1152                field_name: "CHARMM".to_string(),
1153                field_type: ForceFieldType::CHARMM,
1154                parameters: ForceFieldParameters::new(),
1155            },
1156        );
1157        self.force_fields.insert(
1158            "OPLS".to_string(),
1159            ForceField {
1160                field_id: "ff_opls".to_string(),
1161                field_name: "OPLS-AA".to_string(),
1162                field_type: ForceFieldType::OPLS,
1163                parameters: ForceFieldParameters::new(),
1164            },
1165        );
1166        self.force_fields.insert(
1167            "GROMOS".to_string(),
1168            ForceField {
1169                field_id: "ff_gromos".to_string(),
1170                field_name: "GROMOS".to_string(),
1171                field_type: ForceFieldType::GROMOS,
1172                parameters: ForceFieldParameters::new(),
1173            },
1174        );
1175        self.force_fields.insert(
1176            "Universal".to_string(),
1177            ForceField {
1178                field_id: "ff_uff".to_string(),
1179                field_name: "Universal (UFF)".to_string(),
1180                field_type: ForceFieldType::Custom,
1181                parameters: ForceFieldParameters::new(),
1182            },
1183        );
1184    }
1185
1186    /// Look up a registered force field by name.
1187    pub fn get_force_field(&self, name: &str) -> Option<&ForceField> {
1188        self.force_fields.get(name)
1189    }
1190
1191    /// List the names of all registered force fields.
1192    pub fn list_force_fields(&self) -> Vec<String> {
1193        self.force_fields.keys().cloned().collect()
1194    }
1195
1196    /// Register a custom force field under `name`, replacing any existing entry.
1197    pub fn register_force_field(&mut self, name: &str, force_field: ForceField) {
1198        self.force_fields.insert(name.to_string(), force_field);
1199    }
1200}
1201
1202impl ForceField {
1203    pub fn new() -> Self {
1204        Self {
1205            field_id: "ff_1".to_string(),
1206            field_name: "AMBER".to_string(),
1207            field_type: ForceFieldType::AMBER,
1208            parameters: ForceFieldParameters::new(),
1209        }
1210    }
1211}
1212
1213impl ForceFieldParameters {
1214    pub fn new() -> Self {
1215        Self {
1216            bond_parameters: vec![BondParameter::new()],
1217            angle_parameters: vec![AngleParameter::new()],
1218            torsion_parameters: vec![TorsionParameter::new()],
1219            nonbonded_parameters: vec![NonbondedParameter::new()],
1220        }
1221    }
1222}
1223
1224impl BondParameter {
1225    pub fn new() -> Self {
1226        Self {
1227            atom_types: vec!["C".to_string(), "H".to_string()],
1228            equilibrium_length: 1.09,
1229            force_constant: 450.0,
1230        }
1231    }
1232}
1233
1234impl AngleParameter {
1235    pub fn new() -> Self {
1236        Self {
1237            atom_types: vec!["C".to_string(), "H".to_string(), "H".to_string()],
1238            equilibrium_angle: 109.5,
1239            force_constant: 50.0,
1240        }
1241    }
1242}
1243
1244impl TorsionParameter {
1245    pub fn new() -> Self {
1246        Self {
1247            atom_types: vec![
1248                "C".to_string(),
1249                "C".to_string(),
1250                "C".to_string(),
1251                "C".to_string(),
1252            ],
1253            barriers: vec![0.0, 1.0],
1254            phases: vec![0.0, 180.0],
1255            periodicities: vec![1, 2],
1256        }
1257    }
1258}
1259
1260impl NonbondedParameter {
1261    pub fn new() -> Self {
1262        Self {
1263            atom_type: "C".to_string(),
1264            sigma: 3.4,
1265            epsilon: 0.086,
1266            charge: 0.0,
1267        }
1268    }
1269}
1270
1271impl InteractionCalculator {
1272    pub fn new() -> Self {
1273        Self {
1274            bonded_interactions: BondedInteractions::new(),
1275            nonbonded_interactions: NonbondedInteractions::new(),
1276            long_range_interactions: LongRangeInteractions::new(),
1277        }
1278    }
1279
1280    /// Borrow the bonded-interactions calculator.
1281    pub fn bonded_interactions(&self) -> &BondedInteractions {
1282        &self.bonded_interactions
1283    }
1284
1285    /// Mutably borrow the bonded-interactions calculator.
1286    pub fn bonded_interactions_mut(&mut self) -> &mut BondedInteractions {
1287        &mut self.bonded_interactions
1288    }
1289
1290    /// Borrow the nonbonded-interactions calculator.
1291    pub fn nonbonded_interactions(&self) -> &NonbondedInteractions {
1292        &self.nonbonded_interactions
1293    }
1294
1295    /// Mutably borrow the nonbonded-interactions calculator.
1296    pub fn nonbonded_interactions_mut(&mut self) -> &mut NonbondedInteractions {
1297        &mut self.nonbonded_interactions
1298    }
1299
1300    /// Borrow the long-range-interactions calculator.
1301    pub fn long_range_interactions(&self) -> &LongRangeInteractions {
1302        &self.long_range_interactions
1303    }
1304
1305    /// Mutably borrow the long-range-interactions calculator.
1306    pub fn long_range_interactions_mut(&mut self) -> &mut LongRangeInteractions {
1307        &mut self.long_range_interactions
1308    }
1309
1310    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
1311        let _ = self.nonbonded_interactions.lennard_jones();
1312        let _ = self.nonbonded_interactions.coulomb();
1313        let _ = self.nonbonded_interactions.buckingham();
1314        let _ = self.long_range_interactions.ewald_summation();
1315        let _ = self.long_range_interactions.particle_mesh();
1316        let _ = self.long_range_interactions.reaction_field();
1317        Ok(())
1318    }
1319}
1320
1321impl BondedInteractions {
1322    pub fn new() -> Self {
1323        Self {
1324            bond_calculator: BondCalculator::new(),
1325            angle_calculator: AngleCalculator::new(),
1326            torsion_calculator: TorsionCalculator::new(),
1327            improper_calculator: ImproperCalculator::new(),
1328        }
1329    }
1330
1331    /// Borrow the bond-stretch calculator.
1332    pub fn bond_calculator(&self) -> &BondCalculator {
1333        &self.bond_calculator
1334    }
1335
1336    /// Mutably borrow the bond-stretch calculator.
1337    pub fn bond_calculator_mut(&mut self) -> &mut BondCalculator {
1338        &mut self.bond_calculator
1339    }
1340
1341    /// Borrow the angle-bend calculator.
1342    pub fn angle_calculator(&self) -> &AngleCalculator {
1343        &self.angle_calculator
1344    }
1345
1346    /// Mutably borrow the angle-bend calculator.
1347    pub fn angle_calculator_mut(&mut self) -> &mut AngleCalculator {
1348        &mut self.angle_calculator
1349    }
1350
1351    /// Borrow the torsion calculator.
1352    pub fn torsion_calculator(&self) -> &TorsionCalculator {
1353        &self.torsion_calculator
1354    }
1355
1356    /// Mutably borrow the torsion calculator.
1357    pub fn torsion_calculator_mut(&mut self) -> &mut TorsionCalculator {
1358        &mut self.torsion_calculator
1359    }
1360
1361    /// Borrow the improper-torsion calculator.
1362    pub fn improper_calculator(&self) -> &ImproperCalculator {
1363        &self.improper_calculator
1364    }
1365
1366    /// Mutably borrow the improper-torsion calculator.
1367    pub fn improper_calculator_mut(&mut self) -> &mut ImproperCalculator {
1368        &mut self.improper_calculator
1369    }
1370}
1371
1372impl BondCalculator {
1373    pub fn new() -> Self {
1374        Self {
1375            calculator_type: BondCalculatorType::Harmonic,
1376            parameters: BondCalculatorParameters::new(),
1377        }
1378    }
1379}
1380
1381impl BondCalculatorParameters {
1382    pub fn new() -> Self {
1383        Self {
1384            force_constant: 450.0,
1385            equilibrium_length: 1.09,
1386            dissociation_energy: None,
1387        }
1388    }
1389}
1390
1391impl AngleCalculator {
1392    pub fn new() -> Self {
1393        Self {
1394            calculator_type: AngleCalculatorType::Harmonic,
1395            parameters: AngleCalculatorParameters::new(),
1396        }
1397    }
1398}
1399
1400impl AngleCalculatorParameters {
1401    pub fn new() -> Self {
1402        Self {
1403            force_constant: 50.0,
1404            equilibrium_angle: 109.5,
1405            ub_parameters: None,
1406        }
1407    }
1408}
1409
1410impl TorsionCalculator {
1411    pub fn new() -> Self {
1412        Self {
1413            calculator_type: TorsionCalculatorType::Cosine,
1414            parameters: TorsionCalculatorParameters::new(),
1415        }
1416    }
1417}
1418
1419impl TorsionCalculatorParameters {
1420    pub fn new() -> Self {
1421        Self {
1422            barriers: vec![0.0, 1.0],
1423            phases: vec![0.0, 180.0],
1424            periodicities: vec![1, 2],
1425        }
1426    }
1427}
1428
1429impl ImproperCalculator {
1430    pub fn new() -> Self {
1431        Self {
1432            calculator_type: ImproperCalculatorType::Harmonic,
1433            parameters: ImproperCalculatorParameters::new(),
1434        }
1435    }
1436}
1437
1438impl ImproperCalculatorParameters {
1439    pub fn new() -> Self {
1440        Self {
1441            force_constant: 50.0,
1442            equilibrium_angle: 109.5,
1443        }
1444    }
1445}
1446
1447impl NonbondedInteractions {
1448    pub fn new() -> Self {
1449        Self {
1450            lennard_jones: LennardJones::new(),
1451            coulomb: Coulomb::new(),
1452            buckingham: Buckingham::new(),
1453        }
1454    }
1455
1456    /// Borrow the Lennard-Jones potential parameters.
1457    pub fn lennard_jones(&self) -> &LennardJones {
1458        &self.lennard_jones
1459    }
1460
1461    /// Mutably borrow the Lennard-Jones potential parameters.
1462    pub fn lennard_jones_mut(&mut self) -> &mut LennardJones {
1463        &mut self.lennard_jones
1464    }
1465
1466    /// Borrow the Coulomb potential parameters.
1467    pub fn coulomb(&self) -> &Coulomb {
1468        &self.coulomb
1469    }
1470
1471    /// Mutably borrow the Coulomb potential parameters.
1472    pub fn coulomb_mut(&mut self) -> &mut Coulomb {
1473        &mut self.coulomb
1474    }
1475
1476    /// Borrow the Buckingham potential parameters.
1477    pub fn buckingham(&self) -> &Buckingham {
1478        &self.buckingham
1479    }
1480
1481    /// Mutably borrow the Buckingham potential parameters.
1482    pub fn buckingham_mut(&mut self) -> &mut Buckingham {
1483        &mut self.buckingham
1484    }
1485}
1486
1487impl LennardJones {
1488    pub fn new() -> Self {
1489        Self {
1490            epsilon: 0.086,
1491            sigma: 3.4,
1492            cutoff: 12.0,
1493            switching_distance: 10.0,
1494        }
1495    }
1496}
1497
1498impl Coulomb {
1499    pub fn new() -> Self {
1500        Self {
1501            coulomb_constant: 332.06,
1502            dielectric: 1.0,
1503            cutoff: 12.0,
1504            switching_distance: 10.0,
1505        }
1506    }
1507}
1508
1509impl Buckingham {
1510    pub fn new() -> Self {
1511        Self {
1512            a: 1000.0,
1513            b: 3.5,
1514            c: 0.0,
1515            cutoff: 12.0,
1516        }
1517    }
1518}
1519
1520impl LongRangeInteractions {
1521    pub fn new() -> Self {
1522        Self {
1523            ewald_summation: EwaldSummation::new(),
1524            particle_mesh: ParticleMesh::new(),
1525            reaction_field: ReactionField::new(),
1526        }
1527    }
1528
1529    /// Borrow the Ewald-summation parameters.
1530    pub fn ewald_summation(&self) -> &EwaldSummation {
1531        &self.ewald_summation
1532    }
1533
1534    /// Mutably borrow the Ewald-summation parameters.
1535    pub fn ewald_summation_mut(&mut self) -> &mut EwaldSummation {
1536        &mut self.ewald_summation
1537    }
1538
1539    /// Borrow the particle-mesh (PME) parameters.
1540    pub fn particle_mesh(&self) -> &ParticleMesh {
1541        &self.particle_mesh
1542    }
1543
1544    /// Mutably borrow the particle-mesh (PME) parameters.
1545    pub fn particle_mesh_mut(&mut self) -> &mut ParticleMesh {
1546        &mut self.particle_mesh
1547    }
1548
1549    /// Borrow the reaction-field parameters.
1550    pub fn reaction_field(&self) -> &ReactionField {
1551        &self.reaction_field
1552    }
1553
1554    /// Mutably borrow the reaction-field parameters.
1555    pub fn reaction_field_mut(&mut self) -> &mut ReactionField {
1556        &mut self.reaction_field
1557    }
1558}
1559
1560impl EwaldSummation {
1561    pub fn new() -> Self {
1562        Self {
1563            alpha: 0.3,
1564            k_max: 10,
1565            real_cutoff: 12.0,
1566            reciprocal_cutoff: 10.0,
1567        }
1568    }
1569}
1570
1571impl ParticleMesh {
1572    pub fn new() -> Self {
1573        Self {
1574            grid_size: vec![32, 32, 32],
1575            spline_order: 4,
1576            cutoff: 12.0,
1577        }
1578    }
1579}
1580
1581impl ReactionField {
1582    pub fn new() -> Self {
1583        Self {
1584            dielectric_inside: 1.0,
1585            dielectric_outside: 78.5,
1586            cutoff: 12.0,
1587        }
1588    }
1589}
1590
1591impl EnergyCalculator {
1592    pub fn new() -> Self {
1593        Self {
1594            kinetic_energy: KineticEnergy::new(),
1595            potential_energy: PotentialEnergy::new(),
1596            total_energy: TotalEnergy::new(),
1597        }
1598    }
1599
1600    /// Borrow the kinetic-energy state.
1601    pub fn kinetic_energy(&self) -> &KineticEnergy {
1602        &self.kinetic_energy
1603    }
1604
1605    /// Mutably borrow the kinetic-energy state.
1606    pub fn kinetic_energy_mut(&mut self) -> &mut KineticEnergy {
1607        &mut self.kinetic_energy
1608    }
1609
1610    /// Borrow the potential-energy state.
1611    pub fn potential_energy(&self) -> &PotentialEnergy {
1612        &self.potential_energy
1613    }
1614
1615    /// Mutably borrow the potential-energy state.
1616    pub fn potential_energy_mut(&mut self) -> &mut PotentialEnergy {
1617        &mut self.potential_energy
1618    }
1619
1620    /// Borrow the total-energy state.
1621    pub fn total_energy(&self) -> &TotalEnergy {
1622        &self.total_energy
1623    }
1624
1625    /// Mutably borrow the total-energy state.
1626    pub fn total_energy_mut(&mut self) -> &mut TotalEnergy {
1627        &mut self.total_energy
1628    }
1629
1630    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
1631        Ok(())
1632    }
1633}
1634
1635impl KineticEnergy {
1636    pub fn new() -> Self {
1637        Self {
1638            temperature: 300.0,
1639            degrees_of_freedom: 0,
1640            velocities: Vec::new(),
1641        }
1642    }
1643}
1644
1645impl PotentialEnergy {
1646    pub fn new() -> Self {
1647        Self {
1648            bonded_energy: 0.0,
1649            nonbonded_energy: 0.0,
1650            long_range_energy: 0.0,
1651        }
1652    }
1653}
1654
1655impl TotalEnergy {
1656    pub fn new() -> Self {
1657        Self {
1658            kinetic: 0.0,
1659            potential: 0.0,
1660            total: 0.0,
1661            drift: 0.0,
1662        }
1663    }
1664}
1665
1666impl MolecularIntegrator {
1667    pub fn new() -> Self {
1668        Self {
1669            integrator_type: IntegratorType::VelocityVerlet,
1670            integrator_parameters: IntegratorParameters::new(),
1671            constraint_handler: ConstraintHandler::new(),
1672        }
1673    }
1674
1675    /// Borrow the integrator type.
1676    pub fn integrator_type(&self) -> &IntegratorType {
1677        &self.integrator_type
1678    }
1679
1680    /// Set the integrator type.
1681    pub fn set_integrator_type(&mut self, integrator_type: IntegratorType) {
1682        self.integrator_type = integrator_type;
1683    }
1684
1685    /// Borrow the integrator parameters.
1686    pub fn integrator_parameters(&self) -> &IntegratorParameters {
1687        &self.integrator_parameters
1688    }
1689
1690    /// Mutably borrow the integrator parameters.
1691    pub fn integrator_parameters_mut(&mut self) -> &mut IntegratorParameters {
1692        &mut self.integrator_parameters
1693    }
1694
1695    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
1696        self.constraint_handler.initialize()?;
1697        Ok(())
1698    }
1699}
1700
1701impl IntegratorParameters {
1702    pub fn new() -> Self {
1703        Self {
1704            time_step: 0.001,
1705            accuracy: 1e-6,
1706            stability_factor: 0.9,
1707        }
1708    }
1709}
1710
1711impl ConstraintHandler {
1712    pub fn new() -> Self {
1713        Self {
1714            constraint_algorithm: ConstraintAlgorithm::SHAKE,
1715            constraint_parameters: ConstraintParameters::new(),
1716        }
1717    }
1718
1719    /// Borrow the constraint algorithm.
1720    pub fn constraint_algorithm(&self) -> &ConstraintAlgorithm {
1721        &self.constraint_algorithm
1722    }
1723
1724    /// Set the constraint algorithm.
1725    pub fn set_constraint_algorithm(&mut self, algorithm: ConstraintAlgorithm) {
1726        self.constraint_algorithm = algorithm;
1727    }
1728
1729    /// Borrow the constraint parameters.
1730    pub fn constraint_parameters(&self) -> &ConstraintParameters {
1731        &self.constraint_parameters
1732    }
1733
1734    /// Mutably borrow the constraint parameters.
1735    pub fn constraint_parameters_mut(&mut self) -> &mut ConstraintParameters {
1736        &mut self.constraint_parameters
1737    }
1738
1739    pub fn initialize(&mut self) -> Result<(), ChemistryError> {
1740        Ok(())
1741    }
1742}
1743
1744impl ConstraintParameters {
1745    pub fn new() -> Self {
1746        Self {
1747            tolerance: 1e-6,
1748            max_iterations: 100,
1749            relaxation_parameter: 0.1,
1750        }
1751    }
1752}
1753
1754impl BoundaryConditions {
1755    pub fn new() -> Self {
1756        Self {
1757            boundary_type: BoundaryType::Periodic,
1758            box_vectors: vec![
1759                vec![10.0, 0.0, 0.0],
1760                vec![0.0, 10.0, 0.0],
1761                vec![0.0, 0.0, 10.0],
1762            ],
1763            minimum_image: MinimumImage::new(),
1764        }
1765    }
1766
1767    /// Borrow the boundary type.
1768    pub fn boundary_type(&self) -> &BoundaryType {
1769        &self.boundary_type
1770    }
1771
1772    /// Set the boundary type.
1773    pub fn set_boundary_type(&mut self, boundary_type: BoundaryType) {
1774        self.boundary_type = boundary_type;
1775    }
1776
1777    /// Borrow the simulation box vectors.
1778    pub fn box_vectors(&self) -> &Vec<Vec<f64>> {
1779        &self.box_vectors
1780    }
1781
1782    /// Mutably borrow the simulation box vectors.
1783    pub fn box_vectors_mut(&mut self) -> &mut Vec<Vec<f64>> {
1784        &mut self.box_vectors
1785    }
1786
1787    /// Borrow the minimum-image convention state.
1788    pub fn minimum_image(&self) -> &MinimumImage {
1789        &self.minimum_image
1790    }
1791
1792    /// Mutably borrow the minimum-image convention state.
1793    pub fn minimum_image_mut(&mut self) -> &mut MinimumImage {
1794        &mut self.minimum_image
1795    }
1796}
1797
1798impl MinimumImage {
1799    pub fn new() -> Self {
1800        Self {
1801            box_size: vec![10.0, 10.0, 10.0],
1802            periodic: true,
1803        }
1804    }
1805}
qualia_core_db/specialized_libs/chemistry_modeling/simulation.rs

qualia_core_db/specialized_libs/chemistry_modeling/
simulation.rs
