Error with analysis of results in ligandswap tutorial

[Jans Alzate Morales]

Dear Sire developers,

I hope you and your families are doing well and stay healthy.

I am experimenting this error when performing the analysis of result in ligandswap tutorial:

/Users/jalzate/sire.app/lib/python3.7/site-packages/Sire/Tools/ap.py:434: RuntimeWarning: invalid value encountered in double_scalars

  result = int((val - min) / step)

Error plotting graph: cannot convert float NaN to integer# PMFs

Thereby the free energy results are empty for FEP and Bennet methods as you can find the file in this link (https://we.tl/t-G4HUNmKnNB) in WeTransfer.  The previous step, that means monte Carlo calculations of ligands in water and protein, apparently run fine.

Thanks in advance for your support and insight into these first steps with the software.

Best regards and stay safe.

Jans

[Christopher Woods]

Dear Jans,

I've taken a look and the error is because the free energy calculation
has not converged (and is unlikely to converge). FEP has the poorest
convergence property compared to Bennets and TI, and indeed looking at
your results only TI has been able to calculate an actual PMF. If the
free energy doesn't converge then the values are either inf (infinite,
as for your Bennets results) or NaN ("not a number") as is the case
for FEP.

This indicates that ligandswap is not able to work for your system.
This is likely because the two ligands are too different in size, or
the second ligand is overlapping with the protein or water when it is
being swapped. There isn't much you can do to solve this - the method
is limited to swapping very similar ligands. The best you can do is
find a pathway between similarly sized ligands that may link your two
ligands together. However, the cumulative errors over this pathway
will likely make your results poor.

As a guide, the ligandswap free energies should only be +/- 20 kcal
mol-1 maximum. Most values should be in a range of +/- 10 kcal mol-1.
Having values of 80 +/- 160 kcal mol-1, or 320 +/- 250 kcal mol-1 are
red flags that show that the method is not able to work for your
system.

Best wishes,

Christopher

> --
> You received this message because you are subscribed to the Google Groups "Sire Users" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sire-users+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sire-users/38b83530-093b-4ca4-95ef-2f343ef76afd%40googlegroups.com.



--
---------------------------------------------------------
Christopher Woods
+44 (0) 7786 264562
http://chryswoods.com

[Javier García Marín]

Dear Jans,

I am not an expert using ligandswap but I got that kind of error when my simulations

didn't converge. Normally I can fix this problem by using a more equilibrated starting structure

because sometimes the ligand overlaps with protein residues (mostly in docking results with soft potentials).

Best whises

Javier

[Jans Alzate Morales]

Dear Christopher,

Thanks for the quick and detailed response to my questions.

I changed the config file and put bigger numbers in the parameters, something like this: 

nmoves = 10000

nsubmoves = 500000

nequilmoves = 500000 

and now I have free binding energies converged for 3 methods (TI, Bennett and FEP). However, my energies are still 10 kcal/mol higher than those reported in the ligandswap tutorial (around 4.7 kcal/mol).

So, what parameter is better to change in the config file to achieve the millions of evaluations to get those free binding energy values?

On the other hand, I just installed Sire in a 64 core server with cpu (Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz) and gpus (GeForce RTX 2080), but I am experimenting the following error when trying to run the tutorial files of ligandswap:

"[jalzate@c108720 ligandswap]$ $SIRE/bin/ligandswap -t0 rec_fmc.top -c0 fmc.30.crd -l0 FM1 -t1 rec_cti.top -c1 cti.30.crd -l1 CTI -C config

Starting /home/jalzate/sire.app/bin/ligandswap: number of threads equals 64

Loading configuration information from file config

Ligand 0 will be located by finding the first molecule containing residue FM1

Ligand 1 will be located by finding the first molecule containing residue CTI

Running a ligandswap calculation calculating the difference in free energy between

ligands 0 and 1 using files rec_fmc.top|fmc.30.crd and rec_cti.top|cti.30.crd.

Using parameters:

===============

crdfile0 == fmc.30.crd

crdfile1 == cti.30.crd

ligand0 == FM1

ligand1 == CTI

nequilmoves == 100

nmoves == 5

nsubmoves == 10

topfile0 == rec_fmc.top

topfile1 == rec_cti.top

===============

==============================================================

Number of iterations to perform: 5. Number of iterations completed: 5.

Sending anonymous Sire usage statistics to http://siremol.org.

All iterations complete. Simulation complete.

For more information, see http://siremol.org/analytics

To disable, set the environment variable 'SIRE_DONT_PHONEHOME' to 1

Fatal Python error: could not acquire lock for <_io.BufferedWriter name='<stdout>'> at interpreter shutdown, possibly due to daemon threads

Thread 0x00007ff2064aa700 (most recent call first):

  File "/home/jalzate/sire.app/lib/python3.7/site-packages/Sire/__init__.py", line 252 in _uploadUsageData

  File "/home/jalzate/sire.app/lib/python3.7/threading.py", line 865 in run

  File "/home/jalzate/sire.app/lib/python3.7/threading.py", line 917 in _bootstrap_inner

  File "/home/jalzate/sire.app/lib/python3.7/threading.py", line 885 in _bootstrap

Current thread 0x00007ff2aae68780 (most recent call first):

Aborted (core dumped)"

Thanks in advance for your help and support and stay safe.

Jans

> To unsubscribe from this group and stop receiving emails from it, send an email to sire-...@googlegroups.com.

[Jans Alzate Morales]

Dear Javier,

Thanks a lot for your quick response.

I actually change the values for parameters in the config file, and I got free binding energy values more closer (10 kcal/mol higher) to those reported in the ligandswap tutorial.

However, I would like to understand which of those parameters is better to change in order to get the correct free binding energy values.

Any hint is welcomed!

Best regards and stay safe.

Jans

[Christopher Woods]

  Dear Jans,

Sorry that I had not noticed that you were working from the tutorial. This should have converged quickly and without you significantly increasing the length of the simulation (which you are right is one way to drive convergence). This points to a bigger issue that there is likely a bug that has been introduced into the ligandswap code. Ligandswap is rarely used and now better tools for calculating relative binding free energies are available (gromacs, somd etc) such that ligandswap has not been looked after for several years now.

If you still want to use it then I am afraid we don’t have spare resource in the next few weeks available to debug this issue. The error would have been introduced in 2019, so please try one of the 2018 versions of Sire to see if the problem persists. All of the old releases of Sire can be downloaded as binaries from this page - <https://siremol.org/pages/binaries.html>, or you can download a source release from github and compile.

There have been no structural or bug fix changes to Sire in years, so the older versions are fine to use. The bug was likely introduced as an incompatibility with the new file parsing code that was built for BioSimSpace.

  Best wishes,

  Christopher 

--
You received this message because you are subscribed to the Google Groups "Sire Users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to sire-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sire-users/99e796fe-be8a-4a17-99a9-72ad25381622%40googlegroups.com.

--

Sent from Gmail Mobile

[Jans Alzate Morales]

Dear Christopher,

Thanks again for your support.

I will give it a try with an older version of SireMol to see if it works. 

Regarding the other problem with CPUs, I could work it around using the option --ppn=16 (instead of 64 cores) and now I am testing the tutorial in the server.

On the other hand, I am just trying to complete another tutorial about QM/MM, but I am finding this error with the script energy.py:

$SIRE/bin/python energy.py TRP PM3

Traceback (most recent call last):

  File "energy.py", line 139, in <module>

    qmmm_energy = qmmmff.energy()

Boost.Python.ArgumentError: Python argument types in

    QMMMFF.energy(QMMMFF)

did not match C++ signature:

    energy(Squire::QMMMFF {lvalue}, SireFF::EnergyTable {lvalue} energytable, SireCAS::Symbol symbol, double scale_energy=1)

    energy(Squire::QMMMFF {lvalue}, SireFF::EnergyTable {lvalue} energytable, double scale_energy=1)

The version of Sire is 2019, and maybe this is also related to the bug you mentioned.  The other scripts in the tutorial (sample.py and lamba.py) are working properly.

Best regards, stay safe, and thank you in advance for your help.

Jans

To unsubscribe from this group and stop receiving emails from it, send an email to sire-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/sire-users/99e796fe-be8a-4a17-99a9-72ad25381622%40googlegroups.com.

[Christopher Woods]

  Dear Jans,

This is an error caused by the API of the energy function changing in about 2017. The QM/MM tutorial pre-dates that change and has not been updated, and so had become victim to bit rot.

The QM/MM tutorial is designed to only give an idea behind the theory of how it was implemented, and not for actual production use. If you want to run QM/MM free energy simulations then you should use quantomm directly (look at the —help for options). Again, this is now quite old (2012) and so better and newer methods and code may be available and preferred.

Mostly Sire today is used for GPU MD free energy calculations via somd, a little use for absolute binding via waterswap, and mostly as the backend to BioSimSpace for converting file formats, building systems and creating re-usable simulation workflow components.

  Best wishes,

  Christopher 

To unsubscribe from this group and stop receiving emails from it, send an email to sire-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/sire-users/e2240414-7306-4037-b4fe-70e2c53bc529%40googlegroups.com.