OpenBioSim / sire

Sire Molecular Simulations Framework
https://sire.openbiosim.org
GNU General Public License v3.0
39 stars 11 forks source link

Add function to create ions from AMBER templates #210

Closed lohedges closed 2 months ago

lohedges commented 2 months ago

This PR adds simple functions for creating Na+ and Cl- ions using AMBER template files. The idea is to use these in BioSimSpace and SOMD2 to provide templates for alchemical ion pertubations. For example, something like:

import BioSimSpace as BSS
import sire as sr

# Here we have a system called mols that contains a single perturbable molecule (molecule 0) and waters.

# First, find the furthers water from the perturbable molecule.
water = mols["furthest water from molidx 0"]

# Create an ion using the same coordinates as the oxygen (atom 0).
ion = sr.legacy.IO.createChlorineIon(water[0].coordinates(), "tip3p")

# Merge the water and ion. I'd do this directly in Sire, but the default alignment currently fails 
# with:
#     RuntimeError: SireError::program_bug: Attempt to find the alignment matrix failed to produce a
#     valid rotation matrix. Determinant should be 1. It is equal to 0.
# Will need to resolve this.
mapping = BSS.Align.matchAtoms(
    BSS._SireWrappers.Molecule(water),
    BSS._SireWrappers.Molecule(ion),
    scoring_function="rmsd_flex_align"
)
merged = BSS.Align.merge(
    BSS._SireWrappers.Molecule(water),
    BSS._SireWrappers.Molecule(ion),
    mapping=mapping
)

# Replace the water with the perturbable molecule.
mols.remove(water)
mols.add(merged._sire_object)

The ion information is so minimal that we could probably store the template another way if we wished to preserve space. This just mimics the approach used for the water topologies and gives the option of using parameters optimised for the typical water models that we use (can add more in future). I could also consolidate the two functions to call a common _pvt_create_ion function, or similar. However, they were so short that I didn't think it was worth the effort.

Suggested reviewers:

@chryswoods

jmichel80 commented 2 months ago

I think there aren’t enough degrees of freedom in a monoatomic ion for the alignment algorithm to work.

lohedges commented 2 months ago

Yes, exactly. I think we should have some option to support this, though, i.e. detect an ion and just replace the coordinates with those of the mapped atom rather than attempting to align. The issue here is that we use the alignment when scoring the mappings, so it's performed by default.