Wrong intramolecular energy in pose rescoring

[anaome]

Hi !

There seems to be an issue with the calculation of the intramolecular energy for pose rescoring. The term (4) is always different from (2). This discrepancy arises from the terms involving flexible side chains. The system consists in a simple ligand (2 torsions) and a receptor with 15 flexible side chains.

flex = str(sys.argv[1])
ligand = str(sys.argv[2])
maps = str(sys.argv[3])
coords = np.array(str(sys.argv[4]).split(',')).astype(float)

v = Vina(sf_name='vina', cpu=1, verbosity=2, seed=1)
v.set_receptor(maps+'.pdbqt',flex)
v.set_ligand_from_file(ligand)
v.compute_vina_maps(center=coords, box_size=[25, 25, 28])
energy = v.score()
print(energy)

Estimated Free Energy of Binding : -3.583 (kcal/mol) [=(1)+(2)+(3)-(4)]
(1) Final Intermolecular Energy : -2.656 (kcal/mol)
Ligand - Receptor : -0.961 (kcal/mol)
Ligand - Flex side chains : -1.695 (kcal/mol)
(2) Final Total Internal Energy : -20.261 (kcal/mol)
Ligand : 0.000 (kcal/mol)
Flex - Receptor : -12.696 (kcal/mol)
Flex - Flex side chains : -7.565 (kcal/mol)
(3) Torsional Free Energy : 0.314 (kcal/mol)
(4) Unbound System's Energy : -19.019 (kcal/mol)

[diogomart]

Thanks for reporting. I can reproduce. The values should be the same. The discrepancy is caused by the use of grid maps to compute the unbound energy, and direct interactions (no grid maps) to compute the internal energy. If you use the --no_refine flag from command line or v = Vina(no_refine=True) in Python, then the internal energy will be calculated using grid maps and it will be the same as unbound. Alternatively, using finer grids, e.g. v.compute_vina_maps( ..., spacing=0.2) reduces the discrepancy. In any case, we will fix it to have the unbound energy use direct interactions when no_refine=False.

[diogomart]

@anaome could you post the energy decomposition with a spacing of 0.2 for you system? I'm just curious how big of a difference it would make for your 15-sidechain system.

v.compute_vina_maps(center=coords, box_size=[25, 25, 28], spacing=0.2)

the map calculation will take ~7x longer

[anaome]

Thanks for the answer !
Here is the data

- default
Estimated Free Energy of Binding : -3.583 (kcal/mol) [=(1)+(2)+(3)-(4)]
(1) Final Intermolecular Energy : -2.656 (kcal/mol)
Ligand - Receptor : -0.961 (kcal/mol)
Ligand - Flex side chains : -1.695 (kcal/mol)
(2) Final Total Internal Energy : -20.261 (kcal/mol)
Ligand : 0.000 (kcal/mol)
Flex - Receptor : -12.696 (kcal/mol)
Flex - Flex side chains : -7.565 (kcal/mol)
(3) Torsional Free Energy : 0.314 (kcal/mol)
(4) Unbound System's Energy : -19.019 (kcal/mol)

- spacing=0.2
Estimated Free Energy of Binding : -2.802 (kcal/mol) [=(1)+(2)+(3)-(4)]
(4) Unbound System's Energy : -19.868 (kcal/mol)

- spacing=0.1
Estimated Free Energy of Binding : -2.546 (kcal/mol) [=(1)+(2)+(3)-(4)]
(4) Unbound System's Energy : -20.148 (kcal/mol)

- no_refine=True
Estimated Free Energy of Binding : -2.395 (kcal/mol) [=(1)+(2)+(3)-(4)]
(1) Final Intermolecular Energy : -2.605 (kcal/mol)
Ligand - Receptor : -0.910 (kcal/mol)
Ligand - Flex side chains : -1.695 (kcal/mol)
(2) Final Total Internal Energy : -19.019 (kcal/mol)
Ligand : 0.000 (kcal/mol)
Flex - Receptor : -11.454 (kcal/mol)
Flex - Flex side chains : -7.565 (kcal/mol)
(3) Torsional Free Energy : 0.210 (kcal/mol)
(4) Unbound System's Energy : -19.019 (kcal/mol)

[diogomart]

Thank you!