OpenBioSim / sire

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

Add map option to prevent peturbation of the LJ sigma value for ghost atoms #198

Closed lohedges closed 5 months ago

lohedges commented 5 months ago

This PR adds a fix_ghost_sigmas map option that can be used to disable perturbation of the LJ sigma parameter, e.g. when creating a perturbable molecule using a SOMD1 pert file, where both sigma and epsilon will have been zeroed for ghost atoms. This PR supersedes #195.

Suggested reviewers:

@chryswoods

lohedges commented 5 months ago

You're right, I've forgotten to check both end states when creating the list of ghost atoms. Yes, please do take a look. I considered adding this functionality in multiple places and wasn't sure what was best.

lohedges commented 5 months ago

Actually, I think this is okay (now I've remembered why I did it like this). For perturbable molecules, the comstructFromAmber method is called twice, swapping the order of the parameters objects for the end states. Here I'm just only concerned with ghosts atoms in the "reference" each call, where the reference is first lambda = 1, then lambda = 0. If a ghost is present, then I just use sigma from the opposite state. I don't need you detect ghosts in both the states at the same time, since it's processed in two passes. I still think there's probably a better way to do this, though.

chryswoods commented 5 months ago

That makes sense. My thought about moving it to PerturbableOpenMMMoleucle is that it would be an easier loop over the arrays for sigma, and if sigma is zero in one array, then copy the value from another array. The aim is to remove all zero sigmas for perturbing atoms, not just those of ghost atoms.

Also, doing it here would avoid creating and checking the QList you created to hold ghost atoms. I did a bit of optimisation work before to remove containers in these loops as the memory allocation and then O(n) searching per atom ended up being a little slow for large systems (e.g. perturbable proteins). Moving to PeturbableOpenMMMoleucle will be a single O(n) scan per molecule plus no memory allocation.

chryswoods commented 5 months ago

I've made the changes suggested. The LJParameter::isDummy function now only checks epsilon (with updated documentation to say what it is doing). I moved the check for zero sigmas into PerturbableOpenMMMolecule, and added a map option for its constructor so that it can process map options. I changed the name of the option from fix_ghost_sigmas to fix_perturbable_zero_sigmas to better reflect the way I have changed the algorithm (it now prevents any near-zero sigma for any perturbable atom, unless both end states are near zero).

There's a large change because I had to regenerate the wrappers in SireMM and Convert/SireOpenMM. This was the first time since the new copy<T>() function, so all of the C++ wrapper files were updated.