SOMD Gives Incorrect Results with constraints=none

[fjclark]

Hi,

I've recently been running RBFEs for the TYK2 benchmark set with constraints=none, because we want to correct these results using MM -> MM/ML legs with ANI-2x. The results generated using Perses for the same network are fine:

But the SOMD results are very poor:

The largest outlier from the SOMD ddG plot is a H -> methyl transformation (ejm 46 to jmc 28):

The experimental ddG for this edge is -0.3 kcal / mol. Previous results in the group with OFF1.3 were in very good agreement with this. Without constraints, this edge gave an average ddG of ~ - 10 kcal /mol:

Leg	Run1	Run2	Run3
Bound	1.3	1.1	0.6
Free	11.0	11.0	11.4

But rerunning this with constraints=hbonds (I now realise I should have been using hbonds=hbonds-notperturbed) resulted in an average ddG of ~ -1.0 kcal / mol:

Leg	Run1	Run2	Run3
Bound	1.6	1.8	1.9
Free	2.7	2.8	2.7

The main issue seems to be with the free leg. Looking at the free overlap matrix without constraints, there are issues with poor overlap between lambda = 0.9 and 1.0:

#Overlap matrix 
0.6718 0.2828 0.0431 0.0022 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 
0.2828 0.4370 0.2372 0.0404 0.0025 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000 
0.0431 0.2372 0.4320 0.2374 0.0467 0.0035 0.0001 0.0000 0.0000 0.0000 0.0000 
0.0022 0.0404 0.2374 0.4188 0.2461 0.0501 0.0046 0.0004 0.0001 0.0000 0.0000 
0.0000 0.0025 0.0467 0.2461 0.4034 0.2351 0.0557 0.0078 0.0019 0.0007 0.0000 
0.0000 0.0001 0.0035 0.0501 0.2351 0.3998 0.2469 0.0514 0.0097 0.0034 0.0000 
0.0000 0.0000 0.0001 0.0046 0.0557 0.2469 0.4395 0.2298 0.0188 0.0045 0.0001 
0.0000 0.0000 0.0000 0.0004 0.0078 0.0514 0.2298 0.5658 0.1401 0.0046 0.0001 
0.0000 0.0000 0.0000 0.0001 0.0019 0.0097 0.0188 0.1401 0.7699 0.0593 0.0001 
0.0000 0.0000 0.0000 0.0000 0.0007 0.0034 0.0045 0.0046 0.0593 0.9231 0.0044 
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001 0.0044 0.9953

There is a dramatic free energy change from lambda = 0.8 to lambda = 1.0:

#DG from neighbouring lambda in kcal/mol 
0.0000 0.1000 0.0617 0.0281 
0.1000 0.2000 0.0672 0.0291 
0.2000 0.3000 -0.0325 0.0284 
0.3000 0.4000 -0.0788 0.0263 
0.4000 0.5000 -0.1024 0.0262 
0.5000 0.6000 0.0170 0.0245 
0.6000 0.7000 0.7415 0.0300 
0.7000 0.8000 2.1797 0.0509 
0.8000 0.9000 3.9390 0.0921 
0.9000 1.0000 4.2300 0.3859

Adding constraints fixes the overlap issue:

#Overlap matrix 
0.2070 0.1807 0.1584 0.1386 0.1201 0.0995 0.0626 0.0184 0.0089 0.0043 0.0014 
0.1807 0.1685 0.1557 0.1423 0.1281 0.1100 0.0720 0.0231 0.0118 0.0058 0.0019 
0.1584 0.1557 0.1507 0.1437 0.1345 0.1198 0.0821 0.0289 0.0158 0.0080 0.0025 
0.1386 0.1423 0.1437 0.1427 0.1388 0.1286 0.0928 0.0364 0.0213 0.0111 0.0036 
0.1201 0.1281 0.1345 0.1388 0.1404 0.1357 0.1046 0.0470 0.0297 0.0159 0.0052 
0.0995 0.1100 0.1198 0.1286 0.1357 0.1384 0.1197 0.0673 0.0464 0.0260 0.0086 
0.0626 0.0720 0.0821 0.0928 0.1046 0.1197 0.1433 0.1350 0.1044 0.0622 0.0213 
0.0184 0.0231 0.0289 0.0364 0.0470 0.0673 0.1350 0.2341 0.2140 0.1429 0.0528 
0.0089 0.0118 0.0158 0.0213 0.0297 0.0464 0.1044 0.2140 0.2432 0.2088 0.0956 
0.0043 0.0058 0.0080 0.0111 0.0159 0.0260 0.0622 0.1429 0.2088 0.2835 0.2314 
0.0014 0.0019 0.0025 0.0036 0.0052 0.0086 0.0213 0.0528 0.0956 0.2314 0.5758

and prevents the large jumps in ddG near lambda = 1:

#DG from neighbouring lambda in kcal/mol 
0.0000 0.1000 0.1069 0.0007 
0.1000 0.2000 0.1193 0.0006 
0.2000 0.3000 0.1123 0.0005 
0.3000 0.4000 0.1045 0.0005 
0.4000 0.5000 0.1380 0.0007 
0.5000 0.6000 0.3177 0.0019 
0.6000 0.7000 0.2583 0.0031 
0.7000 0.8000 0.0521 0.0017 
0.8000 0.9000 0.4016 0.0026 
0.9000 1.0000 1.1172 0.0050

Example inputs are here.

System Details

• OS: Ubuntu 22.04
• Version of Python: 3.9.16
• Version of sire: 2023.1.3

Additional context
This is not a pressing problem as we will use Perses for these calculations. Julien asked me to raise this as a reminder for when SOMD is redeveloped in future.

Thanks very much.

[chryswoods]

Thanks - this will be a good test case to validate that the redeveloped SOMD is working.

[chryswoods]

Similarly to #47, I'd like to close this as part of spring-cleaning the issues. It is a really good test case that will aid in the development of somd2. I will create a somd2 label so that we can easily retrieve these issues when we ramp up development.