Optimise angles over repeat minimisations. [ref OpenBioSim/cresset#62]

Author: lohedges

src/ghostly/_cli.py

@@ -121,6 +121,17 @@ def run():
         required=False,
     )
 
+    parser.add_argument(
+        "--num-optimise",
+        type=int,
+        help=(
+            "The number of repeats to use when optimising the angle terms involving "
+            "ghost atoms for non-planar triple junctions."
+        ),
+        default=10,
+        required=False,
+    )
+
     parser.add_argument(
         "--output-prefix",
         type=str,
@@ -257,7 +268,13 @@ def run():
     try:
         if args.system is None:
             system = merged.toSystem()._sire_object
-        system = modify(system, k_hard.value(), k_soft.value(), args.optimise_angles)
+        system = modify(
+            system,
+            k_hard.value(),
+            k_soft.value(),
+            args.optimise_angles,
+            args.num_optimise,
+        )
     except Exception as e:
         logger.error(
             f"An error occurred while applying the ghost atom modifications: {e}"

src/ghostly/_ghostly.py

@@ -44,7 +44,7 @@
 del _platform
 
 
-def modify(system, k_hard=100, k_soft=5, optimise_angles=True):
+def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10):
     """
     Apply Boresch modifications to ghost atom bonded terms to avoid non-physical
     coupling between the ghost atoms and the physical region.
@@ -67,6 +67,11 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True):
         Whether to optimise the equilibrium value of the angle terms involving
         ghost atoms for non-planar triple junctions.
 
+    num_optimise : int, optional
+        The number of repeats to average over when optimising the equilibrium
+        value of the angle terms involving ghost atoms for non-planar triple
+        junctions.
+
     Returns
     -------
 
@@ -217,6 +222,7 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True):
                     k_hard=k_hard,
                     k_soft=k_soft,
                     optimise_angles=optimise_angles,
+                    num_optimise=num_optimise,
                 )
 
             # Higher order junction.
@@ -229,6 +235,7 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True):
                     k_hard=k_hard,
                     k_soft=k_soft,
                     optimise_angles=optimise_angles,
+                    num_optimise=num_optimise,
                 )
 
         # Now lambda = 1.
@@ -252,6 +259,7 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True):
                     k_hard=k_hard,
                     k_soft=k_soft,
                     optimise_angles=optimise_angles,
+                    num_optimise=num_optimise,
                     is_lambda1=True,
                 )
 
@@ -265,6 +273,7 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True):
                     k_hard=k_hard,
                     k_soft=k_soft,
                     optimise_angles=optimise_angles,
+                    num_optimise=num_optimise,
                     is_lambda1=True,
                 )
 
@@ -597,6 +606,7 @@ def _triple(
     k_hard=100,
     k_soft=5,
     optimise_angles=True,
+    num_optimise=10,
     is_lambda1=False,
 ):
     r"""
@@ -639,6 +649,10 @@ def _triple(
         Whether to optimise the equilibrium value of the angle terms involving
         ghost atoms for non-planar triple junctions.
 
+    num_optimise : int, optional
+        The number of repeats to use when optimising the angle terms involving
+        ghost atoms for non-planar triple junctions.
+
     is_lambda1 : bool, optional
         Whether the junction is at lambda = 1.
 
@@ -870,32 +884,47 @@ def _triple(
             .commit()
         )
 
-        # Optimise the equilibrium value of theta0 for the softened angle term.
+        # Optimise the equilibrium value of theta0 for the softened angle terms.
         if optimise_angles:
             _logger.debug("  Optimising equilibrium values for softened angles.")
 
             import sire.morph as _morph
             from sire.units import radian as _radian
 
-            # Minimise the molecule.
-            min_mol = _morph.link_to_reference(mol)
-            minimiser = min_mol.minimisation(
-                lambda_value=1.0 if is_lambda1 else 0.0,
-                constraint="none",
-                platform="cpu",
-            )
-            minimiser.run()
-
-            # Commit the changes.
-            min_mol = minimiser.commit()
-
-            # Get the equilibrium angle values.
+            # Initialise the equilibrium angle values.
             theta0s = {}
             for idx in angle_idxs:
-                try:
-                    theta0s[idx] = min_mol.angles(*idx).sizes()[0].to(_radian)
-                except:
-                    raise ValueError(f"Could not find optimised angle term: {idx}")
+                theta0s[idx] = []
+
+            # Perform multiple minimisations to get an average for the theta0 values.
+            for _ in range(num_optimise):
+                # Minimise the molecule.
+                min_mol = _morph.link_to_reference(mol)
+                minimiser = min_mol.minimisation(
+                    lambda_value=1.0 if is_lambda1 else 0.0,
+                    constraint="none",
+                    platform="cpu",
+                )
+                minimiser.run()
+
+                # Commit the changes.
+                min_mol = minimiser.commit()
+
+                # Get the equilibrium angle values.
+                for idx in angle_idxs:
+                    try:
+                        theta0s[idx].append(min_mol.angles(*idx).sizes()[0].to(_radian))
+                    except:
+                        raise ValueError(f"Could not find optimised angle term: {idx}")
+
+            # Compute the mean and standard error.
+            import numpy as _np
+
+            theta0_means = {}
+            theta0_stds = {}
+            for idx in theta0s:
+                theta0_means[idx] = _np.mean(theta0s[idx])
+                theta0_stds[idx] = _np.std(theta0s[idx]) / _np.sqrt(len(theta0s[idx]))
 
             # Get the existing angles.
             angles = mol.property("angle" + suffix)
@@ -913,7 +942,8 @@ def _triple(
                 # This is the softened angle term.
                 if idx in angle_idxs:
                     # Get the optimised equilibrium angle.
-                    theta0 = theta0s[idx]
+                    theta0 = theta0_means[idx]
+                    std = theta0_stds[idx]
 
                     # Create the new angle function.
                     amber_angle = _SireMM.AmberAngle(k_soft, theta0)
@@ -926,7 +956,7 @@ def _triple(
 
                     _logger.debug(
                         f"  Optimising angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], "
-                        f"{p.function()} --> {expression}"
+                        f"{p.function()} --> {expression} (std err: {std:.3f} radian)"
                     )
 
                 else:
@@ -952,6 +982,7 @@ def _higher(
     k_hard=100,
     k_soft=5,
     optimise_angles=True,
+    num_optimise=10,
     is_lambda1=False,
 ):
     r"""
@@ -984,6 +1015,10 @@ def _higher(
         Whether to optimise the equilibrium value of the angle terms involving
         ghost atoms for non-planar triple junctions.
 
+    num_optimise : int, optional
+        The number of repeats to use when optimising the angle terms involving
+        ghost atoms for non-planar triple junctions.
+
     is_lambda1 : bool, optional
         Whether the junction is at lambda = 1.
 
@@ -1072,6 +1107,8 @@ def _higher(
         k_hard=k_hard,
         k_soft=k_soft,
         optimise_angles=optimise_angles,
+        num_optimise=num_optimise,
+        is_lambda1=is_lambda1,
     )