[BUG] More reduced box vector hell

[mark-mackey-cresset]

Describe the bug
Box vectors aren't in reduced form even after reducing them

To Reproduce
Run the following code:

system = BSS.IO.readMolecules(
    [args.topology, args.coordinates], reduce_box=True, rotate_box=True
)
BSS.IO.saveMolecules(args.output_filebase, system, ["prm7", "rst7"])

with the inputs being OrigInput.prm7 and FixWaterInput.rst7. Write the output to "Input.prm7" and "Input.rst7"

Now use the resulting prm7/rst7 in a SOMD calculation, and it fails with "RuntimeError: Periodic box vectors must be in reduced form."

Expected behavior
Reading the system in with reduce_box=True and then writing it should write a box in reduced form

(please complete the following information):

• OS: Linux
• Version of sire: Cresset's Sire fork

Additional context

It looks like the culprit is floating point rounding. The box vectors in the buggy output are:

38.7277946 38.7278023 38.7277973 109.4712324 70.5287914 70.5287857

If I edit this to make the last two numbers the same the SOMD calculation works fine

38.7277946 38.7278023 38.7277973 109.4712324 70.5287914 70.5287914

although using the other value does not work:

38.7277946 38.7278023 38.7277973 109.4712324 70.5287857 70.5287857

acos(1/3) is 70.52877937; presumably 70.5287857 is too far off? It might be worth doing a check to see if an input value is "close" to one of the magic angles like acos(1/3), and clamping it to that value to prevent numerical issues like this?

[lohedges]

Were you meant to post a code snippet and input file?

[mark-mackey-cresset]

Form submitted too early - not sure what I pressed to make that happen! Attachments coming shortly.

[lohedges]

No problem. The post magically updated just as I hit comment :-)

[mark-mackey-cresset]

No description provided.

[mark-mackey-cresset]

Argh, can't seem to upload the files - not clear why. I'll email them to you.

[lohedges]

We are using the exact same reduction scheme as OpenMM. The issue is that their C++ API uses exact numerical precision in their checks, whereas their Python API uses a tolerance of 1e-6. The box vectors are correct when reduced by BioSImSpace, but the issue comes when they are written to fixed-width format, loaded back, and set in OpenMM.

We've previously tried to fix this by adding a bias to the reduction so that we always round in a consistent direction. This worked for OpenMM, but caused issues for NPT simulations with AMBER where the box angles could flip during simulation, which would trigger an exception.

It's one of those annoying issues that seems like it should be trivial to fix, but something that works for one engine but doesn't for another, so there's not a general fix. In this case, I think the easy solution would be to just reduce any box vectors before setting them via the OpenMM C++ API. We already have a workaround like this for the BioSimSpace OpenMM process, so it could just be added to SOMD as well.

[lohedges]

Could you try the attempted "fix" that I've commited here, which is on the fix_305 branch This reapplies the reduction on the triclinic space object just before the box vectors are set via the OpenMM C++ API. This should hopefully solve issues due to rounding in the AMBER RST7 input files for SOMD. Note that I've only applied the fix to the non-PME version of the code, but it is trivial to patch in to all SOMD versions if it works.

Cheers.

[mark-mackey-cresset]

Thanks @lohedges - I'll give that a try.

[lohedges]

Any update on this? If it works, I'll merge across to my other fix branch and get this into devel as soon as possible. If you're unable to test, can you share the pert and config files that go with the AMBER files that you provided so that I can run locally.

[mark-mackey-cresset]

Hi Lester. Not been able to get Sire built with the patches yet (we're re-jigging our build infrastructure at the moment) - hoping to get that done in the next day or two. I'll see if I can dig out the rest of the files for the example test case.

[nigel-palmer-cresset]

Hi @lohedges, will the same change be needed in openmmpmefep.cpp .

Thanks

[lohedges]

Yes, that's right. As mentioned up thread, I only applied this to the regular, non-PME, version of the code. If it works, I'll patch it into all SOMD variants.

[mark-mackey-cresset]

sirebug205.zip

Hi Lester. This doesn't seem to work properly. I've attached a test case: if you run it with our branch of Sire without this patch, then you get a starting energy before minimisation of -13930 kcal/mol, and the calculation goes smoothly. With the patch, the starting energy is 2.969e+20 kcal/mol, the minimisation fails, and then anneal-to-lambda crashes. It looks like munging the box vectors makes the coordinates invalid leading to overlapping atoms?

[lohedges]

Thanks, I'll debug locally to see what's going on.

[lohedges]

It seems to be okay if you only reduce the box, not rotate it. Could you try removing the line:

space.rotate();

I assume that we are only having issues with minor rounding errors in the reduction, so there is no need to do a full box rotation.

With this change I see an initial energy of:

Initial energy: -11308.1 kcal mol-1

And after minimisation:

Energy after the minimization: -13936.8 kcal mol-1

This is the same as what I see with no modification to the box vectors.

[lohedges]

For the rotation, the issue is likely that the center of rotation is wrong given the current box vectors and origin. This can be specified when using the BioSimSpace Python API, but here I was just using the default. In any case, rotation shouldn't be needed if the box vectors are just off by a small delta.

[mark-mackey-cresset]

Hi Lester. Latest patch (just reducing, not rotating) seems to have fixed the problem, thanks!

[lohedges]

Great, thanks for confirming. I'll add this to all SOMD flavours and create a PR.