From bc818c7d0b68227ed971794c0232c321b3dddb02 Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 10 Apr 2026 10:31:23 +0100 Subject: [PATCH 1/2] Add repulsion term to the Morse potential --- .../SireOpenMM/sire_to_openmm_system.cpp | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index 7a9978632..fb494d752 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -346,16 +346,14 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res CODELOC); } - // energy expression of a harmonic bond potential, scaled by rho - // e_restraint = rho*DE*(1-exp(-Bij*dr))**2 - // Bij = sqrt(k/2*DE) - // dr = (r - r0) - const auto energy_expression = QString( - "rho*e_restraint;" - "e_restraint=de*(1-exp(-sqrt(k/(2*de))*delta))^2;" - "delta=(r-r0)") - .toStdString(); + "e_total;" + "e_total = rho * e_morse + e_repulsion;" + "e_morse = de * (1 - exp(-alpha * delta))^2;" + "e_repulsion = e_rep * (r_sigma / r)^r_pow;" + "alpha = sqrt(k / (2 * de));" + "delta = (r - r0)") + .toStdString(); auto *restraintff = new OpenMM::CustomBondForce(energy_expression); @@ -365,7 +363,9 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res restraintff->addPerBondParameter("k"); restraintff->addPerBondParameter("r0"); restraintff->addPerBondParameter("de"); - + restraintff->addPerBondParameter("e_rep"); + restraintff->addPerBondParameter("r_sigma"); + restraintff->addPerBondParameter("r_pow"); restraintff->setUsesPeriodicBoundaryConditions(restraints.usesPbc()); lambda_lever.addRestraintIndex(restraints.name(), @@ -376,7 +376,7 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res const double internal_to_k = (1 * SireUnits::kcal_per_mol / (SireUnits::angstrom2)).to(SireUnits::kJ_per_mol / (SireUnits::nanometer2)); const double internal_to_de = (1 * SireUnits::kcal_per_mol).to(SireUnits::kJ_per_mol); - std::vector custom_params = {1.0, 0.0, 0.0, 0.0}; + std::vector custom_params = {1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; for (const auto &restraint : atom_restraints) { @@ -401,6 +401,9 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res custom_params[1] = restraint.k().value() * internal_to_k; // k custom_params[2] = restraint.r0().value() * internal_to_nm; // r0 custom_params[3] = restraint.de().value() * internal_to_de; // de + custom_params[4] = 41.84; // e_rep (kJ/mol) + custom_params[5] = 0.05; // r_sigma (nm) + custom_params[6] = 12; // r_pow restraintff->addBond(atom0_index, atom1_index, custom_params); } From a3316609313612a5bb0864228e85bed24d05dfec Mon Sep 17 00:00:00 2001 From: Audrius Kalpokas Date: Fri, 17 Apr 2026 11:25:54 +0100 Subject: [PATCH 2/2] Test a different r_sigma parameter for MorserRestraint [ci skip] --- wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index fb494d752..089d9cc26 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -402,7 +402,7 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res custom_params[2] = restraint.r0().value() * internal_to_nm; // r0 custom_params[3] = restraint.de().value() * internal_to_de; // de custom_params[4] = 41.84; // e_rep (kJ/mol) - custom_params[5] = 0.05; // r_sigma (nm) + custom_params[5] = 0.025; // r_sigma (nm) custom_params[6] = 12; // r_pow restraintff->addBond(atom0_index, atom1_index, custom_params);