somd-freenrg crashes when ligand atoms are completely non-interacting

[fjclark]

Describe the bug
I've recently been trying to implement multiple distance restraints for absolute binding free energy calculations in BioSimSpace and SOMD (modified BSS, modified sire, and an example of how to run a complete calculation). The idea is that the ligand intermolecular interactions are removed with many distance restraints active, then there's a final stage where all but one of the restraints are released. The correction for releasing the ligand to the standard state volume can then easily be calculated.

During the final release stage, where the ligand is completely non-interacting I've been getting crashes. The config file I was using with my modified version of sire was:

save coordinates = True
ncycles = 50
nmoves = 25000
ncycles_per_snap = 4
buffered coordinates frequency = 500
timestep = 4 * femtosecond
reaction field dielectric = 78.4
cutoff type = cutoffperiodic
cutoff distance = 12 * angstrom
barostat = True
pressure = 1.00000 atm
thermostat = True
temperature = 300.00 kelvin
constraint = hbonds-notperturbed
energy frequency = 100
lambda array = (0.0, 0.125, 0.25, 0.375, 0.5, 1.0)
lambda_val = 0.0
perturbed residue number = 293
random seed = -1
gpu = 0
center solute = False
minimise = True
hydrogen mass repartitioning factor = 3
use distance restraints = True
use permanent distance restraints = True
distance restraints dictionary = {(1474, 4687): (4.432198728689354, 20.0, 0.5557321930560404), (1472, 4686): (3.164351608202915, 20.0, 0.6751681461642662), (1474, 4690): (4.700974930995326, 20.0, 0.6701990692583388), (1488, 4688): (6.7858681626072315, 20.0, 0.7488970818377805), (2444, 4677): (4.863475388033847, 20.0, 0.7671705708955612), (1472, 4691): (6.367026275053535, 20.0, 0.821252569114769), (2261, 4683): (6.773986513796716, 20.0, 0.8828616026118992), (1490, 4685): (7.720582022157998, 20.0, 0.9347068069890287), (286, 4674): (3.8127377709679453, 20.0, 0.9511268213071902), (374, 4680): (6.181106731821419, 20.0, 0.9637979205880889), (1521, 4692): (3.7598878700775633, 20.0, 0.9830902561404447), (1529, 4684): (7.87209938990775, 20.0, 1.0582870196044247), (374, 4693): (4.975946331072206, 20.0, 1.063773425343931), (346, 4676): (8.470528418491696, 20.0, 1.0492604125738865), (284, 4675): (4.725973194492003, 20.0, 1.1798376555744614), (390, 4682): (8.350568045090748, 20.0, 1.2118321370763052), (2436, 4679): (4.32826597759919, 20.0, 1.2565532487514752), (307, 4678): (7.55090821416431, 20.0, 1.3127014576568437), (1421, 4681): (8.80522670899781, 20.0, 1.3688291245090927), (259, 4673): (9.195983552051707, 20.0, 2.1342422686136686)}
permanent distance restraints dictionary = {(1474, 4689): (3.875183995573175, 20.0, 0.5185270232355128)}
turn on receptor-ligand restraints mode = True

And the pert file just mapped dummy atoms to dummy atoms, e.g.:

version 1
molecule LIG
    atom
        name           C
        initial_type   du
        final_type     du
        initial_LJ     0.00000 0.00000
        final_LJ       0.00000 0.00000
        initial_charge 0.00000
        final_charge   0.00000
    endatom
    atom
        name           C1
        initial_type   du
        final_type     du
        initial_LJ     0.00000 0.00000
        final_LJ       0.00000 0.00000
        initial_charge 0.00000
        final_charge   0.00000
    endatom
...

I still observed crashes at random values of lambda if I dropped the timestep to 1 fs and if I switched to Langevin dynamics.

To check if this was due to changes I had made to sire or the restraints, I ran some sets of calculations with Sire 2023.4.0 (a073701/GitHub/2023-08-07T17:45:41+01:00 | UNCLEAN) but with ligand still non-interacting and with no restraints. Please see the input files. For each test I ran four repeats of each window. I got two crashes at lambda = 0 (Nan) and all simulations at lambda =1 resulted in Segmentation fault (core dumped) somd-freenrg -C somd.cfg -p CUDA -m somd.pert -c somd.rst7 -t somd.prm7.

Have I chosen some configuration parameters poorly, or should this not be happening?

Also, when I watch the trajectories back, the ligand makes large jumps after every 50th frame. There are 600 frames, and 600 / 50 =12 and ncycles / ncycles_per_snap = 12.5, so I assume this is an artifact of how the trajectories are written?

To Reproduce
Please see the input files. Run run.sh. with sire 2023.4.0. (I observed two crashes in four repeats of the stage).

Expected behavior
Stable MD.

(please complete the following information):

• OS: Ubuntu 22.04
• Version of Python: 3.10.12
• Version of sire: 2023.4.0 (a073701/GitHub/2023-08-07T17:45:41+01:00 | UNCLEAN)
• I confirm that I have checked this bug still exists in the latest released version of sire: [no] (As I would rather avoid rerunning too many simulations unless there have been very recent changes which are likely to affect this)

Thank you.

[jmichel80]

hi @fjclark - do you also observe crashes with MIF inputs from the MDR paper repo (that presumably did not cause runtime errors with earlier versions of sire used for the MDR paper) ?

[lohedges]

Can confirm that I see the same segmentation fault at lambda = 1. Will see if I can trace it.

[lohedges]

Just to note that I don't see a mapping from dummy ambertype to dummy ambertype in your example pert files, but this doesn't seem to be the source of the problem:

version 1
molecule LIG
    atom
        name           C
        initial_type   C3
        final_type     du
        initial_LJ     0.00000 0.00000
        final_LJ       0.00000 0.00000
        initial_charge 0.00000
        final_charge   0.00000
    endatom
    atom
        name           C1
        initial_type   C3
        final_type     du
        initial_LJ     0.00000 0.00000
        final_LJ       0.00000 0.00000
        initial_charge 0.00000
        final_charge   0.00000
    endatom
...

[fjclark]

hi @fjclark - do you also observe crashes with MIF inputs from the MDR paper repo (that presumably did not cause runtime errors with earlier versions of sire used for the MDR paper) ?

Hi @jmichel80 I'll check now (yes I did not observe these crashes previously).

[lohedges]

For lambda = 1, I can confirm that the segfault occurs on this line of OpenMMMD.py. I'll add some debugging statements to the C++ code to see where it's blowing up. (Obviously this bit works at lambda = 0, so NaNs there might be completely unrelated.)

[fjclark]

Just to note that I don't see a mapping from dummy ambertype to dummy ambertype in your example pert files, but this doesn't seem to be the source of the problem:

version 1
molecule LIG
    atom
        name           C
        initial_type   C3
        final_type     du
        initial_LJ     0.00000 0.00000
        final_LJ       0.00000 0.00000
        initial_charge 0.00000
        final_charge   0.00000
    endatom
    atom
        name           C1
        initial_type   C3
        final_type     du
        initial_LJ     0.00000 0.00000
        final_LJ       0.00000 0.00000
        initial_charge 0.00000
        final_charge   0.00000
    endatom
...

Sorry, I should have explained - this is because the standard version of sire doesn't allow atoms to start and end as dummy atoms, so I've had to modify the pert file. For my modified version of sire, this is allowed (4513564) and the pert file is as shown in my original post.

[lohedges]

Also, the pert files are missing an endmolecule record at the end, but again this doesn't seem to affect the segfault. Looking at the output of the parser it seems to get to the initial_charge of the final atom, then crashes:

...
"name O"
"initial_type O3"
"final_type du"
"initial_LJ 3.03981 0.21021"
"final_LJ 0.00000 0.00000"
"initial_charge -0.00000"
"final_charge -0.00000"
"endatom"
"atom"
"name O1"
"initial_type O3"
"final_type du"
"initial_LJ 0.00000 0.00000"
"final_LJ 0.00000 0.00000"
"initial_charge0.00000"
[1]    110167 segmentation fault (core dumped)  somd-freenrg -C somd.cfg -t somd.prm7 -c somd.rst7 -m somd.pert -p CUDA

[lohedges]

Ah, I see, there is no space in initial_charge0.00000 for the final record. Printing the line and the words (the records after the line is split) makes it clear what's going wrong. Not sure why I didn't spot this above 🤦‍♂️

("atom")
"name O"
("name", "O")
"initial_type O3"
("initial_type", "O3")
"final_type du"
("final_type", "du")
"initial_LJ 3.03981 0.21021"
("initial_LJ", "3.03981", "0.21021")
"final_LJ 0.00000 0.00000"
("final_LJ", "0.00000", "0.00000")
"initial_charge -0.00000"
("initial_charge", "-0.00000")
"final_charge -0.00000"
("final_charge", "-0.00000")
"endatom"
("endatom")
"atom"
("atom")
"name O1"
("name", "O1")
"initial_type O3"
("initial_type", "O3")
"final_type du"
("final_type", "du")
"initial_LJ 0.00000 0.00000"
("initial_LJ", "0.00000", "0.00000")
"final_LJ 0.00000 0.00000"
("final_LJ", "0.00000", "0.00000")
"initial_charge0.00000"
("initial_charge0.00000")
[1]    114032 segmentation fault (core dumped)  somd-freenrg -C somd.cfg -t somd.prm7 -c somd.rst7 -m somd.pert -p CUDA

After fixing this, the lambda = 1 window runs. Not sure if will NaN at a later date, though.

[fjclark]

Sorry about this - it looks like I've messed up editing a few of the pert files so that initial LJ terms remain for O. The pert file does seem to be correct for lambda = 0 though, where the Nans occur.

[lohedges]

No problem. I'll look at that one to see if I can get it to fail.

[lohedges]

No crashes yet, but I can confirm that I see the ligand "jump" after 50 frames. It almost looks like everything is slightly rotated within the box. It's visible for the protein and ions too, but harder to tell when waters are included. Not sure if this is some buffering error, i.e. there is a larger jump in time between frames 49-50, 99-100, ..., than any other pair of frames, or whether something else is going on.

[fjclark]

It seems like this is due to the use of OFF vs GAFF. When I reparameterise the ligand with GAFF2, keeping the initial coordinates identical, I haven't yet observed any crashes. I think the original ligand FF is OFF2 (set up by a collaborator - I have asked for confirmation). Tests run:

• Applying intermolecular restraints to the decoupled ligand (requires my modified version of Sire - happy to supply inputs if needed). OFF system: 2/ 6 simulations crash with Nan. GAFF2 system: 0 / 18 simulations crash. Using constraint=hbonds-notperturbed.
• Setting constraint=allbonds. OFF system: immediately crashes in all cases. GAFF2 system: runs as expected in all cases.

The input systems are here.

This is another example of the issue discussed here, where constraint=allbonds causes OFF systems to crash.

Thanks.

[chryswoods]

Is this still an open issue? I'm doing some spring cleaning ;-)