From 871e2e2b402dd121a359a132d02116e3182340ae Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 May 2026 10:26:37 +0100 Subject: [PATCH 1/4] Add support for a QM/MM switching function. --- doc/source/changelog.rst | 2 + src/sire/qm/__init__.py | 21 + src/sire/qm/_emle.py | 24 +- .../Convert/SireOpenMM/PyQMEngine.pypp.cpp | 441 ++++++-------- .../Convert/SireOpenMM/TorchQMEngine.pypp.cpp | 424 ++++++------- wrapper/Convert/SireOpenMM/pyqm.cpp | 564 ++++++++++++------ wrapper/Convert/SireOpenMM/pyqm.h | 157 +++-- wrapper/Convert/SireOpenMM/torchqm.cpp | 555 ++++++++++------- wrapper/Convert/SireOpenMM/torchqm.h | 66 +- 9 files changed, 1296 insertions(+), 958 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 69834fabf..c71a60ee9 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -84,6 +84,8 @@ organisation on `GitHub `__. ``reversed.morph(force, lever, initial, final, λ)`` equals ``original.morph(force, lever, final, initial, 1-λ)``. +* Add support for using a switching function for QM/MM simulations. + `2025.4.0 `__ - February 2026 --------------------------------------------------------------------------------------------- diff --git a/src/sire/qm/__init__.py b/src/sire/qm/__init__.py index 6e02245f0..090c31f22 100644 --- a/src/sire/qm/__init__.py +++ b/src/sire/qm/__init__.py @@ -19,6 +19,7 @@ def create_engine( neighbour_list_frequency=0, mechanical_embedding=False, redistribute_charge=True, + switch_width=0.2, map=None, ): """ @@ -74,6 +75,13 @@ def create_engine( molecule, the excess charge is redistributed over the MM atoms within the residues of the QM region. + switch_width : float, optional, default=0.2 + The width of the switching region as a fraction of the cutoff (0 to 1). + A quintic switching function is applied to the MM charges over the last + ``switch_width * cutoff`` angstroms before the cutoff, smoothly scaling + them to zero. Set to 0 or None to disable switching (not recommended + for production MD with charged ML regions). + Returns ------- @@ -133,6 +141,15 @@ def create_engine( if not isinstance(redistribute_charge, bool): raise TypeError("'redistribute_charge' must be of type 'bool'") + if switch_width is None: + switch_width = 0.0 + if not isinstance(switch_width, (int, float)): + raise TypeError("'switch_width' must be of type 'int' or 'float'") + switch_width = float(switch_width) + if switch_width < 0.0 or switch_width > 1.0: + raise ValueError("'switch_width' must be between 0 and 1") + use_switch = switch_width > 0.0 + if map is not None: if not isinstance(map, dict): raise TypeError("'map' must be of type 'dict'") @@ -147,6 +164,10 @@ def create_engine( mechanical_embedding, ) + # Configure the switching function. + engine.set_switch_width(switch_width) + engine.set_use_switch(use_switch) + from ._utils import ( _check_charge, _create_qm_mol_to_atoms, diff --git a/src/sire/qm/_emle.py b/src/sire/qm/_emle.py index 7e6b0eec6..0e619ae46 100644 --- a/src/sire/qm/_emle.py +++ b/src/sire/qm/_emle.py @@ -108,6 +108,7 @@ def emle( cutoff="7.5A", neighbour_list_frequency=0, redistribute_charge=True, + switch_width=0.2, map=None, ): """ @@ -143,6 +144,13 @@ def emle( molecule, the excess charge is redistributed over the MM atoms within the residues of the QM region. + switch_width : float, optional, default=0.2 + The width of the switching region as a fraction of the cutoff (0 to 1). + A quintic switching function is applied to the MM charges over the last + ``switch_width * cutoff`` angstroms before the cutoff, smoothly scaling + them to zero. Set to 0 or None to disable switching (not recommended + for production MD with charged ML regions). + Returns ------- @@ -159,7 +167,6 @@ def emle( try: import torch as _torch - from emle.models import EMLE as _EMLE has_model = True except: @@ -227,6 +234,15 @@ def emle( if not isinstance(redistribute_charge, bool): raise TypeError("'redistribute_charge' must be of type 'bool'") + if switch_width is None: + switch_width = 0.0 + if not isinstance(switch_width, (int, float)): + raise TypeError("'switch_width' must be of type 'int' or 'float'") + switch_width = float(switch_width) + if switch_width < 0.0 or switch_width > 1.0: + raise ValueError("'switch_width' must be between 0 and 1") + use_switch = switch_width > 0.0 + if map is not None: if not isinstance(map, dict): raise TypeError("'map' must be of type 'dict'") @@ -246,7 +262,7 @@ def emle( # Create an engine from an EMLE model. else: try: - from emle.models import EMLE as _EMLE + pass except: raise ImportError( "Could not import emle.models. Please reinstall emle-engine and try again." @@ -276,6 +292,10 @@ def emle( except Exception as e: raise ValueError("Unable to create a TorchEMLEEngine: " + str(e)) + # Configure the switching function. + engine.set_switch_width(switch_width) + engine.set_use_switch(use_switch) + from ._utils import ( _check_charge, _create_qm_mol_to_atoms, diff --git a/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp b/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp index ee197ff51..d86fcb7f1 100644 --- a/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp @@ -2,8 +2,8 @@ // (C) Christopher Woods, GPL >= 3 License -#include "boost/python.hpp" #include "PyQMEngine.pypp.hpp" +#include "boost/python.hpp" namespace bp = boost::python; @@ -51,313 +51,242 @@ namespace bp = boost::python; #include -SireOpenMM::PyQMEngine __copy__(const SireOpenMM::PyQMEngine &other){ return SireOpenMM::PyQMEngine(other); } +SireOpenMM::PyQMEngine __copy__(const SireOpenMM::PyQMEngine &other) { return SireOpenMM::PyQMEngine(other); } #include "Helpers/str.hpp" #include "Helpers/release_gil_policy.hpp" -void register_PyQMEngine_class(){ +void register_PyQMEngine_class() +{ { //::SireOpenMM::PyQMEngine - typedef bp::class_< SireOpenMM::PyQMEngine, bp::bases< SireBase::Property, SireOpenMM::QMEngine > > PyQMEngine_exposer_t; - PyQMEngine_exposer_t PyQMEngine_exposer = PyQMEngine_exposer_t( "PyQMEngine", "", bp::init< >("Default constructor.") ); - bp::scope PyQMEngine_scope( PyQMEngine_exposer ); - PyQMEngine_exposer.def( bp::init< bp::api::object, bp::optional< QString, SireUnits::Dimension::Length, int, bool, double > >(( bp::arg("arg0"), bp::arg("method")="", bp::arg("cutoff")=7.5 * SireUnits::angstrom, bp::arg("neighbour_list_frequency")=(int)(0), bp::arg("is_mechanical")=(bool)(false), bp::arg("lambda")=1. ), "Constructor\nPar:am py_object\nA Python object.\n\nPar:am name\nThe name of the callback method. If empty, then the object is\nassumed to be a callable.\n\nPar:am cutoff\nThe ML cutoff distance.\n\nPar:am neighbour_list_frequency\nThe frequency at which the neighbour list is updated. (Number of steps.)\nIf zero, then no neighbour list is used.\n\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n\nPar:am lambda\nThe lambda weighting factor. This can be used to interpolate between\npotentials for end-state correction calculations.\n") ); - PyQMEngine_exposer.def( bp::init< SireOpenMM::PyQMEngine const & >(( bp::arg("other") ), "Copy constructor.") ); + typedef bp::class_> PyQMEngine_exposer_t; + PyQMEngine_exposer_t PyQMEngine_exposer = PyQMEngine_exposer_t("PyQMEngine", "", bp::init<>("Default constructor.")); + bp::scope PyQMEngine_scope(PyQMEngine_exposer); + PyQMEngine_exposer.def(bp::init>((bp::arg("arg0"), bp::arg("method") = "", bp::arg("cutoff") = 7.5 * SireUnits::angstrom, bp::arg("neighbour_list_frequency") = (int)(0), bp::arg("is_mechanical") = (bool)(false), bp::arg("lambda") = 1.), "Constructor\nPar:am py_object\nA Python object.\n\nPar:am name\nThe name of the callback method. If empty, then the object is\nassumed to be a callable.\n\nPar:am cutoff\nThe ML cutoff distance.\n\nPar:am neighbour_list_frequency\nThe frequency at which the neighbour list is updated. (Number of steps.)\nIf zero, then no neighbour list is used.\n\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n\nPar:am lambda\nThe lambda weighting factor. This can be used to interpolate between\npotentials for end-state correction calculations.\n")); + PyQMEngine_exposer.def(bp::init((bp::arg("other")), "Copy constructor.")); { //::SireOpenMM::PyQMEngine::call - - typedef ::boost::tuples::tuple< double, QVector< QVector< double > >, QVector< QVector< double > >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::PyQMEngine::*call_function_type)( ::QVector< int >,::QVector< double >,::QVector< QVector< double > >,::QVector< QVector< double > >,::QVector< QVector< double > >,::QVector < int > ) const; - call_function_type call_function_value( &::SireOpenMM::PyQMEngine::call ); - - PyQMEngine_exposer.def( - "call" - , call_function_value - , ( bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm") ) - , bp::release_gil_policy() - , "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of the true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n" ); - + + typedef ::boost::tuples::tuple>, QVector>, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::PyQMEngine::*call_function_type)(::QVector, ::QVector, ::QVector>, ::QVector>, ::QVector>, ::QVector) const; + call_function_type call_function_value(&::SireOpenMM::PyQMEngine::call); + + PyQMEngine_exposer.def( + "call", call_function_value, (bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm")), bp::release_gil_policy(), "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of the true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n"); } { //::SireOpenMM::PyQMEngine::getAtoms - - typedef ::QVector< int > ( ::SireOpenMM::PyQMEngine::*getAtoms_function_type)( ) const; - getAtoms_function_type getAtoms_function_value( &::SireOpenMM::PyQMEngine::getAtoms ); - - PyQMEngine_exposer.def( - "getAtoms" - , getAtoms_function_value - , bp::release_gil_policy() - , "Get the indices of the atoms in the QM region.\nReturn:s\nA vector of atom indices for the QM region.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMEngine::*getAtoms_function_type)() const; + getAtoms_function_type getAtoms_function_value(&::SireOpenMM::PyQMEngine::getAtoms); + + PyQMEngine_exposer.def( + "getAtoms", getAtoms_function_value, bp::release_gil_policy(), "Get the indices of the atoms in the QM region.\nReturn:s\nA vector of atom indices for the QM region.\n"); } { //::SireOpenMM::PyQMEngine::getCallback - - typedef ::SireOpenMM::PyQMCallback ( ::SireOpenMM::PyQMEngine::*getCallback_function_type)( ) const; - getCallback_function_type getCallback_function_value( &::SireOpenMM::PyQMEngine::getCallback ); - - PyQMEngine_exposer.def( - "getCallback" - , getCallback_function_value - , bp::release_gil_policy() - , "Get the callback object.\nReturn:s\nA Python object that contains the callback function.\n" ); - + + typedef ::SireOpenMM::PyQMCallback (::SireOpenMM::PyQMEngine::*getCallback_function_type)() const; + getCallback_function_type getCallback_function_value(&::SireOpenMM::PyQMEngine::getCallback); + + PyQMEngine_exposer.def( + "getCallback", getCallback_function_value, bp::release_gil_policy(), "Get the callback object.\nReturn:s\nA Python object that contains the callback function.\n"); } { //::SireOpenMM::PyQMEngine::getCharges - - typedef ::QVector< double > ( ::SireOpenMM::PyQMEngine::*getCharges_function_type)( ) const; - getCharges_function_type getCharges_function_value( &::SireOpenMM::PyQMEngine::getCharges ); - - PyQMEngine_exposer.def( - "getCharges" - , getCharges_function_value - , bp::release_gil_policy() - , "Get the atomic charges of all atoms in the system.\nReturn:s\nA vector of atomic charges for all atoms in the system.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMEngine::*getCharges_function_type)() const; + getCharges_function_type getCharges_function_value(&::SireOpenMM::PyQMEngine::getCharges); + + PyQMEngine_exposer.def( + "getCharges", getCharges_function_value, bp::release_gil_policy(), "Get the atomic charges of all atoms in the system.\nReturn:s\nA vector of atomic charges for all atoms in the system.\n"); } { //::SireOpenMM::PyQMEngine::getCutoff - - typedef ::SireUnits::Dimension::Length ( ::SireOpenMM::PyQMEngine::*getCutoff_function_type)( ) const; - getCutoff_function_type getCutoff_function_value( &::SireOpenMM::PyQMEngine::getCutoff ); - - PyQMEngine_exposer.def( - "getCutoff" - , getCutoff_function_value - , bp::release_gil_policy() - , "Get the QM cutoff distance.\nReturn:s\nThe QM cutoff distance.\n" ); - + + typedef ::SireUnits::Dimension::Length (::SireOpenMM::PyQMEngine::*getCutoff_function_type)() const; + getCutoff_function_type getCutoff_function_value(&::SireOpenMM::PyQMEngine::getCutoff); + + PyQMEngine_exposer.def( + "getCutoff", getCutoff_function_value, bp::release_gil_policy(), "Get the QM cutoff distance.\nReturn:s\nThe QM cutoff distance.\n"); } { //::SireOpenMM::PyQMEngine::getIsMechanical - - typedef bool ( ::SireOpenMM::PyQMEngine::*getIsMechanical_function_type)( ) const; - getIsMechanical_function_type getIsMechanical_function_value( &::SireOpenMM::PyQMEngine::getIsMechanical ); - - PyQMEngine_exposer.def( - "getIsMechanical" - , getIsMechanical_function_value - , bp::release_gil_policy() - , "Get the mechanical embedding flag.\nReturn:s\nA flag to indicate if mechanical embedding is being used.\n" ); - + + typedef bool (::SireOpenMM::PyQMEngine::*getIsMechanical_function_type)() const; + getIsMechanical_function_type getIsMechanical_function_value(&::SireOpenMM::PyQMEngine::getIsMechanical); + + PyQMEngine_exposer.def( + "getIsMechanical", getIsMechanical_function_value, bp::release_gil_policy(), "Get the mechanical embedding flag.\nReturn:s\nA flag to indicate if mechanical embedding is being used.\n"); } { //::SireOpenMM::PyQMEngine::getLambda - - typedef double ( ::SireOpenMM::PyQMEngine::*getLambda_function_type)( ) const; - getLambda_function_type getLambda_function_value( &::SireOpenMM::PyQMEngine::getLambda ); - - PyQMEngine_exposer.def( - "getLambda" - , getLambda_function_value - , bp::release_gil_policy() - , "Get the lambda weighting factor.\nReturn:s\nThe lambda weighting factor.\n" ); - + + typedef double (::SireOpenMM::PyQMEngine::*getLambda_function_type)() const; + getLambda_function_type getLambda_function_value(&::SireOpenMM::PyQMEngine::getLambda); + + PyQMEngine_exposer.def( + "getLambda", getLambda_function_value, bp::release_gil_policy(), "Get the lambda weighting factor.\nReturn:s\nThe lambda weighting factor.\n"); } { //::SireOpenMM::PyQMEngine::getLinkAtoms - - typedef ::boost::tuples::tuple< QMap< int, int >, QMap< int, QVector< int > >, QMap< int, double >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::PyQMEngine::*getLinkAtoms_function_type)( ) const; - getLinkAtoms_function_type getLinkAtoms_function_value( &::SireOpenMM::PyQMEngine::getLinkAtoms ); - - PyQMEngine_exposer.def( - "getLinkAtoms" - , getLinkAtoms_function_value - , bp::release_gil_policy() - , "Get the link atoms associated with each QM atom.\nReturn:s\nA tuple containing:\n\nmm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nmm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nbond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n" ); - + + typedef ::boost::tuples::tuple, QMap>, QMap, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::PyQMEngine::*getLinkAtoms_function_type)() const; + getLinkAtoms_function_type getLinkAtoms_function_value(&::SireOpenMM::PyQMEngine::getLinkAtoms); + + PyQMEngine_exposer.def( + "getLinkAtoms", getLinkAtoms_function_value, bp::release_gil_policy(), "Get the link atoms associated with each QM atom.\nReturn:s\nA tuple containing:\n\nmm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nmm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nbond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n"); } { //::SireOpenMM::PyQMEngine::getMM2Atoms - - typedef ::QVector< int > ( ::SireOpenMM::PyQMEngine::*getMM2Atoms_function_type)( ) const; - getMM2Atoms_function_type getMM2Atoms_function_value( &::SireOpenMM::PyQMEngine::getMM2Atoms ); - - PyQMEngine_exposer.def( - "getMM2Atoms" - , getMM2Atoms_function_value - , bp::release_gil_policy() - , "Get the vector of MM2 atoms.\nReturn:s\nA vector of MM2 atom indices.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMEngine::*getMM2Atoms_function_type)() const; + getMM2Atoms_function_type getMM2Atoms_function_value(&::SireOpenMM::PyQMEngine::getMM2Atoms); + + PyQMEngine_exposer.def( + "getMM2Atoms", getMM2Atoms_function_value, bp::release_gil_policy(), "Get the vector of MM2 atoms.\nReturn:s\nA vector of MM2 atom indices.\n"); } { //::SireOpenMM::PyQMEngine::getNeighbourListFrequency - - typedef int ( ::SireOpenMM::PyQMEngine::*getNeighbourListFrequency_function_type)( ) const; - getNeighbourListFrequency_function_type getNeighbourListFrequency_function_value( &::SireOpenMM::PyQMEngine::getNeighbourListFrequency ); - - PyQMEngine_exposer.def( - "getNeighbourListFrequency" - , getNeighbourListFrequency_function_value - , bp::release_gil_policy() - , "Get the neighbour list frequency.\nReturn:s\nThe neighbour list frequency.\n" ); - + + typedef int (::SireOpenMM::PyQMEngine::*getNeighbourListFrequency_function_type)() const; + getNeighbourListFrequency_function_type getNeighbourListFrequency_function_value(&::SireOpenMM::PyQMEngine::getNeighbourListFrequency); + + PyQMEngine_exposer.def( + "getNeighbourListFrequency", getNeighbourListFrequency_function_value, bp::release_gil_policy(), "Get the neighbour list frequency.\nReturn:s\nThe neighbour list frequency.\n"); } { //::SireOpenMM::PyQMEngine::getNumbers - - typedef ::QVector< int > ( ::SireOpenMM::PyQMEngine::*getNumbers_function_type)( ) const; - getNumbers_function_type getNumbers_function_value( &::SireOpenMM::PyQMEngine::getNumbers ); - - PyQMEngine_exposer.def( - "getNumbers" - , getNumbers_function_value - , bp::release_gil_policy() - , "Get the atomic numbers for the atoms in the QM region.\nReturn:s\nA vector of atomic numbers for the atoms in the QM region.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMEngine::*getNumbers_function_type)() const; + getNumbers_function_type getNumbers_function_value(&::SireOpenMM::PyQMEngine::getNumbers); + + PyQMEngine_exposer.def( + "getNumbers", getNumbers_function_value, bp::release_gil_policy(), "Get the atomic numbers for the atoms in the QM region.\nReturn:s\nA vector of atomic numbers for the atoms in the QM region.\n"); } { //::SireOpenMM::PyQMEngine::operator= - - typedef ::SireOpenMM::PyQMEngine & ( ::SireOpenMM::PyQMEngine::*assign_function_type)( ::SireOpenMM::PyQMEngine const & ) ; - assign_function_type assign_function_value( &::SireOpenMM::PyQMEngine::operator= ); - - PyQMEngine_exposer.def( - "assign" - , assign_function_value - , ( bp::arg("other") ) - , bp::return_self< >() - , "Assignment operator." ); - + + typedef ::SireOpenMM::PyQMEngine &(::SireOpenMM::PyQMEngine::*assign_function_type)(::SireOpenMM::PyQMEngine const &); + assign_function_type assign_function_value(&::SireOpenMM::PyQMEngine::operator=); + + PyQMEngine_exposer.def( + "assign", assign_function_value, (bp::arg("other")), bp::return_self<>(), "Assignment operator."); } { //::SireOpenMM::PyQMEngine::setAtoms - - typedef void ( ::SireOpenMM::PyQMEngine::*setAtoms_function_type)( ::QVector< int > ) ; - setAtoms_function_type setAtoms_function_value( &::SireOpenMM::PyQMEngine::setAtoms ); - - PyQMEngine_exposer.def( - "setAtoms" - , setAtoms_function_value - , ( bp::arg("atoms") ) - , bp::release_gil_policy() - , "Set the list of atom indices for the QM region.\nPar:am atoms\nA vector of atom indices for the QM region.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setAtoms_function_type)(::QVector); + setAtoms_function_type setAtoms_function_value(&::SireOpenMM::PyQMEngine::setAtoms); + + PyQMEngine_exposer.def( + "setAtoms", setAtoms_function_value, (bp::arg("atoms")), bp::release_gil_policy(), "Set the list of atom indices for the QM region.\nPar:am atoms\nA vector of atom indices for the QM region.\n"); } { //::SireOpenMM::PyQMEngine::setCallback - - typedef void ( ::SireOpenMM::PyQMEngine::*setCallback_function_type)( ::SireOpenMM::PyQMCallback ) ; - setCallback_function_type setCallback_function_value( &::SireOpenMM::PyQMEngine::setCallback ); - - PyQMEngine_exposer.def( - "setCallback" - , setCallback_function_value - , ( bp::arg("callback") ) - , bp::release_gil_policy() - , "Set the callback object.\nPar:am callback\nA Python object that contains the callback function.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setCallback_function_type)(::SireOpenMM::PyQMCallback); + setCallback_function_type setCallback_function_value(&::SireOpenMM::PyQMEngine::setCallback); + + PyQMEngine_exposer.def( + "setCallback", setCallback_function_value, (bp::arg("callback")), bp::release_gil_policy(), "Set the callback object.\nPar:am callback\nA Python object that contains the callback function.\n"); } { //::SireOpenMM::PyQMEngine::setCharges - - typedef void ( ::SireOpenMM::PyQMEngine::*setCharges_function_type)( ::QVector< double > ) ; - setCharges_function_type setCharges_function_value( &::SireOpenMM::PyQMEngine::setCharges ); - - PyQMEngine_exposer.def( - "setCharges" - , setCharges_function_value - , ( bp::arg("charges") ) - , bp::release_gil_policy() - , "Set the atomic charges of all atoms in the system.\nPar:am charges\nA vector of atomic charges for all atoms in the system.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setCharges_function_type)(::QVector); + setCharges_function_type setCharges_function_value(&::SireOpenMM::PyQMEngine::setCharges); + + PyQMEngine_exposer.def( + "setCharges", setCharges_function_value, (bp::arg("charges")), bp::release_gil_policy(), "Set the atomic charges of all atoms in the system.\nPar:am charges\nA vector of atomic charges for all atoms in the system.\n"); } { //::SireOpenMM::PyQMEngine::setCutoff - - typedef void ( ::SireOpenMM::PyQMEngine::*setCutoff_function_type)( ::SireUnits::Dimension::Length ) ; - setCutoff_function_type setCutoff_function_value( &::SireOpenMM::PyQMEngine::setCutoff ); - - PyQMEngine_exposer.def( - "setCutoff" - , setCutoff_function_value - , ( bp::arg("cutoff") ) - , bp::release_gil_policy() - , "Set the QM cutoff distance.\nPar:am cutoff\nThe QM cutoff distance.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setCutoff_function_type)(::SireUnits::Dimension::Length); + setCutoff_function_type setCutoff_function_value(&::SireOpenMM::PyQMEngine::setCutoff); + + PyQMEngine_exposer.def( + "setCutoff", setCutoff_function_value, (bp::arg("cutoff")), bp::release_gil_policy(), "Set the QM cutoff distance.\nPar:am cutoff\nThe QM cutoff distance.\n"); } { //::SireOpenMM::PyQMEngine::setIsMechanical - - typedef void ( ::SireOpenMM::PyQMEngine::*setIsMechanical_function_type)( bool ) ; - setIsMechanical_function_type setIsMechanical_function_value( &::SireOpenMM::PyQMEngine::setIsMechanical ); - - PyQMEngine_exposer.def( - "setIsMechanical" - , setIsMechanical_function_value - , ( bp::arg("is_mechanical") ) - , bp::release_gil_policy() - , "Set the mechanical embedding flag.\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setIsMechanical_function_type)(bool); + setIsMechanical_function_type setIsMechanical_function_value(&::SireOpenMM::PyQMEngine::setIsMechanical); + + PyQMEngine_exposer.def( + "setIsMechanical", setIsMechanical_function_value, (bp::arg("is_mechanical")), bp::release_gil_policy(), "Set the mechanical embedding flag.\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n"); } { //::SireOpenMM::PyQMEngine::setLambda - - typedef void ( ::SireOpenMM::PyQMEngine::*setLambda_function_type)( double ) ; - setLambda_function_type setLambda_function_value( &::SireOpenMM::PyQMEngine::setLambda ); - - PyQMEngine_exposer.def( - "setLambda" - , setLambda_function_value - , ( bp::arg("lambda") ) - , bp::release_gil_policy() - , "Set the lambda weighting factor.\nPar:am lambda\nThe lambda weighting factor.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setLambda_function_type)(double); + setLambda_function_type setLambda_function_value(&::SireOpenMM::PyQMEngine::setLambda); + + PyQMEngine_exposer.def( + "setLambda", setLambda_function_value, (bp::arg("lambda")), bp::release_gil_policy(), "Set the lambda weighting factor.\nPar:am lambda\nThe lambda weighting factor.\n"); } { //::SireOpenMM::PyQMEngine::setLinkAtoms - - typedef void ( ::SireOpenMM::PyQMEngine::*setLinkAtoms_function_type)( ::QMap< int, int >,::QMap< int, QVector< int > >,::QMap< int, double > ) ; - setLinkAtoms_function_type setLinkAtoms_function_value( &::SireOpenMM::PyQMEngine::setLinkAtoms ); - - PyQMEngine_exposer.def( - "setLinkAtoms" - , setLinkAtoms_function_value - , ( bp::arg("mm1_to_qm"), bp::arg("mm1_to_mm2"), bp::arg("bond_scale_factors") ) - , bp::release_gil_policy() - , "Set the link atoms associated with each QM atom.\nPar:am mm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nPar:am mm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nPar:am bond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setLinkAtoms_function_type)(::QMap, ::QMap>, ::QMap); + setLinkAtoms_function_type setLinkAtoms_function_value(&::SireOpenMM::PyQMEngine::setLinkAtoms); + + PyQMEngine_exposer.def( + "setLinkAtoms", setLinkAtoms_function_value, (bp::arg("mm1_to_qm"), bp::arg("mm1_to_mm2"), bp::arg("bond_scale_factors")), bp::release_gil_policy(), "Set the link atoms associated with each QM atom.\nPar:am mm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nPar:am mm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nPar:am bond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n"); } { //::SireOpenMM::PyQMEngine::setNeighbourListFrequency - - typedef void ( ::SireOpenMM::PyQMEngine::*setNeighbourListFrequency_function_type)( int ) ; - setNeighbourListFrequency_function_type setNeighbourListFrequency_function_value( &::SireOpenMM::PyQMEngine::setNeighbourListFrequency ); - - PyQMEngine_exposer.def( - "setNeighbourListFrequency" - , setNeighbourListFrequency_function_value - , ( bp::arg("neighbour_list_frequency") ) - , bp::release_gil_policy() - , "Set the neighbour list frequency.\nPar:am neighbour_list_frequency\nThe neighbour list frequency.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setNeighbourListFrequency_function_type)(int); + setNeighbourListFrequency_function_type setNeighbourListFrequency_function_value(&::SireOpenMM::PyQMEngine::setNeighbourListFrequency); + + PyQMEngine_exposer.def( + "setNeighbourListFrequency", setNeighbourListFrequency_function_value, (bp::arg("neighbour_list_frequency")), bp::release_gil_policy(), "Set the neighbour list frequency.\nPar:am neighbour_list_frequency\nThe neighbour list frequency.\n"); } { //::SireOpenMM::PyQMEngine::setNumbers - - typedef void ( ::SireOpenMM::PyQMEngine::*setNumbers_function_type)( ::QVector< int > ) ; - setNumbers_function_type setNumbers_function_value( &::SireOpenMM::PyQMEngine::setNumbers ); - - PyQMEngine_exposer.def( - "setNumbers" - , setNumbers_function_value - , ( bp::arg("numbers") ) - , bp::release_gil_policy() - , "Set the atomic numbers for the atoms in the QM region.\nPar:am numbers\nA vector of atomic numbers for the atoms in the QM region.\n" ); - + + typedef void (::SireOpenMM::PyQMEngine::*setNumbers_function_type)(::QVector); + setNumbers_function_type setNumbers_function_value(&::SireOpenMM::PyQMEngine::setNumbers); + + PyQMEngine_exposer.def( + "setNumbers", setNumbers_function_value, (bp::arg("numbers")), bp::release_gil_policy(), "Set the atomic numbers for the atoms in the QM region.\nPar:am numbers\nA vector of atomic numbers for the atoms in the QM region.\n"); + } + { //::SireOpenMM::PyQMEngine::getSwitchWidth + + typedef double (::SireOpenMM::PyQMEngine::*getSwitchWidth_function_type)() const; + getSwitchWidth_function_type getSwitchWidth_function_value(&::SireOpenMM::PyQMEngine::getSwitchWidth); + + PyQMEngine_exposer.def( + "getSwitchWidth", getSwitchWidth_function_value, bp::release_gil_policy(), "Get the switch width as a fraction of the cutoff.\n"); + } + { //::SireOpenMM::PyQMEngine::setSwitchWidth + + typedef void (::SireOpenMM::PyQMEngine::*setSwitchWidth_function_type)(double); + setSwitchWidth_function_type setSwitchWidth_function_value(&::SireOpenMM::PyQMEngine::setSwitchWidth); + + PyQMEngine_exposer.def( + "setSwitchWidth", setSwitchWidth_function_value, (bp::arg("switch_width")), bp::release_gil_policy(), "Set the switch width as a fraction of the cutoff (0 to 1).\n"); + } + { //::SireOpenMM::PyQMEngine::getUseSwitch + + typedef bool (::SireOpenMM::PyQMEngine::*getUseSwitch_function_type)() const; + getUseSwitch_function_type getUseSwitch_function_value(&::SireOpenMM::PyQMEngine::getUseSwitch); + + PyQMEngine_exposer.def( + "getUseSwitch", getUseSwitch_function_value, bp::release_gil_policy(), "Get whether a switching function is used.\n"); + } + { //::SireOpenMM::PyQMEngine::setUseSwitch + + typedef void (::SireOpenMM::PyQMEngine::*setUseSwitch_function_type)(bool); + setUseSwitch_function_type setUseSwitch_function_value(&::SireOpenMM::PyQMEngine::setUseSwitch); + + PyQMEngine_exposer.def( + "setUseSwitch", setUseSwitch_function_value, (bp::arg("use_switch")), bp::release_gil_policy(), "Set whether to use a switching function.\n"); } { //::SireOpenMM::PyQMEngine::typeName - - typedef char const * ( *typeName_function_type )( ); - typeName_function_type typeName_function_value( &::SireOpenMM::PyQMEngine::typeName ); - - PyQMEngine_exposer.def( - "typeName" - , typeName_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(*typeName_function_type)(); + typeName_function_type typeName_function_value(&::SireOpenMM::PyQMEngine::typeName); + + PyQMEngine_exposer.def( + "typeName", typeName_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } { //::SireOpenMM::PyQMEngine::what - - typedef char const * ( ::SireOpenMM::PyQMEngine::*what_function_type)( ) const; - what_function_type what_function_value( &::SireOpenMM::PyQMEngine::what ); - - PyQMEngine_exposer.def( - "what" - , what_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(::SireOpenMM::PyQMEngine::*what_function_type)() const; + what_function_type what_function_value(&::SireOpenMM::PyQMEngine::what); + + PyQMEngine_exposer.def( + "what", what_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } - PyQMEngine_exposer.staticmethod( "typeName" ); - PyQMEngine_exposer.def( "__copy__", &__copy__); - PyQMEngine_exposer.def( "__deepcopy__", &__copy__); - PyQMEngine_exposer.def( "clone", &__copy__); - PyQMEngine_exposer.def( "__str__", &__str__< ::SireOpenMM::PyQMEngine > ); - PyQMEngine_exposer.def( "__repr__", &__str__< ::SireOpenMM::PyQMEngine > ); + PyQMEngine_exposer.staticmethod("typeName"); + PyQMEngine_exposer.def("__copy__", &__copy__); + PyQMEngine_exposer.def("__deepcopy__", &__copy__); + PyQMEngine_exposer.def("clone", &__copy__); + PyQMEngine_exposer.def("__str__", &__str__<::SireOpenMM::PyQMEngine>); + PyQMEngine_exposer.def("__repr__", &__str__<::SireOpenMM::PyQMEngine>); } - } diff --git a/wrapper/Convert/SireOpenMM/TorchQMEngine.pypp.cpp b/wrapper/Convert/SireOpenMM/TorchQMEngine.pypp.cpp index 29847e98d..b40473582 100644 --- a/wrapper/Convert/SireOpenMM/TorchQMEngine.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/TorchQMEngine.pypp.cpp @@ -2,8 +2,8 @@ // (C) Christopher Woods, GPL >= 3 License -#include "boost/python.hpp" #include "TorchQMEngine.pypp.hpp" +#include "boost/python.hpp" namespace bp = boost::python; @@ -39,300 +39,234 @@ namespace bp = boost::python; #include "torchqm.h" -SireOpenMM::TorchQMEngine __copy__(const SireOpenMM::TorchQMEngine &other){ return SireOpenMM::TorchQMEngine(other); } +SireOpenMM::TorchQMEngine __copy__(const SireOpenMM::TorchQMEngine &other) { return SireOpenMM::TorchQMEngine(other); } #include "Helpers/str.hpp" #include "Helpers/release_gil_policy.hpp" -void register_TorchQMEngine_class(){ +void register_TorchQMEngine_class() +{ { //::SireOpenMM::TorchQMEngine - typedef bp::class_< SireOpenMM::TorchQMEngine, bp::bases< SireBase::Property, SireOpenMM::QMEngine > > TorchQMEngine_exposer_t; - TorchQMEngine_exposer_t TorchQMEngine_exposer = TorchQMEngine_exposer_t( "TorchQMEngine", "", bp::init< >("Default constructor.") ); - bp::scope TorchQMEngine_scope( TorchQMEngine_exposer ); - TorchQMEngine_exposer.def( bp::init< QString, bp::optional< SireUnits::Dimension::Length, int, bool, double > >(( bp::arg("arg0"), bp::arg("cutoff")=7.5 * SireUnits::angstrom, bp::arg("neighbour_list_frequency")=(int)(0), bp::arg("is_mechanical")=(bool)(false), bp::arg("lambda")=1. ), "Constructor\nPar:am module_path\nThe path to the serialised TorchScript module.\n\nPar:am cutoff\nThe ML cutoff distance.\n\nPar:am neighbour_list_frequency\nThe frequency at which the neighbour list is updated. (Number of steps.)\nIf zero, then no neighbour list is used.\n\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n\nPar:am lambda\nThe lambda weighting factor. This can be used to interpolate between\npotentials for end-state correction calculations.\n") ); - TorchQMEngine_exposer.def( bp::init< SireOpenMM::TorchQMEngine const & >(( bp::arg("other") ), "Copy constructor.") ); + typedef bp::class_> TorchQMEngine_exposer_t; + TorchQMEngine_exposer_t TorchQMEngine_exposer = TorchQMEngine_exposer_t("TorchQMEngine", "", bp::init<>("Default constructor.")); + bp::scope TorchQMEngine_scope(TorchQMEngine_exposer); + TorchQMEngine_exposer.def(bp::init>((bp::arg("arg0"), bp::arg("cutoff") = 7.5 * SireUnits::angstrom, bp::arg("neighbour_list_frequency") = (int)(0), bp::arg("is_mechanical") = (bool)(false), bp::arg("lambda") = 1.), "Constructor\nPar:am module_path\nThe path to the serialised TorchScript module.\n\nPar:am cutoff\nThe ML cutoff distance.\n\nPar:am neighbour_list_frequency\nThe frequency at which the neighbour list is updated. (Number of steps.)\nIf zero, then no neighbour list is used.\n\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n\nPar:am lambda\nThe lambda weighting factor. This can be used to interpolate between\npotentials for end-state correction calculations.\n")); + TorchQMEngine_exposer.def(bp::init((bp::arg("other")), "Copy constructor.")); { //::SireOpenMM::TorchQMEngine::getModulePath - - typedef ::QString ( ::SireOpenMM::TorchQMEngine::*getModulePath_function_type)( ) const; - getModulePath_function_type getModulePath_function_value( &::SireOpenMM::TorchQMEngine::getModulePath ); - - TorchQMEngine_exposer.def( - "getModulePath" - , getModulePath_function_value - , bp::release_gil_policy() - , "Get the path to the serialised TorchScript module.\n" ); - + + typedef ::QString (::SireOpenMM::TorchQMEngine::*getModulePath_function_type)() const; + getModulePath_function_type getModulePath_function_value(&::SireOpenMM::TorchQMEngine::getModulePath); + + TorchQMEngine_exposer.def( + "getModulePath", getModulePath_function_value, bp::release_gil_policy(), "Get the path to the serialised TorchScript module.\n"); } { //::SireOpenMM::TorchQMEngine::getAtoms - - typedef ::QVector< int > ( ::SireOpenMM::TorchQMEngine::*getAtoms_function_type)( ) const; - getAtoms_function_type getAtoms_function_value( &::SireOpenMM::TorchQMEngine::getAtoms ); - - TorchQMEngine_exposer.def( - "getAtoms" - , getAtoms_function_value - , bp::release_gil_policy() - , "Get the indices of the atoms in the QM region.\nReturn:s\nA vector of atom indices for the QM region.\n" ); - + + typedef ::QVector (::SireOpenMM::TorchQMEngine::*getAtoms_function_type)() const; + getAtoms_function_type getAtoms_function_value(&::SireOpenMM::TorchQMEngine::getAtoms); + + TorchQMEngine_exposer.def( + "getAtoms", getAtoms_function_value, bp::release_gil_policy(), "Get the indices of the atoms in the QM region.\nReturn:s\nA vector of atom indices for the QM region.\n"); } { //::SireOpenMM::TorchQMEngine::getCharges - - typedef ::QVector< double > ( ::SireOpenMM::TorchQMEngine::*getCharges_function_type)( ) const; - getCharges_function_type getCharges_function_value( &::SireOpenMM::TorchQMEngine::getCharges ); - - TorchQMEngine_exposer.def( - "getCharges" - , getCharges_function_value - , bp::release_gil_policy() - , "Get the atomic charges of all atoms in the system.\nReturn:s\nA vector of atomic charges for all atoms in the system.\n" ); - + + typedef ::QVector (::SireOpenMM::TorchQMEngine::*getCharges_function_type)() const; + getCharges_function_type getCharges_function_value(&::SireOpenMM::TorchQMEngine::getCharges); + + TorchQMEngine_exposer.def( + "getCharges", getCharges_function_value, bp::release_gil_policy(), "Get the atomic charges of all atoms in the system.\nReturn:s\nA vector of atomic charges for all atoms in the system.\n"); } { //::SireOpenMM::TorchQMEngine::getCutoff - - typedef ::SireUnits::Dimension::Length ( ::SireOpenMM::TorchQMEngine::*getCutoff_function_type)( ) const; - getCutoff_function_type getCutoff_function_value( &::SireOpenMM::TorchQMEngine::getCutoff ); - - TorchQMEngine_exposer.def( - "getCutoff" - , getCutoff_function_value - , bp::release_gil_policy() - , "Get the QM cutoff distance.\nReturn:s\nThe QM cutoff distance.\n" ); - + + typedef ::SireUnits::Dimension::Length (::SireOpenMM::TorchQMEngine::*getCutoff_function_type)() const; + getCutoff_function_type getCutoff_function_value(&::SireOpenMM::TorchQMEngine::getCutoff); + + TorchQMEngine_exposer.def( + "getCutoff", getCutoff_function_value, bp::release_gil_policy(), "Get the QM cutoff distance.\nReturn:s\nThe QM cutoff distance.\n"); } { //::SireOpenMM::TorchQMEngine::getIsMechanical - - typedef bool ( ::SireOpenMM::TorchQMEngine::*getIsMechanical_function_type)( ) const; - getIsMechanical_function_type getIsMechanical_function_value( &::SireOpenMM::TorchQMEngine::getIsMechanical ); - - TorchQMEngine_exposer.def( - "getIsMechanical" - , getIsMechanical_function_value - , bp::release_gil_policy() - , "Get the mechanical embedding flag.\nReturn:s\nA flag to indicate if mechanical embedding is being used.\n" ); - + + typedef bool (::SireOpenMM::TorchQMEngine::*getIsMechanical_function_type)() const; + getIsMechanical_function_type getIsMechanical_function_value(&::SireOpenMM::TorchQMEngine::getIsMechanical); + + TorchQMEngine_exposer.def( + "getIsMechanical", getIsMechanical_function_value, bp::release_gil_policy(), "Get the mechanical embedding flag.\nReturn:s\nA flag to indicate if mechanical embedding is being used.\n"); } { //::SireOpenMM::TorchQMEngine::getLambda - - typedef double ( ::SireOpenMM::TorchQMEngine::*getLambda_function_type)( ) const; - getLambda_function_type getLambda_function_value( &::SireOpenMM::TorchQMEngine::getLambda ); - - TorchQMEngine_exposer.def( - "getLambda" - , getLambda_function_value - , bp::release_gil_policy() - , "Get the lambda weighting factor.\nReturn:s\nThe lambda weighting factor.\n" ); - + + typedef double (::SireOpenMM::TorchQMEngine::*getLambda_function_type)() const; + getLambda_function_type getLambda_function_value(&::SireOpenMM::TorchQMEngine::getLambda); + + TorchQMEngine_exposer.def( + "getLambda", getLambda_function_value, bp::release_gil_policy(), "Get the lambda weighting factor.\nReturn:s\nThe lambda weighting factor.\n"); } { //::SireOpenMM::TorchQMEngine::getLinkAtoms - - typedef ::boost::tuples::tuple< QMap< int, int >, QMap< int, QVector< int > >, QMap< int, double >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::TorchQMEngine::*getLinkAtoms_function_type)( ) const; - getLinkAtoms_function_type getLinkAtoms_function_value( &::SireOpenMM::TorchQMEngine::getLinkAtoms ); - - TorchQMEngine_exposer.def( - "getLinkAtoms" - , getLinkAtoms_function_value - , bp::release_gil_policy() - , "Get the link atoms associated with each QM atom.\nReturn:s\nA tuple containing:\n\nmm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nmm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nbond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n" ); - + + typedef ::boost::tuples::tuple, QMap>, QMap, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::TorchQMEngine::*getLinkAtoms_function_type)() const; + getLinkAtoms_function_type getLinkAtoms_function_value(&::SireOpenMM::TorchQMEngine::getLinkAtoms); + + TorchQMEngine_exposer.def( + "getLinkAtoms", getLinkAtoms_function_value, bp::release_gil_policy(), "Get the link atoms associated with each QM atom.\nReturn:s\nA tuple containing:\n\nmm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nmm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nbond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n"); } { //::SireOpenMM::TorchQMEngine::getMM2Atoms - - typedef ::QVector< int > ( ::SireOpenMM::TorchQMEngine::*getMM2Atoms_function_type)( ) const; - getMM2Atoms_function_type getMM2Atoms_function_value( &::SireOpenMM::TorchQMEngine::getMM2Atoms ); - - TorchQMEngine_exposer.def( - "getMM2Atoms" - , getMM2Atoms_function_value - , bp::release_gil_policy() - , "Get the vector of MM2 atoms.\nReturn:s\nA vector of MM2 atom indices.\n" ); - + + typedef ::QVector (::SireOpenMM::TorchQMEngine::*getMM2Atoms_function_type)() const; + getMM2Atoms_function_type getMM2Atoms_function_value(&::SireOpenMM::TorchQMEngine::getMM2Atoms); + + TorchQMEngine_exposer.def( + "getMM2Atoms", getMM2Atoms_function_value, bp::release_gil_policy(), "Get the vector of MM2 atoms.\nReturn:s\nA vector of MM2 atom indices.\n"); } { //::SireOpenMM::TorchQMEngine::getNeighbourListFrequency - - typedef int ( ::SireOpenMM::TorchQMEngine::*getNeighbourListFrequency_function_type)( ) const; - getNeighbourListFrequency_function_type getNeighbourListFrequency_function_value( &::SireOpenMM::TorchQMEngine::getNeighbourListFrequency ); - - TorchQMEngine_exposer.def( - "getNeighbourListFrequency" - , getNeighbourListFrequency_function_value - , bp::release_gil_policy() - , "Get the neighbour list frequency.\nReturn:s\nThe neighbour list frequency.\n" ); - + + typedef int (::SireOpenMM::TorchQMEngine::*getNeighbourListFrequency_function_type)() const; + getNeighbourListFrequency_function_type getNeighbourListFrequency_function_value(&::SireOpenMM::TorchQMEngine::getNeighbourListFrequency); + + TorchQMEngine_exposer.def( + "getNeighbourListFrequency", getNeighbourListFrequency_function_value, bp::release_gil_policy(), "Get the neighbour list frequency.\nReturn:s\nThe neighbour list frequency.\n"); } { //::SireOpenMM::TorchQMEngine::getNumbers - - typedef ::QVector< int > ( ::SireOpenMM::TorchQMEngine::*getNumbers_function_type)( ) const; - getNumbers_function_type getNumbers_function_value( &::SireOpenMM::TorchQMEngine::getNumbers ); - - TorchQMEngine_exposer.def( - "getNumbers" - , getNumbers_function_value - , bp::release_gil_policy() - , "Get the atomic numbers for the atoms in the QM region.\nReturn:s\nA vector of atomic numbers for the atoms in the QM region.\n" ); - + + typedef ::QVector (::SireOpenMM::TorchQMEngine::*getNumbers_function_type)() const; + getNumbers_function_type getNumbers_function_value(&::SireOpenMM::TorchQMEngine::getNumbers); + + TorchQMEngine_exposer.def( + "getNumbers", getNumbers_function_value, bp::release_gil_policy(), "Get the atomic numbers for the atoms in the QM region.\nReturn:s\nA vector of atomic numbers for the atoms in the QM region.\n"); } { //::SireOpenMM::TorchQMEngine::operator= - - typedef ::SireOpenMM::TorchQMEngine & ( ::SireOpenMM::TorchQMEngine::*assign_function_type)( ::SireOpenMM::TorchQMEngine const & ) ; - assign_function_type assign_function_value( &::SireOpenMM::TorchQMEngine::operator= ); - - TorchQMEngine_exposer.def( - "assign" - , assign_function_value - , ( bp::arg("other") ) - , bp::return_self< >() - , "Assignment operator." ); - + + typedef ::SireOpenMM::TorchQMEngine &(::SireOpenMM::TorchQMEngine::*assign_function_type)(::SireOpenMM::TorchQMEngine const &); + assign_function_type assign_function_value(&::SireOpenMM::TorchQMEngine::operator=); + + TorchQMEngine_exposer.def( + "assign", assign_function_value, (bp::arg("other")), bp::return_self<>(), "Assignment operator."); } { //::SireOpenMM::TorchQMEngine::setModulePath - - typedef void ( ::SireOpenMM::TorchQMEngine::*setModulePath_function_type)( ::QString ) ; - setModulePath_function_type setModulePath_function_value( &::SireOpenMM::TorchQMEngine::setModulePath ); - - TorchQMEngine_exposer.def( - "setModulePath" - , setModulePath_function_value - , ( bp::arg("module_path") ) - , bp::release_gil_policy() - , "Set the path to the serialised TorchScript module.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setModulePath_function_type)(::QString); + setModulePath_function_type setModulePath_function_value(&::SireOpenMM::TorchQMEngine::setModulePath); + + TorchQMEngine_exposer.def( + "setModulePath", setModulePath_function_value, (bp::arg("module_path")), bp::release_gil_policy(), "Set the path to the serialised TorchScript module.\n"); } { //::SireOpenMM::TorchQMEngine::setAtoms - - typedef void ( ::SireOpenMM::TorchQMEngine::*setAtoms_function_type)( ::QVector< int > ) ; - setAtoms_function_type setAtoms_function_value( &::SireOpenMM::TorchQMEngine::setAtoms ); - - TorchQMEngine_exposer.def( - "setAtoms" - , setAtoms_function_value - , ( bp::arg("atoms") ) - , bp::release_gil_policy() - , "Set the list of atom indices for the QM region.\nPar:am atoms\nA vector of atom indices for the QM region.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setAtoms_function_type)(::QVector); + setAtoms_function_type setAtoms_function_value(&::SireOpenMM::TorchQMEngine::setAtoms); + + TorchQMEngine_exposer.def( + "setAtoms", setAtoms_function_value, (bp::arg("atoms")), bp::release_gil_policy(), "Set the list of atom indices for the QM region.\nPar:am atoms\nA vector of atom indices for the QM region.\n"); } { //::SireOpenMM::TorchQMEngine::setCharges - - typedef void ( ::SireOpenMM::TorchQMEngine::*setCharges_function_type)( ::QVector< double > ) ; - setCharges_function_type setCharges_function_value( &::SireOpenMM::TorchQMEngine::setCharges ); - - TorchQMEngine_exposer.def( - "setCharges" - , setCharges_function_value - , ( bp::arg("charges") ) - , bp::release_gil_policy() - , "Set the atomic charges of all atoms in the system.\nPar:am charges\nA vector of atomic charges for all atoms in the system.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setCharges_function_type)(::QVector); + setCharges_function_type setCharges_function_value(&::SireOpenMM::TorchQMEngine::setCharges); + + TorchQMEngine_exposer.def( + "setCharges", setCharges_function_value, (bp::arg("charges")), bp::release_gil_policy(), "Set the atomic charges of all atoms in the system.\nPar:am charges\nA vector of atomic charges for all atoms in the system.\n"); } { //::SireOpenMM::TorchQMEngine::setCutoff - - typedef void ( ::SireOpenMM::TorchQMEngine::*setCutoff_function_type)( ::SireUnits::Dimension::Length ) ; - setCutoff_function_type setCutoff_function_value( &::SireOpenMM::TorchQMEngine::setCutoff ); - - TorchQMEngine_exposer.def( - "setCutoff" - , setCutoff_function_value - , ( bp::arg("cutoff") ) - , bp::release_gil_policy() - , "Set the QM cutoff distance.\nPar:am cutoff\nThe QM cutoff distance.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setCutoff_function_type)(::SireUnits::Dimension::Length); + setCutoff_function_type setCutoff_function_value(&::SireOpenMM::TorchQMEngine::setCutoff); + + TorchQMEngine_exposer.def( + "setCutoff", setCutoff_function_value, (bp::arg("cutoff")), bp::release_gil_policy(), "Set the QM cutoff distance.\nPar:am cutoff\nThe QM cutoff distance.\n"); } { //::SireOpenMM::TorchQMEngine::setIsMechanical - - typedef void ( ::SireOpenMM::TorchQMEngine::*setIsMechanical_function_type)( bool ) ; - setIsMechanical_function_type setIsMechanical_function_value( &::SireOpenMM::TorchQMEngine::setIsMechanical ); - - TorchQMEngine_exposer.def( - "setIsMechanical" - , setIsMechanical_function_value - , ( bp::arg("is_mechanical") ) - , bp::release_gil_policy() - , "Set the mechanical embedding flag.\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setIsMechanical_function_type)(bool); + setIsMechanical_function_type setIsMechanical_function_value(&::SireOpenMM::TorchQMEngine::setIsMechanical); + + TorchQMEngine_exposer.def( + "setIsMechanical", setIsMechanical_function_value, (bp::arg("is_mechanical")), bp::release_gil_policy(), "Set the mechanical embedding flag.\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n"); } { //::SireOpenMM::TorchQMEngine::setLambda - - typedef void ( ::SireOpenMM::TorchQMEngine::*setLambda_function_type)( double ) ; - setLambda_function_type setLambda_function_value( &::SireOpenMM::TorchQMEngine::setLambda ); - - TorchQMEngine_exposer.def( - "setLambda" - , setLambda_function_value - , ( bp::arg("lambda") ) - , bp::release_gil_policy() - , "Set the lambda weighting factor.\nPar:am lambda\nThe lambda weighting factor.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setLambda_function_type)(double); + setLambda_function_type setLambda_function_value(&::SireOpenMM::TorchQMEngine::setLambda); + + TorchQMEngine_exposer.def( + "setLambda", setLambda_function_value, (bp::arg("lambda")), bp::release_gil_policy(), "Set the lambda weighting factor.\nPar:am lambda\nThe lambda weighting factor.\n"); } { //::SireOpenMM::TorchQMEngine::setLinkAtoms - - typedef void ( ::SireOpenMM::TorchQMEngine::*setLinkAtoms_function_type)( ::QMap< int, int >,::QMap< int, QVector< int > >,::QMap< int, double > ) ; - setLinkAtoms_function_type setLinkAtoms_function_value( &::SireOpenMM::TorchQMEngine::setLinkAtoms ); - - TorchQMEngine_exposer.def( - "setLinkAtoms" - , setLinkAtoms_function_value - , ( bp::arg("mm1_to_qm"), bp::arg("mm1_to_mm2"), bp::arg("bond_scale_factors") ) - , bp::release_gil_policy() - , "Set the link atoms associated with each QM atom.\nPar:am mm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nPar:am mm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nPar:am bond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setLinkAtoms_function_type)(::QMap, ::QMap>, ::QMap); + setLinkAtoms_function_type setLinkAtoms_function_value(&::SireOpenMM::TorchQMEngine::setLinkAtoms); + + TorchQMEngine_exposer.def( + "setLinkAtoms", setLinkAtoms_function_value, (bp::arg("mm1_to_qm"), bp::arg("mm1_to_mm2"), bp::arg("bond_scale_factors")), bp::release_gil_policy(), "Set the link atoms associated with each QM atom.\nPar:am mm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nPar:am mm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nPar:am bond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n"); } { //::SireOpenMM::TorchQMEngine::setNeighbourListFrequency - - typedef void ( ::SireOpenMM::TorchQMEngine::*setNeighbourListFrequency_function_type)( int ) ; - setNeighbourListFrequency_function_type setNeighbourListFrequency_function_value( &::SireOpenMM::TorchQMEngine::setNeighbourListFrequency ); - - TorchQMEngine_exposer.def( - "setNeighbourListFrequency" - , setNeighbourListFrequency_function_value - , ( bp::arg("neighbour_list_frequency") ) - , bp::release_gil_policy() - , "Set the neighbour list frequency.\nPar:am neighbour_list_frequency\nThe neighbour list frequency.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setNeighbourListFrequency_function_type)(int); + setNeighbourListFrequency_function_type setNeighbourListFrequency_function_value(&::SireOpenMM::TorchQMEngine::setNeighbourListFrequency); + + TorchQMEngine_exposer.def( + "setNeighbourListFrequency", setNeighbourListFrequency_function_value, (bp::arg("neighbour_list_frequency")), bp::release_gil_policy(), "Set the neighbour list frequency.\nPar:am neighbour_list_frequency\nThe neighbour list frequency.\n"); } { //::SireOpenMM::TorchQMEngine::setNumbers - - typedef void ( ::SireOpenMM::TorchQMEngine::*setNumbers_function_type)( ::QVector< int > ) ; - setNumbers_function_type setNumbers_function_value( &::SireOpenMM::TorchQMEngine::setNumbers ); - - TorchQMEngine_exposer.def( - "setNumbers" - , setNumbers_function_value - , ( bp::arg("numbers") ) - , bp::release_gil_policy() - , "Set the atomic numbers for the atoms in the QM region.\nPar:am numbers\nA vector of atomic numbers for the atoms in the QM region.\n" ); - + + typedef void (::SireOpenMM::TorchQMEngine::*setNumbers_function_type)(::QVector); + setNumbers_function_type setNumbers_function_value(&::SireOpenMM::TorchQMEngine::setNumbers); + + TorchQMEngine_exposer.def( + "setNumbers", setNumbers_function_value, (bp::arg("numbers")), bp::release_gil_policy(), "Set the atomic numbers for the atoms in the QM region.\nPar:am numbers\nA vector of atomic numbers for the atoms in the QM region.\n"); + } + { //::SireOpenMM::TorchQMEngine::getSwitchWidth + + typedef double (::SireOpenMM::TorchQMEngine::*getSwitchWidth_function_type)() const; + getSwitchWidth_function_type getSwitchWidth_function_value(&::SireOpenMM::TorchQMEngine::getSwitchWidth); + + TorchQMEngine_exposer.def( + "getSwitchWidth", getSwitchWidth_function_value, bp::release_gil_policy(), "Get the switch width as a fraction of the cutoff.\n"); + } + { //::SireOpenMM::TorchQMEngine::setSwitchWidth + + typedef void (::SireOpenMM::TorchQMEngine::*setSwitchWidth_function_type)(double); + setSwitchWidth_function_type setSwitchWidth_function_value(&::SireOpenMM::TorchQMEngine::setSwitchWidth); + + TorchQMEngine_exposer.def( + "setSwitchWidth", setSwitchWidth_function_value, (bp::arg("switch_width")), bp::release_gil_policy(), "Set the switch width as a fraction of the cutoff (0 to 1).\n"); + } + { //::SireOpenMM::TorchQMEngine::getUseSwitch + + typedef bool (::SireOpenMM::TorchQMEngine::*getUseSwitch_function_type)() const; + getUseSwitch_function_type getUseSwitch_function_value(&::SireOpenMM::TorchQMEngine::getUseSwitch); + + TorchQMEngine_exposer.def( + "getUseSwitch", getUseSwitch_function_value, bp::release_gil_policy(), "Get whether a switching function is used.\n"); + } + { //::SireOpenMM::TorchQMEngine::setUseSwitch + + typedef void (::SireOpenMM::TorchQMEngine::*setUseSwitch_function_type)(bool); + setUseSwitch_function_type setUseSwitch_function_value(&::SireOpenMM::TorchQMEngine::setUseSwitch); + + TorchQMEngine_exposer.def( + "setUseSwitch", setUseSwitch_function_value, (bp::arg("use_switch")), bp::release_gil_policy(), "Set whether to use a switching function.\n"); } { //::SireOpenMM::TorchQMEngine::typeName - - typedef char const * ( *typeName_function_type )( ); - typeName_function_type typeName_function_value( &::SireOpenMM::TorchQMEngine::typeName ); - - TorchQMEngine_exposer.def( - "typeName" - , typeName_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(*typeName_function_type)(); + typeName_function_type typeName_function_value(&::SireOpenMM::TorchQMEngine::typeName); + + TorchQMEngine_exposer.def( + "typeName", typeName_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } { //::SireOpenMM::TorchQMEngine::what - - typedef char const * ( ::SireOpenMM::TorchQMEngine::*what_function_type)( ) const; - what_function_type what_function_value( &::SireOpenMM::TorchQMEngine::what ); - - TorchQMEngine_exposer.def( - "what" - , what_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(::SireOpenMM::TorchQMEngine::*what_function_type)() const; + what_function_type what_function_value(&::SireOpenMM::TorchQMEngine::what); + + TorchQMEngine_exposer.def( + "what", what_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } - TorchQMEngine_exposer.staticmethod( "typeName" ); - TorchQMEngine_exposer.def( "__copy__", &__copy__); - TorchQMEngine_exposer.def( "__deepcopy__", &__copy__); - TorchQMEngine_exposer.def( "clone", &__copy__); - TorchQMEngine_exposer.def( "__str__", &__str__< ::SireOpenMM::TorchQMEngine > ); - TorchQMEngine_exposer.def( "__repr__", &__str__< ::SireOpenMM::TorchQMEngine > ); + TorchQMEngine_exposer.staticmethod("typeName"); + TorchQMEngine_exposer.def("__copy__", &__copy__); + TorchQMEngine_exposer.def("__deepcopy__", &__copy__); + TorchQMEngine_exposer.def("clone", &__copy__); + TorchQMEngine_exposer.def("__str__", &__str__<::SireOpenMM::TorchQMEngine>); + TorchQMEngine_exposer.def("__repr__", &__str__<::SireOpenMM::TorchQMEngine>); } - } diff --git a/wrapper/Convert/SireOpenMM/pyqm.cpp b/wrapper/Convert/SireOpenMM/pyqm.cpp index c8e0f1334..de703b1bf 100644 --- a/wrapper/Convert/SireOpenMM/pyqm.cpp +++ b/wrapper/Convert/SireOpenMM/pyqm.cpp @@ -26,6 +26,7 @@ * \*********************************************/ +#include #include #include @@ -54,8 +55,9 @@ static const double VIRTUAL_PC_DELTA = 0.01; class GILLock { public: - GILLock() { state_ = PyGILState_Ensure(); } - ~GILLock() { PyGILState_Release(state_); } + GILLock() { state_ = PyGILState_Ensure(); } + ~GILLock() { PyGILState_Release(state_); } + private: PyGILState_STATE state_; }; @@ -85,8 +87,9 @@ bp::object getPythonObject(QString uuid) if (not py_object_registry.contains(uuid)) { throw SireError::invalid_key(QObject::tr( - "Unable to find UUID %1 in the PyQMForce callback registry.").arg(uuid), - CODELOC); + "Unable to find UUID %1 in the PyQMForce callback registry.") + .arg(uuid), + CODELOC); } return py_object_registry[uuid]; @@ -137,8 +140,7 @@ PyQMCallback::PyQMCallback() { } -PyQMCallback::PyQMCallback(bp::object py_object, QString name) : - py_object(py_object), name(name) +PyQMCallback::PyQMCallback(bp::object py_object, QString name) : py_object(py_object), name(name) { // Is this a method or free function. if (name.isEmpty()) @@ -172,15 +174,14 @@ PyQMCallback::call( xyz_qm, xyz_mm, cell, - idx_mm - ); + idx_mm); } catch (const bp::error_already_set &) { PyErr_Print(); throw SireError::process_error(QObject::tr( - "An error occurred when calling the QM Python callback method"), - CODELOC); + "An error occurred when calling the QM Python callback method"), + CODELOC); } } else @@ -194,19 +195,76 @@ PyQMCallback::call( xyz_qm, xyz_mm, cell, - idx_mm - ); + idx_mm); } catch (const bp::error_already_set &) { PyErr_Print(); throw SireError::process_error(QObject::tr( - "An error occurred when calling the QM Python callback method"), - CODELOC); + "An error occurred when calling the QM Python callback method"), + CODELOC); } } } +boost::tuple>, QVector>, QVector> +PyQMCallback::call4( + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const +{ + // Acquire GIL before calling Python code. + GILLock lock; + + bp::object result; + try + { + if (this->is_method) + { + result = bp::call_method( + this->py_object.ptr(), + this->name.toStdString().c_str(), + numbers_qm, charges_mm, xyz_qm, xyz_mm, cell, idx_mm); + } + else + { + result = bp::call( + this->py_object.ptr(), + numbers_qm, charges_mm, xyz_qm, xyz_mm, cell, idx_mm); + } + } + catch (const bp::error_already_set &) + { + PyErr_Print(); + throw SireError::process_error(QObject::tr( + "An error occurred when calling the QM Python callback method"), + CODELOC); + } + + // Extract mandatory first three elements while GIL is held. + const auto energy = bp::extract(result[0])(); + const auto forces_qm = bp::extract>>(result[1])(); + const auto forces_mm = bp::extract>>(result[2])(); + + // Extract optional fourth element: dE/dcharges_mm (empty if not returned). + QVector dE_dcharges_mm; + if (bp::len(result) > 3) + { + try + { + dE_dcharges_mm = bp::extract>(result[3])(); + } + catch (...) + { + } + } + + return boost::make_tuple(energy, forces_qm, forces_mm, dE_dcharges_mm); +} + const char *PyQMCallback::typeName() { return QMetaType::typeName(qMetaTypeId()); @@ -225,14 +283,15 @@ static const RegisterMetaType r_pyqmforce(NO_ROOT); QDataStream &operator<<(QDataStream &ds, const PyQMForce &pyqmforce) { - writeHeader(ds, r_pyqmforce, 1); + writeHeader(ds, r_pyqmforce, 2); SharedDataStream sds(ds); sds << pyqmforce.callback << pyqmforce.cutoff << pyqmforce.neighbour_list_frequency - << pyqmforce.is_mechanical << pyqmforce.lambda << pyqmforce.atoms + << pyqmforce.is_mechanical << pyqmforce.lambda << pyqmforce.atoms << pyqmforce.mm1_to_qm << pyqmforce.mm1_to_mm2 << pyqmforce.bond_scale_factors - << pyqmforce.mm2_atoms << pyqmforce.numbers << pyqmforce.charges; + << pyqmforce.mm2_atoms << pyqmforce.numbers << pyqmforce.charges + << pyqmforce.switch_width << pyqmforce.use_switch; return ds; } @@ -241,17 +300,23 @@ QDataStream &operator>>(QDataStream &ds, PyQMForce &pyqmforce) { VersionID v = readHeader(ds, r_pyqmforce); - if (v == 1) + if (v == 2) + { + SharedDataStream sds(ds); + + sds >> pyqmforce.callback >> pyqmforce.cutoff >> pyqmforce.neighbour_list_frequency >> pyqmforce.is_mechanical >> pyqmforce.lambda >> pyqmforce.atoms >> pyqmforce.mm1_to_qm >> pyqmforce.mm1_to_mm2 >> pyqmforce.bond_scale_factors >> pyqmforce.mm2_atoms >> pyqmforce.numbers >> pyqmforce.charges >> pyqmforce.switch_width >> pyqmforce.use_switch; + } + else if (v == 1) { SharedDataStream sds(ds); - sds >> pyqmforce.callback >> pyqmforce.cutoff >> pyqmforce.neighbour_list_frequency - >> pyqmforce.is_mechanical >> pyqmforce.lambda >> pyqmforce.atoms - >> pyqmforce.mm1_to_qm >> pyqmforce.mm1_to_mm2 >> pyqmforce.bond_scale_factors - >> pyqmforce.mm2_atoms >> pyqmforce.numbers >> pyqmforce.charges; + sds >> pyqmforce.callback >> pyqmforce.cutoff >> pyqmforce.neighbour_list_frequency >> pyqmforce.is_mechanical >> pyqmforce.lambda >> pyqmforce.atoms >> pyqmforce.mm1_to_qm >> pyqmforce.mm1_to_mm2 >> pyqmforce.bond_scale_factors >> pyqmforce.mm2_atoms >> pyqmforce.numbers >> pyqmforce.charges; + + pyqmforce.switch_width = 0.2; + pyqmforce.use_switch = true; } else - throw version_error(v, "1", r_pyqmforce, CODELOC); + throw version_error(v, "2", r_pyqmforce, CODELOC); return ds; } @@ -272,35 +337,39 @@ PyQMForce::PyQMForce( QMap bond_scale_factors, QVector mm2_atoms, QVector numbers, - QVector charges) : - callback(callback), - cutoff(cutoff), - neighbour_list_frequency(neighbour_list_frequency), - is_mechanical(is_mechanical), - lambda(lambda), - atoms(atoms), - mm1_to_qm(mm1_to_qm), - mm1_to_mm2(mm1_to_mm2), - bond_scale_factors(bond_scale_factors), - mm2_atoms(mm2_atoms), - numbers(numbers), - charges(charges) -{ -} - -PyQMForce::PyQMForce(const PyQMForce &other) : - callback(other.callback), - cutoff(other.cutoff), - neighbour_list_frequency(other.neighbour_list_frequency), - is_mechanical(other.is_mechanical), - lambda(other.lambda), - atoms(other.atoms), - mm1_to_qm(other.mm1_to_qm), - mm1_to_mm2(other.mm1_to_mm2), - mm2_atoms(other.mm2_atoms), - bond_scale_factors(other.bond_scale_factors), - numbers(other.numbers), - charges(other.charges) + QVector charges, + double switch_width, + bool use_switch) : callback(callback), + cutoff(cutoff), + neighbour_list_frequency(neighbour_list_frequency), + is_mechanical(is_mechanical), + lambda(lambda), + atoms(atoms), + mm1_to_qm(mm1_to_qm), + mm1_to_mm2(mm1_to_mm2), + bond_scale_factors(bond_scale_factors), + mm2_atoms(mm2_atoms), + numbers(numbers), + charges(charges), + switch_width(switch_width), + use_switch(use_switch) +{ +} + +PyQMForce::PyQMForce(const PyQMForce &other) : callback(other.callback), + cutoff(other.cutoff), + neighbour_list_frequency(other.neighbour_list_frequency), + is_mechanical(other.is_mechanical), + lambda(other.lambda), + atoms(other.atoms), + mm1_to_qm(other.mm1_to_qm), + mm1_to_mm2(other.mm1_to_mm2), + mm2_atoms(other.mm2_atoms), + bond_scale_factors(other.bond_scale_factors), + numbers(other.numbers), + charges(other.charges), + switch_width(other.switch_width), + use_switch(other.use_switch) { } @@ -318,6 +387,8 @@ PyQMForce &PyQMForce::operator=(const PyQMForce &other) this->bond_scale_factors = other.bond_scale_factors; this->numbers = other.numbers; this->charges = other.charges; + this->switch_width = other.switch_width; + this->use_switch = other.use_switch; return *this; } @@ -362,7 +433,8 @@ int PyQMForce::getNeighbourListFrequency() const bool PyQMForce::getIsMechanical() const { - return this->is_mechanical;; + return this->is_mechanical; + ; } QVector PyQMForce::getAtoms() const @@ -390,6 +462,16 @@ QVector PyQMForce::getCharges() const return this->charges; } +double PyQMForce::getSwitchWidth() const +{ + return this->switch_width; +} + +bool PyQMForce::getUseSwitch() const +{ + return this->use_switch; +} + const char *PyQMForce::typeName() { return QMetaType::typeName(qMetaTypeId()); @@ -412,75 +494,88 @@ PyQMForce::call( return this->callback.call(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell, idx_mm); } +boost::tuple>, QVector>, QVector> +PyQMForce::call4( + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const +{ + return this->callback.call4(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell, idx_mm); +} + ///////// ///////// OpenMM Serialization ///////// namespace OpenMM { - class PyQMForceProxy : public SerializationProxy { - public: - PyQMForceProxy() : SerializationProxy("PyQMForce") - { - }; + class PyQMForceProxy : public SerializationProxy + { + public: + PyQMForceProxy() : SerializationProxy("PyQMForce") { + }; - void serialize(const void* object, SerializationNode& node) const - { - // Serialize the object. - QByteArray data; - QDataStream ds(&data, QIODevice::WriteOnly); - PyQMForce pyqmforce = *static_cast(object); - ds << pyqmforce; + void serialize(const void *object, SerializationNode &node) const + { + // Serialize the object. + QByteArray data; + QDataStream ds(&data, QIODevice::WriteOnly); + PyQMForce pyqmforce = *static_cast(object); + ds << pyqmforce; - // Set the version. - node.setIntProperty("version", 0); + // Set the version. + node.setIntProperty("version", 0); - // Set the note attribute. - node.setStringProperty("note", - "This force only supports partial serialization, so can only be used " - "within the same session and memory space."); + // Set the note attribute. + node.setStringProperty("note", + "This force only supports partial serialization, so can only be used " + "within the same session and memory space."); - // Set the data by converting the QByteArray to a hexidecimal string. - node.setStringProperty("data", data.toHex().data()); - }; + // Set the data by converting the QByteArray to a hexidecimal string. + node.setStringProperty("data", data.toHex().data()); + }; - void* deserialize(const SerializationNode& node) const + void *deserialize(const SerializationNode &node) const + { + // Check the version. + int version = node.getIntProperty("version"); + if (version != 0) { - // Check the version. - int version = node.getIntProperty("version"); - if (version != 0) - { - throw OpenMM::OpenMMException("Unsupported version number"); - } + throw OpenMM::OpenMMException("Unsupported version number"); + } - // Get the data as a std::string. - auto string = node.getStringProperty("data"); + // Get the data as a std::string. + auto string = node.getStringProperty("data"); - // Convert to hexidecimal. - auto hex = QByteArray::fromRawData(string.data(), string.size()); + // Convert to hexidecimal. + auto hex = QByteArray::fromRawData(string.data(), string.size()); - // Convert to a QByteArray. - auto data = QByteArray::fromHex(hex); + // Convert to a QByteArray. + auto data = QByteArray::fromHex(hex); - // Deserialize the object. - QDataStream ds(data); - PyQMForce pyqmforce; + // Deserialize the object. + QDataStream ds(data); + PyQMForce pyqmforce; - try - { - ds >> pyqmforce; - } - catch (...) - { - throw OpenMM::OpenMMException("Unable to find UUID in the PyQMForce callback registry."); - } + try + { + ds >> pyqmforce; + } + catch (...) + { + throw OpenMM::OpenMMException("Unable to find UUID in the PyQMForce callback registry."); + } - return new PyQMForce(pyqmforce); - }; + return new PyQMForce(pyqmforce); + }; }; // Register the PyQMForce serialization proxy. - extern "C" void registerPyQMSerializationProxies() { + extern "C" void registerPyQMSerializationProxies() + { SerializationProxy::registerProxy(typeid(PyQMForce), new PyQMForceProxy()); } }; @@ -503,9 +598,8 @@ OpenMM::ForceImpl *PyQMForce::createImpl() const } #ifdef SIRE_USE_CUSTOMCPPFORCE -PyQMForceImpl::PyQMForceImpl(const PyQMForce &owner) : - OpenMM::CustomCPPForceImpl(owner), - owner(owner) +PyQMForceImpl::PyQMForceImpl(const PyQMForce &owner) : OpenMM::CustomCPPForceImpl(owner), + owner(owner) { } @@ -530,13 +624,17 @@ double PyQMForceImpl::computeForce( this->cutoff = this->owner.getCutoff().value(); // The neighbour list cutoff is 20% larger than the cutoff. - this->neighbour_list_cutoff = 1.2*this->cutoff; + this->neighbour_list_cutoff = 1.2 * this->cutoff; // Store the neighbour list update frequency. this->neighbour_list_frequency = this->owner.getNeighbourListFrequency(); // Flag whether a neighbour list is used. this->is_neighbour_list = this->neighbour_list_frequency > 0; + + // Cache switching function parameters. + this->use_switch = this->owner.getUseSwitch(); + this->r_switch = (1.0 - this->owner.getSwitchWidth()) * this->cutoff; } // Get the current box vectors in nanometers. @@ -545,17 +643,15 @@ double PyQMForceImpl::computeForce( // Create a triclinic space, converting to Angstrom. TriclinicBox space( - Vector(10*box_x[0], 10*box_x[1], 10*box_x[2]), - Vector(10*box_y[0], 10*box_y[1], 10*box_y[2]), - Vector(10*box_z[0], 10*box_z[1], 10*box_z[2]) - ); + Vector(10 * box_x[0], 10 * box_x[1], 10 * box_x[2]), + Vector(10 * box_y[0], 10 * box_y[1], 10 * box_y[2]), + Vector(10 * box_z[0], 10 * box_z[1], 10 * box_z[2])); // Store the cell vectors in Angstrom. QVector> cell = { - {10*box_x[0], 10*box_x[1], 10*box_x[2]}, - {10*box_y[0], 10*box_y[1], 10*box_y[2]}, - {10*box_z[0], 10*box_z[1], 10*box_z[2]} - }; + {10 * box_x[0], 10 * box_x[1], 10 * box_x[2]}, + {10 * box_y[0], 10 * box_y[1], 10 * box_y[2]}, + {10 * box_z[0], 10 * box_z[1], 10 * box_z[2]}}; // Store the QM atomic indices and numbers. auto qm_atoms = this->owner.getAtoms(); @@ -578,12 +674,12 @@ double PyQMForceImpl::computeForce( for (const auto &idx : qm_atoms) { const auto &pos = positions[idx]; - Vector qm_vec(10*pos[0], 10*pos[1], 10*pos[2]); + Vector qm_vec(10 * pos[0], 10 * pos[1], 10 * pos[2]); xyz_qm_vec[i] = qm_vec; i++; } - // Next sure that the QM atoms are whole (unwrapped). + // Make sure that the QM atoms are whole (unwrapped). xyz_qm_vec = space.makeWhole(xyz_qm_vec); // Get the center of the QM atoms. We will use this as a reference when @@ -608,6 +704,23 @@ double PyQMForceImpl::computeForce( // Store the current number of MM atoms. unsigned int num_mm = 0; + // ds/dr of the quintic switching function (zero outside the switching region). + // Defined at outer scope so the chain-rule correction loop can use it after + // the electrostatic-embedding block closes. + auto switch_deriv = [&](double r) -> double + { + if (not this->use_switch or r <= this->r_switch or r >= this->cutoff) + return 0.0; + const double x = (r - this->r_switch) / (this->cutoff - this->r_switch); + return -30.0 * x * x * (x - 1.0) * (x - 1.0) / (this->cutoff - this->r_switch); + }; + + // Unscaled charges and chain-rule data for accepted MM atoms. + // Declared at outer scope for the same reason as switch_deriv above. + QVector charges_unscaled; + QVector min_dists; + QVector nearest_qm_vecs; + // If we are using electrostatic embedding, the work out the MM point charges and // build the neighbour list. if (not this->owner.getIsMechanical()) @@ -617,7 +730,19 @@ double PyQMForceImpl::computeForce( QVector> xyz_virtual; QVector charges_virtual; - // Manually work out the MM point charges and build the neigbour list. + // Quintic switching function: scales charges smoothly to zero at the cutoff. + // Continuous through second derivative; r_switch and use_switch cached at step 0. + auto switching_function = [&](double r) -> double + { + if (not this->use_switch or r <= this->r_switch) + return 1.0; + if (r >= this->cutoff) + return 0.0; + const double x = (r - this->r_switch) / (this->cutoff - this->r_switch); + return 1.0 - x * x * x * (6.0 * x * x - 15.0 * x + 10.0); + }; + + // Manually work out the MM point charges and build the neighbour list. if (not this->is_neighbour_list or this->step_count % this->neighbour_list_frequency == 0) { // Clear the neighbour list. @@ -636,36 +761,48 @@ double PyQMForceImpl::computeForce( not mm2_atoms.contains(i)) { // Store the MM atom position in Sire Vector format. - Vector mm_vec(10*pos[0], 10*pos[1], 10*pos[2]); + Vector mm_vec(10 * pos[0], 10 * pos[1], 10 * pos[2]); + + // Find the minimum distance to any QM atom. + double min_dist = std::numeric_limits::max(); + Vector nearest_qm_vec; // Loop over all of the QM atoms. for (const auto &qm_vec : xyz_qm_vec) { - // Work out the distance between the current MM atom and QM atoms. + // Work out the distance between the current MM atom and QM atom. const auto dist = space.calcDist(mm_vec, qm_vec); + if (dist < min_dist) + { + min_dist = dist; + nearest_qm_vec = qm_vec; + } + // The current MM atom is within the neighbour list cutoff. if (this->is_neighbour_list and dist < this->neighbour_list_cutoff) { // Insert the MM atom index into the neighbour list. this->neighbour_list.insert(i); } + } - // The current MM atom is within the cutoff, add it. - if (dist < cutoff) - { - // Work out the minimum image position with respect to the - // reference position and add to the vector. - mm_vec = space.getMinimumImage(mm_vec, center); - xyz_mm.append(QVector({mm_vec[0], mm_vec[1], mm_vec[2]})); + // The current MM atom is within the cutoff: add it. + if (min_dist < cutoff) + { + // Work out the minimum image position with respect to the + // reference position and add to the vector. + mm_vec = space.getMinimumImage(mm_vec, center); + xyz_mm.append(QVector({mm_vec[0], mm_vec[1], mm_vec[2]})); - // Add the charge and index. - charges_mm.append(this->owner.getCharges()[i]); - idx_mm.append(i); + const double q = this->owner.getCharges()[i]; + charges_unscaled.append(q); + min_dists.append(min_dist); + nearest_qm_vecs.append(nearest_qm_vec); - // Exit the inner loop. - break; - } + // Scale charge by switching function. + charges_mm.append(q * switching_function(min_dist)); + idx_mm.append(i); } } @@ -680,41 +817,51 @@ double PyQMForceImpl::computeForce( for (const auto &idx : this->neighbour_list) { // Store the MM atom position in Sire Vector format. - Vector mm_vec(10*positions[idx][0], 10*positions[idx][1], 10*positions[idx][2]); + Vector mm_vec(10 * positions[idx][0], 10 * positions[idx][1], 10 * positions[idx][2]); + + // Find the minimum distance to any QM atom. + double min_dist = std::numeric_limits::max(); + Vector nearest_qm_vec; - // Loop over all of the QM atoms. for (const auto &qm_vec : xyz_qm_vec) { - // The current MM atom is within the cutoff, add it. - if (space.calcDist(mm_vec, qm_vec) < cutoff) + const auto dist = space.calcDist(mm_vec, qm_vec); + if (dist < min_dist) { - // Work out the minimum image position with respect to the - // reference position and add to the vector. - mm_vec = space.getMinimumImage(mm_vec, center); - xyz_mm.append(QVector({mm_vec[0], mm_vec[1], mm_vec[2]})); + min_dist = dist; + nearest_qm_vec = qm_vec; + } + } - // Add the charge and index. - charges_mm.append(this->owner.getCharges()[idx]); - idx_mm.append(idx); + // The current MM atom is within the cutoff: add it. + if (min_dist < cutoff) + { + mm_vec = space.getMinimumImage(mm_vec, center); + xyz_mm.append(QVector({mm_vec[0], mm_vec[1], mm_vec[2]})); - // Exit the inner loop. - break; - } + const double q = this->owner.getCharges()[idx]; + charges_unscaled.append(q); + min_dists.append(min_dist); + nearest_qm_vecs.append(nearest_qm_vec); + + // Scale charge by switching function. + charges_mm.append(q * switching_function(min_dist)); + idx_mm.append(idx); } } } // Handle link atoms via the Charge Shift method. // See: https://www.ks.uiuc.edu/Research/qmmm - for (const auto &idx: mm1_to_mm2.keys()) + for (const auto &idx : mm1_to_mm2.keys()) { // Get the QM atom to which the current MM atom is bonded. const auto qm_idx = mm1_to_qm[idx]; // Store the MM1 position in Sire Vector format, along with the // position of the QM atom to which it is bonded. - Vector mm1_vec(10*positions[idx][0], 10*positions[idx][1], 10*positions[idx][2]); - Vector qm_vec(10*positions[qm_idx][0], 10*positions[qm_idx][1], 10*positions[qm_idx][2]); + Vector mm1_vec(10 * positions[idx][0], 10 * positions[idx][1], 10 * positions[idx][2]); + Vector qm_vec(10 * positions[qm_idx][0], 10 * positions[qm_idx][1], 10 * positions[qm_idx][2]); // Work out the minimum image positions with respect to the reference position. mm1_vec = space.getMinimumImage(mm1_vec, center); @@ -725,7 +872,7 @@ double PyQMForceImpl::computeForce( // where R0(QM-L) is the equilibrium bond length for the QM and link (L) // elements, and R0(QM-MM1) is the equilibrium bond length for the QM // and MM1 elements. - const auto link_vec = qm_vec + bond_scale_factors[idx]*(mm1_vec - qm_vec); + const auto link_vec = qm_vec + bond_scale_factors[idx] * (mm1_vec - qm_vec); // Add to the QM positions. xyz_qm.append(QVector({link_vec[0], link_vec[1], link_vec[2]})); @@ -747,10 +894,10 @@ double PyQMForceImpl::computeForce( // charge is redistributed over the MM2 atoms and two virtual point // charges are added either side of the MM2 atoms in order to preserve // the MM1-MM2 dipole. - for (const auto& mm2_idx : mm1_to_mm2[idx]) + for (const auto &mm2_idx : mm1_to_mm2[idx]) { // Store the MM2 position in Sire Vector format. - Vector mm2_vec(10*positions[mm2_idx][0], 10*positions[mm2_idx][1], 10*positions[mm2_idx][2]); + Vector mm2_vec(10 * positions[mm2_idx][0], 10 * positions[mm2_idx][1], 10 * positions[mm2_idx][2]); // Work out the minimum image position with respect to the reference position. mm2_vec = space.getMinimumImage(mm2_vec, center); @@ -768,12 +915,12 @@ double PyQMForceImpl::computeForce( const auto normal = (mm2_vec - mm1_vec).normalise(); // Positive direction. (Away from MM1 atom.) - auto xyz = mm2_vec + VIRTUAL_PC_DELTA*normal; + auto xyz = mm2_vec + VIRTUAL_PC_DELTA * normal; xyz_virtual.append(QVector({xyz[0], xyz[1], xyz[2]})); charges_virtual.append(-frac_charge); // Negative direction (Towards MM1 atom.) - xyz = mm2_vec - VIRTUAL_PC_DELTA*normal; + xyz = mm2_vec - VIRTUAL_PC_DELTA * normal; xyz_virtual.append(QVector({xyz[0], xyz[1], xyz[2]})); charges_virtual.append(frac_charge); } @@ -791,20 +938,20 @@ double PyQMForceImpl::computeForce( } } - // Call the callback. - auto result = this->owner.call( + // Call the callback, requesting the optional dE/dcharges_mm fourth element. + auto result = this->owner.call4( numbers, charges_mm, xyz_qm, xyz_mm, cell, - idx_mm - ); + idx_mm); // Extract the results. These will automatically be returned in OpenMM units. auto energy = result.get<0>(); auto forces_qm = result.get<1>(); auto forces_mm = result.get<2>(); + const auto dE_dcharges_mm = result.get<3>(); // The current interpolation (weighting) parameter. double lambda; @@ -841,7 +988,7 @@ double PyQMForceImpl::computeForce( // Now update the force vector. // First the QM atoms. - for (int i=0; i(), - callback(py_object, name), - cutoff(cutoff), - neighbour_list_frequency(neighbour_list_frequency), - is_mechanical(is_mechanical), - lambda(lambda) + double lambda, + double switch_width, + bool use_switch) : ConcreteProperty(), + callback(py_object, name), + cutoff(cutoff), + neighbour_list_frequency(neighbour_list_frequency), + is_mechanical(is_mechanical), + lambda(lambda), + switch_width(switch_width), + use_switch(use_switch) { // Register the serialization proxies. OpenMM::registerPyQMSerializationProxies(); @@ -915,19 +1087,20 @@ PyQMEngine::PyQMEngine( } } -PyQMEngine::PyQMEngine(const PyQMEngine &other) : - callback(other.callback), - cutoff(other.cutoff), - neighbour_list_frequency(other.neighbour_list_frequency), - is_mechanical(other.is_mechanical), - lambda(other.lambda), - atoms(other.atoms), - mm1_to_qm(other.mm1_to_qm), - mm1_to_mm2(other.mm1_to_mm2), - mm2_atoms(other.mm2_atoms), - bond_scale_factors(other.bond_scale_factors), - numbers(other.numbers), - charges(other.charges) +PyQMEngine::PyQMEngine(const PyQMEngine &other) : callback(other.callback), + cutoff(other.cutoff), + neighbour_list_frequency(other.neighbour_list_frequency), + is_mechanical(other.is_mechanical), + lambda(other.lambda), + switch_width(other.switch_width), + use_switch(other.use_switch), + atoms(other.atoms), + mm1_to_qm(other.mm1_to_qm), + mm1_to_mm2(other.mm1_to_mm2), + mm2_atoms(other.mm2_atoms), + bond_scale_factors(other.bond_scale_factors), + numbers(other.numbers), + charges(other.charges) { } @@ -938,6 +1111,8 @@ PyQMEngine &PyQMEngine::operator=(const PyQMEngine &other) this->neighbour_list_frequency = other.neighbour_list_frequency; this->is_mechanical = other.is_mechanical; this->lambda = other.lambda; + this->switch_width = other.switch_width; + this->use_switch = other.use_switch; this->atoms = other.atoms; this->mm1_to_qm = other.mm1_to_qm; this->mm1_to_mm2 = other.mm1_to_mm2; @@ -1069,6 +1244,30 @@ void PyQMEngine::setCharges(QVector charges) this->charges = charges; } +double PyQMEngine::getSwitchWidth() const +{ + return this->switch_width; +} + +void PyQMEngine::setSwitchWidth(double switch_width) +{ + if (switch_width < 0.0) + switch_width = 0.0; + else if (switch_width > 1.0) + switch_width = 1.0; + this->switch_width = switch_width; +} + +bool PyQMEngine::getUseSwitch() const +{ + return this->use_switch; +} + +void PyQMEngine::setUseSwitch(bool use_switch) +{ + this->use_switch = use_switch; +} + const char *PyQMEngine::typeName() { return QMetaType::typeName(qMetaTypeId()); @@ -1091,7 +1290,7 @@ PyQMEngine::call( return this->callback.call(numbers_qm, charges_mm, xyz_qm, xyz_mm, cell, idx_mm); } -QMForce* PyQMEngine::createForce() const +QMForce *PyQMEngine::createForce() const { return new PyQMForce( this->callback, @@ -1105,6 +1304,7 @@ QMForce* PyQMEngine::createForce() const this->bond_scale_factors, this->mm2_atoms, this->numbers, - this->charges - ); + this->charges, + this->switch_width, + this->use_switch); } diff --git a/wrapper/Convert/SireOpenMM/pyqm.h b/wrapper/Convert/SireOpenMM/pyqm.h index 5b149add6..d26777885 100644 --- a/wrapper/Convert/SireOpenMM/pyqm.h +++ b/wrapper/Convert/SireOpenMM/pyqm.h @@ -96,7 +96,7 @@ namespace SireOpenMM - A list of forces for the MM atoms in kJ/mol/nm. If empty, then the object is assumed to be a callable. */ - PyQMCallback(bp::object, QString name=""); + PyQMCallback(bp::object, QString name = ""); //! Call the callback function. /*! \param numbers_qm @@ -126,13 +126,27 @@ namespace SireOpenMM - A vector of forces for the MM atoms in kJ/mol/nm. */ boost::tuple>, QVector>> call( - QVector numbers_qm, - QVector charges_mm, - QVector> xyz_qm, - QVector> xyz_mm, - QVector> cell, - QVector idx_mm - ) const; + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const; + + //! Call the callback and return an optional dE/dcharges_mm fourth element. + /*! Same arguments as call(). The 4th element of the returned tuple is the + gradient of the energy w.r.t. the effective MM charges (dE/dq_eff). + If the callback returns only 3 elements the 4th is an empty vector, + providing backward compatibility with existing callbacks. + */ + boost::tuple>, QVector>, QVector> + call4( + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const; //! Return the C++ name for this class. static const char *typeName(); @@ -212,8 +226,9 @@ namespace SireOpenMM QMap bond_scale_factors, QVector mm2_atoms, QVector numbers, - QVector charges - ); + QVector charges, + double switch_width = 0.2, + bool use_switch = true); //! Copy constructor. PyQMForce(const PyQMForce &other); @@ -309,6 +324,18 @@ namespace SireOpenMM */ QVector getCharges() const; + //! Get the switch width. + /*! \returns + The switch width as a fraction of the cutoff. + */ + double getSwitchWidth() const; + + //! Get whether a switching function is used. + /*! \returns + Whether a switching function is used. + */ + bool getUseSwitch() const; + //! Return the C++ name for this class. static const char *typeName(); @@ -316,40 +343,24 @@ namespace SireOpenMM const char *what() const; //! Call the callback function. - /*! \param numbers_qm - A vector of atomic numbers for the atoms in the ML region. - - \param charges_mm - A vector of the charges on the MM atoms in mod electron charge. - - \param xyz_qm - A vector of positions for the atoms in the ML region in Angstrom. - - \param xyz_mm - A vector of positions for the atoms in the MM region in Angstrom. - - \param cell - A vector of the 3 cell vectors in Angstrom. - - \param idx_mm - A vector of MM atom indices. Note that len(idx_mm) <= len(charges_mm) - since it only contains the indices of true MM atoms, not link atoms - or virtual charges. - - \returns - A tuple containing: - - The energy in kJ/mol. - - A vector of forces for the QM atoms in kJ/mol/nm. - - A vector of forces for the MM atoms in kJ/mol/nm. - */ + /*! \returns A tuple: (energy, forces_qm, forces_mm). */ boost::tuple>, QVector>> call( - QVector numbers_qm, - QVector charges_mm, - QVector> xyz_qm, - QVector> xyz_mm, - QVector> cell, - QVector idx_mm - ) const; + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const; + + //! Call the callback, returning optional dE/dcharges_mm as fourth element. + boost::tuple>, QVector>, QVector> + call4( + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const; protected: OpenMM::ForceImpl *createImpl() const; @@ -360,6 +371,8 @@ namespace SireOpenMM int neighbour_list_frequency; bool is_mechanical; double lambda; + double switch_width; + bool use_switch; QVector atoms; QMap mm1_to_qm; QMap> mm1_to_mm2; @@ -385,8 +398,10 @@ namespace SireOpenMM private: const PyQMForce &owner; - unsigned long long step_count=0; + unsigned long long step_count = 0; double cutoff; + double r_switch; + bool use_switch; bool is_neighbour_list; int neighbour_list_frequency; double neighbour_list_cutoff; @@ -424,12 +439,13 @@ namespace SireOpenMM */ PyQMEngine( bp::object, - QString method="", - SireUnits::Dimension::Length cutoff=7.5*SireUnits::angstrom, - int neighbour_list_frequency=0, - bool is_mechanical=false, - double lambda=1.0 - ); + QString method = "", + SireUnits::Dimension::Length cutoff = 7.5 * SireUnits::angstrom, + int neighbour_list_frequency = 0, + bool is_mechanical = false, + double lambda = 1.0, + double switch_width = 0.2, + bool use_switch = true); //! Copy constructor. PyQMEngine(const PyQMEngine &other); @@ -580,6 +596,30 @@ namespace SireOpenMM */ void setCharges(QVector charges); + //! Get the switch width. + /*! \returns + The switch width as a fraction of the cutoff. + */ + double getSwitchWidth() const; + + //! Set the switch width. + /*! \param switch_width + The switch width as a fraction of the cutoff (0 to 1). + */ + void setSwitchWidth(double switch_width); + + //! Get whether a switching function is used. + /*! \returns + Whether a switching function is used. + */ + bool getUseSwitch() const; + + //! Set whether a switching function is used. + /*! \param use_switch + Whether to use a switching function. + */ + void setUseSwitch(bool use_switch); + //! Return the C++ name for this class. static const char *typeName(); @@ -614,16 +654,15 @@ namespace SireOpenMM - A vector of forces for the MM atoms in kJ/mol/nm. */ boost::tuple>, QVector>> call( - QVector numbers_qm, - QVector charges_mm, - QVector> xyz_qm, - QVector> xyz_mm, - QVector> cell, - QVector idx_mm - ) const; + QVector numbers_qm, + QVector charges_mm, + QVector> xyz_qm, + QVector> xyz_mm, + QVector> cell, + QVector idx_mm) const; //! Create an EMLE force object. - QMForce* createForce() const; + QMForce *createForce() const; private: PyQMCallback callback; @@ -631,6 +670,8 @@ namespace SireOpenMM int neighbour_list_frequency; bool is_mechanical; double lambda; + double switch_width; + bool use_switch; QVector atoms; QMap mm1_to_qm; QMap> mm1_to_mm2; diff --git a/wrapper/Convert/SireOpenMM/torchqm.cpp b/wrapper/Convert/SireOpenMM/torchqm.cpp index 8b36db7c8..6f055d412 100644 --- a/wrapper/Convert/SireOpenMM/torchqm.cpp +++ b/wrapper/Convert/SireOpenMM/torchqm.cpp @@ -33,6 +33,8 @@ #include #endif +#include + #include "SireError/errors.h" #include "SireMaths/vector.h" #include "SireStream/datastream.h" @@ -61,14 +63,15 @@ static const RegisterMetaType r_torchqmforce(NO_ROOT); QDataStream &operator<<(QDataStream &ds, const TorchQMForce &torchqmforce) { - writeHeader(ds, r_torchqmforce, 1); + writeHeader(ds, r_torchqmforce, 2); SharedDataStream sds(ds); sds << torchqmforce.module_path << torchqmforce.cutoff << torchqmforce.neighbour_list_frequency << torchqmforce.is_mechanical << torchqmforce.lambda << torchqmforce.atoms << torchqmforce.mm1_to_qm << torchqmforce.mm1_to_mm2 << torchqmforce.bond_scale_factors - << torchqmforce.mm2_atoms << torchqmforce.numbers << torchqmforce.charges; + << torchqmforce.mm2_atoms << torchqmforce.numbers << torchqmforce.charges + << torchqmforce.switch_width << torchqmforce.use_switch; return ds; } @@ -77,20 +80,29 @@ QDataStream &operator>>(QDataStream &ds, TorchQMForce &torchqmforce) { VersionID v = readHeader(ds, r_torchqmforce); - if (v == 1) + if (v == 2) + { + SharedDataStream sds(ds); + + sds >> torchqmforce.module_path >> torchqmforce.cutoff >> torchqmforce.neighbour_list_frequency >> torchqmforce.is_mechanical >> torchqmforce.lambda >> torchqmforce.atoms >> torchqmforce.mm1_to_qm >> torchqmforce.mm1_to_mm2 >> torchqmforce.bond_scale_factors >> torchqmforce.mm2_atoms >> torchqmforce.numbers >> torchqmforce.charges >> torchqmforce.switch_width >> torchqmforce.use_switch; + + // Re-load the Torch module. + torchqmforce.setModulePath(torchqmforce.getModulePath()); + } + else if (v == 1) { SharedDataStream sds(ds); - sds >> torchqmforce.module_path >> torchqmforce.cutoff >> torchqmforce.neighbour_list_frequency - >> torchqmforce.is_mechanical >> torchqmforce.lambda >> torchqmforce.atoms - >> torchqmforce.mm1_to_qm >> torchqmforce.mm1_to_mm2 >> torchqmforce.bond_scale_factors - >> torchqmforce.mm2_atoms >> torchqmforce.numbers >> torchqmforce.charges; + sds >> torchqmforce.module_path >> torchqmforce.cutoff >> torchqmforce.neighbour_list_frequency >> torchqmforce.is_mechanical >> torchqmforce.lambda >> torchqmforce.atoms >> torchqmforce.mm1_to_qm >> torchqmforce.mm1_to_mm2 >> torchqmforce.bond_scale_factors >> torchqmforce.mm2_atoms >> torchqmforce.numbers >> torchqmforce.charges; + + torchqmforce.switch_width = 0.2; + torchqmforce.use_switch = true; // Re-load the Torch module. torchqmforce.setModulePath(torchqmforce.getModulePath()); } else - throw version_error(v, "1", r_torchqmforce, CODELOC); + throw version_error(v, "2", r_torchqmforce, CODELOC); return ds; } @@ -111,37 +123,41 @@ TorchQMForce::TorchQMForce( QMap bond_scale_factors, QVector mm2_atoms, QVector numbers, - QVector charges) : - cutoff(cutoff), - neighbour_list_frequency(neighbour_list_frequency), - is_mechanical(is_mechanical), - lambda(lambda), - atoms(atoms), - mm1_to_qm(mm1_to_qm), - mm1_to_mm2(mm1_to_mm2), - bond_scale_factors(bond_scale_factors), - mm2_atoms(mm2_atoms), - numbers(numbers), - charges(charges) + QVector charges, + double switch_width, + bool use_switch) : cutoff(cutoff), + neighbour_list_frequency(neighbour_list_frequency), + is_mechanical(is_mechanical), + lambda(lambda), + atoms(atoms), + mm1_to_qm(mm1_to_qm), + mm1_to_mm2(mm1_to_mm2), + bond_scale_factors(bond_scale_factors), + mm2_atoms(mm2_atoms), + numbers(numbers), + charges(charges), + switch_width(switch_width), + use_switch(use_switch) { // Try to load the Torch module. this->setModulePath(module_path); } -TorchQMForce::TorchQMForce(const TorchQMForce &other) : - module_path(other.module_path), - torch_module(other.torch_module), - cutoff(other.cutoff), - neighbour_list_frequency(other.neighbour_list_frequency), - is_mechanical(other.is_mechanical), - lambda(other.lambda), - atoms(other.atoms), - mm1_to_qm(other.mm1_to_qm), - mm1_to_mm2(other.mm1_to_mm2), - mm2_atoms(other.mm2_atoms), - bond_scale_factors(other.bond_scale_factors), - numbers(other.numbers), - charges(other.charges) +TorchQMForce::TorchQMForce(const TorchQMForce &other) : module_path(other.module_path), + torch_module(other.torch_module), + cutoff(other.cutoff), + neighbour_list_frequency(other.neighbour_list_frequency), + is_mechanical(other.is_mechanical), + lambda(other.lambda), + atoms(other.atoms), + mm1_to_qm(other.mm1_to_qm), + mm1_to_mm2(other.mm1_to_mm2), + mm2_atoms(other.mm2_atoms), + bond_scale_factors(other.bond_scale_factors), + numbers(other.numbers), + charges(other.charges), + switch_width(other.switch_width), + use_switch(other.use_switch) { } @@ -160,6 +176,8 @@ TorchQMForce &TorchQMForce::operator=(const TorchQMForce &other) this->bond_scale_factors = other.bond_scale_factors; this->numbers = other.numbers; this->charges = other.charges; + this->switch_width = other.switch_width; + this->use_switch = other.use_switch; return *this; } @@ -173,13 +191,14 @@ void TorchQMForce::setModulePath(QString module_path) this->torch_module = torch::jit::load(module_path.toStdString()); this->torch_module.eval(); } - catch (const c10::Error& e) + catch (const c10::Error &e) { throw SireError::io_error( - QObject::tr( - "Unable to load the TorchScript module '%1'. The error was '%2'.") - .arg(module_path).arg(e.what()), - CODELOC); + QObject::tr( + "Unable to load the TorchScript module '%1'. The error was '%2'.") + .arg(module_path) + .arg(e.what()), + CODELOC); } #endif @@ -194,7 +213,7 @@ QString TorchQMForce::getModulePath() const #ifdef SIRE_USE_TORCH torch::jit::script::Module TorchQMForce::getTorchModule() const #else -void* TorchQMForce::getTorchModule() const +void *TorchQMForce::getTorchModule() const #endif { return this->torch_module; @@ -231,7 +250,8 @@ int TorchQMForce::getNeighbourListFrequency() const bool TorchQMForce::getIsMechanical() const { - return this->is_mechanical;; + return this->is_mechanical; + ; } QVector TorchQMForce::getAtoms() const @@ -259,6 +279,16 @@ QVector TorchQMForce::getCharges() const return this->charges; } +double TorchQMForce::getSwitchWidth() const +{ + return this->switch_width; +} + +bool TorchQMForce::getUseSwitch() const +{ + return this->use_switch; +} + const char *TorchQMForce::typeName() { return QMetaType::typeName(qMetaTypeId()); @@ -275,69 +305,70 @@ const char *TorchQMForce::what() const namespace OpenMM { - class TorchQMForceProxy : public SerializationProxy { - public: - TorchQMForceProxy() : SerializationProxy("TorchQMForce") - { - }; + class TorchQMForceProxy : public SerializationProxy + { + public: + TorchQMForceProxy() : SerializationProxy("TorchQMForce") { + }; - void serialize(const void* object, SerializationNode& node) const - { - // Serialize the object. - QByteArray data; - QDataStream ds(&data, QIODevice::WriteOnly); - TorchQMForce torchqmforce = *static_cast(object); - ds << torchqmforce; + void serialize(const void *object, SerializationNode &node) const + { + // Serialize the object. + QByteArray data; + QDataStream ds(&data, QIODevice::WriteOnly); + TorchQMForce torchqmforce = *static_cast(object); + ds << torchqmforce; - // Set the version. - node.setIntProperty("version", 0); + // Set the version. + node.setIntProperty("version", 0); - // Set the note attribute. - node.setStringProperty("note", - "This force only supports partial serialization, so can only be used " - "within the same session and memory space."); + // Set the note attribute. + node.setStringProperty("note", + "This force only supports partial serialization, so can only be used " + "within the same session and memory space."); - // Set the data by converting the QByteArray to a hexidecimal string. - node.setStringProperty("data", data.toHex().data()); - }; + // Set the data by converting the QByteArray to a hexidecimal string. + node.setStringProperty("data", data.toHex().data()); + }; - void* deserialize(const SerializationNode& node) const + void *deserialize(const SerializationNode &node) const + { + // Check the version. + int version = node.getIntProperty("version"); + if (version != 0) { - // Check the version. - int version = node.getIntProperty("version"); - if (version != 0) - { - throw OpenMM::OpenMMException("Unsupported version number"); - } + throw OpenMM::OpenMMException("Unsupported version number"); + } - // Get the data as a std::string. - auto string = node.getStringProperty("data"); + // Get the data as a std::string. + auto string = node.getStringProperty("data"); - // Convert to hexidecimal. - auto hex = QByteArray::fromRawData(string.data(), string.size()); + // Convert to hexidecimal. + auto hex = QByteArray::fromRawData(string.data(), string.size()); - // Convert to a QByteArray. - auto data = QByteArray::fromHex(hex); + // Convert to a QByteArray. + auto data = QByteArray::fromHex(hex); - // Deserialize the object. - QDataStream ds(data); - TorchQMForce torchqmforce; + // Deserialize the object. + QDataStream ds(data); + TorchQMForce torchqmforce; - try - { - ds >> torchqmforce; - } - catch (...) - { - throw OpenMM::OpenMMException("Unable deserialize TorchQMForce"); - } + try + { + ds >> torchqmforce; + } + catch (...) + { + throw OpenMM::OpenMMException("Unable deserialize TorchQMForce"); + } - return new TorchQMForce(torchqmforce); - }; + return new TorchQMForce(torchqmforce); + }; }; // Register the TorchQMForce serialization proxy. - extern "C" void registerTorchQMSerializationProxies() { + extern "C" void registerTorchQMSerializationProxies() + { SerializationProxy::registerProxy(typeid(TorchQMForce), new TorchQMForceProxy()); } }; @@ -360,9 +391,8 @@ OpenMM::ForceImpl *TorchQMForce::createImpl() const } #if defined(SIRE_USE_CUSTOMCPPFORCE) and defined(SIRE_USE_TORCH) -TorchQMForceImpl::TorchQMForceImpl(const TorchQMForce &owner) : - OpenMM::CustomCPPForceImpl(owner), - owner(owner) +TorchQMForceImpl::TorchQMForceImpl(const TorchQMForce &owner) : OpenMM::CustomCPPForceImpl(owner), + owner(owner) { this->torch_module = owner.getTorchModule(); } @@ -409,13 +439,17 @@ double TorchQMForceImpl::computeForce( this->cutoff = this->owner.getCutoff().value(); // The neighbour list cutoff is 20% larger than the cutoff. - this->neighbour_list_cutoff = 1.2*this->cutoff; + this->neighbour_list_cutoff = 1.2 * this->cutoff; // Store the neighbour list update frequency. this->neighbour_list_frequency = this->owner.getNeighbourListFrequency(); // Flag whether a neighbour list is used. this->is_neighbour_list = this->neighbour_list_frequency > 0; + + // Cache switching function parameters. + this->use_switch = this->owner.getUseSwitch(); + this->r_switch = (1.0 - this->owner.getSwitchWidth()) * this->cutoff; } // Get the current box vectors in nanometers. @@ -424,17 +458,15 @@ double TorchQMForceImpl::computeForce( // Create a triclinic space, converting to Angstrom. TriclinicBox space( - Vector(10*box_x[0], 10*box_x[1], 10*box_x[2]), - Vector(10*box_y[0], 10*box_y[1], 10*box_y[2]), - Vector(10*box_z[0], 10*box_z[1], 10*box_z[2]) - ); + Vector(10 * box_x[0], 10 * box_x[1], 10 * box_x[2]), + Vector(10 * box_y[0], 10 * box_y[1], 10 * box_y[2]), + Vector(10 * box_z[0], 10 * box_z[1], 10 * box_z[2])); // Store the cell vectors in Angstrom. QVector> cell = { - {10*box_x[0], 10*box_x[1], 10*box_x[2]}, - {10*box_y[0], 10*box_y[1], 10*box_y[2]}, - {10*box_z[0], 10*box_z[1], 10*box_z[2]} - }; + {10 * box_x[0], 10 * box_x[1], 10 * box_x[2]}, + {10 * box_y[0], 10 * box_y[1], 10 * box_y[2]}, + {10 * box_z[0], 10 * box_z[1], 10 * box_z[2]}}; // Store the QM atomic indices and numbers. auto qm_atoms = this->owner.getAtoms(); @@ -450,19 +482,19 @@ double TorchQMForceImpl::computeForce( // Initialise a vector to hold the current positions for the QM atoms. QVector xyz_qm_vec(qm_atoms.size()); - std::vector xyz_qm(3*qm_atoms.size()); + std::vector xyz_qm(3 * qm_atoms.size()); // First loop over all QM atoms and store the positions. int i = 0; for (const auto &idx : qm_atoms) { const auto &pos = positions[idx]; - Vector qm_vec(10*pos[0], 10*pos[1], 10*pos[2]); + Vector qm_vec(10 * pos[0], 10 * pos[1], 10 * pos[2]); xyz_qm_vec[i] = qm_vec; i++; } - // Next sure that the QM atoms are whole (unwrapped). + // Make sure that the QM atoms are whole (unwrapped). xyz_qm_vec = space.makeWhole(xyz_qm_vec); // Get the center of the QM atoms. We will use this as a reference when @@ -471,9 +503,9 @@ double TorchQMForceImpl::computeForce( i = 0; for (const auto &qm_vec : xyz_qm_vec) { - xyz_qm[3*i] = qm_vec[0]; - xyz_qm[3*i+1] = qm_vec[1]; - xyz_qm[3*i+2] = qm_vec[2]; + xyz_qm[3 * i] = qm_vec[0]; + xyz_qm[3 * i + 1] = qm_vec[1]; + xyz_qm[3 * i + 2] = qm_vec[2]; center += qm_vec; i++; } @@ -489,6 +521,23 @@ double TorchQMForceImpl::computeForce( // Store the current number of MM atoms. unsigned int num_mm = 0; + // ds/dr of the quintic switching function (zero outside the switching region). + // Defined at outer scope so the chain-rule correction loop can use it after + // the electrostatic-embedding block closes. + auto switch_deriv = [&](double r) -> double + { + if (not this->use_switch or r <= this->r_switch or r >= this->cutoff) + return 0.0; + const double x = (r - this->r_switch) / (this->cutoff - this->r_switch); + return -30.0 * x * x * (x - 1.0) * (x - 1.0) / (this->cutoff - this->r_switch); + }; + + // Unscaled charges and chain-rule data for accepted MM atoms. + // Declared at outer scope for the same reason as switch_deriv above. + QVector charges_unscaled; + QVector min_dists; + QVector nearest_qm_vecs; + // If we are using electrostatic embedding, the work out the MM point charges and // build the neighbour list. if (not this->owner.getIsMechanical()) @@ -498,7 +547,19 @@ double TorchQMForceImpl::computeForce( std::vector xyz_virtual; QVector charges_virtual; - // Manually work out the MM point charges and build the neigbour list. + // Quintic switching function: scales charges smoothly to zero at the cutoff. + // Continuous through second derivative; r_switch and use_switch cached at step 0. + auto switching_function = [&](double r) -> double + { + if (not this->use_switch or r <= this->r_switch) + return 1.0; + if (r >= this->cutoff) + return 0.0; + const double x = (r - this->r_switch) / (this->cutoff - this->r_switch); + return 1.0 - x * x * x * (6.0 * x * x - 15.0 * x + 10.0); + }; + + // Manually work out the MM point charges and build the neighbour list. if (not this->is_neighbour_list or this->step_count % this->neighbour_list_frequency == 0) { // Clear the neighbour list. @@ -517,38 +578,50 @@ double TorchQMForceImpl::computeForce( not mm2_atoms.contains(i)) { // Store the MM atom position in Sire Vector format. - Vector mm_vec(10*pos[0], 10*pos[1], 10*pos[2]); + Vector mm_vec(10 * pos[0], 10 * pos[1], 10 * pos[2]); + + // Find the minimum distance to any QM atom. + double min_dist = std::numeric_limits::max(); + Vector nearest_qm_vec; // Loop over all of the QM atoms. for (const auto &qm_vec : xyz_qm_vec) { - // Work out the distance between the current MM atom and QM atoms. + // Work out the distance between the current MM atom and QM atom. const auto dist = space.calcDist(mm_vec, qm_vec); + if (dist < min_dist) + { + min_dist = dist; + nearest_qm_vec = qm_vec; + } + // The current MM atom is within the neighbour list cutoff. if (this->is_neighbour_list and dist < this->neighbour_list_cutoff) { // Insert the MM atom index into the neighbour list. this->neighbour_list.insert(i); } + } - // The current MM atom is within the cutoff, add it. - if (dist < cutoff) - { - // Work out the minimum image position with respect to the - // reference position and add to the vector. - mm_vec = space.getMinimumImage(mm_vec, center); - xyz_mm.push_back(mm_vec[0]); - xyz_mm.push_back(mm_vec[1]); - xyz_mm.push_back(mm_vec[2]); - - // Add the charge and index. - charges_mm.append(this->owner.getCharges()[i]); - idx_mm.append(i); - - // Exit the inner loop. - break; - } + // The current MM atom is within the cutoff: add it. + if (min_dist < cutoff) + { + // Work out the minimum image position with respect to the + // reference position and add to the vector. + mm_vec = space.getMinimumImage(mm_vec, center); + xyz_mm.push_back(mm_vec[0]); + xyz_mm.push_back(mm_vec[1]); + xyz_mm.push_back(mm_vec[2]); + + const double q = this->owner.getCharges()[i]; + charges_unscaled.append(q); + min_dists.append(min_dist); + nearest_qm_vecs.append(nearest_qm_vec); + + // Scale charge by switching function. + charges_mm.append(q * switching_function(min_dist)); + idx_mm.append(i); } } @@ -563,43 +636,53 @@ double TorchQMForceImpl::computeForce( for (const auto &idx : this->neighbour_list) { // Store the MM atom position in Sire Vector format. - Vector mm_vec(10*positions[idx][0], 10*positions[idx][1], 10*positions[idx][2]); + Vector mm_vec(10 * positions[idx][0], 10 * positions[idx][1], 10 * positions[idx][2]); + + // Find the minimum distance to any QM atom. + double min_dist = std::numeric_limits::max(); + Vector nearest_qm_vec; - // Loop over all of the QM atoms. for (const auto &qm_vec : xyz_qm_vec) { - // The current MM atom is within the cutoff, add it. - if (space.calcDist(mm_vec, qm_vec) < cutoff) + const auto dist = space.calcDist(mm_vec, qm_vec); + if (dist < min_dist) { - // Work out the minimum image position with respect to the - // reference position and add to the vector. - mm_vec = space.getMinimumImage(mm_vec, center); - xyz_mm.push_back(mm_vec[0]); - xyz_mm.push_back(mm_vec[1]); - xyz_mm.push_back(mm_vec[2]); - - // Add the charge and index. - charges_mm.append(this->owner.getCharges()[idx]); - idx_mm.append(idx); - - // Exit the inner loop. - break; + min_dist = dist; + nearest_qm_vec = qm_vec; } } + + // The current MM atom is within the cutoff: add it. + if (min_dist < cutoff) + { + mm_vec = space.getMinimumImage(mm_vec, center); + xyz_mm.push_back(mm_vec[0]); + xyz_mm.push_back(mm_vec[1]); + xyz_mm.push_back(mm_vec[2]); + + const double q = this->owner.getCharges()[idx]; + charges_unscaled.append(q); + min_dists.append(min_dist); + nearest_qm_vecs.append(nearest_qm_vec); + + // Scale charge by switching function. + charges_mm.append(q * switching_function(min_dist)); + idx_mm.append(idx); + } } } // Handle link atoms via the Charge Shift method. // See: https://www.ks.uiuc.edu/Research/qmmm - for (const auto &idx: mm1_to_mm2.keys()) + for (const auto &idx : mm1_to_mm2.keys()) { // Get the QM atom to which the current MM atom is bonded. const auto qm_idx = mm1_to_qm[idx]; // Store the MM1 position in Sire Vector format, along with the // position of the QM atom to which it is bonded. - Vector mm1_vec(10*positions[idx][0], 10*positions[idx][1], 10*positions[idx][2]); - Vector qm_vec(10*positions[qm_idx][0], 10*positions[qm_idx][1], 10*positions[qm_idx][2]); + Vector mm1_vec(10 * positions[idx][0], 10 * positions[idx][1], 10 * positions[idx][2]); + Vector qm_vec(10 * positions[qm_idx][0], 10 * positions[qm_idx][1], 10 * positions[qm_idx][2]); // Work out the minimum image positions with respect to the reference position. mm1_vec = space.getMinimumImage(mm1_vec, center); @@ -610,7 +693,7 @@ double TorchQMForceImpl::computeForce( // where R0(QM-L) is the equilibrium bond length for the QM and link (L) // elements, and R0(QM-MM1) is the equilibrium bond length for the QM // and MM1 elements. - const auto link_vec = qm_vec + bond_scale_factors[idx]*(mm1_vec - qm_vec); + const auto link_vec = qm_vec + bond_scale_factors[idx] * (mm1_vec - qm_vec); // Add to the QM positions. xyz_qm.push_back(link_vec[0]); @@ -634,10 +717,10 @@ double TorchQMForceImpl::computeForce( // charge is redistributed over the MM2 atoms and two virtual point // charges are added either side of the MM2 atoms in order to preserve // the MM1-MM2 dipole. - for (const auto& mm2_idx : mm1_to_mm2[idx]) + for (const auto &mm2_idx : mm1_to_mm2[idx]) { // Store the MM2 position in Sire Vector format. - Vector mm2_vec(10*positions[mm2_idx][0], 10*positions[mm2_idx][1], 10*positions[mm2_idx][2]); + Vector mm2_vec(10 * positions[mm2_idx][0], 10 * positions[mm2_idx][1], 10 * positions[mm2_idx][2]); // Work out the minimum image position with respect to the reference position. mm2_vec = space.getMinimumImage(mm2_vec, center); @@ -657,14 +740,14 @@ double TorchQMForceImpl::computeForce( const auto normal = (mm2_vec - mm1_vec).normalise(); // Positive direction. (Away from MM1 atom.) - auto xyz = mm2_vec + VIRTUAL_PC_DELTA*normal; + auto xyz = mm2_vec + VIRTUAL_PC_DELTA * normal; xyz_virtual.push_back(xyz[0]); xyz_virtual.push_back(xyz[1]); xyz_virtual.push_back(xyz[2]); charges_virtual.append(-frac_charge); // Negative direction (Towards MM1 atom.) - xyz = mm2_vec - VIRTUAL_PC_DELTA*normal; + xyz = mm2_vec - VIRTUAL_PC_DELTA * normal; xyz_virtual.push_back(xyz[0]); xyz_virtual.push_back(xyz[1]); xyz_virtual.push_back(xyz[2]); @@ -694,38 +777,44 @@ double TorchQMForceImpl::computeForce( // Resize the charges and positions vectors to the maximum number of MM atoms. // This is to try to preserve a static compute graph to avoid re-jitting. charges_mm.resize(this->max_num_mm); - xyz_mm.resize(3*this->max_num_mm); + xyz_mm.resize(3 * this->max_num_mm); } } // Convert input to Torch tensors. - // MM charges. + // MM charges. requires_grad_(true) must be set before the forward pass so + // that the computation graph tracks the dependency on charges, enabling the + // chain-rule correction for the switching function. torch::Tensor charges_mm_torch = torch::from_blob(charges_mm.data(), {charges_mm.size()}, - torch::TensorOptions().dtype(torch::kFloat64)) - .to(torch::kFloat32).to(device); + torch::TensorOptions().dtype(torch::kFloat64)) + .to(torch::kFloat32) + .to(device); + charges_mm_torch.requires_grad_(true); // Atomic numbers. torch::Tensor atomic_numbers_torch = torch::from_blob(numbers.data(), {numbers.size()}, - torch::TensorOptions().dtype(torch::kInt32)) - .to(torch::kInt64).to(device); + torch::TensorOptions().dtype(torch::kInt32)) + .to(torch::kInt64) + .to(device); // QM positions. torch::Tensor xyz_qm_torch = torch::from_blob(xyz_qm.data(), {numbers.size(), 3}, - torch::TensorOptions().dtype(torch::kFloat32)) - .to(device); + torch::TensorOptions().dtype(torch::kFloat32)) + .to(device); xyz_qm_torch.requires_grad_(true); // MM positions. torch::Tensor xyz_mm_torch = torch::from_blob(xyz_mm.data(), {charges_mm.size(), 3}, - torch::TensorOptions().dtype(torch::kFloat32)) - .to(device); + torch::TensorOptions().dtype(torch::kFloat32)) + .to(device); xyz_mm_torch.requires_grad_(true); // Cell vectors. torch::Tensor cell_torch = torch::from_blob(cell.data(), {3, 3}, - torch::TensorOptions().dtype(torch::kFloat64)) - .to(torch::kFloat32).to(device); + torch::TensorOptions().dtype(torch::kFloat64)) + .to(torch::kFloat32) + .to(device); cell_torch.requires_grad_(false); // Create the input vector. @@ -734,8 +823,7 @@ double TorchQMForceImpl::computeForce( charges_mm_torch, xyz_qm_torch, xyz_mm_torch, - cell_torch - }; + cell_torch}; // Compute the energies. auto energies = this->torch_module.forward(input).toTensor(); @@ -743,19 +831,27 @@ double TorchQMForceImpl::computeForce( // Store the sum of the energy in kJ. const auto energy = energies.sum().item() * HARTREE_TO_KJ_MOL; - // If there are no MM atoms, then we need to allow unused tensors. - bool allow_unused = num_mm == 0; + // Always allow unused: xyz_mm/charges_mm may not be connected when there are + // no MM atoms, and the TorchScript model may not propagate gradients through + // charges even when they are in the graph. We guard each gradient before use. + const bool allow_unused = true; - // Compute the gradients. + // Compute the gradients w.r.t. QM positions, MM positions, and MM charges. const auto gradients = torch::autograd::grad( - {energies.sum()}, {xyz_qm_torch, xyz_mm_torch}, {}, c10::nullopt, false, allow_unused); + {energies.sum()}, {xyz_qm_torch, xyz_mm_torch, charges_mm_torch}, + {}, c10::nullopt, false, allow_unused); // Compute the forces, converting from Hatree/Anstrom to kJ/mol/nm. const auto forces_qm = -(gradients[0] * HARTREE_TO_KJ_MOL * 10).detach().cpu(); torch::Tensor forces_mm; + torch::Tensor dE_dq; if (num_mm > 0) { forces_mm = -(gradients[1] * HARTREE_TO_KJ_MOL * 10).detach().cpu(); + if (gradients[2].defined()) + { + dE_dq = gradients[2].detach().cpu(); + } } else { @@ -805,35 +901,55 @@ double TorchQMForceImpl::computeForce( forces_mm.data_ptr() + forces_mm.numel()); // First the QM atoms. - for (int i=0; i() * HARTREE_TO_KJ_MOL * 10.0 * charges_unscaled[i] * dsdr; + forces[idx] += lambda * OpenMM::Vec3( + correction * r_hat[0], + correction * r_hat[1], + correction * r_hat[2]); + } + } } // Update the step count. @@ -859,19 +975,22 @@ TorchQMEngine::TorchQMEngine( SireUnits::Dimension::Length cutoff, int neighbour_list_frequency, bool is_mechanical, - double lambda) : - ConcreteProperty(), - module_path(module_path), - cutoff(cutoff), - neighbour_list_frequency(neighbour_list_frequency), - is_mechanical(is_mechanical), - lambda(lambda) + double lambda, + double switch_width, + bool use_switch) : ConcreteProperty(), + module_path(module_path), + cutoff(cutoff), + neighbour_list_frequency(neighbour_list_frequency), + is_mechanical(is_mechanical), + lambda(lambda), + switch_width(switch_width), + use_switch(use_switch) { #ifndef SIRE_USE_TORCH throw SireError::unsupported(QObject::tr( - "Unable to create an TorchQMEngine because Sire has been compiled " - "without Torch support."), - CODELOC); + "Unable to create an TorchQMEngine because Sire has been compiled " + "without Torch support."), + CODELOC); #endif // Register the serialization proxies. @@ -891,19 +1010,20 @@ TorchQMEngine::TorchQMEngine( } } -TorchQMEngine::TorchQMEngine(const TorchQMEngine &other) : - module_path(other.module_path), - cutoff(other.cutoff), - neighbour_list_frequency(other.neighbour_list_frequency), - is_mechanical(other.is_mechanical), - lambda(other.lambda), - atoms(other.atoms), - mm1_to_qm(other.mm1_to_qm), - mm1_to_mm2(other.mm1_to_mm2), - mm2_atoms(other.mm2_atoms), - bond_scale_factors(other.bond_scale_factors), - numbers(other.numbers), - charges(other.charges) +TorchQMEngine::TorchQMEngine(const TorchQMEngine &other) : module_path(other.module_path), + cutoff(other.cutoff), + neighbour_list_frequency(other.neighbour_list_frequency), + is_mechanical(other.is_mechanical), + lambda(other.lambda), + switch_width(other.switch_width), + use_switch(other.use_switch), + atoms(other.atoms), + mm1_to_qm(other.mm1_to_qm), + mm1_to_mm2(other.mm1_to_mm2), + mm2_atoms(other.mm2_atoms), + bond_scale_factors(other.bond_scale_factors), + numbers(other.numbers), + charges(other.charges) { } @@ -914,6 +1034,8 @@ TorchQMEngine &TorchQMEngine::operator=(const TorchQMEngine &other) this->neighbour_list_frequency = other.neighbour_list_frequency; this->is_mechanical = other.is_mechanical; this->lambda = other.lambda; + this->switch_width = other.switch_width; + this->use_switch = other.use_switch; this->atoms = other.atoms; this->mm1_to_qm = other.mm1_to_qm; this->mm1_to_mm2 = other.mm1_to_mm2; @@ -1045,6 +1167,30 @@ void TorchQMEngine::setCharges(QVector charges) this->charges = charges; } +double TorchQMEngine::getSwitchWidth() const +{ + return this->switch_width; +} + +void TorchQMEngine::setSwitchWidth(double switch_width) +{ + if (switch_width < 0.0) + switch_width = 0.0; + else if (switch_width > 1.0) + switch_width = 1.0; + this->switch_width = switch_width; +} + +bool TorchQMEngine::getUseSwitch() const +{ + return this->use_switch; +} + +void TorchQMEngine::setUseSwitch(bool use_switch) +{ + this->use_switch = use_switch; +} + const char *TorchQMEngine::typeName() { return QMetaType::typeName(qMetaTypeId()); @@ -1055,7 +1201,7 @@ const char *TorchQMEngine::what() const return TorchQMEngine::typeName(); } -QMForce* TorchQMEngine::createForce() const +QMForce *TorchQMEngine::createForce() const { return new TorchQMForce( this->module_path, @@ -1069,6 +1215,7 @@ QMForce* TorchQMEngine::createForce() const this->bond_scale_factors, this->mm2_atoms, this->numbers, - this->charges - ); + this->charges, + this->switch_width, + this->use_switch); } diff --git a/wrapper/Convert/SireOpenMM/torchqm.h b/wrapper/Convert/SireOpenMM/torchqm.h index 97969d75c..a6b8e4ca2 100644 --- a/wrapper/Convert/SireOpenMM/torchqm.h +++ b/wrapper/Convert/SireOpenMM/torchqm.h @@ -136,8 +136,9 @@ namespace SireOpenMM QMap bond_scale_factors, QVector mm2_atoms, QVector numbers, - QVector charges - ); + QVector charges, + double switch_width = 0.2, + bool use_switch = true); //! Copy constructor. TorchQMForce(const TorchQMForce &other); @@ -164,7 +165,7 @@ namespace SireOpenMM #ifdef SIRE_USE_TORCH torch::jit::script::Module getTorchModule() const; #else - void* getTorchModule() const; + void *getTorchModule() const; #endif //! Get the lambda weighting factor. @@ -243,6 +244,18 @@ namespace SireOpenMM */ QVector getCharges() const; + //! Get the switch width. + /*! \returns + The switch width as a fraction of the cutoff. + */ + double getSwitchWidth() const; + + //! Get whether a switching function is used. + /*! \returns + Whether a switching function is used. + */ + bool getUseSwitch() const; + //! Return the C++ name for this class. static const char *typeName(); @@ -270,6 +283,8 @@ namespace SireOpenMM QVector mm2_atoms; QVector numbers; QVector charges; + double switch_width; + bool use_switch; }; #if defined(SIRE_USE_CUSTOMCPPFORCE) && defined(SIRE_USE_TORCH) @@ -289,13 +304,15 @@ namespace SireOpenMM private: const TorchQMForce &owner; torch::jit::script::Module torch_module; - unsigned long long step_count=0; + unsigned long long step_count = 0; double cutoff; + double r_switch; + bool use_switch; bool is_neighbour_list; int neighbour_list_frequency; double neighbour_list_cutoff; QSet neighbour_list; - int max_num_mm=0; + int max_num_mm = 0; c10::DeviceType device; }; #endif @@ -326,11 +343,12 @@ namespace SireOpenMM */ TorchQMEngine( QString module_path, - SireUnits::Dimension::Length cutoff=7.5*SireUnits::angstrom, - int neighbour_list_frequency=0, - bool is_mechanical=false, - double lambda=1.0 - ); + SireUnits::Dimension::Length cutoff = 7.5 * SireUnits::angstrom, + int neighbour_list_frequency = 0, + bool is_mechanical = false, + double lambda = 1.0, + double switch_width = 0.2, + bool use_switch = true); //! Copy constructor. TorchQMEngine(const TorchQMEngine &other); @@ -481,6 +499,30 @@ namespace SireOpenMM */ void setCharges(QVector charges); + //! Get the switch width. + /*! \returns + The switch width as a fraction of the cutoff. + */ + double getSwitchWidth() const; + + //! Set the switch width. + /*! \param switch_width + The switch width as a fraction of the cutoff (0 to 1). + */ + void setSwitchWidth(double switch_width); + + //! Get whether a switching function is used. + /*! \returns + Whether a switching function is used. + */ + bool getUseSwitch() const; + + //! Set whether a switching function is used. + /*! \param use_switch + Whether to use a switching function. + */ + void setUseSwitch(bool use_switch); + //! Return the C++ name for this class. static const char *typeName(); @@ -488,7 +530,7 @@ namespace SireOpenMM const char *what() const; //! Create an EMLE force object. - QMForce* createForce() const; + QMForce *createForce() const; private: QString module_path; @@ -496,6 +538,8 @@ namespace SireOpenMM int neighbour_list_frequency; bool is_mechanical; double lambda; + double switch_width; + bool use_switch; QVector atoms; QMap mm1_to_qm; QMap> mm1_to_mm2; From 276e8cd78c6dced39a74c8a8e426ed20cb851e8d Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 May 2026 13:17:08 +0100 Subject: [PATCH 2/4] Add CUDA as emle system requirement. --- pixi.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pixi.toml b/pixi.toml index 19092abb9..f0b424f23 100644 --- a/pixi.toml +++ b/pixi.toml @@ -133,6 +133,9 @@ loguru = "*" pygit2 = "*" pyyaml = "*" +[feature.emle.system-requirements] +cuda = "12" + [feature.emle.target.linux-64.dependencies] ambertools = ">=22" deepmd-kit = "*" From 6eb4c815219ef9f3210c630497578de37075992f Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 May 2026 13:37:52 +0100 Subject: [PATCH 3/4] Fix docstrings. --- .../Convert/SireOpenMM/PyQMCallback.pypp.cpp | 89 +++-- .../Convert/SireOpenMM/PyQMEngine.pypp.cpp | 2 +- wrapper/Convert/SireOpenMM/PyQMForce.pypp.cpp | 304 +++++++----------- wrapper/Convert/SireOpenMM/pyqm.h | 4 + 4 files changed, 161 insertions(+), 238 deletions(-) diff --git a/wrapper/Convert/SireOpenMM/PyQMCallback.pypp.cpp b/wrapper/Convert/SireOpenMM/PyQMCallback.pypp.cpp index 885a15bda..f5fa83fd9 100644 --- a/wrapper/Convert/SireOpenMM/PyQMCallback.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PyQMCallback.pypp.cpp @@ -2,8 +2,8 @@ // (C) Christopher Woods, GPL >= 3 License -#include "boost/python.hpp" #include "PyQMCallback.pypp.hpp" +#include "boost/python.hpp" namespace bp = boost::python; @@ -51,69 +51,56 @@ namespace bp = boost::python; #include -SireOpenMM::PyQMCallback __copy__(const SireOpenMM::PyQMCallback &other){ return SireOpenMM::PyQMCallback(other); } +SireOpenMM::PyQMCallback __copy__(const SireOpenMM::PyQMCallback &other) { return SireOpenMM::PyQMCallback(other); } #include "Qt/qdatastream.hpp" -const char* pvt_get_name(const SireOpenMM::PyQMCallback&){ return "SireOpenMM::PyQMCallback";} +const char *pvt_get_name(const SireOpenMM::PyQMCallback &) { return "SireOpenMM::PyQMCallback"; } #include "Helpers/release_gil_policy.hpp" -void register_PyQMCallback_class(){ +void register_PyQMCallback_class() +{ { //::SireOpenMM::PyQMCallback - typedef bp::class_< SireOpenMM::PyQMCallback > PyQMCallback_exposer_t; - PyQMCallback_exposer_t PyQMCallback_exposer = PyQMCallback_exposer_t( "PyQMCallback", "A callback wrapper class to interface with external QM engines\nvia the CustomCPPForceImpl.", bp::init< >("Default constructor.") ); - bp::scope PyQMCallback_scope( PyQMCallback_exposer ); - PyQMCallback_exposer.def( bp::init< bp::api::object, bp::optional< QString > >(( bp::arg("arg0"), bp::arg("name")="" ), "Constructor\nPar:am py_object\nA Python object that contains the callback function.\n\nPar:am name\nThe name of a callback method that take the following arguments:\n- numbers_qm: A list of atomic numbers for the atoms in the ML region.\n- charges_mm: A list of the MM charges in mod electron charge.\n- xyz_qm: A list of positions for the atoms in the ML region in Angstrom.\n- xyz_mm: A list of positions for the atoms in the MM region in Angstrom.\n- cell: A list of cell vectors in Angstrom.\n- idx_mm: A list of indices for the MM atoms in the QM/MM region.\nThe callback should return a tuple containing:\n- The energy in kJmol.\n- A list of forces for the QM atoms in kJmolnm.\n- A list of forces for the MM atoms in kJmolnm.\nIf empty, then the object is assumed to be a callable.\n") ); + typedef bp::class_ PyQMCallback_exposer_t; + PyQMCallback_exposer_t PyQMCallback_exposer = PyQMCallback_exposer_t("PyQMCallback", "A callback wrapper class to interface with external QM engines\nvia the CustomCPPForceImpl.", bp::init<>("Default constructor.")); + bp::scope PyQMCallback_scope(PyQMCallback_exposer); + PyQMCallback_exposer.def(bp::init>((bp::arg("arg0"), bp::arg("name") = ""), "Constructor\nPar:am py_object\nA Python object that contains the callback function.\n\nPar:am name\nThe name of a callback method that take the following arguments:\n- numbers_qm: A list of atomic numbers for the atoms in the ML region.\n- charges_mm: A list of the MM charges in mod electron charge.\n- xyz_qm: A list of positions for the atoms in the ML region in Angstrom.\n- xyz_mm: A list of positions for the atoms in the MM region in Angstrom.\n- cell: A list of cell vectors in Angstrom.\n- idx_mm: A list of indices for the MM atoms in the QM/MM region.\nThe callback should return a tuple containing:\n- The energy in kJmol.\n- A list of forces for the QM atoms in kJmolnm.\n- A list of forces for the MM atoms in kJmolnm.\n- (Optional) The gradient of the energy w.r.t. the effective MM charges in kJmol/e, used for the chain-rule switching correction.\nIf empty, then the object is assumed to be a callable.\n")); { //::SireOpenMM::PyQMCallback::call - - typedef ::boost::tuples::tuple< double, QVector< QVector< double > >, QVector< QVector< double > >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::PyQMCallback::*call_function_type)( ::QVector< int >,::QVector< double >,::QVector< QVector< double > >,::QVector< QVector< double > >,::QVector >, ::QVector< int > ) const; - call_function_type call_function_value( &::SireOpenMM::PyQMCallback::call ); - - PyQMCallback_exposer.def( - "call" - , call_function_value - , ( bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm") ) - , bp::release_gil_policy() - , "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n" ); - + + typedef ::boost::tuples::tuple>, QVector>, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::PyQMCallback::*call_function_type)(::QVector, ::QVector, ::QVector>, ::QVector>, ::QVector>, ::QVector) const; + call_function_type call_function_value(&::SireOpenMM::PyQMCallback::call); + + PyQMCallback_exposer.def( + "call", call_function_value, (bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm")), bp::release_gil_policy(), "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n- (Optional) The gradient of the energy w.r.t. the effective MM charges in kJmol/e, used for the chain-rule switching correction.\n"); } { //::SireOpenMM::PyQMCallback::typeName - - typedef char const * ( *typeName_function_type )( ); - typeName_function_type typeName_function_value( &::SireOpenMM::PyQMCallback::typeName ); - - PyQMCallback_exposer.def( - "typeName" - , typeName_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(*typeName_function_type)(); + typeName_function_type typeName_function_value(&::SireOpenMM::PyQMCallback::typeName); + + PyQMCallback_exposer.def( + "typeName", typeName_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } { //::SireOpenMM::PyQMCallback::what - - typedef char const * ( ::SireOpenMM::PyQMCallback::*what_function_type)( ) const; - what_function_type what_function_value( &::SireOpenMM::PyQMCallback::what ); - - PyQMCallback_exposer.def( - "what" - , what_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(::SireOpenMM::PyQMCallback::*what_function_type)() const; + what_function_type what_function_value(&::SireOpenMM::PyQMCallback::what); + + PyQMCallback_exposer.def( + "what", what_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } - PyQMCallback_exposer.staticmethod( "typeName" ); - PyQMCallback_exposer.def( "__copy__", &__copy__); - PyQMCallback_exposer.def( "__deepcopy__", &__copy__); - PyQMCallback_exposer.def( "clone", &__copy__); - PyQMCallback_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireOpenMM::PyQMCallback >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); - PyQMCallback_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireOpenMM::PyQMCallback >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); - PyQMCallback_exposer.def_pickle(sire_pickle_suite< ::SireOpenMM::PyQMCallback >()); - PyQMCallback_exposer.def( "__str__", &pvt_get_name); - PyQMCallback_exposer.def( "__repr__", &pvt_get_name); + PyQMCallback_exposer.staticmethod("typeName"); + PyQMCallback_exposer.def("__copy__", &__copy__); + PyQMCallback_exposer.def("__deepcopy__", &__copy__); + PyQMCallback_exposer.def("clone", &__copy__); + PyQMCallback_exposer.def("__rlshift__", &__rlshift__QDataStream<::SireOpenMM::PyQMCallback>, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>()); + PyQMCallback_exposer.def("__rrshift__", &__rrshift__QDataStream<::SireOpenMM::PyQMCallback>, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>()); + PyQMCallback_exposer.def_pickle(sire_pickle_suite<::SireOpenMM::PyQMCallback>()); + PyQMCallback_exposer.def("__str__", &pvt_get_name); + PyQMCallback_exposer.def("__repr__", &pvt_get_name); } - } diff --git a/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp b/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp index d86fcb7f1..37308fa66 100644 --- a/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PyQMEngine.pypp.cpp @@ -72,7 +72,7 @@ void register_PyQMEngine_class() call_function_type call_function_value(&::SireOpenMM::PyQMEngine::call); PyQMEngine_exposer.def( - "call", call_function_value, (bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm")), bp::release_gil_policy(), "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of the true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n"); + "call", call_function_value, (bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm")), bp::release_gil_policy(), "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of the true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n- (Optional) The gradient of the energy w.r.t. the effective MM charges in kJmol/e, used for the chain-rule switching correction.\n"); } { //::SireOpenMM::PyQMEngine::getAtoms diff --git a/wrapper/Convert/SireOpenMM/PyQMForce.pypp.cpp b/wrapper/Convert/SireOpenMM/PyQMForce.pypp.cpp index b8b8bebd3..31a6c82d2 100644 --- a/wrapper/Convert/SireOpenMM/PyQMForce.pypp.cpp +++ b/wrapper/Convert/SireOpenMM/PyQMForce.pypp.cpp @@ -2,8 +2,8 @@ // (C) Christopher Woods, GPL >= 3 License -#include "boost/python.hpp" #include "PyQMForce.pypp.hpp" +#include "boost/python.hpp" namespace bp = boost::python; @@ -51,229 +51,161 @@ namespace bp = boost::python; #include -SireOpenMM::PyQMForce __copy__(const SireOpenMM::PyQMForce &other){ return SireOpenMM::PyQMForce(other); } +SireOpenMM::PyQMForce __copy__(const SireOpenMM::PyQMForce &other) { return SireOpenMM::PyQMForce(other); } #include "Qt/qdatastream.hpp" -const char* pvt_get_name(const SireOpenMM::PyQMForce&){ return "SireOpenMM::PyQMForce";} +const char *pvt_get_name(const SireOpenMM::PyQMForce &) { return "SireOpenMM::PyQMForce"; } #include "Helpers/release_gil_policy.hpp" -void register_PyQMForce_class(){ +void register_PyQMForce_class() +{ { //::SireOpenMM::PyQMForce - typedef bp::class_< SireOpenMM::PyQMForce, bp::bases< SireOpenMM::QMForce > > PyQMForce_exposer_t; - PyQMForce_exposer_t PyQMForce_exposer = PyQMForce_exposer_t( "PyQMForce", "", bp::init< >("Default constructor.") ); - bp::scope PyQMForce_scope( PyQMForce_exposer ); - PyQMForce_exposer.def( bp::init< SireOpenMM::PyQMCallback, SireUnits::Dimension::Length, int, bool, double, QVector< int >, QMap< int, int >, QMap< int, QVector< int > >, QMap< int, double >, QVector< int >, QVector< int >, QVector< double > >(( bp::arg("callback"), bp::arg("cutoff"), bp::arg("neighbour_list_frequency"), bp::arg("is_mechanical"), bp::arg("lambda"), bp::arg("atoms"), bp::arg("mm1_to_qm"), bp::arg("mm1_to_mm2"), bp::arg("bond_scale_factors"), bp::arg("mm2_atoms"), bp::arg("numbers"), bp::arg("charges") ), "Constructor.\nPar:am callback\nThe PyQMCallback object.\n\nPar:am cutoff\nThe ML cutoff distance.\n\nPar:am neighbour_list_frequency\nThe frequency at which the neighbour list is updated. (Number of steps.)\nIf zero, then no neighbour list is used.\n\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n\nPar:am lambda\nThe lambda weighting factor. This can be used to interpolate between\npotentials for end-state correction calculations.\n\nPar:am atoms\nA vector of atom indices for the QM region.\n\nPar:am mm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nPar:am mm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nPar:am bond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\nPar:am mm2_atoms\nA vector of MM2 atom indices.\n\nPar:am numbers\nA vector of atomic charges for all atoms in the system.\n\nPar:am charges\nA vector of atomic charges for all atoms in the system.\n") ); - PyQMForce_exposer.def( bp::init< SireOpenMM::PyQMForce const & >(( bp::arg("other") ), "Copy constructor.") ); + typedef bp::class_> PyQMForce_exposer_t; + PyQMForce_exposer_t PyQMForce_exposer = PyQMForce_exposer_t("PyQMForce", "", bp::init<>("Default constructor.")); + bp::scope PyQMForce_scope(PyQMForce_exposer); + PyQMForce_exposer.def(bp::init, QMap, QMap>, QMap, QVector, QVector, QVector>((bp::arg("callback"), bp::arg("cutoff"), bp::arg("neighbour_list_frequency"), bp::arg("is_mechanical"), bp::arg("lambda"), bp::arg("atoms"), bp::arg("mm1_to_qm"), bp::arg("mm1_to_mm2"), bp::arg("bond_scale_factors"), bp::arg("mm2_atoms"), bp::arg("numbers"), bp::arg("charges")), "Constructor.\nPar:am callback\nThe PyQMCallback object.\n\nPar:am cutoff\nThe ML cutoff distance.\n\nPar:am neighbour_list_frequency\nThe frequency at which the neighbour list is updated. (Number of steps.)\nIf zero, then no neighbour list is used.\n\nPar:am is_mechanical\nA flag to indicate if mechanical embedding is being used.\n\nPar:am lambda\nThe lambda weighting factor. This can be used to interpolate between\npotentials for end-state correction calculations.\n\nPar:am atoms\nA vector of atom indices for the QM region.\n\nPar:am mm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nPar:am mm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nPar:am bond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\nPar:am mm2_atoms\nA vector of MM2 atom indices.\n\nPar:am numbers\nA vector of atomic charges for all atoms in the system.\n\nPar:am charges\nA vector of atomic charges for all atoms in the system.\n")); + PyQMForce_exposer.def(bp::init((bp::arg("other")), "Copy constructor.")); { //::SireOpenMM::PyQMForce::call - - typedef ::boost::tuples::tuple< double, QVector< QVector< double > >, QVector< QVector< double > >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::PyQMForce::*call_function_type)( ::QVector< int >,::QVector< double >,::QVector< QVector< double > >,::QVector< QVector< double > >, ::QVector< QVector< double > >, ::QVector < int > ) const; - call_function_type call_function_value( &::SireOpenMM::PyQMForce::call ); - - PyQMForce_exposer.def( - "call" - , call_function_value - , ( bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm") ) - , bp::release_gil_policy() - , "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n" ); - + + typedef ::boost::tuples::tuple>, QVector>, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::PyQMForce::*call_function_type)(::QVector, ::QVector, ::QVector>, ::QVector>, ::QVector>, ::QVector) const; + call_function_type call_function_value(&::SireOpenMM::PyQMForce::call); + + PyQMForce_exposer.def( + "call", call_function_value, (bp::arg("numbers_qm"), bp::arg("charges_mm"), bp::arg("xyz_qm"), bp::arg("xyz_mm"), bp::arg("cell"), bp::arg("idx_mm")), bp::release_gil_policy(), "Call the callback function.\nPar:am numbers_qm\nA vector of atomic numbers for the atoms in the ML region.\n\nPar:am charges_mm\nA vector of the charges on the MM atoms in mod electron charge.\n\nPar:am xyz_qm\nA vector of positions for the atoms in the ML region in Angstrom.\n\nPar:am xyz_mm\nA vector of positions for the atoms in the MM region in Angstrom.\n\nPar:am cell A list of cell vectors in Angstrom.\n\nPar:am idx_mm A vector of indices for the MM atoms in the QM/MM region. Note that len(idx_mm) <= len(charges_mm) since it only contains the indices of true MM atoms, not link atoms or virtual charges.\n\nReturn:s\nA tuple containing:\n- The energy in kJmol.\n- A vector of forces for the QM atoms in kJmolnm.\n- A vector of forces for the MM atoms in kJmolnm.\n- (Optional) The gradient of the energy w.r.t. the effective MM charges in kJmol/e, used for the chain-rule switching correction.\n"); } { //::SireOpenMM::PyQMForce::getAtoms - - typedef ::QVector< int > ( ::SireOpenMM::PyQMForce::*getAtoms_function_type)( ) const; - getAtoms_function_type getAtoms_function_value( &::SireOpenMM::PyQMForce::getAtoms ); - - PyQMForce_exposer.def( - "getAtoms" - , getAtoms_function_value - , bp::release_gil_policy() - , "Get the indices of the atoms in the QM region.\nReturn:s\nA vector of atom indices for the QM region.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMForce::*getAtoms_function_type)() const; + getAtoms_function_type getAtoms_function_value(&::SireOpenMM::PyQMForce::getAtoms); + + PyQMForce_exposer.def( + "getAtoms", getAtoms_function_value, bp::release_gil_policy(), "Get the indices of the atoms in the QM region.\nReturn:s\nA vector of atom indices for the QM region.\n"); } { //::SireOpenMM::PyQMForce::getCallback - - typedef ::SireOpenMM::PyQMCallback ( ::SireOpenMM::PyQMForce::*getCallback_function_type)( ) const; - getCallback_function_type getCallback_function_value( &::SireOpenMM::PyQMForce::getCallback ); - - PyQMForce_exposer.def( - "getCallback" - , getCallback_function_value - , bp::release_gil_policy() - , "Get the callback object.\nReturn:s\nA Python object that contains the callback function.\n" ); - + + typedef ::SireOpenMM::PyQMCallback (::SireOpenMM::PyQMForce::*getCallback_function_type)() const; + getCallback_function_type getCallback_function_value(&::SireOpenMM::PyQMForce::getCallback); + + PyQMForce_exposer.def( + "getCallback", getCallback_function_value, bp::release_gil_policy(), "Get the callback object.\nReturn:s\nA Python object that contains the callback function.\n"); } { //::SireOpenMM::PyQMForce::getCharges - - typedef ::QVector< double > ( ::SireOpenMM::PyQMForce::*getCharges_function_type)( ) const; - getCharges_function_type getCharges_function_value( &::SireOpenMM::PyQMForce::getCharges ); - - PyQMForce_exposer.def( - "getCharges" - , getCharges_function_value - , bp::release_gil_policy() - , "Get the atomic charges of all atoms in the system.\nReturn:s\nA vector of atomic charges for all atoms in the system.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMForce::*getCharges_function_type)() const; + getCharges_function_type getCharges_function_value(&::SireOpenMM::PyQMForce::getCharges); + + PyQMForce_exposer.def( + "getCharges", getCharges_function_value, bp::release_gil_policy(), "Get the atomic charges of all atoms in the system.\nReturn:s\nA vector of atomic charges for all atoms in the system.\n"); } { //::SireOpenMM::PyQMForce::getCutoff - - typedef ::SireUnits::Dimension::Length ( ::SireOpenMM::PyQMForce::*getCutoff_function_type)( ) const; - getCutoff_function_type getCutoff_function_value( &::SireOpenMM::PyQMForce::getCutoff ); - - PyQMForce_exposer.def( - "getCutoff" - , getCutoff_function_value - , bp::release_gil_policy() - , "Get the QM cutoff distance.\nReturn:s\nThe QM cutoff distance.\n" ); - + + typedef ::SireUnits::Dimension::Length (::SireOpenMM::PyQMForce::*getCutoff_function_type)() const; + getCutoff_function_type getCutoff_function_value(&::SireOpenMM::PyQMForce::getCutoff); + + PyQMForce_exposer.def( + "getCutoff", getCutoff_function_value, bp::release_gil_policy(), "Get the QM cutoff distance.\nReturn:s\nThe QM cutoff distance.\n"); } { //::SireOpenMM::PyQMForce::getIsMechanical - - typedef bool ( ::SireOpenMM::PyQMForce::*getIsMechanical_function_type)( ) const; - getIsMechanical_function_type getIsMechanical_function_value( &::SireOpenMM::PyQMForce::getIsMechanical ); - - PyQMForce_exposer.def( - "getIsMechanical" - , getIsMechanical_function_value - , bp::release_gil_policy() - , "Get the mechanical embedding flag.\nReturn:s\nA flag to indicate if mechanical embedding is being used.\n" ); - + + typedef bool (::SireOpenMM::PyQMForce::*getIsMechanical_function_type)() const; + getIsMechanical_function_type getIsMechanical_function_value(&::SireOpenMM::PyQMForce::getIsMechanical); + + PyQMForce_exposer.def( + "getIsMechanical", getIsMechanical_function_value, bp::release_gil_policy(), "Get the mechanical embedding flag.\nReturn:s\nA flag to indicate if mechanical embedding is being used.\n"); } { //::SireOpenMM::PyQMForce::getLambda - - typedef double ( ::SireOpenMM::PyQMForce::*getLambda_function_type)( ) const; - getLambda_function_type getLambda_function_value( &::SireOpenMM::PyQMForce::getLambda ); - - PyQMForce_exposer.def( - "getLambda" - , getLambda_function_value - , bp::release_gil_policy() - , "Get the lambda weighting factor.\nReturn:s\nThe lambda weighting factor.\n" ); - + + typedef double (::SireOpenMM::PyQMForce::*getLambda_function_type)() const; + getLambda_function_type getLambda_function_value(&::SireOpenMM::PyQMForce::getLambda); + + PyQMForce_exposer.def( + "getLambda", getLambda_function_value, bp::release_gil_policy(), "Get the lambda weighting factor.\nReturn:s\nThe lambda weighting factor.\n"); } { //::SireOpenMM::PyQMForce::getLinkAtoms - - typedef ::boost::tuples::tuple< QMap< int, int >, QMap< int, QVector< int > >, QMap< int, double >, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type > ( ::SireOpenMM::PyQMForce::*getLinkAtoms_function_type)( ) const; - getLinkAtoms_function_type getLinkAtoms_function_value( &::SireOpenMM::PyQMForce::getLinkAtoms ); - - PyQMForce_exposer.def( - "getLinkAtoms" - , getLinkAtoms_function_value - , bp::release_gil_policy() - , "Get the link atoms associated with each QM atom.\nReturn:s\nA tuple containing:\n\nmm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nmm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nbond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n" ); - + + typedef ::boost::tuples::tuple, QMap>, QMap, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type, boost::tuples::null_type> (::SireOpenMM::PyQMForce::*getLinkAtoms_function_type)() const; + getLinkAtoms_function_type getLinkAtoms_function_value(&::SireOpenMM::PyQMForce::getLinkAtoms); + + PyQMForce_exposer.def( + "getLinkAtoms", getLinkAtoms_function_value, bp::release_gil_policy(), "Get the link atoms associated with each QM atom.\nReturn:s\nA tuple containing:\n\nmm1_to_qm\nA dictionary mapping link atom (MM1) indices to the QM atoms to\nwhich they are bonded.\n\nmm1_to_mm2\nA dictionary of link atoms indices (MM1) to a list of the MM\natoms to which they are bonded (MM2).\n\nbond_scale_factors\nA dictionary of link atom indices (MM1) to a list of the bond\nlength scale factors between the QM and MM1 atoms. The scale\nfactors are the ratio of the equilibrium bond lengths for the\nQM-L (QM-link) atom and QM-MM1 atom, i.e. R0(QM-L) R0(QM-MM1),\ntaken from the MM force field parameters for the molecule.\n\n"); } { //::SireOpenMM::PyQMForce::getMM2Atoms - - typedef ::QVector< int > ( ::SireOpenMM::PyQMForce::*getMM2Atoms_function_type)( ) const; - getMM2Atoms_function_type getMM2Atoms_function_value( &::SireOpenMM::PyQMForce::getMM2Atoms ); - - PyQMForce_exposer.def( - "getMM2Atoms" - , getMM2Atoms_function_value - , bp::release_gil_policy() - , "Get the vector of MM2 atoms.\nReturn:s\nA vector of MM2 atom indices.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMForce::*getMM2Atoms_function_type)() const; + getMM2Atoms_function_type getMM2Atoms_function_value(&::SireOpenMM::PyQMForce::getMM2Atoms); + + PyQMForce_exposer.def( + "getMM2Atoms", getMM2Atoms_function_value, bp::release_gil_policy(), "Get the vector of MM2 atoms.\nReturn:s\nA vector of MM2 atom indices.\n"); } { //::SireOpenMM::PyQMForce::getNeighbourListFrequency - - typedef int ( ::SireOpenMM::PyQMForce::*getNeighbourListFrequency_function_type)( ) const; - getNeighbourListFrequency_function_type getNeighbourListFrequency_function_value( &::SireOpenMM::PyQMForce::getNeighbourListFrequency ); - - PyQMForce_exposer.def( - "getNeighbourListFrequency" - , getNeighbourListFrequency_function_value - , bp::release_gil_policy() - , "Get the neighbour list frequency.\nReturn:s\nThe neighbour list frequency.\n" ); - + + typedef int (::SireOpenMM::PyQMForce::*getNeighbourListFrequency_function_type)() const; + getNeighbourListFrequency_function_type getNeighbourListFrequency_function_value(&::SireOpenMM::PyQMForce::getNeighbourListFrequency); + + PyQMForce_exposer.def( + "getNeighbourListFrequency", getNeighbourListFrequency_function_value, bp::release_gil_policy(), "Get the neighbour list frequency.\nReturn:s\nThe neighbour list frequency.\n"); } { //::SireOpenMM::PyQMForce::getNumbers - - typedef ::QVector< int > ( ::SireOpenMM::PyQMForce::*getNumbers_function_type)( ) const; - getNumbers_function_type getNumbers_function_value( &::SireOpenMM::PyQMForce::getNumbers ); - - PyQMForce_exposer.def( - "getNumbers" - , getNumbers_function_value - , bp::release_gil_policy() - , "Get the atomic numbers for the atoms in the QM region.\nReturn:s\nA vector of atomic numbers for the atoms in the QM region.\n" ); - + + typedef ::QVector (::SireOpenMM::PyQMForce::*getNumbers_function_type)() const; + getNumbers_function_type getNumbers_function_value(&::SireOpenMM::PyQMForce::getNumbers); + + PyQMForce_exposer.def( + "getNumbers", getNumbers_function_value, bp::release_gil_policy(), "Get the atomic numbers for the atoms in the QM region.\nReturn:s\nA vector of atomic numbers for the atoms in the QM region.\n"); } { //::SireOpenMM::PyQMForce::operator= - - typedef ::SireOpenMM::PyQMForce & ( ::SireOpenMM::PyQMForce::*assign_function_type)( ::SireOpenMM::PyQMForce const & ) ; - assign_function_type assign_function_value( &::SireOpenMM::PyQMForce::operator= ); - - PyQMForce_exposer.def( - "assign" - , assign_function_value - , ( bp::arg("other") ) - , bp::return_self< >() - , "Assignment operator." ); - + + typedef ::SireOpenMM::PyQMForce &(::SireOpenMM::PyQMForce::*assign_function_type)(::SireOpenMM::PyQMForce const &); + assign_function_type assign_function_value(&::SireOpenMM::PyQMForce::operator=); + + PyQMForce_exposer.def( + "assign", assign_function_value, (bp::arg("other")), bp::return_self<>(), "Assignment operator."); } { //::SireOpenMM::PyQMForce::setCallback - - typedef void ( ::SireOpenMM::PyQMForce::*setCallback_function_type)( ::SireOpenMM::PyQMCallback ) ; - setCallback_function_type setCallback_function_value( &::SireOpenMM::PyQMForce::setCallback ); - - PyQMForce_exposer.def( - "setCallback" - , setCallback_function_value - , ( bp::arg("callback") ) - , bp::release_gil_policy() - , "Set the callback object.\nPar:am callback\nA Python object that contains the callback function.\n" ); - + + typedef void (::SireOpenMM::PyQMForce::*setCallback_function_type)(::SireOpenMM::PyQMCallback); + setCallback_function_type setCallback_function_value(&::SireOpenMM::PyQMForce::setCallback); + + PyQMForce_exposer.def( + "setCallback", setCallback_function_value, (bp::arg("callback")), bp::release_gil_policy(), "Set the callback object.\nPar:am callback\nA Python object that contains the callback function.\n"); } { //::SireOpenMM::PyQMForce::setLambda - - typedef void ( ::SireOpenMM::PyQMForce::*setLambda_function_type)( double ) ; - setLambda_function_type setLambda_function_value( &::SireOpenMM::PyQMForce::setLambda ); - - PyQMForce_exposer.def( - "setLambda" - , setLambda_function_value - , ( bp::arg("lambda") ) - , bp::release_gil_policy() - , "Set the lambda weighting factor\nPar:am lambda\nThe lambda weighting factor.\n" ); - + + typedef void (::SireOpenMM::PyQMForce::*setLambda_function_type)(double); + setLambda_function_type setLambda_function_value(&::SireOpenMM::PyQMForce::setLambda); + + PyQMForce_exposer.def( + "setLambda", setLambda_function_value, (bp::arg("lambda")), bp::release_gil_policy(), "Set the lambda weighting factor\nPar:am lambda\nThe lambda weighting factor.\n"); } { //::SireOpenMM::PyQMForce::typeName - - typedef char const * ( *typeName_function_type )( ); - typeName_function_type typeName_function_value( &::SireOpenMM::PyQMForce::typeName ); - - PyQMForce_exposer.def( - "typeName" - , typeName_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - + + typedef char const *(*typeName_function_type)(); + typeName_function_type typeName_function_value(&::SireOpenMM::PyQMForce::typeName); + + PyQMForce_exposer.def( + "typeName", typeName_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); } { //::SireOpenMM::PyQMForce::what - - typedef char const * ( ::SireOpenMM::PyQMForce::*what_function_type)( ) const; - what_function_type what_function_value( &::SireOpenMM::PyQMForce::what ); - - PyQMForce_exposer.def( - "what" - , what_function_value - , bp::release_gil_policy() - , "Return the C++ name for this class." ); - - } - PyQMForce_exposer.staticmethod( "typeName" ); - PyQMForce_exposer.def( "__copy__", &__copy__); - PyQMForce_exposer.def( "__deepcopy__", &__copy__); - PyQMForce_exposer.def( "clone", &__copy__); - PyQMForce_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireOpenMM::PyQMForce >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); - PyQMForce_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireOpenMM::PyQMForce >, - bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() ); - PyQMForce_exposer.def_pickle(sire_pickle_suite< ::SireOpenMM::PyQMForce >()); - PyQMForce_exposer.def( "__str__", &pvt_get_name); - PyQMForce_exposer.def( "__repr__", &pvt_get_name); - } + typedef char const *(::SireOpenMM::PyQMForce::*what_function_type)() const; + what_function_type what_function_value(&::SireOpenMM::PyQMForce::what); + + PyQMForce_exposer.def( + "what", what_function_value, bp::release_gil_policy(), "Return the C++ name for this class."); + } + PyQMForce_exposer.staticmethod("typeName"); + PyQMForce_exposer.def("__copy__", &__copy__); + PyQMForce_exposer.def("__deepcopy__", &__copy__); + PyQMForce_exposer.def("clone", &__copy__); + PyQMForce_exposer.def("__rlshift__", &__rlshift__QDataStream<::SireOpenMM::PyQMForce>, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>()); + PyQMForce_exposer.def("__rrshift__", &__rrshift__QDataStream<::SireOpenMM::PyQMForce>, + bp::return_internal_reference<1, bp::with_custodian_and_ward<1, 2>>()); + PyQMForce_exposer.def_pickle(sire_pickle_suite<::SireOpenMM::PyQMForce>()); + PyQMForce_exposer.def("__str__", &pvt_get_name); + PyQMForce_exposer.def("__repr__", &pvt_get_name); + } } diff --git a/wrapper/Convert/SireOpenMM/pyqm.h b/wrapper/Convert/SireOpenMM/pyqm.h index d26777885..2acb32b34 100644 --- a/wrapper/Convert/SireOpenMM/pyqm.h +++ b/wrapper/Convert/SireOpenMM/pyqm.h @@ -124,6 +124,8 @@ namespace SireOpenMM - The energy in kJ/mol. - A vector of forces for the QM atoms in kJ/mol/nm. - A vector of forces for the MM atoms in kJ/mol/nm. + - (Optional) The gradient of the energy w.r.t. the effective MM + charges in kJ/mol/e, used for the chain-rule switching correction. */ boost::tuple>, QVector>> call( QVector numbers_qm, @@ -652,6 +654,8 @@ namespace SireOpenMM - The energy in kJ/mol. - A vector of forces for the QM atoms in kJ/mol/nm. - A vector of forces for the MM atoms in kJ/mol/nm. + - (Optional) The gradient of the energy w.r.t. the effective MM + charges in kJ/mol/e, used for the chain-rule switching correction. */ boost::tuple>, QVector>> call( QVector numbers_qm, From 5af3b0bc65a8c3c7ecba5d44e33470e3de4c07af Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 May 2026 14:24:30 +0100 Subject: [PATCH 4/4] Document new switch_width option. --- doc/source/tutorial/part08/01_intro.rst | 19 ++++++++++++++----- doc/source/tutorial/part08/02_emle.rst | 18 ++++++++++++------ 2 files changed, 26 insertions(+), 11 deletions(-) diff --git a/doc/source/tutorial/part08/01_intro.rst b/doc/source/tutorial/part08/01_intro.rst index f68beef45..f9559715c 100644 --- a/doc/source/tutorial/part08/01_intro.rst +++ b/doc/source/tutorial/part08/01_intro.rst @@ -37,6 +37,7 @@ an engine to perform the calculation: ... cutoff="7.5A", ... neighbour_list_frequency=20, ... mechanical_embedding=False, +... switch_width=0.2, ... ) Here the first argument is the molecules that we are simulating, the second @@ -65,7 +66,7 @@ signature: xyz_mm: List[List[float]], cell: Optional[List[List[float]]] = None, idx_mm: Optional[List[int]] = None, - ) -> Tuple[float, List[List[float]], List[List[float]]]: + ) -> Tuple[float, List[List[float]], List[List[float]], Optional[List[float]]]: The function takes the atomic numbers of the QM atoms, the charges of the MM atoms in mod electron charge, the coordinates of the QM atoms in Angstrom, and @@ -75,10 +76,18 @@ QM/MM region. This is useful for obtaining any additional atomic properties that may be required by the callback. (Note that link atoms and virtual charges are always placed last in the list of MM charges and positions.) The function should return the calculated energy in kJ/mol, the forces on the QM atoms in -kJ/mol/nm, and the forces on the MM atoms in kJ/mol/nm. The remaining arguments -are optional and specify the QM cutoff distance, the neighbour list update -frequency, and whether the electrostatics should be treated with mechanical -embedding. When mechanical embedding is used, the electrostatics are treated +kJ/mol/nm, and the forces on the MM atoms in kJ/mol/nm. Optionally, it may also +return a fourth element: the gradient of the energy with respect to the effective +MM charges (dE/dq) in kJ/mol/e. When present, this is used to apply a chain-rule +force correction that accounts for the positional dependence of the charge switching +function (see below), which is important for systems with a charged ML region. The +remaining arguments are optional and specify the QM cutoff distance, the neighbour +list update frequency, whether the electrostatics should be treated with mechanical +embedding, and the width of the switching region as a fraction of the cutoff +(``switch_width``, default 0.2). The switching function smoothly scales the MM +charges to zero between ``(1 - switch_width) * cutoff`` and ``cutoff``, which +avoids force discontinuities as MM atoms enter or leave the cutoff sphere. Setting +``switch_width=0`` disables switching (not recommended for charged ML regions). When mechanical embedding is used, the electrostatics are treated at the MM level by ``OpenMM``. Note that this doesn't change the signature of the callback function, i.e. it will be passed empty lists for the MM specific arguments and should return an empty list for the MM forces. Atomic positions diff --git a/doc/source/tutorial/part08/02_emle.rst b/doc/source/tutorial/part08/02_emle.rst index b171d6292..d75662225 100644 --- a/doc/source/tutorial/part08/02_emle.rst +++ b/doc/source/tutorial/part08/02_emle.rst @@ -55,14 +55,20 @@ an engine to perform the calculation: ... mols[0], ... calculator, ... cutoff="7.5A", -... neighbour_list_frequency=20 +... neighbour_list_frequency=20, +... switch_width=0.2, ... ) Here the first argument is the molecules that we are simulating, the second selection coresponding to the QM region (here this is the first molecule), and -the third is calculator that was created above. The fourth and fifth arguments -are optional, and specify the QM cutoff distance and the neighbour list update -frequency respectively. (Shown are the default values.) The function returns a +the third is calculator that was created above. The remaining arguments are +optional and specify the QM cutoff distance, the neighbour list update frequency, +and the width of the switching region as a fraction of the cutoff (``switch_width``, +default 0.2). (Shown are the default values.) The switching function smoothly +scales the MM charges to zero between ``(1 - switch_width) * cutoff`` and +``cutoff``, which avoids force discontinuities as MM atoms enter or leave the +cutoff sphere. This is particularly important for systems with a charged ML region. +Setting ``switch_width=0`` disables switching. The function returns a modified version of the molecules containing a "merged" dipeptide that can be interpolated between MM and QM levels of theory, along with an engine. The engine registers a Python callback that uses ``emle-engine`` to perform the QM @@ -182,7 +188,7 @@ computes the electrostatic embedding: Next we create a new engine bound to the calculator: >>> _, engine = sr.qm.emle( ->>> ... mols, mols[0], calculator, cutoff="7.5A", neighbour_list_frequency=20 +>>> ... mols, mols[0], calculator, cutoff="7.5A", neighbour_list_frequency=20, switch_width=0.2 >>> ... ) .. note:: @@ -423,7 +429,7 @@ The model is serialisable, so can be saved and loaded using the standard It is also possible to use the model with Sire when performing QM/MM dynamics: >>> qm_mols, engine = sr.qm.emle( -... mols, mols[0], model, cutoff="7.5A", neighbour_list_frequency=20 +... mols, mols[0], model, cutoff="7.5A", neighbour_list_frequency=20, switch_width=0.2 ... ) The model will be serialised and loaded into a C++ ``TorchQMEngine`` object,