Merge pull request #198 from OpenBioSim/fix_ghost_sigmas

Add map option to prevent peturbation of the LJ sigma value for ghost atoms
OpenBioSim · Jun 2, 2024 · 74743b5 · 74743b5
2 parents 8cff51b + 002ada4
commit 74743b5
Show file tree

Hide file tree

Showing 199 changed files with 1,040 additions and 587 deletions.
diff --git a/corelib/src/libs/SireIO/grotop.cpp b/corelib/src/libs/SireIO/grotop.cpp
@@ -5399,6 +5399,12 @@ QVector<QString> GroTop::preprocess(const QVector<QString> &lines, QHash<QString
             }
         }
 
+        // skip BioSimSpace position restraint includes
+        if (line.contains("#include \"posre"))
+        {
+            continue;
+        }
+
         // now look for #include lines
         if (line.startsWith("#include"))
         {

diff --git a/corelib/src/libs/SireMM/ljparameter.h b/corelib/src/libs/SireMM/ljparameter.h
@@ -141,12 +141,14 @@ namespace SireMM
         return not operator==(other);
     }
 
-    /** Return whether or not this is a dummy LJ parameter */
+    /** Return whether or not this is a dummy LJ parameter. A dummy LJ
+     *  parameter is one with a zero epsilon value (i.e. it will have a
+     *  a zero LJ interaction energy with all other particles, regardless
+     *  of their LJ parameters)
+     */
     SIRE_ALWAYS_INLINE bool LJParameter::isDummy() const
     {
-        // we only need to compare sqrtsig as this will be set to zero if
-        // sqrteps is zero
-        return sqrtsig == 0.0;
+        return sqrteps == 0.0;
     }
 
     /** Return whether or not this parameter has non-zero LJ parameters */

diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst
@@ -42,6 +42,12 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
   standard trajectory save functions, e.g.
   ``sire.save(mols.trajectory(), "output", format=["PRMTOP", "RST"])``.
 
+* Ignore BioSimSpace format position restraint include directives when
+  parsing GROMACS topology files.
+
+* Add a map option to prevent perturbation of the Lennard-Jones sigma
+  parameter for ghost atoms during alchemical free energy simulations.
+
 * Please add an item to this changelog when you create your PR
 
 `2024.1.0 <https://github.com/openbiosim/sire/compare/2023.5.2...2024.1.0>`__ - April 2024

diff --git a/tests/io/test_grotop.py b/tests/io/test_grotop.py
@@ -33,3 +33,9 @@ def test_grospace(tmpdir, kigaki_mols):
     mols = sr.load(f, show_warnings=False)
 
     assert energy.value() == pytest.approx(mols.energy().value())
+
+
+def test_posre():
+    # Make sure we can parse a file with BioSimSpace position restraint include
+    # directives.
+    mols = sr.load_test_files("posre.top")
diff --git a/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp b/wrapper/Convert/SireOpenMM/LambdaLever.pypp.cpp
@@ -29,6 +29,8 @@ namespace bp = boost::python;
 
 SireOpenMM::LambdaLever __copy__(const SireOpenMM::LambdaLever &other){ return SireOpenMM::LambdaLever(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Helpers/str.hpp"
 
 #include "Helpers/release_gil_policy.hpp"
@@ -55,14 +57,13 @@ void register_LambdaLever_class(){
         }
         { //::SireOpenMM::LambdaLever::addPerturbableMolecule
 
-            typedef int ( ::SireOpenMM::LambdaLever::*addPerturbableMolecule_function_type)( ::SireOpenMM::OpenMMMolecule const &,::QHash< QString, int > const & ) ;
+            typedef int ( ::SireOpenMM::LambdaLever::*addPerturbableMolecule_function_type)( ::SireOpenMM::OpenMMMolecule const &,::QHash< QString, int > const &,::SireBase::PropertyMap const & ) ;
             addPerturbableMolecule_function_type addPerturbableMolecule_function_value( &::SireOpenMM::LambdaLever::addPerturbableMolecule );
 
             LambdaLever_exposer.def( 
                 "addPerturbableMolecule"
                 , addPerturbableMolecule_function_value
-                , ( bp::arg("molecule"), bp::arg("start_indicies") )
-                , bp::release_gil_policy()
+                , ( bp::arg("molecule"), bp::arg("start_indicies"), bp::arg("map")=SireBase::PropertyMap() )
                 , "Add info for the passed perturbable OpenMMMolecule, returning\n  its index in the list of perturbable molecules\n" );
 
         }
@@ -115,7 +116,7 @@ void register_LambdaLever_class(){
                 , getLeverValues_function_value
                 , ( bp::arg("lambda_values"), bp::arg("mol") )
                 , bp::release_gil_policy()
-                , "" );
+                , "Get all of the lever values that would be set for the passed\n  lambda values using the current context. This returns a PropertyList\n  of columns, where each column is a PropertyMap with the column name\n  and either double or QString array property of values.\n\n  This is designed to be used by a higher-level python function that\n  will convert this output into, e.g. a pandas DataFrame\n" );
 
         }
         { //::SireOpenMM::LambdaLever::getPerturbableMoleculeMaps
@@ -272,9 +273,9 @@ void register_LambdaLever_class(){
 
         }
         LambdaLever_exposer.staticmethod( "typeName" );
-        LambdaLever_exposer.def( "__copy__", &__copy__);
-        LambdaLever_exposer.def( "__deepcopy__", &__copy__);
-        LambdaLever_exposer.def( "clone", &__copy__);
+        LambdaLever_exposer.def( "__copy__", &__copy__<SireOpenMM::LambdaLever>);
+        LambdaLever_exposer.def( "__deepcopy__", &__copy__<SireOpenMM::LambdaLever>);
+        LambdaLever_exposer.def( "clone", &__copy__<SireOpenMM::LambdaLever>);
         LambdaLever_exposer.def( "__str__", &__str__< ::SireOpenMM::LambdaLever > );
         LambdaLever_exposer.def( "__repr__", &__str__< ::SireOpenMM::LambdaLever > );
     }

diff --git a/wrapper/Convert/SireOpenMM/OpenMMMetaData.pypp.cpp b/wrapper/Convert/SireOpenMM/OpenMMMetaData.pypp.cpp
@@ -129,6 +129,8 @@ namespace bp = boost::python;
 
 SireOpenMM::OpenMMMetaData __copy__(const SireOpenMM::OpenMMMetaData &other){ return SireOpenMM::OpenMMMetaData(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Helpers/str.hpp"
 
 #include "Helpers/release_gil_policy.hpp"
@@ -273,9 +275,9 @@ void register_OpenMMMetaData_class(){
 
         }
         OpenMMMetaData_exposer.staticmethod( "typeName" );
-        OpenMMMetaData_exposer.def( "__copy__", &__copy__);
-        OpenMMMetaData_exposer.def( "__deepcopy__", &__copy__);
-        OpenMMMetaData_exposer.def( "clone", &__copy__);
+        OpenMMMetaData_exposer.def( "__copy__", &__copy__<SireOpenMM::OpenMMMetaData>);
+        OpenMMMetaData_exposer.def( "__deepcopy__", &__copy__<SireOpenMM::OpenMMMetaData>);
+        OpenMMMetaData_exposer.def( "clone", &__copy__<SireOpenMM::OpenMMMetaData>);
         OpenMMMetaData_exposer.def( "__str__", &__str__< ::SireOpenMM::OpenMMMetaData > );
         OpenMMMetaData_exposer.def( "__repr__", &__str__< ::SireOpenMM::OpenMMMetaData > );
     }

diff --git a/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp b/wrapper/Convert/SireOpenMM/PerturbableOpenMMMolecule.pypp.cpp
@@ -113,6 +113,8 @@ namespace bp = boost::python;
 
 SireOpenMM::PerturbableOpenMMMolecule __copy__(const SireOpenMM::PerturbableOpenMMMolecule &other){ return SireOpenMM::PerturbableOpenMMMolecule(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Helpers/str.hpp"
 
 #include "Helpers/release_gil_policy.hpp"
@@ -123,8 +125,8 @@ void register_PerturbableOpenMMMolecule_class(){
         typedef bp::class_< SireOpenMM::PerturbableOpenMMMolecule > PerturbableOpenMMMolecule_exposer_t;
         PerturbableOpenMMMolecule_exposer_t PerturbableOpenMMMolecule_exposer = PerturbableOpenMMMolecule_exposer_t( "PerturbableOpenMMMolecule", "This class holds all of the information of an OpenMM molecule\nthat can be perturbed using a LambdaSchedule. The data is held\nin easy-to-access arrays, with guarantees that the arrays are\ncompatible and the data is aligned.\n", bp::init< >("Null constructor") );
         bp::scope PerturbableOpenMMMolecule_scope( PerturbableOpenMMMolecule_exposer );
-        PerturbableOpenMMMolecule_exposer.def( bp::init< SireOpenMM::OpenMMMolecule const & >(( bp::arg("mol") ), "Construct from the passed OpenMMMolecule") );
-        PerturbableOpenMMMolecule_exposer.def( bp::init< SireMol::Molecule const &, SireBase::PropertyMap const & >(( bp::arg("mol"), bp::arg("map") ), "Construct from a passed molecule and map") );
+        PerturbableOpenMMMolecule_exposer.def( bp::init< SireOpenMM::OpenMMMolecule const &, bp::optional< SireBase::PropertyMap const & > >(( bp::arg("mol"), bp::arg("map")=SireBase::PropertyMap() ), "Construct from the passed OpenMMMolecule") );
+        PerturbableOpenMMMolecule_exposer.def( bp::init< SireMol::Molecule const &, bp::optional< SireBase::PropertyMap const & > >(( bp::arg("mol"), bp::arg("map")=SireBase::PropertyMap() ), "Construct from a passed molecule and map") );
         PerturbableOpenMMMolecule_exposer.def( bp::init< SireOpenMM::PerturbableOpenMMMolecule const & >(( bp::arg("other") ), "Copy constructor") );
         { //::SireOpenMM::PerturbableOpenMMMolecule::angles
 
@@ -595,6 +597,18 @@ void register_PerturbableOpenMMMolecule_class(){
                 , bp::release_gil_policy()
                 , "Return true if the atom is a ghost atom in the\n  referenece or perturbed states" );
 
+        }
+        { //::SireOpenMM::PerturbableOpenMMMolecule::isNull
+
+            typedef bool ( ::SireOpenMM::PerturbableOpenMMMolecule::*isNull_function_type)(  ) const;
+            isNull_function_type isNull_function_value( &::SireOpenMM::PerturbableOpenMMMolecule::isNull );
+
+            PerturbableOpenMMMolecule_exposer.def( 
+                "isNull"
+                , isNull_function_value
+                , bp::release_gil_policy()
+                , "Return whether or not this is null" );
+
         }
         PerturbableOpenMMMolecule_exposer.def( bp::self != bp::self );
         { //::SireOpenMM::PerturbableOpenMMMolecule::operator=
@@ -686,9 +700,9 @@ void register_PerturbableOpenMMMolecule_class(){
 
         }
         PerturbableOpenMMMolecule_exposer.staticmethod( "typeName" );
-        PerturbableOpenMMMolecule_exposer.def( "__copy__", &__copy__);
-        PerturbableOpenMMMolecule_exposer.def( "__deepcopy__", &__copy__);
-        PerturbableOpenMMMolecule_exposer.def( "clone", &__copy__);
+        PerturbableOpenMMMolecule_exposer.def( "__copy__", &__copy__<SireOpenMM::PerturbableOpenMMMolecule>);
+        PerturbableOpenMMMolecule_exposer.def( "__deepcopy__", &__copy__<SireOpenMM::PerturbableOpenMMMolecule>);
+        PerturbableOpenMMMolecule_exposer.def( "clone", &__copy__<SireOpenMM::PerturbableOpenMMMolecule>);
         PerturbableOpenMMMolecule_exposer.def( "__str__", &__str__< ::SireOpenMM::PerturbableOpenMMMolecule > );
         PerturbableOpenMMMolecule_exposer.def( "__repr__", &__str__< ::SireOpenMM::PerturbableOpenMMMolecule > );
     }

diff --git a/wrapper/Convert/SireOpenMM/lambdalever.cpp b/wrapper/Convert/SireOpenMM/lambdalever.cpp
@@ -1964,10 +1964,11 @@ QList<OpenMM::Force *> LambdaLever::getRestraints(const QString &name,
  *  its index in the list of perturbable molecules
  */
 int LambdaLever::addPerturbableMolecule(const OpenMMMolecule &molecule,
-                                        const QHash<QString, qint32> &starts)
+                                        const QHash<QString, qint32> &starts,
+                                        const PropertyMap &map)
 {
     // should add in some sanity checks for these inputs
-    this->perturbable_mols.append(PerturbableOpenMMMolecule(molecule));
+    this->perturbable_mols.append(PerturbableOpenMMMolecule(molecule, map));
     this->start_indicies.append(starts);
     this->perturbable_maps.insert(molecule.number, molecule.perturtable_map);
     this->lambda_cache.clear();

diff --git a/wrapper/Convert/SireOpenMM/lambdalever.h b/wrapper/Convert/SireOpenMM/lambdalever.h
@@ -123,7 +123,8 @@ namespace SireOpenMM
         void addRestraintIndex(const QString &force, int index);
 
         int addPerturbableMolecule(const OpenMMMolecule &molecule,
-                                   const QHash<QString, qint32> &start_indicies);
+                                   const QHash<QString, qint32> &start_indicies,
+                                   const SireBase::PropertyMap &map = SireBase::PropertyMap());
 
         void setExceptionIndicies(int idx, const QString &ff,
                                   const QVector<boost::tuple<int, int>> &exception_idxs);

diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.cpp b/wrapper/Convert/SireOpenMM/openmmmolecule.cpp
@@ -2073,7 +2073,7 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const Molecule &mol,
                                                      const PropertyMap &map)
     : ConcreteProperty<PerturbableOpenMMMolecule, Property>()
 {
-    this->operator=(PerturbableOpenMMMolecule(OpenMMMolecule(mol, map)));
+    this->operator=(PerturbableOpenMMMolecule(OpenMMMolecule(mol, map), map));
 }
 
 /** Return whether or not this is null */
@@ -2085,7 +2085,8 @@ bool PerturbableOpenMMMolecule::isNull() const
 }
 
 /** Construct from the passed OpenMMMolecule */
-PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol)
+PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol,
+                                                     const PropertyMap &map)
     : ConcreteProperty<PerturbableOpenMMMolecule, Property>()
 {
     if (mol.perturbed.get() == 0)
@@ -2175,6 +2176,45 @@ PerturbableOpenMMMolecule::PerturbableOpenMMMolecule(const OpenMMMolecule &mol)
     from_ghost_idxs = mol.from_ghost_idxs;
 
     perturbable_constraints = mol.perturbable_constraints;
+
+    bool fix_perturbable_zero_sigmas = false;
+
+    if (map.specified("fix_perturbable_zero_sigmas"))
+    {
+        fix_perturbable_zero_sigmas = map["fix_perturbable_zero_sigmas"].value().asABoolean();
+    }
+
+    if (fix_perturbable_zero_sigmas)
+    {
+        const int nats = sig0.count();
+
+        if (sig1.count() != nats)
+        {
+            throw SireError::program_bug(QObject::tr(
+                                             "The number of atoms in the reference (%1) and perturbed (%2) "
+                                             "states do not match.")
+                                             .arg(nats)
+                                             .arg(sig1.count()),
+                                         CODELOC);
+        }
+
+        const auto *sig0_data = sig0.constData();
+        const auto *sig1_data = sig1.constData();
+
+        for (int i = 0; i < nats; ++i)
+        {
+            if (std::abs(sig0_data[i]) < 1e-9)
+            {
+                sig0[i] = sig1_data[i];
+                sig0_data = sig0.constData();
+            }
+            else if (std::abs(sig1_data[i] < 1e-9))
+            {
+                sig1[i] = sig0_data[i];
+                sig1_data = sig1.constData();
+            }
+        }
+    }
 }
 
 /** Copy constructor */

diff --git a/wrapper/Convert/SireOpenMM/openmmmolecule.h b/wrapper/Convert/SireOpenMM/openmmmolecule.h
@@ -204,10 +204,11 @@ namespace SireOpenMM
     public:
         PerturbableOpenMMMolecule();
 
-        PerturbableOpenMMMolecule(const OpenMMMolecule &mol);
+        PerturbableOpenMMMolecule(const OpenMMMolecule &mol,
+                                  const SireBase::PropertyMap &map = SireBase::PropertyMap());
 
         PerturbableOpenMMMolecule(const SireMol::Molecule &mol,
-                                  const SireBase::PropertyMap &map);
+                                  const SireBase::PropertyMap &map = SireBase::PropertyMap());
 
         PerturbableOpenMMMolecule(const PerturbableOpenMMMolecule &other);
 

diff --git a/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp b/wrapper/Convert/SireOpenMM/sire_to_openmm_system.cpp
@@ -1144,7 +1144,8 @@ OpenMMMetaData SireOpenMM::sire_to_openmm_system(OpenMM::System &system,
             // of perturbable molecules (e.g. the first perturbable
             // molecule we find has index 0)
             auto pert_idx = lambda_lever.addPerturbableMolecule(mol,
-                                                                start_indicies);
+                                                                start_indicies,
+                                                                map);
 
             // and we can record the map from the molecule index
             // to the perturbable molecule index

diff --git a/wrapper/Convert/__init__.py b/wrapper/Convert/__init__.py
@@ -398,7 +398,7 @@ def sire_to_openmm(mols, map):
                     platform = platforms["cpu"]
 
                 elif len(platforms) > 0:
-                    platform = platforms[platforms.keys()[0]]
+                    platform = platforms[list(platforms.keys())[0]]
 
             if platform is None:
                 raise ValueError(

diff --git a/wrapper/MM/AmberAngle.pypp.cpp b/wrapper/MM/AmberAngle.pypp.cpp
@@ -73,6 +73,8 @@ namespace bp = boost::python;
 
 SireMM::AmberAngle __copy__(const SireMM::AmberAngle &other){ return SireMM::AmberAngle(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Qt/qdatastream.hpp"
 
 #include "Helpers/str.hpp"
@@ -217,9 +219,9 @@ void register_AmberAngle_class(){
 
         }
         AmberAngle_exposer.staticmethod( "typeName" );
-        AmberAngle_exposer.def( "__copy__", &__copy__);
-        AmberAngle_exposer.def( "__deepcopy__", &__copy__);
-        AmberAngle_exposer.def( "clone", &__copy__);
+        AmberAngle_exposer.def( "__copy__", &__copy__<SireMM::AmberAngle>);
+        AmberAngle_exposer.def( "__deepcopy__", &__copy__<SireMM::AmberAngle>);
+        AmberAngle_exposer.def( "clone", &__copy__<SireMM::AmberAngle>);
         AmberAngle_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AmberAngle >,
                             bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() );
         AmberAngle_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AmberAngle >,

diff --git a/wrapper/MM/AmberBond.pypp.cpp b/wrapper/MM/AmberBond.pypp.cpp
@@ -73,6 +73,8 @@ namespace bp = boost::python;
 
 SireMM::AmberBond __copy__(const SireMM::AmberBond &other){ return SireMM::AmberBond(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Qt/qdatastream.hpp"
 
 #include "Helpers/str.hpp"
@@ -217,9 +219,9 @@ void register_AmberBond_class(){
 
         }
         AmberBond_exposer.staticmethod( "typeName" );
-        AmberBond_exposer.def( "__copy__", &__copy__);
-        AmberBond_exposer.def( "__deepcopy__", &__copy__);
-        AmberBond_exposer.def( "clone", &__copy__);
+        AmberBond_exposer.def( "__copy__", &__copy__<SireMM::AmberBond>);
+        AmberBond_exposer.def( "__deepcopy__", &__copy__<SireMM::AmberBond>);
+        AmberBond_exposer.def( "clone", &__copy__<SireMM::AmberBond>);
         AmberBond_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AmberBond >,
                             bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() );
         AmberBond_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AmberBond >,

diff --git a/wrapper/MM/AmberDihPart.pypp.cpp b/wrapper/MM/AmberDihPart.pypp.cpp
@@ -73,6 +73,8 @@ namespace bp = boost::python;
 
 SireMM::AmberDihPart __copy__(const SireMM::AmberDihPart &other){ return SireMM::AmberDihPart(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Qt/qdatastream.hpp"
 
 #include "Helpers/str.hpp"
@@ -215,9 +217,9 @@ void register_AmberDihPart_class(){
 
         }
         AmberDihPart_exposer.staticmethod( "typeName" );
-        AmberDihPart_exposer.def( "__copy__", &__copy__);
-        AmberDihPart_exposer.def( "__deepcopy__", &__copy__);
-        AmberDihPart_exposer.def( "clone", &__copy__);
+        AmberDihPart_exposer.def( "__copy__", &__copy__<SireMM::AmberDihPart>);
+        AmberDihPart_exposer.def( "__deepcopy__", &__copy__<SireMM::AmberDihPart>);
+        AmberDihPart_exposer.def( "clone", &__copy__<SireMM::AmberDihPart>);
         AmberDihPart_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AmberDihPart >,
                             bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() );
         AmberDihPart_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AmberDihPart >,

diff --git a/wrapper/MM/AmberDihedral.pypp.cpp b/wrapper/MM/AmberDihedral.pypp.cpp
@@ -73,6 +73,8 @@ namespace bp = boost::python;
 
 SireMM::AmberDihedral __copy__(const SireMM::AmberDihedral &other){ return SireMM::AmberDihedral(other); }
 
+#include "Helpers/copy.hpp"
+
 #include "Qt/qdatastream.hpp"
 
 #include "Helpers/str.hpp"
@@ -203,9 +205,9 @@ void register_AmberDihedral_class(){
 
         }
         AmberDihedral_exposer.staticmethod( "typeName" );
-        AmberDihedral_exposer.def( "__copy__", &__copy__);
-        AmberDihedral_exposer.def( "__deepcopy__", &__copy__);
-        AmberDihedral_exposer.def( "clone", &__copy__);
+        AmberDihedral_exposer.def( "__copy__", &__copy__<SireMM::AmberDihedral>);
+        AmberDihedral_exposer.def( "__deepcopy__", &__copy__<SireMM::AmberDihedral>);
+        AmberDihedral_exposer.def( "clone", &__copy__<SireMM::AmberDihedral>);
         AmberDihedral_exposer.def( "__rlshift__", &__rlshift__QDataStream< ::SireMM::AmberDihedral >,
                             bp::return_internal_reference<1, bp::with_custodian_and_ward<1,2> >() );
         AmberDihedral_exposer.def( "__rrshift__", &__rrshift__QDataStream< ::SireMM::AmberDihedral >,