TI with gaussian quadrature

[hl2500]

Added an optional gaussian quadrature functionality to estimate the TI integral.

This functionality could be triggered when initializing: TI(gauss_qua=True)

The lambda values and weights were referred the Amber manual (page 493): http://ambermd.org/doc12/Amber22.pdf

The formula used to estimate the free energy and standard errors were referred the the Amber manual (eq. 25.2) and this literature (eq. 3): https://pubs.acs.org/doi/10.1021/acs.jcim.2c01052?ref=pdf, respectively.

The final two dataframes with free energy values and errors are upper triangluar matrix, because the values are cumulative sum upto each lambda state rather than differences. The diagonal values are the estimation at each lambda.

I have tested 12 lambdas and 9 lambdas with Amber one-step method to estimate the free energy and error of each leg in the thermodynamic cycle. The results are comparable to the ones estimated with trapezoidal rule. The other number of lambdas (1, 2, 3, 5, 7) should work as well.

[codecov[bot]]

Codecov Report

Merging #304 (13764cb) into master (f5cf43a) will increase coverage by 0.05%. The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master     #304      +/-   ##
==========================================
+ Coverage   98.76%   98.81%   +0.05%     
==========================================
  Files          27       28       +1     
  Lines        1778     1860      +82     
  Branches      391      402      +11     
==========================================
+ Hits         1756     1838      +82     
  Misses          2        2              
  Partials       20       20              

Impacted Files	Coverage Δ
src/alchemlyb/estimators/__init__.py	100.00% <100.00%> (ø)
...rc/alchemlyb/estimators/ti_gaussian_quadrature_.py	100.00% <100.00%> (ø)

[orbeckst]

Sorry, I don't have much time this/next week to review. I hope @xiki-tempula can have a look. I agree with the suggestion to make it a separate estimator instead of changing standard TI.

[hl2500]

Thank you for the comments @orbeckst @xiki-tempula I may need some help to add tests to cover all code. I didn't find a testset in alchemtest for Amber gaussian quadrature (with suggested lambda values). Do I also need to upload an example to alchemtest for the purpose of testing and adding an example on TI documentation page? Thanks!

[hl2500]

Thank you so much for the suggestions. I have added two testsets to alchemtest (PR: Add gaussian quadrature test) and made some changes. @orbeckst @xiki-tempula

[xiki-tempula]

@hl2500 Thanks for the new test files updated in alchemtest. Do you mind test them in this PR? So I could know what are the test for and how they have been used.

[hl2500]

@xiki-tempula Okay. Could I ask how to implement test in this PR? Thanks.

[xiki-tempula]

I think if you just change this line https://github.com/alchemistry/alchemlyb/blob/master/.github/workflows/ci.yaml#L58 to python -m pip install git+https://github.com/hl2500/alchemtest.git@Add_gaussian_quadrature_test You should be able to get this in the CI.

[hl2500]

@xiki-tempula I have changed this line and some of the ubuntu tests failed. Could I ask what should I do? Thanks.

[xiki-tempula]

Just Codecov acting weird. It is fine now.

[hl2500]

@xiki-tempula Thanks, could I ask what other tests I need to do?

[xiki-tempula]

@hl2500 Thanks for the contribution, so the test should be in a format that is similar to https://github.com/alchemistry/alchemlyb/blob/master/src/alchemlyb/tests/test_ti_estimators.py, which is for the TI estimator. Which shall test if the result is numerically good. The other thing that you need to test is if the test is covering all of the code, you could check the coverage here https://github.com/alchemistry/alchemlyb/pull/304/checks?check_run_id=12878246768, we would like to pump it as high as possible.

[hl2500]

@xiki-tempula Thank you for the information! Some tests have been added to improve the coverage.

[hl2500]

@xiki-tempula Thank you so much for the suggestions! Please feel free to let me know if any other changes are needed.

[xiki-tempula]

@hl2500 Sorry for this late comment. I was playing with this estimator on some datasets and found a problem. The current layout for TI, BAR, MBAR is that for deltaf of three lambda [0, 0.5, 1.0]. The [0,0] is zero, and the [0,1] is the difference between the zeroth element in the lambda list [0] and the first element in the lambda list [0.5], the [0,2] is the difference between the zeroth element in the lambda list [0] and the first element in the lambda list [1.0].

Do you mind making this estimator follow this convention as well? Thank you.

I think I'm quite happy with the state of the code now, I will do some further tests to see if it is behaving as I would like it to be and we are good to go.

[xiki-tempula]

@orbeckst Do you mind having a look at this PR as well? I think the code part is good and we just need to polish the doc a bit.

[hl2500]

Thanks very much for all the effort. I think I'm quite happy with this PR now. However, apologise that we cannot merge this in now as it is a significant and important addition. So I want to make sure that no one would have strong opposition to this PR. I purpose that we leave the PR here for one month and if no one is feeling strongly about it, I will merge this in (aim for 2.2.0). What do you think @orbeckst ?

Okay, thank you so much for your guidance and time into this.

[hl2500]

@xiki-tempula This part has been removed. Thank you!

[xiki-tempula]

I have merge this PR thanks for all the contribution and the patience.

[hl2500]

@xiki-tempula Thank you again for your guidance in this work!

[xiki-tempula]

@hl2500 No worries. Looking forward to further contribution from you.