Merge remote-tracking branch 'exscientia/devel'

Author: lohedges

.github/workflows/Sandpit_exs.yml

@@ -19,7 +19,7 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest", ]
-        python-version: ["3.10",]
+        python-version: ["3.11",]
 
     steps:
       - uses: actions/checkout@v2
@@ -33,7 +33,7 @@ jobs:
 
       - name: Install dependency
         run: |
-          conda install -c conda-forge -c openbiosim/label/main biosimspace python=3.10 ambertools gromacs "sire=2024.2.0" "alchemlyb>=2.1" pytest openff-interchange pint=0.21 rdkit "jaxlib>0.3.7" tqdm
+          conda install -c conda-forge -c openbiosim/label/main biosimspace python=3.11 ambertools gromacs "sire=2024.3.0" "alchemlyb>=2.1" pytest openff-interchange pint=0.21 rdkit "jaxlib>0.3.7" tqdm
           python -m pip install git+https://github.com/Exscientia/MDRestraintsGenerator.git
           # For the testing of BSS.FreeEnergy.AlchemicalFreeEnergy.analysis
           python -m pip install https://github.com/alchemistry/alchemtest/archive/master.zip

python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint.py

@@ -37,7 +37,7 @@
 )
 from sire.units import GeneralUnit as _sire_GeneralUnit
 
-from BioSimSpace.Sandpit.Exscientia.Types._general_unit import (
+from ..Types._general_unit import (
     GeneralUnit as _GeneralUnit,
 )
 from ..Types import Angle as _Angle, Length as _Length, Temperature as _Temperature
@@ -450,7 +450,7 @@ def format_bond(equilibrium_values, force_constants):
             )
 
         # Format the parameters for the angles and dihedrals
-        def format_angle(equilibrium_values, force_constants, restraint_lambda):
+        def format_angle(equilibrium_values, force_constants, restraint_lambda, perturbed=True):
             """
             Format the angle equilibrium values and force constant
             in into the Gromacs topology format.
@@ -465,20 +465,30 @@ def format_angle(equilibrium_values, force_constants, restraint_lambda):
 
             When restraint_lambda is True, the dihedrals will be stored in the dihedral_restraints.
             """
-            converted_equ_val = (
+            if isinstance(equilibrium_values, _Angle):
+                converted_equ_val = equilibrium_values / _degree
+            else:
+                converted_equ_val = (
                 self._restraint_dict["equilibrium_values"][equilibrium_values] / _degree
             )
-            converted_fc = self._restraint_dict["force_constants"][force_constants] / (
-                _kj_per_mol / (_radian * _radian)
-            )
+
+            if isinstance(force_constants, _GeneralUnit):
+                converted_fc = force_constants / (
+                    _kj_per_mol / (_radian * _radian)
+                )
+            else:
+                converted_fc = self._restraint_dict["force_constants"][force_constants] / (
+                    _kj_per_mol / (_radian * _radian)
+                )
+
             par_string = (
                 dihedral_restraints_parameters_string
                 if restraint_lambda
                 else parameters_string
             )
             return par_string.format(
                 eq0="{:.3f}".format(converted_equ_val),
-                fc0="{:.2f}".format(0),
+                fc0="{:.2f}".format(0 if perturbed else converted_fc),
                 eq1="{:.3f}".format(converted_equ_val),
                 fc1="{:.2f}".format(converted_fc),
             )
@@ -493,12 +503,12 @@ def write_bond(key_list, equilibrium_values, force_constants):
                 parameters=format_bond(equilibrium_values, force_constants),
             )
 
-        def write_angle(key_list, equilibrium_values, force_constants):
+        def write_angle(key_list, equilibrium_values, force_constants, func_type=1, perturbed=True):
             return master_string.format(
                 index=format_index(key_list),
-                func_type=1,
+                func_type=func_type,
                 parameters=format_angle(
-                    equilibrium_values, force_constants, restraint_lambda=False
+                    equilibrium_values, force_constants, restraint_lambda=False, perturbed=perturbed
                 ),
             )
 
@@ -539,6 +549,22 @@ def write_dihedral(
         output.append(write_angle(("r2", "r1", "l1"), "thetaA0", "kthetaA"))
         # Angles: r1-l1-l2 (thetaB0, kthetaB)
         output.append(write_angle(("r1", "l1", "l2"), "thetaB0", "kthetaB"))
+        # Bent angle: r2-r1-l1
+        # Center is 90 degree
+        # force constant is set by evaluating the free energy of adding a harmonic angle potential with fc at kcal/mol/rad2 at 135 degree
+        #| fc (kcal/mol/rad2) | FE (kcal/mol) |
+        #|--------------------|---------------|
+        #| control            | 1.05          |
+        #| 0                  | 1.05          |
+        #| 0.1                | 0.99          |
+        #| 1                  | 1.17          |
+        #| 5                  | 2.14          |
+        #| 10                 | 3.46          |
+        #| 100                | 14.92         |
+        # Thus 1 kcal/mol/rad2 is choosen
+        output.append(write_angle(("r2", "r1", "l1"), 90 * _degree, 1 * _kcal_per_mol / _radian**2, func_type=10, perturbed=False))
+        # Bent angle: r2-r1-l1
+        output.append(write_angle(("r1", "l1", "l2"), 90 * _degree, 1 * _kcal_per_mol / _radian**2, func_type=10, perturbed=False))
 
         if restraint_lambda:
             output.append("[ dihedral_restraints ]")

tests/Sandpit/Exscientia/FreeEnergy/test_restraint.py

@@ -14,7 +14,6 @@
 from BioSimSpace.Sandpit.Exscientia.Units.Energy import kcal_per_mol
 from BioSimSpace.Sandpit.Exscientia.Units.Length import angstrom
 from BioSimSpace.Sandpit.Exscientia.Units.Temperature import kelvin
-
 from tests.conftest import has_gromacs
 
 # Store the tutorial URL.
@@ -201,8 +200,28 @@ def test_angle(self, Topology):
         assert aj == "1496"
         assert ak == "1497"
 
+    def test_bented_angle(self, Topology):
+        ai, aj, ak, type, thA, kA, thB, kB = Topology[8].split()
+        assert ai == "2"
+        assert aj == "1"
+        assert ak == "1496"
+        assert type == "10"
+        assert thA == "90.000"
+        assert kA == "4.18"
+        assert thB == "90.000"
+        assert kB == "4.18"
+        ai, aj, ak, type, thA, kA, thB, kB = Topology[9].split()
+        assert ai == "1"
+        assert aj == "1496"
+        assert ak == "1497"
+        assert type == "10"
+        assert thA == "90.000"
+        assert kA == "4.18"
+        assert thB == "90.000"
+        assert kB == "4.18"
+
     def test_dihedral(self, Topology):
-        ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[10].split()
+        ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[12].split()
         assert ai == "3"
         assert aj == "2"
         assert ak == "1"
@@ -212,12 +231,12 @@ def test_dihedral(self, Topology):
         assert kA == "0.00"
         assert phiB == "148.396"
         assert kB == "41.84"
-        ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[11].split()
+        ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[13].split()
         assert ai == "2"
         assert aj == "1"
         assert ak == "1496"
         assert al == "1497"
-        ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[12].split()
+        ai, aj, ak, al, type, phiA, kA, phiB, kB = Topology[14].split()
         assert ai == "1"
         assert aj == "1496"
         assert ak == "1497"
@@ -233,8 +252,8 @@ def Topology(boresch_restraint):
         ).split("\n")
 
     def test_dihedral(self, Topology):
-        assert "dihedral_restraints" in Topology[8]
-        ai, aj, ak, al, type, phiA, dphiA, kA, phiB, dphiB, kB = Topology[10].split()
+        assert "dihedral_restraints" in Topology[10]
+        ai, aj, ak, al, type, phiA, dphiA, kA, phiB, dphiB, kB = Topology[12].split()
         assert ai == "3"
         assert aj == "2"
         assert ak == "1"