Positional restraints error in runFreeNrg: property format mismatch between setupPositionalRestraints and atomNumVectorListToProperty

[Roy-Haolin-Du]

Hi Lester, hope you are well and sorry to bother you!
One of our collaborators would like to use position restraints in somd-freenrg.
And I just have a quick question regarding this.

Summary:
runFreeNrg() with position restraints fail because the properties in a different format than expected by the integrator initialization.

runFreeNrg() use setupPositionalRestraints() via PositionalRestraintsToProperty creates properties like: Anchor(0), Atom(0), k(0), nrestrainedatoms.
However, this error showed : AtomNum(0), was looked seems by atomNumVectorListToProperty and propertyToAtomNumVectorList expected, which is not what we created.
So raise an error:

Traceback (most recent call last):
  File "/home/fredpowell/anaconda3/envs/a3fe/share/Sire/scripts/somd-freenrg.py", line 224, in <module>
    OpenMMMD.runFreeNrg(params)
  File "/home/fredpowell/anaconda3/envs/a3fe/lib/python3.12/site-packages/sire/legacy/Tools/__init__.py", line 175, in inner
    retval = func()
             ^^^^^^
  File "/home/fredpowell/anaconda3/envs/a3fe/lib/python3.12/site-packages/sire/legacy/Tools/OpenMMMD.py", line 2834, in runFreeNrg
    moves = setupMovesFreeEnergy(system, debug_seed.val, gpu.val, lambda_val.val)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/fredpowell/anaconda3/envs/a3fe/lib/python3.12/site-packages/sire/legacy/Tools/OpenMMMD.py", line 2158, in setupMovesFreeEnergy
    Integrator_OpenMM.initialise()
KeyError: 'SireBase::missing_property: There is no property with name "AtomNum(0)". Available properties are [ Atom(15), Atom(5), k(8), Anchor(8), k(9), k(2), Anchor(10), k(14), Anchor(4), k(7), Atom(6), Anchor(9), k(16), Atom(2), Anchor(0), Atom(14), Anchor(16), nrestrainedatoms, Atom(11), Anchor(1), Anchor(2), Anchor(12), Anchor(6), k(0), Anchor(13), Atom(1), k(10), Atom(10), Atom(4), Atom(13), k(1), k(11), Atom(7), Atom(8), k(15), k(6), Anchor(5), Atom(16), k(4), Anchor(14), k(13), Atom(12), Anchor(11), Anchor(3), k(3), Anchor(7), k(12), Atom(3), Atom(0), k(5), Anchor(15), Atom(9) ]. (call sire.error.get_last_error_details() for more info)'
srun: error: BarkerGroup: task 0: Exited with exit code 1

The restrain setting in somd.cfg :

restrained atoms = {38387: 454, 38388: 810, 38389: 1942, 38390: 2123, 38391: 2198, 38392: 2287, 38393: 2302, 38394: 2319, 38395: 2326, 38396: 2434, 38397: 2809, 38398: 2860, 38399: 2879, 38400: 2932, 38401: 2951, 38402: 2995, 38403: 3012}
restraint force constant = 2750

Is something related Integrator_OpenMM.initialise() to use propertyToAtomNumVectorList ?
I Try to clear myself:

runFreeNrg()
├── setupPositionalRestraints()     # ✅ create new properties
├── setupForceFieldsFreeEnergy()    # ✅ setup FF
├── setupMovesFreeEnergy()          # ← error 
    └── Integrator_OpenMM.initialise()

Probably C++:

Integrator_OpenMM.initialise()
  ↓ (somehow)
Python: propertyToAtomNumVectorList()
  ↓
Python: prop["AtomNum(0)"]  # Error!

Sorry for my messy explanation. I’m still trying to figure it out, but it’s quite difficult for me. Hope you can help me if you have some time.

Thank you very very much!

Best Wishes,

Roy

[lohedges]

Hi @Roy-Haolin-Du. Are you able to share a somd working directory so that I can debug locally? I tried adding some restraints to a local test file that I use, but it didn't trigger the exception that you show.

[Roy-Haolin-Du]

Hi Lester,

Thanks a lot. But I currently only have these two files. I’ll ask Fred for the rest so we can fully reproduce the issue.

Also, I’d like to double-check my understanding of the restraint setup:
restrained_atoms = {38387: 454}
Does this mean that atom 38387 is being restrained to the position of atom 454?

If my goal is to restrain atom 38387 to its initial position instead, wouldn’t it be more appropriate to set:
restrained_atoms = {38387: 38387}

Just want to make sure I understand this correctly.

somd.cfg.txt
somd-array-gpu-35992.4294967294.out.txt

Thanks a lot for your help!

Best,

[lohedges]

Does this mean that atom 38387 is being restrained to the position of atom 454?

No, the dictionary keys specify the indices of additional dummy atoms that are harmonically bound to real atoms (the values). This functionality was added in this PR to be consistent with the way that position restraints are used by BioSImSpace with OpenMM. This is the recommended approach for position restraints that still works correctly when running simulations in the NPT ensemble.

However, it looks like he code doesn't automatically add the required dummy atoms for you, which is something that is done by BioSImSpace when running with OpenMM. I'll check with @mb2055 and get back to you with the recommended setup.

[lohedges]

The functionality was intended to be combined with the addanchors.py file contained in this archive. Apparently @jmichel80 wrote the code and @mb2055 pushed it, so it's probably worth chatting with @jmichel80 when he's back if things still don't make sense. Essentially the script should take the input form the somd1 working directory, then write back an updated coordinate and topology file containing the additional dummy atoms. You can pass a CSV with the atoms that you want to restrain, or instead the script will find heavy atoms for you, although the logic doesn't look the must reliable.

[Roy-Haolin-Du]

Hi Lester,

Thanks so much for letting me know more about the position restraints, that’s really helpful!
Actually, we do have one working case where somd.cfg uses:

cutoff type = cutoffperiodic  
cutoff distance = 12 * angstrom  
reaction field dielectric = 78.3 

and it runs fine.
However, when we switch to:

cutoff type = PME  
cutoff distance = 10 * angstrom  

it fails for charged ligands.

Do you have any suggestions on what might be causing this, or if there’s something specific we should check when using PME with charged systems? You can reproduce it by using attach simple files and change to PME or RF to see differents.

positionRestrain_PME.zip

Any advice would be much appreciated!

Best regards,

[lohedges]

I can only assume that the required code (what was in the PR linked to above) was never added to the PME version of SOMD.

[Roy-Haolin-Du]

Thanks so much for your kind help, I really appreciate it.
I completely agree with you, especially considering the logic:

if cutoff_type.val == "PME":
    fep_cls = Sire.Move.OpenMMPMEFEP
else:  # no cutoff and RF
    fep_cls = Sire.Move.OpenMMFrEnergyST

We’ll wait for Mathew and Julien to get back so we can understand more about this.

Thanks again for your support!

Have a good afternoon!

[lohedges]

No problem. From what I recall this was a prototype implementation intended to be tested by industrial collaborators, hence why it not fully featured. It probably wasn't intended to be widely used at this stage.