Add functionality to fix LJ sigmas for ghost atoms.

src/somd2/runner/_runner.py

@@ -151,12 +151,14 @@ def __init__(self, system, config):
 
             # Only check for light atoms by the maxium end state mass if running
             # in SOMD1 compatibility mode. Ghost atoms are considered light when
-            # adding bond constraints.
+            # adding bond constraints. Also fix the LJ sigma for ghost atoms so
+            # it isn't scaled to zero.
             self._config._extra_args["ghosts_are_light"] = True
             self._config._extra_args["check_for_h_by_max_mass"] = True
             self._config._extra_args["check_for_h_by_mass"] = False
             self._config._extra_args["check_for_h_by_element"] = False
             self._config._extra_args["check_for_h_by_ambertype"] = False
+            self._config._extra_args["fix_ghost_sigmas"] = True
 
         # Check for a periodic space.
         self._check_space()

src/somd2/runner/_somd1.py

@@ -59,8 +59,9 @@ def _make_compatible(system):
     except KeyError:
         raise KeyError("No perturbable molecules in the system")
 
-    # Store a dummy element.
-    dummy = _SireMol.Element("Xx")
+    # Store a dummy element and ambertype.
+    element_dummy = _SireMol.Element("Xx")
+    amber_dummy = "du"
 
     for mol in pert_mols:
         # Store the molecule info.
@@ -69,9 +70,43 @@ def _make_compatible(system):
         # Get an editable version of the molecule.
         edit_mol = mol.edit()
 
-        ##########################
-        # First process the bonds.
-        ##########################
+        #####################################
+        # First fix the ghost atom LJ sigmas.
+        #####################################
+
+        for atom in mol.atoms():
+            # Lambda = 0 state is a dummy, use sigma from the lambda = 1 state.
+            if (
+                atom.property("element0") == element_dummy
+                or atom.property("ambertype0") == amber_dummy
+            ):
+                lj0 = atom.property("LJ0")
+                lj1 = atom.property("LJ1")
+                edit_mol = (
+                    edit_mol.atom(atom.index())
+                    .set_property(
+                        "LJ0", _SireMM.LJParameter(lj1.sigma(), lj0.epsilon())
+                    )
+                    .molecule()
+                )
+            # Lambda = 1 state is a dummy, use sigma from the lambda = 0 state.
+            elif (
+                atom.property("element1") == element_dummy
+                or atom.property("ambertype1") == amber_dummy
+            ):
+                lj0 = atom.property("LJ0")
+                lj1 = atom.property("LJ1")
+                edit_mol = (
+                    edit_mol.atom(atom.index())
+                    .set_property(
+                        "LJ1", _SireMM.LJParameter(lj0.sigma(), lj1.epsilon())
+                    )
+                    .molecule()
+                )
+
+        ########################
+        # Now process the bonds.
+        ########################
 
         new_bonds0 = _SireMM.TwoAtomFunctions(mol.info())
         new_bonds1 = _SireMM.TwoAtomFunctions(mol.info())
@@ -534,17 +569,26 @@ def _has_dummy(mol, idxs, is_lambda1=False):
         Whether a dummy atom is present.
     """
 
-    # Set the element property associated with the end state.
+    # We need to check by ambertype too since this molecule may have been
+    # created via sire.morph.create_from_pertfile, in which case the element
+    # property will have been set to the end state with the largest mass, i.e.
+    # may no longer by a dummy.
     if is_lambda1:
-        prop = "element1"
+        element_prop = "element1"
+        ambertype_prop = "ambertype1"
     else:
-        prop = "element0"
+        element_prop = "element0"
+        ambertype_prop = "ambertype0"
 
-    dummy = _SireMol.Element(0)
+    element_dummy = _SireMol.Element(0)
+    ambertype_dummy = "du"
 
     # Check whether an of the atoms is a dummy.
     for idx in idxs:
-        if mol.atom(idx).property(prop) == dummy:
+        if (
+            mol.atom(idx).property(element_prop) == element_dummy
+            or mol.atom(idx).property(ambertype_prop) == ambertype_dummy
+        ):
             return True
 
     return False
@@ -573,21 +617,36 @@ def _is_dummy(mol, idxs, is_lambda1=False):
         Whether each atom is a dummy.
     """
 
-    # Set the element property associated with the end state.
+    # We need to check by ambertype too since this molecule may have been
+    # created via sire.morph.create_from_pertfile, in which case the element
+    # property will have been set to the end state with the largest mass, i.e.
+    # may no longer by a dummy.
     if is_lambda1:
-        prop = "element1"
+        element_prop = "element1"
+        ambertype_prop = "ambertype1"
     else:
-        prop = "element0"
+        element_prop = "element0"
+        ambertype_prop = "ambertype0"
 
-    # Store a dummy element.
-    dummy = _SireMol.Element(0)
+    if is_lambda1:
+        element_prop = "element1"
+        ambertype_prop = "ambertype1"
+    else:
+        element_prop = "element0"
+        ambertype_prop = "ambertype0"
+
+    element_dummy = _SireMol.Element(0)
+    ambertype_dummy = "du"
 
     # Initialise a list to store the state of each atom.
     is_dummy = []
 
     # Check whether each of the atoms is a dummy.
     for idx in idxs:
-        is_dummy.append(mol.atom(idx).property(prop) == dummy)
+        is_dummy.append(
+            mol.atom(idx).property(element_prop) == element_dummy
+            or mol.atom(idx).property(ambertype_prop) == ambertype_dummy
+        )
 
     return is_dummy