[WIP] Allow GROMACS water molecules to be flagged as crystal waters

corelib/src/libs/SireIO/biosimspace.cpp

@@ -271,7 +271,7 @@ namespace SireIO
         // Copy across all properties that are unique to the original molecule.
         for (const auto &prop : molecule.propertyKeys())
         {
-            if (not water.hasProperty(prop))
+            if (not water.hasProperty(prop) and not prop.compare("is_non_searchable_water"))
             {
                 edit_mol = edit_mol.setProperty(prop, molecule.property(prop)).molecule();
             }
@@ -485,15 +485,15 @@ namespace SireIO
     }
 
     Molecule _pvt_setGromacsWater(Molecule &molecule, const Molecule &water, const QString &model, bool has_virtual,
-                                  const PropertyMap &map)
+                                  const PropertyMap &map, bool is_crystal)
     {
         // Make the template water molecule editable and renumber it.
         auto edit_mol = water.edit().renumber(molecule.number());
 
         // Copy across all properties that are unique to the original molecule.
         for (const auto &prop : molecule.propertyKeys())
         {
-            if (not water.hasProperty(prop))
+            if (not water.hasProperty(prop) and not prop.compare("is_non_searchable_water"))
             {
                 edit_mol = edit_mol.setProperty(prop, molecule.property(prop)).molecule();
             }
@@ -539,6 +539,13 @@ namespace SireIO
                 coord_oxygen = molecule.atom(idx).property<Vector>(map["coordinates"]);
         }
 
+        // Rename the molecule and residue to XTL if this is a crystal water.
+        if (is_crystal)
+        {
+            edit_mol = edit_mol.rename(MolName("XTL")).molecule();
+            edit_mol = edit_mol.residue(ResIdx(0)).rename(ResName("XTL")).molecule();
+        }
+
         // Replace the atomic coordinates in the template.
         edit_mol = edit_mol.atom(AtomIdx(0)).setProperty(map["coordinates"], coord_oxygen).molecule();
         edit_mol = edit_mol.atom(AtomIdx(1)).setProperty(map["coordinates"], coord_hydrogen0).molecule();
@@ -652,7 +659,7 @@ namespace SireIO
         return new_system;
     }
 
-    System setGromacsWater(const System &system, const QString &model, const PropertyMap &map)
+    System setGromacsWater(const System &system, const QString &model, const PropertyMap &map, bool is_crystal)
     {
         // Create a new system object.
         System new_system;
@@ -705,7 +712,7 @@ namespace SireIO
 
             if (isWater(molecule))
             {
-                molecule = _pvt_setGromacsWater(molecule, template_molecule, _model, has_virtual, map);
+                molecule = _pvt_setGromacsWater(molecule, template_molecule, _model, has_virtual, map, is_crystal);
             }
 
             molgroup.add(molecule);

corelib/src/libs/SireIO/biosimspace.h

@@ -97,7 +97,7 @@ namespace SireIO
                                 const PropertyMap &map = PropertyMap());
 
     Molecule _pvt_setGromacsWater(Molecule &molecule, const Molecule &water, const QString &model, bool has_virtual,
-                                  const PropertyMap &map = PropertyMap());
+                                  const PropertyMap &map = PropertyMap(), bool is_crystal = false);
 
     //! Set all water molecules in the passed system to the appropriate AMBER
     /*! format topology.
@@ -128,11 +128,15 @@ namespace SireIO
         \param map
             A dictionary of user-defined molecular property names.
 
+        \param is_crystal
+            Whether this is a crystal water molecule. If true, then the molecule
+            and residue name will be set to XTL rather than SOL.
+
         \retval system
             The system with updated water topology.
      */
     SIREIO_EXPORT SireSystem::System setGromacsWater(const SireSystem::System &system, const QString &model,
-                                                     const PropertyMap &map = PropertyMap());
+                                                     const PropertyMap &map = PropertyMap(), bool is_crystal = false);
 
     //! Set all water molecules in the passed system to the appropriate AMBER
     /*! format topology.

corelib/src/libs/SireIO/grotop.cpp

@@ -4341,17 +4341,40 @@ GroTop::GroTop(const SireSystem::System &system, const PropertyMap &map)
         isSorted = map["sort"].value().asA<BooleanProperty>().value();
     }
 
-    // Search the system for water molecules.
-    auto waters = system.search("water");
+    // Search for waters and crystal waters. The user can speficy the residue name
+    // for crystal waters using the "crystal_water" property in the map. If the user
+    // wishes to preserve a custom water topology naming, then they can use "skip_water".
+    SelectResult waters;
+    SelectResult xtal_waters;
+    if (map.specified("crystal_water"))
+    {
+        auto xtal_water_resname = map["crystal_water"].source();
+        xtal_waters = system.search(QString("resname %1").arg(xtal_water_resname));
+
+        if (not map.specified("skip_water"))
+        {
+            waters =
+                system.search(
+                    QString("(not mols with property is_non_searchable_water) and (water and not resname %1")
+                    .arg(xtal_water_resname));
+        }
+    }
+    else
+    {
+        waters = system.search("(not mols with property is_non_searchable_water) and water");
+    }
 
     // Extract the molecule numbers of the water molecules.
     auto water_nums = waters.molNums();
 
+    // Extract the molecule numbers of the crystal water molecules.
+    auto xtal_water_nums = xtal_waters.molNums();
+
     // Loop over the molecules to find the non-water molecules.
     QList<MolNum> non_water_nums;
     for (const auto &num : molnums)
     {
-        if (not water_nums.contains(num))
+        if (not water_nums.contains(num) and not xtal_water_nums.contains(num))
             non_water_nums.append(num);
     }
 
@@ -4456,6 +4479,32 @@ GroTop::GroTop(const SireSystem::System &system, const PropertyMap &map)
         }
     }
 
+    // Now add the crystal waters.
+    if (xtal_waters.count() > 0)
+    {
+        // Extract the GroMolType of the first water molecule.
+        auto water_type = GroMolType(system[xtal_water_nums[0]].molecule(), map);
+        auto name = water_type.name();
+        auto molnum = xtal_water_nums[0];
+        auto idx = molnum_to_idx[molnum];
+
+        // Populate the mappings.
+        name_to_mtyp.insert(name, water_type);
+        idx_name_to_mtyp.insert(QPair<int, QString>(idx, water_type.name()), water_type);
+        idx_name_to_example.insert(QPair<int, QString>(idx, name), system[molnum].molecule());
+
+        for (int i = 0; i < xtal_water_nums.count(); ++i)
+        {
+            // Extract the molecule number of the molecule and work out
+            // the index in the system.
+            auto molnum = xtal_water_nums[i];
+            auto idx = molnum_to_idx[molnum];
+
+            // Store the name of the molecule type.
+            mol_to_moltype[idx] = name;
+        }
+    }
+
     QStringList errors;
 
     // first, we need to extract the common forcefield from the molecules

doc/source/changelog.rst

@@ -34,6 +34,10 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 * Forced all new-style modules to import when `sr.use_new_api()` is called.
   This will make it easier to use sire with multiprocessing.
 
+* Added option to allow GROMACS water molecules to be flagged as crystal waters.
+  This means that they will be ignored by ``gmx genion`` when choosing water
+  molecules to replace with ions.
+
 * Please add the changelog entry for your PR here. We will add the link to your PR
   during the code review :-)
 

wrapper/IO/AmberRst.pypp.cpp

@@ -9,6 +9,8 @@ namespace bp = boost::python;
 
 #include "SireBase/generalunitproperty.h"
 
+#include "SireBase/getinstalldir.h"
+
 #include "SireBase/parallel.h"
 
 #include "SireBase/progressbar.h"

wrapper/IO/CMakeAutogenFile.txt

@@ -1,44 +1,44 @@
 # WARNING - AUTOGENERATED FILE - CONTENTS WILL BE OVERWRITTEN!
 set ( PYPP_SOURCES
-       Amber.pypp.cpp
-       AmberPrm.pypp.cpp
+       PDB.pypp.cpp
+       ZmatrixMaker.pypp.cpp
        AmberRst.pypp.cpp
-       AmberRst7.pypp.cpp
-       AmberTraj.pypp.cpp
-       BrokenParser.pypp.cpp
-       CharmmPSF.pypp.cpp
-       Cube.pypp.cpp
-       DCD.pypp.cpp
-       FileTrajectoryParser.pypp.cpp
-       FlexibilityLibrary.pypp.cpp
-       FlexibilityTemplate.pypp.cpp
-       Gro87.pypp.cpp
+       _IO_free_functions.pypp.cpp
        GroAtom.pypp.cpp
-       GroMolType.pypp.cpp
-       GroSystem.pypp.cpp
-       GroTop.pypp.cpp
-       IOBase.pypp.cpp
+       PDBParameters.pypp.cpp
        IOParametersBase.pypp.cpp
-       Mol2.pypp.cpp
-       MoleculeParser.pypp.cpp
-       NullIO.pypp.cpp
-       NullParser.pypp.cpp
-       PDB.pypp.cpp
+       TinkerParameters.pypp.cpp
+       Supplementary.pypp.cpp
        PDB2.pypp.cpp
-       PDBParameters.pypp.cpp
+       AmberPrm.pypp.cpp
+       GroSystem.pypp.cpp
+       FlexibilityTemplate.pypp.cpp
+       TrajectoryMonitor.pypp.cpp
+       IOBase.pypp.cpp
        PerturbationsLibrary.pypp.cpp
+       AmberRst7.pypp.cpp
+       Tinker.pypp.cpp
+       Cube.pypp.cpp
        PerturbationsTemplate.pypp.cpp
+       MoleculeParser.pypp.cpp
+       GroMolType.pypp.cpp
+       Amber.pypp.cpp
+       Gro87.pypp.cpp
+       FlexibilityLibrary.pypp.cpp
+       CharmmPSF.pypp.cpp
+       Mol2.pypp.cpp
        ProtoMS.pypp.cpp
        ProtoMSParameters.pypp.cpp
+       NullIO.pypp.cpp
+       GroTop.pypp.cpp
+       NullParser.pypp.cpp
+       FileTrajectoryParser.pypp.cpp
        SDF.pypp.cpp
-       Supplementary.pypp.cpp
        TRR.pypp.cpp
-       Tinker.pypp.cpp
-       TinkerParameters.pypp.cpp
-       TrajectoryMonitor.pypp.cpp
+       BrokenParser.pypp.cpp
        XTC.pypp.cpp
-       ZmatrixMaker.pypp.cpp
-       _IO_free_functions.pypp.cpp
+       DCD.pypp.cpp
+       AmberTraj.pypp.cpp
        SireIO_containers.cpp
        SireIO_properties.cpp
        SireIO_registrars.cpp

wrapper/IO/DCD.pypp.cpp

@@ -65,10 +65,14 @@ namespace bp = boost::python;
 
 #include "dcd.h"
 
+#include "sire_version.h"
+
 #include "tostring.h"
 
 #include <QDataStream>
 
+#include <QDateTime>
+
 #include <QDebug>
 
 #include <QFile>

wrapper/IO/SireIO_properties.cpp

@@ -4,16 +4,6 @@
 #include "Base/convertproperty.hpp"
 #include "SireIO_properties.h"
 
-#include "SireError/errors.h"
-#include "SireMol/cuttingfunction.h"
-#include "SireMol/molecule.h"
-#include "SireMol/molidx.h"
-#include "SireMol/mover.hpp"
-#include "SireStream/datastream.h"
-#include "iobase.h"
-#include <QDebug>
-#include <QFile>
-#include "iobase.h"
 #include "SireBase/booleanproperty.h"
 #include "SireBase/parallel.h"
 #include "SireBase/progressbar.h"
@@ -49,8 +39,18 @@
 #include <QMutex>
 #include <QTextStream>
 #include "moleculeparser.h"
+#include "SireError/errors.h"
+#include "SireMol/cuttingfunction.h"
+#include "SireMol/molecule.h"
+#include "SireMol/molidx.h"
+#include "SireMol/mover.hpp"
+#include "SireStream/datastream.h"
+#include "iobase.h"
+#include <QDebug>
+#include <QFile>
+#include "iobase.h"
 void register_SireIO_properties()
 {
-    register_property_container< SireIO::IOPtr, SireIO::IOBase >();
     register_property_container< SireIO::MoleculeParserPtr, SireIO::MoleculeParser >();
+    register_property_container< SireIO::IOPtr, SireIO::IOBase >();
 }