TypeError: 'restraint' must be of type 'BioSimSpace.FreeEnergy.Restraint' for free leg of simulation.

[yasirkhanqu]

Hi,

I am working with a system with a protein, an acceptor substrate of (three) carbohydrate subunits parameterised with GLYCAM, MG2+ ion, and a donor substrate (LIG), parameterised with gaff forcefield.

Everything works fine for the bound leg until this error erupts (which is for the free leg as I can see from the error). I have attached the input files and the overall output.log file. I have also attached the restraint configuration file from bound/restraint/input/somd.cfg. Additionally, I have attached the param files and the restrain folder from bound/restrain.

Scoring candidate Boresch anchor points. Anchor set no:  96%|█████████▌| 48/50 [24:44<00:59, 29.84s/it]
Scoring candidate Boresch anchor points. Anchor set no: 100%|██████████| 50/50 [25:46<00:00, 30.92s/it]
INFO - 2025-02-18 02:01:45,283 - Leg (type = BOUND)_306 - Writing input files for BOUND leg RESTRAIN stage
INFO - 2025-02-18 02:01:45,307 - Leg (type = BOUND)_306 - Perturbation type: restraint
╭───────────────────── Traceback (most recent call last) ──────────────────────╮
│ /home/yasir/a3fe/glft2_param/abfe.py:3 in <module>                           │
│                                                                              │
│    1 import a3fe as a3                                                       │
│    2 calc = a3.Calculation(ensemble_size = 5)                                │
│ ❱  3 calc.setup()                                                            │
│    4 # Get optimised lambda schedule with thermodynamic speed                │
│    5 # of 2 kcal mol-1                                                       │
│    6 calc.get_optimal_lam_vals(delta_er = 2)                                 │
│                                                                              │
│ /home/yasir/.conda/envs/a3fe/lib/python3.12/site-packages/a3fe/run/calculati │
│ on.py:211 in setup                                                           │
│                                                                              │
│   208 │   │   │   │   stream_log_level=self.stream_log_level,                │
│   209 │   │   │   )                                                          │
│   210 │   │   │   self.legs.append(leg)                                      │
│ ❱ 211 │   │   │   leg.setup(configs[leg_type])                               │
│   212 │   │                                                                  │
│   213 │   │   # Save the state                                               │
│   214 │   │   self.setup_complete = True                                     │
│                                                                              │
│ /home/yasir/.conda/envs/a3fe/lib/python3.12/site-packages/a3fe/run/leg.py:24 │
│ 2 in setup                                                                   │
│                                                                              │
│    239 │   │   │   system = self.run_ensemble_equilibration(sysprep_config=c │
│    240 │   │                                                                 │
│    241 │   │   # Write input files                                           │
│ ❱  242 │   │   self.write_input_files(system, config=cfg)                    │
│    243 │   │                                                                 │
│    244 │   │   # Make sure the stored restraints reflect the restraints used │
│    245 │   │   # make this more robust my using the SOMD functionality to ex │
│                                                                              │
│ /home/yasir/.conda/envs/a3fe/lib/python3.12/site-packages/a3fe/run/leg.py:76 │
│ 2 in write_input_files                                                       │
│                                                                              │
│    759 │   │   │   )                                                         │
│    760 │   │   │   self._logger.info(f"Perturbation type: {stage_type.bss_pe │
│    761 │   │   │   # Ensure we remove the velocites to avoid RST7 file writi │
│ ❱  762 │   │   │   _BSS.FreeEnergy.AlchemicalFreeEnergy(                     │
│    763 │   │   │   │   pre_equilibrated_system,                              │
│    764 │   │   │   │   protocol,                                             │
│    765 │   │   │   │   engine="SOMD",                                        │
│                                                                              │
│ /home/yasir/.conda/envs/a3fe/lib/python3.12/site-packages/BioSimSpace/Sandpi │
│ t/Exscientia/FreeEnergy/_alchemical_free_energy.py:361 in __init__           │
│                                                                              │
│    358 │   │   # For free leg simulations, the restraint will be None.       │
│    359 │   │   if restraint is not None:                                     │
│    360 │   │   │   if not isinstance(restraint, _Restraint):                 │
│ ❱  361 │   │   │   │   raise TypeError(                                      │
│    362 │   │   │   │   │   "'restraint' must be of type 'BioSimSpace.FreeEne │
│    363 │   │   │   │   )                                                     │
│    364 │   │   │   else:                                                     │
╰──────────────────────────────────────────────────────────────────────────────╯
TypeError: 'restraint' must be of type 'BioSimSpace.FreeEnergy.Restraint'.

Here is bound/restraint/input/somd.cfg file.

### For information on options and defaults, run `somd-freenrg --help-config`

### Integrator - ncycles modified as required by a3fe ###
nmoves = 25000
ncycles = 60
timestep = 4 * femtosecond
constraint = hbonds
hydrogen mass repartitioning factor = 3.0
integrator = langevinmiddle
inverse friction = 1 * picosecond
temperature = 25 * celsius
# Thermostatting already handled by langevin integrator
thermostat = False

### Barostat ###
barostat = True
pressure = 1 * atm

### Non-Bonded Interactions ###
cutoff type = PME
cutoff distance = 10 * angstrom
reaction field dielectric = 78.3

### Trajectory ###
buffered coordinates frequency = 5000
center solute = True

### Minimisation ###
minimise = True

### Alchemistry - restraints added by a3fe ###
perturbed residue number = 1
energy frequency = 200

perturbed_residue number = 1
use boresch restraints = True
turn on receptor-ligand restraints mode = True
charge difference = 0
lambda array = 0.0, 0.125, 0.25, 0.375, 0.5, 1.0

param_log_restrain.zip

[Roy-Haolin-Du]

Hi, thanks so much for your inofrmation.

I guess you resubmitted the run.
But in your first run, you had already created discharge, vanish, and restrain directories under bound. When rerun, it would load again, it detects these existing files and doesn’t recreate them. However, the boresch restrain file must be in a specific format, so it throws an error.

So what you need to do before rerunning is delete Calculation.log and Calculation.pkl. Also, delete bound/discharge, vanish, restrain, Leg.log, Leg.pkl, and virtual_quee.log and free/discharge, vanish, Leg.log, Leg.pkl, and virtual_quee.log" Only Only keep bound/ensemble_quil_12345, free/ensemble_quil_12345, and the input folder.😊

Give it a try and let me know if the progress is going well.

Cheers,

[Roy-Haolin-Du]

By the way, for you somd.cfg

### Non-Bonded Interactions ###
cutoff type = PME
cutoff distance = 10 * angstrom
reaction field dielectric = 78.3

when using PME, it's best to remove reaction field dielectric = 78.3 under the ### Non-Bonded Interactions ### section.

Thanks~

[yasirkhanqu]

Hi,

I have removed reaction field dielectric = 78.3 and run the calculations.

The results are here, and it seems like the restrain and vanish stages are fine, but I have problems in the discharge state; can you please check and let me know what could possibly be wrong?

The files I used as input are the same as above.

Here is the comparison of results for MIF from the paper and my target.

MIFvsGlfT2.pdf

[Roy-Haolin-Du]

Hi,

I have removed reaction field dielectric = 78.3 and run the calculations.

The results are here, and it seems like the restrain and vanish stages are fine, but I have problems in the discharge state; can you please check and let me know what could possibly be wrong?

The files I used as input are the same as above.

Here is the comparison of results for MIF from the paper and my target.

MIFvsGlfT2.pdf

Thank you for the information and very great that you got the results!

From your discharge plot, it seems like the issue might be related to somd.pert. Do you use a specific somd.pert file for different stages? Since the discharge somd.pert file controls the electrostatic change, it could be the source of the problem.

Would you mind having a quick look? Thank you!

[fjclark]

Thanks for the details @yasirkhanqu and thanks for looking into this @Roy-Haolin-Du.

Your PMFs look reasonable to me! Are your final dGs reasonable? The reason your discharge PMFs look very different to ours for MIF is that we used reaction field electrostatics whereas you're using PME. While the PMFs look very different, the difference between the bound and free discharge dGs should be very similar with PME or RF (with a neutral ligand).

Also, while it's a good idea to remove reaction field dielectric = 78 for clarity, as Roy says, it shouldn't affect the calculation as it will be ignored when you set cutofd type = PME.

[yasirkhanqu]

The problem is that my energies are positive:

###################################### Free Energies ########################################
Mean free energy: 11.979 + /- 2.930 kcal/mol
Free energy from run 1: 14.505 +/- 0.133 kcal/mol
Free energy from run 2: 13.935 +/- 0.133 kcal/mol
Free energy from run 3: 8.984 +/- 0.132 kcal/mol
Free energy from run 4: 10.239 +/- 0.133 kcal/mol
Free energy from run 5: 12.233 +/- 0.134 kcal/mol
Errors are 95 % C.I.s based on the assumption of a Gaussian distribution of free energies

Also, I used the same somd.pert file, as generated automatically by the workflow. I understand the higher energies for the discharge leg as the active site has a Mn2+ ion. I don't know why the overall dGs are positive. The ligand is charged with -2e, but i have neutralised the whole system.

[fjclark]

Ok, thanks for the details. Could you please share all inputs/ outputs to make it easier for us to debug? The easiest way to do this is likely:

• Copy the entire calculation to another directory
• Load the calculation and delete all large files with:

import a3fe as a3
calc = a3.Calculation()
calc.lighten()

• Compress and share the lightened calculation here

We haven't run any calculations with a -2 (or +2) charge - this is a challenging case but I would hope that we could do better than > 10 kcal mol-1 estimate, especially as the variance isn't crazy. Thanks.

[yasirkhanqu]

Thank you so much, @fjclark; I have done as you suggested. The size of the directory is still reasonably high, and it fails to bind here, so I uploaded it to Google Drive.

Here is the link: https://drive.google.com/file/d/1gC0VcwNN_rni-QFjEMrJK1pk3uAiWe_c/view?usp=drive_link

Please send a request, and I will grant you permissions.

Looking forward to seeing how it works out.

[fjclark]

Thanks @yasirkhanqu - I've asked for access.

I'll take a look when I can, but I'm busy and might not get to it this week. I'm not sure if @Roy-Haolin-Du has time? Either way, we'll get to it as soon as we cann

[Roy-Haolin-Du]

Thanks @yasirkhanqu I’ve also requested access and just got approved. Apologies for the late reply and I’ll review everything in the next couple of days.

[Roy-Haolin-Du]

Hi @yasirkhanqu , hope you are well.

For the unexpected data, I noticed that in your bound system, the calculation was simulated without using Using co-alchemical ion approach, whereas in the free system, it was used.

You can see it in the leg.log file:
free/Leg.log
INFO - 2025-02-23 10:00:15,757 - Leg (type = FREE)_305 - Ligand has charge -2. Using co-alchemical ion approach to maintain neutrality.
But no this sentence in bound/Leg.log.

Then I use to check your ligand charge:

import BioSimSpace as BSS
system = BSS.IO.readMolecules(["bound_preequil.prm7", "bound_preequil.rst7"])
lig = system.search("resname LIG").molecules()[0]
lig_charge = round(lig.charge().value())
print(lig_charge)
0

system_lig = BSS.IO.readMolecules(["free_preequil.prm7", "free_preequil.rst7"])
lig = system_lig.search("resname LIG").molecules()[0]
lig_charge = round(lig.charge().value())
print(lig_charge)
-2

So for your bound system, when you parameterised your complex system, the charged ligand lost its charge, which caused 0 charge recognized.
I dont know how you parameterised your system, maybe try to solvate within a3fe/biosimspace and see if this problem still exists.

Thanks!

[yasirkhanqu]

Hi, @Roy-Haolin-Du. Thank you. I'm OK, hope you are too.

I have reprepared the systems. Here is how it looks right now

`In [2]: system = BSS.IO.readMolecules(["bound_param.prm7", "bound_param.rst7"])
...: lig = system.search("resname LIG").molecules()[0]
...: lig_charge = round(lig.charge().value())
...: print(lig_charge)
-2

In [3]: system_lig = BSS.IO.readMolecules(["free_param.prm7", "free_param.rst7"])
...: lig = system_lig.search("resname LIG").molecules()[0]
...: lig_charge = round(lig.charge().value())
...: print(lig_charge)
-2`

I will let you know about the outcomes...

btw, I am doing parameterization in LEAP. I have tried several times adding ions and solvents too (in leap), but then a3fe would not recognise the box with an error something like there is no box found or something related... I tried using biosimspace, too, but couldn't get it right. I would be glad if you would let me know. It would help me add a suitable amount of salt as well.

[Roy-Haolin-Du]

Good!, looking forward your new system result.

For the box issue, OpenMM (like some other engines) has additional constraints on the orientation of the lattice vectors.
If you run this would get right box.

import BioSImSpace as BSS

system = BSS.IO.readMolecules("*7")
system.rotateBoxVectors()
system.reduceBoxVectors()

BSS.IO.saveMolecules("somd", system, "rst7")

If you use a3fe/biosimspace, the box should be fine and could be recognised well.

There are two related issues, you can have a look to give more detailed explaination.
OpenBioSim/biosimspace#238 (comment).
OpenBioSim/sire#285.

Thanks.