Implement iterative average structure (Issue#2039)

[leonwehrhan]

Fixes #2039

Changes made in this Pull Request:

• Implemented average structures with iterative algorithm from DOI 10.1021/acs.jpcb.7b11988. (Issue Iterative Average structures #2039)

PR Checklist

• Tests?
• Docs?
• CHANGELOG updated?
• Issue raised/referenced?

Developers certificate of origin

• I certify that this contribution is covered by the LGPLv2.1+ license as defined in our LICENSE and adheres to the Developer Certificate of Origin.

---

📚 Documentation preview 📚: https://mdanalysis--4524.org.readthedocs.build/en/4524/

[pep8speaks]

Hello @leonwehrhan! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2024-03-29 18:38:13 UTC

[github-actions]

Hello there first time contributor! Welcome to the MDAnalysis community! We ask that all contributors abide by our Code of Conduct and that first time contributors introduce themselves on GitHub Discussions so we can get to know you. You can learn more about participating here. Please also add yourself to package/AUTHORS as part of this PR.

[github-actions]

Linter Bot Results:

Hi @leonwehrhan! Thanks for making this PR. We linted your code and found the following:

Some issues were found with the formatting of your code.

Code Location	Outcome
main package	⚠️ Possible failure
testsuite	⚠️ Possible failure

Please have a look at the darker-main-code and darker-test-code steps here for more details: https://github.com/MDAnalysis/mdanalysis/actions/runs/8484369067/job/23247178106

---

Please note: The black linter is purely informational, you can safely ignore these outcomes if there are no flake8 failures!

[yuxuanzhuang]

Hello @leonwehrhan! Welcome to MDAnalysis! Please add #2039 directly after "Fixes" so the issue will be successfully closed after this PR is merged.

Feel free to ping us when you feel it's ready for review :)

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.47%. Comparing base (dfcd06f) to head (0499f98).
Report is 1 commits behind head on develop.

Additional details and impacted files

@@             Coverage Diff             @@
##           develop    #4524      +/-   ##
===========================================
- Coverage    93.65%   93.47%   -0.18%     
===========================================
  Files          168      180      +12     
  Lines        21216    22319    +1103     
  Branches      3908     3913       +5     
===========================================
+ Hits         19870    20863     +993     
- Misses         888      996     +108     
- Partials       458      460       +2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[leonwehrhan]

Hello @yuxuanzhuang, I think the PR is ready for a review.

Could you maybe give some guidance on how to resolve the failed automatic checks? I am not sure how to tackle the issue with the mypy linter, as my local mypy does not see any issues with the files I have touched. Also, my tests should cover every line of new code (codecov/patch 100%) so I do not know why the codecov/project test failed.

Thank you very much!

[yuxuanzhuang]

Hi, there are a few aspects up for discussion (also pinging @lilyminium, since you added the AverageStructure):

1. Should it be a method or an Analysis class? It also makes sense to include it as an option to run AverageStructure (e.g., iterative=True). For now, I think it's acceptable to leave it as is.

2. In scenarios where you have a large system and want to obtain the average structure of, for example, a protein, the current implementation does not offer the option to do subselections, which AverageStructure can.

3. The code you provided doesn't actually accomplish the task. If you check for verbosity, you can see that the RMSD is not changing at all. This is because ref.positions = res.positions does not actually modify the position of atoms in the reference universe. What you can do is use ref = res.universe, which would represent the average structure from the last round of iteration. Please double-check the output of your code.

4. The test also fails to achieve its goal (except for checking for regression). One problem is that u1 = mda.Merge(u.select_atoms("bynum 1:10")) turns the universe into a single-frame trajectory (check u1.trajectory).

Don't worry about the failing of linter and main_tests (ubuntu-latest, 3.12, true, true), they are unrelated issues.

[leonwehrhan]

Hi @yuxuanzhuang,

thank you very much for the review! I have updated the code and adressed your concerns as follows:

1. While I also think an implementation as part of the AverageStructure class makes sense, I have left the method as it is. I could open another issue and work on it once this PR is merged?
2. I have implemented the optional select keyword. As the reference structure is updated throughout the iterative calculation, the selection for the reference is set in the beginning, by transforming the reference universe into a one-frame universe of the selected atoms of the current frame of the reference. The selection for the mobile universe is passed to the AverageStructure class in each iteration.
3. I have implented the code you suggested. I checked the output using verbose=True to see the RMSD is changing. I ran iterative_average on u=mda.Universe(PSF, DCD) and visually checked that the resulting structure looked reasonable.
4. Given the implementation of the select keyword, it is now not necessary to use an atom-sliced universe as input for the iterative workflow, so the original universe u is now returned and processed in the tests. I included some regression tests and some tests checking that the output of the first 10 atoms is the same that I have visually checked.

Other things I have changed:

• The weights keyword could not be actually used before, because calling rms.rmsd would produce an error. This is fixed by calling get_weights and is included in the tests.
• The default eps=1e-8 was not really useful to get to a converged structure (even after 1000 iterations), so I set it to 1e-6.
• Updated documentation.

[yuxuanzhuang]

Looks great! Just a few comments on the documentation side.

And please fix pep8 issue reported above.

[yuxuanzhuang]

I think this function should actually be put inside the analysis.align module since what it does is alignment. Could you move the code (and tests) there?

The other issue is the documentation is not showing up on the webpage (there's a preview you can find above https://mdanalysis--4524.org.readthedocs.build/en/4524/documentation_pages/analysis/rms.html). You need to add explicitly .. autofunction:: xxx to the documentation at the beginning of the file.

Additionally, please leave the review unresolved so that the reviewer can determine whether the issue has indeed been addressed.

[leonwehrhan]

I have moved the function to analysis.align and added the remark to the documentation in the beginning of the file. The documentation of iterative_average now shows up on https://mdanalysis--4524.org.readthedocs.build/en/4524/documentation_pages/analysis/align.html.

[orbeckst]

@yuxuanzhuang I added you as the person responsible for the PR.

@leonwehrhan I un-assigned you because we use the "Assignee" field to indicate who from the developer team manages review and merge.

[yuxuanzhuang]

lgtm!

[orbeckst]

Looking good but we need a convergence check and we need to handle text output via logging.

[leonwehrhan]

@yuxuanzhuang @orbeckst I have made the requested changes.

[orbeckst]

Thanks for making the changes.

My main comment is for the tests: please split your tests so that you always just test one thing as this will help tremendously in debugging if anything ever fails. Please see the comments inline.

[leonwehrhan]

@orbeckst I have changed the tests, so that every testing method only contains one call of iterative average and changed the documentation as you suggested.

[orbeckst]

Thanks for changing the tests. Just found two small things. Almost there!!

[leonwehrhan]

@orbeckst @yuxuanzhuang the requested changes are pushed, but it looks like auto-merge has been disabled by my commit because I do not have write access. Could one of you do the merge?

[orbeckst]

Auto-merge just means that the merge will be performed when all requirements are met

• all reviews are approving
• all required status checks passing

It's just a short cut for maintainers to essentially "queue up" the merge.

Let me have a quick look.

[orbeckst]

Please check the imports that I flagged.

In order to remove them you can simply use my suggestions in the GH interface.

[leonwehrhan]

Thanks for the clarification! I thought I might have unintentionally halted the merge process :)

I will adress your changes now.

[leonwehrhan]

@orbeckst I have comitted the changes via your comment. The reason for these imports being there was that I originally implemented the function in analysis.rms instead of analysis.align and later moved the code (see comments above).

[orbeckst]

ok. the other test failure is not due to your stuff, see #4540

[orbeckst]

Thank you for your contribution! 🎉

[leonwehrhan]

Thanks a lot for guiding me through the PR!

[yuxuanzhuang]

Congrats! @leonwehrhan