Don't perturb LJ sigma for ghost atoms

[lohedges]

This PR switches to not perturbing the LJ sigma parameter for ghost atoms so that merged molecules created by BioSimSpace are compatible with SOMD2 by default. We still preserve the existing behaviour (zero sigma and epsilon) for SOMD1. The PR also includes a few modifications to harden the checking of dummy atoms during the pert file writing stage. (This isn't needed during merging, since we are explicitly setting the appropriate molecule and atom properties.)

• I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): [y]
• I confirm that I have permission to release this code under the GPL3 license: [y]

Suggested reviewers:

@chryswoods

[lohedges]

Looks like ambertype isn't guaranteed to be present in all molecules, e.g. the alchemical ions from the Exs sandpit. I'll add a check for the property before accessing it tomorrow.

[lohedges]

I've fixed the issue with the missing ambertype property. There still appears to be an issue using non-zero LJ sigma values for dummy atoms with GROMACS that I am debugging. With the merge code updated to fix the sigma to the opposite end state value, gmx grompp now fails when trying to create a binary run file, e.g.:

E               RuntimeError: Unable to generate GROMACS binary run input file.
E
E               'gmx grompp' reported the following errors:
E               ERROR 1 [file gromacs.top, line 69]:
E                 Atomtype hcx_du not found

I'll check the topology file generated during the test to see what's gone wrong.

[lohedges]

This is because the GroTop parser uses this logic to detect dummy atoms:

            if (elem.nProtons() == 0 and lj.isDummy())
            {
                if (is_perturbable)
                    atomtype += "_du";

                // Only label dummies for regular simulations.
                else if (not was_perturbable)
                    particle_type = "D";
            }

Clearly lj.isDummy() is no longer valid, since we support dummy atoms with non-zero sigma. Perhaps I should just check the element, or the element and the LJ epsilon value. @chryswoods: we'll need to add a fix to whatever PR supersedes this one.

[chryswoods]

Will do - I will sort that this weekend. I'll change isDummy over to testing only the epsilon value (as described in the other thread)

[lohedges]

Brilliant, many thanks.

[lohedges]

Everything is now working following the latest Sire update.

[lohedges]

Just to confirm that I have now checked that the three different approaches for fixing zero LJ sigma values are giving the same result. We now have:

1. The fix_perturbable_zero_sigmas map option in sire that applies the fix when constructing the perturbable OpenMMMolecule object.
2. A fix in the somd1 compatibility layer within somd2.
3. The updated code in this PR, which ensures that no zero LJ sigma values are present following a merge in BioSimSpace.

The following plot shows a minimised lambda potential energy scan for the ethane --> methanol perturbation. The sire and somd2 results were generated using an old BioSimSpace merged molecule, where zero LJ sigma values were present. The BioSimSpace result just uses the new merging code here, with the fixes from 1) and 2) not applied. I have also generated XML files from OpenMM systems for each data point and have confirmed that (after sorting) they are identical.

[lohedges]

I'm just wondering about how to preserve the existing somd1 behaviour with the new approach. If we only fix LJ sigmas for dummy atoms, then everything is easy since it's trivial to detect dummy atoms when the pert file is written. When fixing all zero LJ sigma values, as is done here, then things are more difficult. As an example, for ethane --> methanol I find the following atom record is different when I diff the old and new pert files:

Old:

    atom
        name           H4
        initial_type   hc
        final_type     ho
        initial_LJ     2.64953 0.01570
        final_LJ       0.00000 0.00000
        initial_charge 0.03145
        final_charge   0.39600
    endatom

New:

    atom
        name           H4
        initial_type   hc
        final_type     ho
        initial_LJ     2.64953 0.01570
        final_LJ       2.64953 0.00000
        initial_charge 0.03145
        final_charge   0.39600
    endatom

@jmichel80: Are there any issues with this approach? If so, then I could always add another atom property, e.g. LJ0_somd1 and LJ1_somd1 that stores specific LJ parameters to use with somd1 only.

[jmichel80]

This would change the behaviour of the softcore I think and could cause changes in the efficiency of transformations with somd1, it would be best to preserve the somd1 protocol.

[lohedges]

Here's a comparison of the ethane <--> methanol PMFs for somd1 with and without the LJ sigma modification, i.e. where the only difference is the change to the non-dummy zero LJ sigma term (for hc to ho) shown above. These are all short runs and each have unique starting inputs. The overal dG is comparable (within error) and both are reversible,

[chryswoods]

Looks good :-)