diff --git a/corelib/src/libs/SireIO/biosimspace.cpp b/corelib/src/libs/SireIO/biosimspace.cpp index ee3604d36..475fd6a2f 100644 --- a/corelib/src/libs/SireIO/biosimspace.cpp +++ b/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,7 +485,7 @@ 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()); @@ -493,7 +493,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(); } @@ -539,6 +539,13 @@ namespace SireIO coord_oxygen = molecule.atom(idx).property(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); diff --git a/corelib/src/libs/SireIO/biosimspace.h b/corelib/src/libs/SireIO/biosimspace.h index a0b18d7b2..5f7541c5a 100644 --- a/corelib/src/libs/SireIO/biosimspace.h +++ b/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. diff --git a/corelib/src/libs/SireIO/grotop.cpp b/corelib/src/libs/SireIO/grotop.cpp index 00cd17a02..653c543cf 100644 --- a/corelib/src/libs/SireIO/grotop.cpp +++ b/corelib/src/libs/SireIO/grotop.cpp @@ -4341,17 +4341,40 @@ GroTop::GroTop(const SireSystem::System &system, const PropertyMap &map) isSorted = map["sort"].value().asA().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 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(idx, water_type.name()), water_type); + idx_name_to_example.insert(QPair(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 diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 2f76be628..d21c3c84e 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -34,6 +34,10 @@ organisation on `GitHub `__. * 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 :-) diff --git a/wrapper/IO/AmberRst.pypp.cpp b/wrapper/IO/AmberRst.pypp.cpp index 1b503dd36..70b9cdf5b 100644 --- a/wrapper/IO/AmberRst.pypp.cpp +++ b/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" diff --git a/wrapper/IO/CMakeAutogenFile.txt b/wrapper/IO/CMakeAutogenFile.txt index ef4a7028c..8d36de296 100644 --- a/wrapper/IO/CMakeAutogenFile.txt +++ b/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 diff --git a/wrapper/IO/DCD.pypp.cpp b/wrapper/IO/DCD.pypp.cpp index ecd65690c..3c6c65ef5 100644 --- a/wrapper/IO/DCD.pypp.cpp +++ b/wrapper/IO/DCD.pypp.cpp @@ -65,10 +65,14 @@ namespace bp = boost::python; #include "dcd.h" +#include "sire_version.h" + #include "tostring.h" #include +#include + #include #include diff --git a/wrapper/IO/SireIO_properties.cpp b/wrapper/IO/SireIO_properties.cpp index 9f26f536d..9315b796c 100644 --- a/wrapper/IO/SireIO_properties.cpp +++ b/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 -#include -#include "iobase.h" #include "SireBase/booleanproperty.h" #include "SireBase/parallel.h" #include "SireBase/progressbar.h" @@ -49,8 +39,18 @@ #include #include #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 +#include +#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 >(); } diff --git a/wrapper/IO/SireIO_registrars.cpp b/wrapper/IO/SireIO_registrars.cpp index b80000e3c..c65429d55 100644 --- a/wrapper/IO/SireIO_registrars.cpp +++ b/wrapper/IO/SireIO_registrars.cpp @@ -3,77 +3,77 @@ #include "SireIO_registrars.h" -#include "amber.h" +#include "moleculeparser.h" +#include "pdb2.h" +#include "trajectorymonitor.h" #include "amberprm.h" -#include "amberrst.h" #include "amberrst7.h" -#include "ambertraj.h" #include "charmmpsf.h" -#include "dcd.h" -#include "filetrajectory.h" -#include "filetrajectoryparser.h" -#include "flexibilitylibrary.h" -#include "gro87.h" -#include "grotop.h" +#include "tinker.h" #include "iobase.h" -#include "mol2.h" -#include "moleculeparser.h" #include "pdb.h" -#include "pdb2.h" -#include "perturbationslibrary.h" #include "protoms.h" -#include "sdf.h" +#include "flexibilitylibrary.h" +#include "perturbationslibrary.h" +#include "mol2.h" +#include "amberrst.h" +#include "grotop.h" +#include "zmatrixmaker.h" #include "supplementary.h" -#include "tinker.h" -#include "trajectorymonitor.h" -#include "trr.h" +#include "amber.h" +#include "gro87.h" +#include "filetrajectory.h" +#include "sdf.h" +#include "dcd.h" #include "xtc.h" -#include "zmatrixmaker.h" +#include "ambertraj.h" +#include "filetrajectoryparser.h" +#include "trr.h" #include "Helpers/objectregistry.hpp" void register_SireIO_objects() { - ObjectRegistry::registerConverterFor< SireIO::Amber >(); + ObjectRegistry::registerConverterFor< SireIO::NullParser >(); + ObjectRegistry::registerConverterFor< SireIO::BrokenParser >(); + ObjectRegistry::registerConverterFor< SireIO::PDBAtom >(); + ObjectRegistry::registerConverterFor< SireIO::PDB2 >(); + ObjectRegistry::registerConverterFor< SireIO::TrajectoryMonitor >(); ObjectRegistry::registerConverterFor< SireIO::AmberPrm >(); - ObjectRegistry::registerConverterFor< SireIO::AmberRst >(); ObjectRegistry::registerConverterFor< SireIO::AmberRst7 >(); - ObjectRegistry::registerConverterFor< SireIO::AmberTraj >(); ObjectRegistry::registerConverterFor< SireIO::PSFAtom >(); ObjectRegistry::registerConverterFor< SireIO::CharmmParam >(); ObjectRegistry::registerConverterFor< SireIO::CharmmPSF >(); - ObjectRegistry::registerConverterFor< SireIO::DCD >(); - ObjectRegistry::registerConverterFor< SireIO::FileTrajectory >(); - ObjectRegistry::registerConverterFor< SireIO::FileTrajectoryParser >(); + ObjectRegistry::registerConverterFor< SireIO::Tinker >(); + ObjectRegistry::registerConverterFor< SireIO::NullIO >(); + ObjectRegistry::registerConverterFor< SireIO::PDB >(); + ObjectRegistry::registerConverterFor< SireIO::ProtoMS >(); ObjectRegistry::registerConverterFor< SireIO::FlexibilityLibrary >(); ObjectRegistry::registerConverterFor< SireIO::FlexibilityTemplate >(); - ObjectRegistry::registerConverterFor< SireIO::Gro87 >(); - ObjectRegistry::registerConverterFor< SireIO::GroTop >(); - ObjectRegistry::registerConverterFor< SireIO::GroMolType >(); - ObjectRegistry::registerConverterFor< SireIO::GroAtom >(); - ObjectRegistry::registerConverterFor< SireIO::GroSystem >(); - ObjectRegistry::registerConverterFor< SireIO::NullIO >(); + ObjectRegistry::registerConverterFor< SireIO::PerturbationsLibrary >(); + ObjectRegistry::registerConverterFor< SireIO::PerturbationsTemplate >(); ObjectRegistry::registerConverterFor< SireIO::Mol2Atom >(); ObjectRegistry::registerConverterFor< SireIO::Mol2Bond >(); ObjectRegistry::registerConverterFor< SireIO::Mol2Molecule >(); ObjectRegistry::registerConverterFor< SireIO::Mol2Substructure >(); ObjectRegistry::registerConverterFor< SireIO::Mol2 >(); - ObjectRegistry::registerConverterFor< SireIO::NullParser >(); - ObjectRegistry::registerConverterFor< SireIO::BrokenParser >(); - ObjectRegistry::registerConverterFor< SireIO::PDB >(); - ObjectRegistry::registerConverterFor< SireIO::PDBAtom >(); - ObjectRegistry::registerConverterFor< SireIO::PDB2 >(); - ObjectRegistry::registerConverterFor< SireIO::PerturbationsLibrary >(); - ObjectRegistry::registerConverterFor< SireIO::PerturbationsTemplate >(); - ObjectRegistry::registerConverterFor< SireIO::ProtoMS >(); - ObjectRegistry::registerConverterFor< SireIO::SDF >(); + ObjectRegistry::registerConverterFor< SireIO::AmberRst >(); + ObjectRegistry::registerConverterFor< SireIO::GroTop >(); + ObjectRegistry::registerConverterFor< SireIO::GroMolType >(); + ObjectRegistry::registerConverterFor< SireIO::GroAtom >(); + ObjectRegistry::registerConverterFor< SireIO::GroSystem >(); + ObjectRegistry::registerConverterFor< SireIO::ZmatrixMaker >(); ObjectRegistry::registerConverterFor< SireIO::Supplementary >(); - ObjectRegistry::registerConverterFor< SireIO::Tinker >(); - ObjectRegistry::registerConverterFor< SireIO::TrajectoryMonitor >(); - ObjectRegistry::registerConverterFor< SireIO::TRR >(); + ObjectRegistry::registerConverterFor< SireIO::Amber >(); + ObjectRegistry::registerConverterFor< SireIO::Gro87 >(); + ObjectRegistry::registerConverterFor< SireIO::FileTrajectory >(); + ObjectRegistry::registerConverterFor< SireIO::SDF >(); + ObjectRegistry::registerConverterFor< SireIO::DCD >(); ObjectRegistry::registerConverterFor< SireIO::XTC >(); - ObjectRegistry::registerConverterFor< SireIO::ZmatrixMaker >(); + ObjectRegistry::registerConverterFor< SireIO::AmberTraj >(); + ObjectRegistry::registerConverterFor< SireIO::FileTrajectoryParser >(); + ObjectRegistry::registerConverterFor< SireIO::TRR >(); } diff --git a/wrapper/IO/_IO_free_functions.pypp.cpp b/wrapper/IO/_IO_free_functions.pypp.cpp index 17ff3dda0..4787605a1 100644 --- a/wrapper/IO/_IO_free_functions.pypp.cpp +++ b/wrapper/IO/_IO_free_functions.pypp.cpp @@ -651,14 +651,14 @@ void register_free_functions(){ { //::SireIO::setGromacsWater - typedef ::SireSystem::System ( *setGromacsWater_function_type )( ::SireSystem::System const &,::QString const &,::SireBase::PropertyMap const & ); + typedef ::SireSystem::System ( *setGromacsWater_function_type )( ::SireSystem::System const &,::QString const &,::SireBase::PropertyMap const &,bool ); setGromacsWater_function_type setGromacsWater_function_value( &::SireIO::setGromacsWater ); bp::def( "setGromacsWater" , setGromacsWater_function_value - , ( bp::arg("system"), bp::arg("model"), bp::arg("map")=SireBase::PropertyMap() ) - , "Set all water molecules in the passed system to the appropriate GROMACS\nformat topology.\n\nPar:am system\nThe molecular system of interest.\n\nPar:am model\nThe name of the water model.\n\nPar:am map\nA dictionary of user-defined molecular property names.\n\nRetval: system\nThe system with updated water topology.\n" ); + , ( bp::arg("system"), bp::arg("model"), bp::arg("map")=SireBase::PropertyMap(), bp::arg("is_crystal")=(bool)(false) ) + , "Set all water molecules in the passed system to the appropriate GROMACS\nformat topology.\n\nPar:am system\nThe molecular system of interest.\n\nPar:am model\nThe name of the water model.\n\nPar:am map\nA dictionary of user-defined molecular property names.\n\nPar:am is_crystal\nWhether this is a crystal water molecule. If true, then the molecule\nand residue name will be set to XTL rather than SOL.\n\nRetval: system\nThe system with updated water topology.\n" ); }