diff --git a/src/BioSimSpace/Process/_amber.py b/src/BioSimSpace/Process/_amber.py index e7c94c254..f8dd37342 100644 --- a/src/BioSimSpace/Process/_amber.py +++ b/src/BioSimSpace/Process/_amber.py @@ -1413,6 +1413,74 @@ def getCurrentDihedralEnergy(self, time_series=False, region=0, soft_core=False) time_series=time_series, region=region, soft_core=soft_core, block=False ) + def getCMAPEnergy(self, time_series=False, region=0, soft_core=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord( + "CMAP", + time_series=time_series, + unit=_Units.Energy.kcal_per_mol, + region=region, + soft_core=soft_core, + block=block, + ) + + def getCurrentCMAPEnergy(self, time_series=False, region=0, soft_core=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy( + time_series=time_series, region=region, soft_core=soft_core, block=False + ) + def getElectrostaticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): diff --git a/src/BioSimSpace/Process/_gromacs.py b/src/BioSimSpace/Process/_gromacs.py index bdcb49f0b..16c8e6f68 100644 --- a/src/BioSimSpace/Process/_gromacs.py +++ b/src/BioSimSpace/Process/_gromacs.py @@ -1382,6 +1382,47 @@ def getCurrentImproperEnergy(self, time_series=False): """ return self.getImproperEnergy(time_series, block=False) + def getCMAPEnergy(self, time_series=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord("CMAPDIH", time_series, _Units.Energy.kj_per_mol, block) + + def getCurrentCMAPEnergy(self, time_series=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy(time_series, block=False) + def getLennardJones14(self, time_series=False, block="AUTO"): """ Get the Lennard-Jones energy between atoms 1 and 4. diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py index add46a4fb..36a52762b 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py @@ -1497,6 +1497,74 @@ def getCurrentDihedralEnergy(self, time_series=False, region=0, soft_core=False) time_series=time_series, region=region, soft_core=soft_core, block=False ) + def getCMAPEnergy(self, time_series=False, region=0, soft_core=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord( + "CMAP", + time_series=time_series, + unit=_Units.Energy.kcal_per_mol, + region=region, + soft_core=soft_core, + block=block, + ) + + def getCurrentCMAPEnergy(self, time_series=False, region=0, soft_core=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy( + time_series=time_series, region=region, soft_core=soft_core, block=False + ) + def getElectrostaticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py index c90c47914..b7125b955 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py @@ -1440,6 +1440,47 @@ def getCurrentImproperEnergy(self, time_series=False): """ return self.getImproperEnergy(time_series, block=False) + def getCMAPEnergy(self, time_series=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord("CMAPDIH", time_series, _Units.Energy.kj_per_mol, block) + + def getCurrentCMAPEnergy(self, time_series=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy(time_series, block=False) + def getLennardJones14(self, time_series=False, block="AUTO"): """ Get the Lennard-Jones energy between atoms 1 and 4. diff --git a/tests/Process/test_single_point_energy.py b/tests/Process/test_single_point_energy.py index 99c6e501f..b86dda785 100644 --- a/tests/Process/test_single_point_energy.py +++ b/tests/Process/test_single_point_energy.py @@ -12,6 +12,16 @@ def ubiquitin_system(): ) +@pytest.fixture(scope="module") +def cmap_system(): + """An ff19SB system, which carries CMAP backbone correction terms.""" + import sire as sr + + return BSS._SireWrappers.System( + sr.load_test_files("zero_k_torsions.gro", "zero_k_torsions.top")._system + ) + + @pytest.mark.skipif( has_amber is False or has_gromacs is False, reason="Requires that both AMBER and GROMACS are installed.", @@ -100,3 +110,42 @@ def test_amber_gromacs_triclinic(ubiquitin_system): nrg_amb = process_amb.getDihedralEnergy().kj_per_mol().value() nrg_gmx = process_gmx.getDihedralEnergy().kj_per_mol().value() assert nrg_amb == pytest.approx(nrg_gmx, rel=1e-2) + + +@pytest.mark.skipif( + has_amber is False or has_gromacs is False, + reason="Requires that both AMBER and GROMACS are installed.", +) +def test_amber_gromacs_cmap(cmap_system): + """Single point CMAP energy comparison between AMBER and GROMACS.""" + + # Create a single-step minimisation protocol. + protocol = BSS.Protocol.Minimisation(steps=1) + + # Create a process to run with AMBER. + process_amb = BSS.Process.Amber(cmap_system, protocol) + + # Create a process to run with GROMACS. + process_gmx = BSS.Process.Gromacs( + cmap_system, protocol, extra_options={"nsteps": 0} + ) + + # Run the AMBER process and wait for it to finish. + process_amb.start() + process_amb.wait() + + # Run the GROMACS process and wait for it to finish. + process_gmx.start() + process_gmx.wait() + + # Compare CMAP energies. (In kJ / mol) + nrg_amb = process_amb.getCMAPEnergy() + nrg_gmx = process_gmx.getCMAPEnergy() + + # The comparison is meaningless if either engine didn't report the term. + assert nrg_amb is not None + assert nrg_gmx is not None + + assert nrg_amb.kj_per_mol().value() == pytest.approx( + nrg_gmx.kj_per_mol().value(), rel=1e-2 + ) diff --git a/tests/Sandpit/Exscientia/Process/test_single_point_energy.py b/tests/Sandpit/Exscientia/Process/test_single_point_energy.py index 3fe986c0d..5fb6e089e 100644 --- a/tests/Sandpit/Exscientia/Process/test_single_point_energy.py +++ b/tests/Sandpit/Exscientia/Process/test_single_point_energy.py @@ -12,6 +12,16 @@ def system(): ) +@pytest.fixture(scope="session") +def cmap_system(): + """An ff19SB system, which carries CMAP backbone correction terms.""" + import sire as sr + + return BSS._SireWrappers.System( + sr.load_test_files("zero_k_torsions.gro", "zero_k_torsions.top")._system + ) + + @pytest.mark.skipif( has_amber is False or has_gromacs is False or has_pyarrow is False, reason="Requires that AMBER, GROMACS, and pyarrow are installed.", @@ -96,3 +106,42 @@ def test_amber_gromacs_triclinic(system): nrg_amb = process_amb.getDihedralEnergy().kj_per_mol().value() nrg_gmx = process_gmx.getDihedralEnergy().kj_per_mol().value() assert nrg_amb == pytest.approx(nrg_gmx, rel=1e-2) + + +@pytest.mark.skipif( + has_amber is False or has_gromacs is False or has_pyarrow is False, + reason="Requires that AMBER, GROMACS, and pyarrow are installed.", +) +def test_amber_gromacs_cmap(cmap_system): + """Single point CMAP energy comparison between AMBER and GROMACS.""" + + # Create a single-step minimisation protocol. + protocol = BSS.Protocol.Minimisation(steps=1) + + # Create a process to run with AMBER. + process_amb = BSS.Process.Amber(cmap_system, protocol) + + # Create a process to run with GROMACS. + process_gmx = BSS.Process.Gromacs( + cmap_system, protocol, extra_options={"nsteps": 0} + ) + + # Run the AMBER process and wait for it to finish. + process_amb.start() + process_amb.wait() + + # Run the GROMACS process and wait for it to finish. + process_gmx.start() + process_gmx.wait() + + # Compare CMAP energies. (In kJ / mol) + nrg_amb = process_amb.getCMAPEnergy() + nrg_gmx = process_gmx.getCMAPEnergy() + + # The comparison is meaningless if either engine didn't report the term. + assert nrg_amb is not None + assert nrg_gmx is not None + + assert nrg_amb.kj_per_mol().value() == pytest.approx( + nrg_gmx.kj_per_mol().value(), rel=1e-2 + )