diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index a2eca4a2c..b4f29d783 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -52,6 +52,8 @@ organisation on `GitHub `__. * Print warning when ``sire.legacy.IO.AmberPrm`` parser re-orders residues due to covalent bonds between molecules. +* Added internal implementation of virtual sites in the Sire-to-OpenMM conversion layer + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/tests/convert/test_openmm_vsites.py b/tests/convert/test_openmm_vsites.py new file mode 100644 index 000000000..fd1208c58 --- /dev/null +++ b/tests/convert/test_openmm_vsites.py @@ -0,0 +1,300 @@ +import sire as sr +import pytest + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_vsite_params(ethane_12dichloroethane, openmm_platform): + # Can we create an openmm system with the correct vsite parameters and charges + mols = ethane_12dichloroethane + + # Just dichloroethane + mol = mols[0].property("molecule1") + + # Set vsite properties + vsite_dict = { + "0": { + "vs_indices": [0, 1, 2], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + "1": { + "vs_indices": [3, 2, 1], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + } + + parents_dict = {str(atom_i): [] for atom_i in range(mol.num_atoms())} + for v, vs in enumerate(vsite_dict): + parent = vsite_dict[vs]["vs_indices"][0] + parents_dict[str(parent)].append(v) + + n_virtual_sites = len(vsite_dict) + vs_charges = [0.2, 0.2] + + cursor = mol.cursor() + cursor.set("n_virtual_sites", n_virtual_sites) + cursor.set("vs_charges", vs_charges) + cursor.set("virtual_sites", vsite_dict) + cursor.set("parents", parents_dict) + mol = cursor.commit() + + map = {"platform": openmm_platform} + + # Create openmm system + omm = sr.convert.to(mol, "openmm", map=map) + + omm_system = omm.getSystem() + + # Check the right number of virtual sites are added + assert omm_system.getNumParticles() == mol.num_atoms() + n_virtual_sites + + # Check the VS parameters and charges are correct + + from openmm import LocalCoordinatesSite, Vec3, unit + + nb_force = next( + force for force in omm_system.getForces() if force.getName() == "NonbondedForce" + ) + + for vs_index in range(n_virtual_sites): + omm_vs = omm_system.getVirtualSite(mol.num_atoms() + vs_index) + assert isinstance(omm_vs, LocalCoordinatesSite) + + origin_weights = omm_vs.getOriginWeights() + assert list(origin_weights) == vsite_dict[str(vs_index)]["vs_ows"] + + x_weights = omm_vs.getXWeights() + assert list(x_weights) == vsite_dict[str(vs_index)]["vs_xs"] + + y_weights = omm_vs.getYWeights() + assert list(y_weights) == vsite_dict[str(vs_index)]["vs_ys"] + + local_position = omm_vs.getLocalPosition() + assert local_position.value_in_unit(unit.nanometer) == Vec3( + *[x for x in vsite_dict[str(vs_index)]["vs_local"]] + ) + + nb_params = nb_force.getParticleParameters(mol.num_atoms() + vs_index) + assert vs_charges[vs_index] == nb_params[0].value_in_unit( + unit.elementary_charge + ) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_vsite_pertubation(ethane_12dichloroethane, openmm_platform): + # Are vsite parameters scaled correctly by lambda + mols = ethane_12dichloroethane + + # Just dichloroethane + mol = mols[0] + + # Set vsite properties + vsite_dict = { + "0": { + "vs_indices": [0, 1, 2], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + "1": { + "vs_indices": [3, 2, 1], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + }, + } + + parents_dict = {str(atom_i): [] for atom_i in range(mol.num_atoms())} + for v, vs in enumerate(vsite_dict): + parent = vsite_dict[vs]["vs_indices"][0] + parents_dict[str(parent)].append(v) + + n_virtual_sites = len(vsite_dict) + vs_charges0 = [0.1, 0.1] + vs_charges1 = [0.2, 0.2] + + cursor = mol.cursor() + cursor.set("n_virtual_sites", n_virtual_sites) + cursor.set("vs_charges0", vs_charges0) + cursor.set("vs_charges1", vs_charges1) + cursor.set("virtual_sites", vsite_dict) + cursor.set("parents", parents_dict) + mol = cursor.commit() + + mol = sr.morph.link_to_reference(mol) + + from openmm import unit + + for lam in [0.0, 0.5, 1.0]: + d = mol.dynamics(lambda_value=lam, platform=openmm_platform) + system = d.context().getSystem() + nb_force = next( + force for force in system.getForces() if force.getName() == "NonbondedForce" + ) + for vs_index in range(n_virtual_sites): + expected_charge = (1 - lam) * vs_charges0[vs_index] + lam * vs_charges1[ + vs_index + ] + nb_charge = nb_force.getParticleParameters(mol.num_atoms() + vs_index)[0] + assert expected_charge == nb_charge.value_in_unit(unit.elementary_charge) + + +@pytest.mark.skipif( + "openmm" not in sr.convert.supported_formats(), + reason="openmm support is not available", +) +def test_vsite_restraints(solvated_neopentane_methane, openmm_platform): + # Are restraints added correctly to vsite systems + mols = solvated_neopentane_methane + + mol0 = mols[0] + + # Arbitrary vsite definition, we just want to check that the restraints + # are correctly mapped to the new system with vsites + vsite_dict = { + "0": { + "vs_indices": [0, 1, 2], + "vs_ows": [1, 0, 0], + "vs_xs": [1, -1, 0], + "vs_ys": [0, 1, -1], + "vs_local": [0.03, 0, 0], + } + } + + parents_dict = {str(atom_i): [] for atom_i in range(mol0.num_atoms())} + for v, vs in enumerate(vsite_dict): + parent = vsite_dict[vs]["vs_indices"][0] + parents_dict[str(parent)].append(v) + + n_virtual_sites = len(vsite_dict) + vs_charges0 = [0.1, 0.1] + vs_charges1 = [0.2, 0.2] + + cursor = mol0.cursor() + cursor.set("n_virtual_sites", n_virtual_sites) + cursor.set("vs_charges0", vs_charges0) + cursor.set("vs_charges1", vs_charges1) + cursor.set("virtual_sites", vsite_dict) + cursor.set("parents", parents_dict) + mol0 = cursor.commit() + mols.update(mol0) + + mols = sr.morph.link_to_reference(mols) + + start_mol1 = mol0.num_atoms() + n_virtual_sites + base_particles = sum(m.num_atoms() for m in mols) + n_virtual_sites + + # Applying restraints to the first N water atoms (not realistic but we + # just want to check the mapping is correct) + restr_atoms = mols[1:].atoms() + + restraint_cases = [ + ( + "bond", + lambda system: sr.restraints.bond( + system, atoms0=restr_atoms[0], atoms1=restr_atoms[1] + ), + "BondRestraintForce", + lambda force: ( + force.getBondParameters(0)[:2] == [start_mol1 + 0, start_mol1 + 1] + ), + ), + ( + "inverse_bond", + lambda system: sr.restraints.inverse_bond( + system, atoms0=restr_atoms[1], atoms1=restr_atoms[2] + ), + "InverseBondRestraintForce", + lambda force: ( + force.getBondParameters(0)[:2] == [start_mol1 + 1, start_mol1 + 2] + ), + ), + ( + "morse_potential", + lambda system: sr.restraints.morse_potential( + system, + atoms0=restr_atoms[0], + atoms1=restr_atoms[1], + r0="1.5A", + k="100 kcal mol-1 A-2", + de="25 kcal mol-1", + auto_parametrise=False, + )[0], + "MorsePotentialRestraintForce", + lambda force: ( + force.getBondParameters(0)[:2] == [start_mol1 + 0, start_mol1 + 1] + ), + ), + ( + "positional", + lambda system: sr.restraints.positional(system, atoms=restr_atoms[0]), + "PositionalRestraintForce", + lambda force: ( + force.getBondParameters(0)[0] == start_mol1 + 0 + and force.getBondParameters(0)[1] >= base_particles + or force.getBondParameters(0)[1] == start_mol1 + 0 + and force.getBondParameters(0)[0] >= base_particles + ), + ), + ( + "angle", + lambda system: sr.restraints.angle(system, atoms=restr_atoms[:3]), + "AngleRestraintForce", + lambda force: ( + force.getAngleParameters(0)[:3] + == [start_mol1 + 0, start_mol1 + 1, start_mol1 + 2] + ), + ), + ( + "dihedral", + lambda system: sr.restraints.dihedral(system, atoms=restr_atoms[:4]), + "TorsionRestraintForce", + lambda force: ( + force.getTorsionParameters(0)[:4] + == [start_mol1 + 0, start_mol1 + 1, start_mol1 + 2, start_mol1 + 3] + ), + ), + ( + "boresch", + lambda system: sr.restraints.boresch( + system, + receptor=restr_atoms[0:3], + ligand=restr_atoms[3:6], + ), + "BoreschRestraintForce", + lambda force: ( + list(force.getBondParameters(0)[0]) + == [ + start_mol1 + 0, + start_mol1 + 1, + start_mol1 + 2, + start_mol1 + 3, + start_mol1 + 4, + start_mol1 + 5, + ] + ), + ), + ] + + for restraint_name, build_restraints, force_name, validate in restraint_cases: + restraints = build_restraints(mols) + d = mols.dynamics(restraints=restraints, platform=openmm_platform) + system = d.context().getSystem() + force = next( + force for force in system.getForces() if force.getName() == force_name + ) + + assert validate(force), f"Incorrect atom mapping for {restraint_name}" diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index f87fb6580..b5d8233cf 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -69,6 +69,9 @@ def __init__( # place the coordinates and velocities into the context set_openmm_coordinates_and_velocities(self, metadata) + # set the positions of any virtual sites + self.computeVirtualSites() + # Check for a REST2 scaling factor. if map.specified("rest2_scale"): try: diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp index 23d3f77fc..fb3d4fc35 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp @@ -1,4 +1,3 @@ - #include "openmmmolecule.h" #include "SireMol/atomcharges.h" @@ -67,6 +66,9 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol, return; } + + // Set up virtual site properties + bool is_perturbable = false; bool swap_end_states = false; @@ -89,6 +91,26 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol, ffinfo = mol.property(map["forcefield"]).asA(); } + if (mol.hasProperty("n_virtual_sites") and mol.property("n_virtual_sites").asAnInteger() > 0) + { + this->has_vs = true; + this->vs_parents = mol.property("parents").asA(); + this->vs_properties = mol.property("virtual_sites").asA(); + this->n_vs = mol.property("n_virtual_sites").asAnInteger(); + if (is_perturbable) + { + this->vs_charges = mol.property("vs_charges0").asAnArray(); + } + else + { + this->vs_charges = mol.property("vs_charges").asAnArray(); + } + } + else + { + this->has_vs = false; + } + if (map.specified("constraint")) { const auto c = map["constraint"].source().toLower().simplified().replace("_", "-"); @@ -261,7 +283,7 @@ OpenMMMolecule::OpenMMMolecule(const Molecule &mol, "dihedral", "element", "forcefield", "gb_radii", "gb_screening", "improper", "intrascale", "mass", "name", - "parameters", "treechain"}; + "parameters", "treechain", "vs_charges"}; // we can't specialise these globally in case other molecules // are not of amber type @@ -701,7 +723,13 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, const auto params_charges = params.charges(); const auto params_ljs = params.ljs(); - this->cljs = QVector>(nats, boost::make_tuple(0.0, 0.0, 0.0)); + int n_params = nats; + if (this->has_vs) + { + n_params += this->n_vs; + } + this->cljs = QVector>(n_params, boost::make_tuple(0.0, 0.0, 0.0)); + auto cljs_data = cljs.data(); for (int i = 0; i < nats; ++i) @@ -730,6 +758,19 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol, cljs_data[i] = boost::make_tuple(chg, sig, eps); } + if (this->has_vs) + { + auto vs_charges = mol.property(map["vs_charges"]).asAnArray(); + for (int vs = 0; vs < this->n_vs; ++vs) + { + double chg = vs_charges.at(vs).asADouble(); + double sig = 1e-9; + double eps = 0.0; + + cljs_data[nats + vs] = boost::make_tuple(chg, sig, eps); + } + } + this->bond_params.clear(); this->constraints.clear(); this->perturbable_constraints.clear(); @@ -1454,7 +1495,7 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) this->perturbed->alphas = this->alphas; this->perturbed->kappas = this->kappas; - for (int i = 0; i < cljs.count(); ++i) + for (int i = 0; i < this->nAtoms(); ++i) { const auto &clj0 = cljs.at(i); const auto &clj1 = perturbed->cljs.at(i); @@ -1478,6 +1519,19 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) // kappa is 1 for both end states for ghost atoms this->kappas[i] = 1.0; this->perturbed->kappas[i] = 1.0; + + if (this->has_vs) + { + auto atom_vs = this->vs_parents.property(std::to_string(i).c_str()).asAnArray(); + for (int vs = 0; vs < atom_vs.size(); ++vs) + { + int vs_index = this->nAtoms() + atom_vs.at(vs).asAnInteger(); + from_ghost_idxs.insert(vs_index); + this->alphas[vs_index] = 1.0; + this->kappas[vs_index] = 1.0; + this->perturbed->kappas[vs_index] = 1.0; + } + } } } else if (is_ghost(clj1)) @@ -1491,6 +1545,19 @@ void OpenMMMolecule::alignInternals(const PropertyMap &map) // kappa is 1 for both end states for ghost atoms this->kappas[i] = 1.0; this->perturbed->kappas[i] = 1.0; + + if (this->has_vs) + { + auto atom_vs = this->vs_parents.property(std::to_string(i).c_str()).asAnArray(); + for (int vs = 0; vs < atom_vs.size(); ++vs) + { + int vs_index = this->nAtoms() + atom_vs.at(vs).asAnInteger(); + to_ghost_idxs.insert(vs_index); + this->perturbed->alphas[vs_index] = 1.0; + this->kappas[vs_index] = 1.0; + this->perturbed->kappas[vs_index] = 1.0; + } + } } } } @@ -1901,7 +1968,7 @@ void OpenMMMolecule::buildExceptions(const Molecule &mol, if (unbonded_atoms.isEmpty()) unbonded_atoms = QSet(); - const int nats = this->cljs.count(); + const int nats = this->atoms.count(); const auto &nbpairs = mol.property(map["intrascale"]).asA(); @@ -2135,6 +2202,47 @@ void OpenMMMolecule::buildExceptions(const Molecule &mol, } } } + + // Add virtual site exceptions + if (this->has_vs) + { + int n_exceptions = exception_params.count(); + for (int exc = 0; exc < n_exceptions; ++exc) + { + const auto ¶m = exception_params[exc]; + const auto atom0 = boost::get<0>(param); + const auto atom1 = boost::get<1>(param); + const auto coul_14_scale = boost::get<2>(param); + const auto lj_14_scale = boost::get<3>(param); + + int n_atoms = this->coords.count(); + + auto vs_on_0 = vs_parents.property(std::to_string(atom0).c_str()).asAnArray(); + auto vs_on_1 = vs_parents.property(std::to_string(atom1).c_str()).asAnArray(); + + // Not add_exception because we don't want to add constraints between virtual sites + // Scaling factors are inherited from the parent atom exception + for (int vs0 = 0; vs0 < vs_on_0.size(); ++vs0) + { + int vs0_index = vs_on_0.at(vs0).asAnInteger() + n_atoms; + exception_params.append(boost::make_tuple(vs0_index, atom1, coul_14_scale, lj_14_scale)); + for (int vs1 = 0; vs1 < vs_on_1.size(); ++vs1) + { + int vs1_index = vs_on_1.at(vs1).asAnInteger() + n_atoms; + exception_params.append(boost::make_tuple(vs0_index, vs1_index, coul_14_scale, lj_14_scale)); + } + } + + // We need an additional loop for the case where atom 0 has no virtual sites + for (int vs1 = 0; vs1 < vs_on_1.size(); ++vs1) + { + int vs1_index = vs_on_1.at(vs1).asAnInteger() + n_atoms; + exception_params.append(boost::make_tuple(atom0, vs1_index, coul_14_scale, lj_14_scale)); + } + } + } + + } void OpenMMMolecule::copyInCoordsAndVelocities(OpenMM::Vec3 *c, OpenMM::Vec3 *v) const diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h index dd5481efb..79af7a268 100644 --- a/wrapper/Convert/SireOpenMM/openmmmolecule.h +++ b/wrapper/Convert/SireOpenMM/openmmmolecule.h @@ -224,6 +224,13 @@ namespace SireOpenMM /** The starting index of the first OpenMM atom in the original Sire system. */ int start_atom_idx; + /** Virtual site properties */ + bool has_vs; + int n_vs; + SireBase::Properties vs_parents; + SireBase::PropertyList vs_charges; + SireBase::Properties vs_properties; + private: void constructFromAmber(const SireMol::Molecule &mol, const SireMM::AmberParams ¶ms, diff --git a/wrapper/Convert/SireOpenMM/sire_openmm.cpp b/wrapper/Convert/SireOpenMM/sire_openmm.cpp index 36794bb13..efe5293b8 100644 --- a/wrapper/Convert/SireOpenMM/sire_openmm.cpp +++ b/wrapper/Convert/SireOpenMM/sire_openmm.cpp @@ -495,6 +495,10 @@ namespace SireOpenMM { offsets[i] = offset; offset += mols[i].nAtoms(); + if (mols[i].hasProperty("n_virtual_sites")) + { + offset += mols[i].property("n_virtual_sites").asAnInteger(); + } } const auto offsets_data = offsets.constData(); @@ -822,6 +826,10 @@ namespace SireOpenMM { offsets[i] = offset; offset += mols[i].nAtoms(); + if (mols[i].hasProperty("n_virtual_sites")) + { + offset += mols[i].property("n_virtual_sites").asAnInteger(); + } } const auto offsets_data = offsets.constData(); diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp index fdb09c3e7..2cd400f13 100644 --- a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp +++ b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp @@ -1,4 +1,3 @@ - #include "sire_openmm.h" #include @@ -21,13 +20,13 @@ #include "SireMol/moleditor.h" #include "SireMM/amberparams.h" -#include "SireMM/cmapparameter.h" #include "SireMM/anglerestraints.h" #include "SireMM/atomljs.h" #include "SireMM/bondrestraints.h" -#include "SireMM/inversebondrestraints.h" #include "SireMM/boreschrestraints.h" +#include "SireMM/cmapparameter.h" #include "SireMM/dihedralrestraints.h" +#include "SireMM/inversebondrestraints.h" #include "SireMM/morsepotentialrestraints.h" #include "SireMM/positionalrestraints.h" #include "SireMM/rmsdrestraints.h" @@ -72,7 +71,7 @@ using namespace SireOpenMM; */ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -159,8 +158,8 @@ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, for (int i = 0; i < 3; ++i) { - particles[i] = restraint.receptorAtoms()[i]; - particles[i + 3] = restraint.ligandAtoms()[i]; + particles[i] = real_atoms[restraint.receptorAtoms()[i]]; + particles[i + 3] = real_atoms[restraint.ligandAtoms()[i]]; if (particles[i] < 0 or particles[i] >= natoms or particles[i + 3] < 0 or particles[i + 3] >= natoms) @@ -199,7 +198,7 @@ void _add_boresch_restraints(const SireMM::BoreschRestraints &restraints, */ void _add_bond_restraints(const SireMM::BondRestraints &restraints, OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -240,8 +239,8 @@ void _add_bond_restraints(const SireMM::BondRestraints &restraints, for (const auto &restraint : atom_restraints) { - int atom0_index = restraint.atom0(); - int atom1_index = restraint.atom1(); + int atom0_index = real_atoms[restraint.atom0()]; + int atom1_index = real_atoms[restraint.atom1()]; if (atom0_index < 0 or atom0_index >= natoms) throw SireError::invalid_index(QObject::tr( @@ -266,8 +265,8 @@ void _add_bond_restraints(const SireMM::BondRestraints &restraints, } void _add_inverse_bond_restraints(const SireMM::InverseBondRestraints &restraints, - OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -280,10 +279,10 @@ void _add_inverse_bond_restraints(const SireMM::InverseBondRestraints &restraint } const auto energy_expression = QString( - "rho*k*delta*delta*step;" - "delta=(r-r0);" - "step=max(0,min(1,(r0 - r)))") - .toStdString(); + "rho*k*delta*delta*step;" + "delta=(r-r0);" + "step=max(0,min(1,(r0 - r)))") + .toStdString(); auto *restraintff = new OpenMM::CustomBondForce(energy_expression); restraintff->setName("InverseBondRestraintForce"); @@ -296,7 +295,7 @@ void _add_inverse_bond_restraints(const SireMM::InverseBondRestraints &restraint restraintff->setForceGroup(force_group_counter); lambda_lever.addRestraintIndex(restraints.name(), - system.addForce(restraintff)); + system.addForce(restraintff)); lambda_lever.setRestraintForceGroup(restraints.name(), force_group_counter++); const auto atom_restraints = restraints.atomRestraints(); @@ -310,22 +309,22 @@ void _add_inverse_bond_restraints(const SireMM::InverseBondRestraints &restraint for (const auto &restraint : atom_restraints) { - int atom0_index = restraint.atom0(); - int atom1_index = restraint.atom1(); + int atom0_index = real_atoms[restraint.atom0()]; + int atom1_index = real_atoms[restraint.atom1()]; if (atom0_index < 0 or atom0_index >= natoms) - throw SireError::invalid_index(QObject::tr( - "Invalid particle index! %1 from %2") - .arg(atom0_index) - .arg(natoms), - CODELOC); + throw SireError::invalid_index(QObject::tr( + "Invalid particle index! %1 from %2") + .arg(atom0_index) + .arg(natoms), + CODELOC); if (atom1_index < 0 or atom1_index >= natoms) - throw SireError::invalid_index(QObject::tr( - "Invalid particle index! %1 from %2") - .arg(atom1_index) - .arg(natoms), - CODELOC); + throw SireError::invalid_index(QObject::tr( + "Invalid particle index! %1 from %2") + .arg(atom1_index) + .arg(natoms), + CODELOC); custom_params[0] = 1.0; // rho - always equal to 1 (scaled by lever) custom_params[1] = restraint.k().value() * internal_to_k; // k @@ -339,17 +338,17 @@ void _add_inverse_bond_restraints(const SireMM::InverseBondRestraints &restraint * system, which is acted on by the passed LambdaLever. */ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &restraints, - OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) - return; + return; if (restraints.hasCentroidRestraints()) { throw SireError::unsupported(QObject::tr( - "Centroid bond restraints aren't yet supported..."), - CODELOC); + "Centroid bond restraints aren't yet supported..."), + CODELOC); } // energy expression of a harmonic bond potential, scaled by rho @@ -358,11 +357,10 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res // 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(); - + "rho*e_restraint;" + "e_restraint=de*(1-exp(-sqrt(k/(2*de))*delta))^2;" + "delta=(r-r0)") + .toStdString(); auto *restraintff = new OpenMM::CustomBondForce(energy_expression); restraintff->setName("MorsePotentialRestraintForce"); @@ -376,7 +374,7 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res restraintff->setForceGroup(force_group_counter); lambda_lever.addRestraintIndex(restraints.name(), - system.addForce(restraintff)); + system.addForce(restraintff)); lambda_lever.setRestraintForceGroup(restraints.name(), force_group_counter++); const auto atom_restraints = restraints.atomRestraints(); @@ -388,22 +386,22 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res for (const auto &restraint : atom_restraints) { - int atom0_index = restraint.atom0(); - int atom1_index = restraint.atom1(); + int atom0_index = real_atoms[restraint.atom0()]; + int atom1_index = real_atoms[restraint.atom1()]; if (atom0_index < 0 or atom0_index >= natoms) - throw SireError::invalid_index(QObject::tr( - "Invalid particle index! %1 from %2") - .arg(atom0_index) - .arg(natoms), - CODELOC); + throw SireError::invalid_index(QObject::tr( + "Invalid particle index! %1 from %2") + .arg(atom0_index) + .arg(natoms), + CODELOC); if (atom1_index < 0 or atom1_index >= natoms) - throw SireError::invalid_index(QObject::tr( - "Invalid particle index! %1 from %2") - .arg(atom1_index) - .arg(natoms), - CODELOC); + throw SireError::invalid_index(QObject::tr( + "Invalid particle index! %1 from %2") + .arg(atom1_index) + .arg(natoms), + CODELOC); custom_params[0] = 1.0; // rho - always equal to 1 (scaled by lever) custom_params[1] = restraint.k().value() * internal_to_k; // k @@ -423,7 +421,7 @@ void _add_morse_potential_restraints(const SireMM::MorsePotentialRestraints &res void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, OpenMM::System &system, LambdaLever &lambda_lever, std::vector &anchor_coords, - int natoms, int &force_group_counter) + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -484,7 +482,7 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, // we need to add all of the positions as anchor particles for (const auto &restraint : atom_restraints) { - int atom_index = restraint.atom(); + int atom_index = real_atoms[restraint.atom()]; if (atom_index < 0 or atom_index >= natoms) throw SireError::invalid_index(QObject::tr( @@ -558,8 +556,8 @@ void _add_positional_restraints(const SireMM::PositionalRestraints &restraints, } void _add_rmsd_restraints(const SireMM::RMSDRestraints &restraints, - OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + OpenMM::System &system, LambdaLever &lambda_lever, + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -579,16 +577,19 @@ void _add_rmsd_restraints(const SireMM::RMSDRestraints &restraints, }; // Count the number of existing RMSD forces in the system - std::vector forces; + std::vector forces; - for (int i = 0; i < system.getNumForces(); ++i) { + for (int i = 0; i < system.getNumForces(); ++i) + { forces.push_back(&system.getForce(i)); } int n_CVForces = 0; - for (auto* force : forces) { - if (dynamic_cast(force)) { + for (auto *force : forces) + { + if (dynamic_cast(force)) + { n_CVForces++; } } @@ -606,7 +607,7 @@ void _add_rmsd_restraints(const SireMM::RMSDRestraints &restraints, // energy expression of a flat-bottom well potential, scaled by rho const auto energy_expression = rho_unique + "*" + k_unique + "*step(delta)*delta*delta;" + - "delta=(" + rmsd_unique + "-" + rmsd_b_unique + ")"; + "delta=(" + rmsd_unique + "-" + rmsd_b_unique + ")"; double k = restraint.k().value() * internal_to_k; double r0 = restraint.r0().value() * internal_to_nm; @@ -618,7 +619,7 @@ void _add_rmsd_restraints(const SireMM::RMSDRestraints &restraints, for (int i = 0; i < n_particles; ++i) { - particles[i] = restraint.atoms()[i]; + particles[i] = real_atoms[restraint.atoms()[i]]; } // Extract reference positions and convert to correct units @@ -650,7 +651,7 @@ void _add_rmsd_restraints(const SireMM::RMSDRestraints &restraints, } restraintff->setForceGroup(grp); lambda_lever.addRestraintIndex(restraints.name(), - system.addForce(restraintff)); + system.addForce(restraintff)); lambda_lever.setRestraintForceGroup(restraints.name(), grp); // Update the counter for number of CustomCVForces @@ -664,7 +665,7 @@ void _add_rmsd_restraints(const SireMM::RMSDRestraints &restraints, */ void _add_angle_restraints(const SireMM::AngleRestraints &restraints, OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -707,7 +708,7 @@ void _add_angle_restraints(const SireMM::AngleRestraints &restraints, for (int i = 0; i < 3; ++i) { - particles[i] = restraint.atoms()[i]; + particles[i] = real_atoms[restraint.atoms()[i]]; } std::vector parameters; @@ -724,7 +725,7 @@ void _add_angle_restraints(const SireMM::AngleRestraints &restraints, void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, OpenMM::System &system, LambdaLever &lambda_lever, - int natoms, int &force_group_counter) + int natoms, QVector &real_atoms, int &force_group_counter) { if (restraints.isEmpty()) return; @@ -772,7 +773,7 @@ void _add_dihedral_restraints(const SireMM::DihedralRestraints &restraints, for (int i = 0; i < 4; ++i) { - particles[i] = restraint.atoms()[i]; + particles[i] = real_atoms[restraint.atoms()[i]]; } std::vector parameters; @@ -1100,7 +1101,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, start_atom_index[0] = 0; for (int i = 1; i < nmols; ++i) { - start_atom_index[i] = start_atom_index[i-1] + mols[i-1].nAtoms(); + start_atom_index[i] = start_atom_index[i - 1] + mols[i - 1].nAtoms(); } if (SireBase::should_run_in_parallel(nmols, map)) @@ -1608,29 +1609,39 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, QVector ghost_atoms; QVector non_ghost_atoms; + // indices for real atoms (i.e. not virtual sites) + // required so that restraints are applied to the correct particles + QVector real_atoms; + // count the number of atoms and ghost atoms int n_atoms = 0; int n_ghost_atoms = 0; + int n_vs = 0; for (int i = 0; i < nmols; ++i) { const auto &mol = openmm_mols_data[i]; n_atoms += mol.nAtoms(); n_ghost_atoms += mol.nGhostAtoms(); + if (mol.has_vs) + { + n_vs += mol.n_vs; + } } // there's probably lots of optimisations we can make if the // number of ghost atoms is zero... - ghost_atoms.reserve(n_ghost_atoms); - non_ghost_atoms.reserve(n_atoms - n_ghost_atoms); + // Making sure there is space in all arrays for virtual sites to be ghosts + ghost_atoms.reserve(n_ghost_atoms + n_vs); + non_ghost_atoms.reserve(n_atoms - n_ghost_atoms + n_vs); // the set of all ghost atoms, with the value // indicating if this is a from-ghost (true) or // a to-ghost (false) QVector from_ghost_idxs; QVector to_ghost_idxs; - from_ghost_idxs.reserve(n_ghost_atoms); - to_ghost_idxs.reserve(n_ghost_atoms); + from_ghost_idxs.reserve(n_ghost_atoms + n_vs); + to_ghost_idxs.reserve(n_ghost_atoms + n_vs); // map from CMAPParameter to OpenMM map index for non-perturbable molecules // (enables sharing identical grids across molecules) @@ -1770,6 +1781,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ghost_ghostff->addParticle(custom_params); ghost_nonghostff->addParticle(custom_params); + real_atoms.append(atom_index); + if (is_from_ghost or is_to_ghost) { // this is a ghost atom! We need to record this @@ -1829,6 +1842,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, // Add the particle to the standard CLJ forcefield cljff->addParticle(boost::get<0>(clj), boost::get<1>(clj), boost::get<2>(clj)); + real_atoms.append(atom_index); + // We need to add this molecule to the ghost and ghost // forcefields if there are any perturbable molecules if (any_perturbable) @@ -1849,6 +1864,101 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } } } + if (mol.has_vs) + { + int start_vs = start_index + mol.molinfo.nAtoms(); + + for (int k = 0; k < mol.n_vs; ++k) + { + SireBase::Properties vs_params = mol.vs_properties.property(std::to_string(k).c_str()).asA(); + // Add parameters to system + // Virtual sites with non-zero LJ interactions are not supported + const int atom_index = start_vs + k; + system.addParticle(0.0); + + // Calculate virtual site parameters + SireBase::PropertyList indices = vs_params.property("vs_indices").asAnArray(); + std::vector indices_vec = {}; + for (int a = 0; a < indices.size(); ++a) + { + indices_vec.push_back(indices.at(a).asAnInteger() + start_index); + } + int parent_idx = indices.at(0).asAnInteger(); + + SireBase::PropertyList ows = vs_params.property("vs_ows").asAnArray(); + std::vector ows_vec = {}; + for (int a = 0; a < ows.size(); ++a) + { + ows_vec.push_back(ows.at(a).asADouble()); + } + + SireBase::PropertyList xs = vs_params.property("vs_xs").asAnArray(); + std::vector xs_vec = {}; + for (int a = 0; a < xs.size(); ++a) + { + xs_vec.push_back(xs.at(a).asADouble()); + } + + SireBase::PropertyList ys = vs_params.property("vs_ys").asAnArray(); + std::vector ys_vec = {}; + for (int a = 0; a < ys.size(); ++a) + { + ys_vec.push_back(ys.at(a).asADouble()); + } + + SireBase::PropertyList local = vs_params.property("vs_local").asAnArray(); + OpenMM::Vec3 local_vec = {local.at(0).asADouble(), local.at(1).asADouble(), local.at(2).asADouble()}; + + OpenMM::LocalCoordinatesSite *new_vs = new OpenMM::LocalCoordinatesSite(indices_vec, ows_vec, xs_vec, ys_vec, local_vec); + system.setVirtualSite(atom_index, new_vs); + + // Add to forcefield, depending on whether the system is perturbable + // Note that VS with LJ parameters are currently not supported, so epsilon and sigma are hard-coded to 0 in all cases + double vs_charge = mol.vs_charges.at(k).asADouble(); + cljff->addParticle(vs_charge, 1.0, 0.0); + + if (any_perturbable and mol.isPerturbable()) + { + // reduced_q + custom_params[0] = vs_charge; + // half_sigma + custom_params[1] = 1.0; + // two_sqrt_epsilon + custom_params[2] = 0.0; + // alpha + custom_params[3] = alphas_data[mol.molinfo.nAtoms() + k]; + // kappa + custom_params[4] = kappas_data[mol.molinfo.nAtoms() + k]; + + ghost_ghostff->addParticle(custom_params); + ghost_nonghostff->addParticle(custom_params); + + const bool vs_to_ghost = mol.to_ghost_idxs.contains(parent_idx); + const bool vs_from_ghost = mol.from_ghost_idxs.contains(parent_idx); + + // Append virtual sites to ghost atom list + + if (vs_to_ghost) + { + ghost_atoms.append(atom_index); + to_ghost_idxs.append(atom_index); + } + else if (vs_from_ghost) + { + ghost_atoms.append(atom_index); + from_ghost_idxs.append(atom_index); + } + } + else if (any_perturbable) + { + // Add to ghost FFs if necessary + custom_params = {vs_charge, 1.0, 0.0, 0.0, 0.0}; + ghost_ghostff->addParticle(custom_params); + ghost_nonghostff->addParticle(custom_params); + non_ghost_atoms.append(atom_index); + } + } + } // Register virtual sites (OPC, TIP4P, TIP5P, …). // setVirtualSite() must be called after all particles have been added @@ -1856,9 +1966,9 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, for (const auto &vs : mol.virtual_sites) { const int vsite_atom = vs.vsite_idx + start_index; - const int p1 = vs.p1_idx + start_index; - const int p2 = vs.p2_idx + start_index; - const int p3 = vs.p3_idx + start_index; + const int p1 = vs.p1_idx + start_index; + const int p2 = vs.p2_idx + start_index; + const int p3 = vs.p3_idx + start_index; if (vs.type == VirtualSiteInfo::ThreeParticleAverage) { @@ -2011,6 +2121,10 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, } start_index += mol.masses.count(); + if (mol.has_vs) + { + start_index += mol.n_vs; + } } /// Finally tell the ghost forcefields about the ghost and non-ghost @@ -2080,6 +2194,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, { int start_index = start_indexes[i]; const auto &mol = openmm_mols_data[i]; + auto cljs_data = mol.cljs.constData(); QVector> exception_idxs; QVector constraint_idxs; @@ -2213,10 +2328,47 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, { auto pert_idx = idx_to_pert_idx.value(i, openmm_mols.count() + 1); lambda_lever.setExceptionIndicies(pert_idx, - "clj", exception_idxs); + "clj", exception_idxs); lambda_lever.setConstraintIndicies(pert_idx, constraint_idxs); } + + // Exclusions/exceptions between virtual sites on the same atom, and with the parent atom + if (mol.has_vs) + { + int n_atom_vs; + int vs_start = start_index + mol.nAtoms(); + for (int a = 0; a < mol.nAtoms(); ++a) + { + SireBase::PropertyList atom_vs = mol.vs_parents.property(std::to_string(a).c_str()).asAnArray(); + n_atom_vs = atom_vs.size(); + for (int v0 = 0; v0 < atom_vs.size(); ++v0) + { + int vs0_index = vs_start + atom_vs.at(v0).asAnInteger(); + cljff->addException(vs0_index, start_index + a, + 0.0, 1, + 0, false); + if (ghost_ghostff != 0) + { + ghost_ghostff->addExclusion(vs0_index, start_index + a); + ghost_nonghostff->addExclusion(vs0_index, start_index + a); + } + + for (int v1 = v0 + 1; v1 < atom_vs.size(); ++v1) + { + int vs1_index = vs_start + atom_vs.at(v1).asAnInteger(); + cljff->addException(vs0_index, vs1_index, + 0.0, 1, + 0, false); + if (ghost_ghostff != 0) + { + ghost_ghostff->addExclusion(vs0_index, vs1_index); + ghost_nonghostff->addExclusion(vs0_index, vs1_index); + } + } + } + } + } } // go through all of the ghost atoms and exclude interactions @@ -2244,7 +2396,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, ghost_ghostff->addExclusion(from_ghost_idx, to_ghost_idx); ghost_nonghostff->addExclusion(from_ghost_idx, to_ghost_idx); cljff->addException(from_ghost_idx, to_ghost_idx, - 0.0, 1e-9, 1e-9, true); + 0.0, 1e-9, 1e-9, false); } } } @@ -2290,43 +2442,43 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, if (prop.read().isA()) { _add_dihedral_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_angle_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_positional_restraints(prop.read().asA(), - system, lambda_lever, anchor_coords, start_index, + system, lambda_lever, anchor_coords, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_morse_potential_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_rmsd_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_bond_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_inverse_bond_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } else if (prop.read().isA()) { _add_boresch_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } } } @@ -2351,7 +2503,7 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, if (prop.read().isA()) { _add_inverse_bond_restraints(prop.read().asA(), - system, lambda_lever, start_index, force_group_counter); + system, lambda_lever, start_index, real_atoms, force_group_counter); } } } @@ -2390,6 +2542,16 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, const auto &mol = openmm_mols_data[i]; mol.copyInCoordsAndVelocities(coords_data + start_index, vels_data + start_index); + if (mol.has_vs) + { + // Initiate all VS with zero coords, as they will need to be + // calculated in the openmm context anyway + for (int vs = 0; vs < mol.n_vs; ++vs) + { + coords_data[start_index+mol.nAtoms()+vs] = OpenMM::Vec3(0,0,0); + vels_data[start_index+mol.nAtoms()+vs] = OpenMM::Vec3(0,0,0); + } + } } }); } else @@ -2400,6 +2562,14 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system, const auto &mol = openmm_mols_data[i]; mol.copyInCoordsAndVelocities(coords_data + start_index, vels_data + start_index); + if (mol.has_vs) + { + for (int vs = 0; vs < mol.n_vs; ++vs) + { + coords_data[start_index + mol.nAtoms() + vs] = OpenMM::Vec3(0, 0, 0); + vels_data[start_index + mol.nAtoms() + vs] = OpenMM::Vec3(0, 0, 0); + } + } } }