Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unwrap with default settings fails for virtual sites #3168

Open
mnmelo opened this issue Mar 15, 2021 · 26 comments
Open

Unwrap with default settings fails for virtual sites #3168

mnmelo opened this issue Mar 15, 2021 · 26 comments

Comments

@mnmelo
Copy link
Member

mnmelo commented Mar 15, 2021

Expected behavior

For an ag containing virtual-sites (which are always massless in GROMACS topologies) ag.unwrap() should unwrap ag atoms.

Actual behavior

Traceback (most recent call last):
  File "unwrap_err.py", line 4, in <module>
    u.atoms.unwrap()
  File "/.../MDAnalysis/core/groups.py", line 1726, in unwrap
    "the {} is zero.".format(comp))
ValueError: Cannot perform unwrap with reference='com' because the total mass of at least one of the fragments is zero.

Code to reproduce the behavior

import MDAnalysis as mda
from MDAnalysisTests.datafiles import GRO, TPR
u = mda.Universe(TPR, GRO)
u.atoms.unwrap()

Cause / possible solutions

  1. Virtual-sites share no bonds with the rest of the molecule. This is an issue on its own. We should probably add bonds between virtual sites and all constructing atoms. If these were in place, virtual sites would no longer be each in a fragment by themselves.
  2. The unwrap code could check for compound size when com shifting and skip unwrap+weighting for single-particle compounds (for these, a simple apply_PBC suffices). This also should speedup the unwrapping of systems with any single-particle compound (ions, CG waters, etc.). This solution sidesteps the error but may allow virtual sites to be unwrapped away from their constructing atoms.

Current version of MDAnalysis

  • Which version are you using? 2.0.0-dev0 at commit 6282385
  • Which version of Python (python -V)? 3.6.9
  • Which operating system? Ubuntu 18.04

(EDIT: numbered solutions)

@mnmelo
Copy link
Member Author

mnmelo commented Mar 15, 2021

This bug is a manifestation of what @orbeckst already stated in #1954, and follows the discussion there and in #588.
I am preparing a PR that addresses solution 2 above, but the most elegant would be what @jbarnoud proposes in #1954.

@BartBruininks
Copy link
Contributor

I was wondering if this is till an active topic. It seems to still be broken.

@jbarnoud
Copy link
Contributor

jbarnoud commented Oct 9, 2023

To my knowledge, nobody is working on this. The issue is that MDAnalysis does not know how virtual sites are related to the rest of the atoms.

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 9, 2023 via email

@jbarnoud
Copy link
Contributor

jbarnoud commented Oct 9, 2023

The latest solution we discussed was to encode the virtual sites as a list of duplets (virtual site, underlying atom) and to use consider these the same as bonds when building fragments. This would go in the TPR parser with some changes in the code about fragments.

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 9, 2023 via email

@orbeckst
Copy link
Member

orbeckst commented Oct 9, 2023

@BartBruininks contributions are very welcome! Start with a PR and then it's a back-and-forth on the PR with multiple revisions until it has the right shape.

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 9, 2023 via email

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@orbeckst orbeckst linked a pull request Oct 10, 2023 that will close this issue
5 tasks
@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@orbeckst
Copy link
Member

Your PR #4318 came through just fine. I had wrongly linked it as fixing this issue; it would be nice if it did (although you said in #3168 (comment) that it doesn't, right?). Your PR does partially address #1954, though.

@BartBruininks
Copy link
Contributor

Sorry for my unclarity. I think the issue vsites -> bonds is fixed in the ITPReader. However, the problem I had with actually installing the git clone is not.

@orbeckst
Copy link
Member

Not sure what's happening with your installation. When you import MDAnalysis, make sure to not be in the checked out sources (at least that's often an issue that I encounter).

We also have https://userguide.mdanalysis.org/stable/contributing_code.html with typical installation and set-up instructions.

The discord server is a bit better suited for quick back and forth about installation issues so I'd move further discussion about installation to #developers (there are also more people who can help than in this issue). (To join the Discord server see https://www.mdanalysis.org/#participating ).

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@BartBruininks
Copy link
Contributor

BartBruininks commented Oct 10, 2023 via email

@orbeckst
Copy link
Member

I'll continue discussion on PR #4318 :-)

@jbarnoud
Copy link
Contributor

Hi Bart. What you describe is the general idea. However, I would not add the virtual sites to the bonds themselves. Indeed, while they are connections, they are not bonds despite the very loose definition of the term. Think about how you would draw a molecule or how you would expect to look in vmd, you would not expect some sort of star coming out of your virtual sites. Same goes with bond length distributions.

What I would rather have is to add the virtual site to a new VirtualSites topology attribute (to define in mda.core.topologyattrs) that would look very similar to the Bonds topology attribute. Then, I would alter the function that builds the fragments (that is a method of the Bonds topology attribute defined in mda.core.topologyattrs iirc) to account forvthe virtual sites in addition to the bonds.

One issue I foresee but you don't have to care about just now is that the fragment method is transplanted from the Bonds attribute and I don't know from the top of my mind how to transplant it from the VirtualSites attribute as well. Perhaps @richardjgowers would have an idea how to handle that. In practice, all the format that define virtual sites also define bonds, so this may be nitpicking.

@BartBruininks
Copy link
Contributor

I see your point, although I would love to see that star and calculate the vsite lengths as bonds.

However, I agree that one way or another you would want to have control on including them in the bonds or not (without re-reading the topology).

@hejamu
Copy link
Contributor

hejamu commented Oct 11, 2023

To add to @jbarnoud's point, some MD softwares treat virtual sites fully implicitly, meaning your trajectory will not contain them. With a TopologyAttribute that one can add to the universe that would make working with explicit and implicit virtual sites consistent.

@mattwthompson
Copy link
Contributor

As a user on the hunt for a tools that provide good visualizations of trajectories with virtual sites, I'd love to see better support for them somewhere. I recognize, and can say with some first-hand experience, that it's a difficult problem. But for both intuitive and practical reasons, I don't think storing them as normal bonds is the way to go. OpenMM's and GROMACS's implementations do not do this and theirs are the two best implementation of virtual sites I've seen in engines so far. (In general, they're dummy atoms that are massless, barely exist in the topology, and are subject to a weight average of the forces on the atoms that define them.) Amber adds them in as explicit bonded neighbors, which works fine for some cases but makes things more complicated than 4- and 5-site water models too hard for me to decipher at the moment.

@orbeckst
Copy link
Member

I summarized the discussion above at #1954 (comment) .

It sounds to me that we should focus on first defining properly how vsites ought to be handled and then proceed with PR #4318.

If you can provide additional technical insights on how to do this or — even better — want to get involved then I suggest to have the initial technical discussion on the vsites/dummies data structures on #1954 .

@kimdn
Copy link

kimdn commented May 29, 2024

Btw I seem to have a problem with installing MDAnalysis using pip -e . in the package folder: It install just fine in my venv python3.10.13 but then when I try to import it (python3.10 -c import MDAnalysis) I get: Traceback (most recent call last): File "/media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/lib/util.py", line 220, in from ._cutil import unique_int_1d ImportError: /media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/lib/_cutil.cpython-310-x86_64-linux-gnu.so: failed to map segment from shared object During handling of the above exception, another exception occurred: Traceback (most recent call last): File "", line 1, in File "/media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/init.py", line 193, in from .lib import log File "/media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/lib/init.py", line 34, in from . import transformations File "/media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/lib/transformations.py", line 172, in from .mdamath import angle as vecangle File "/media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/lib/mdamath.py", line 63, in from . import util File "/media/windows/Users/Bart/projects/mda_development/mdanalysis/package/MDAnalysis/lib/util.py", line 222, in raise ImportError("MDAnalysis not installed properly. " ImportError: MDAnalysis not installed properly. This can happen if your C extensions have not been built. Op di 10 okt 2023 om 16:20 schreef Bart Bruininks @.>:

I think/hope that I just send a PR to MDA. Op di 10 okt 2023 om 15:50 schreef Bart Bruininks @.
>: > > Ok, i think I have added the virtual sites and everything is working in short: > > 1) add the virtual site parsers for each of the virtual_sites[1,2,3,4,n]. > 2) write the respective parsers in Molecule and add each base-vsite as a bond. > 3) fragments seem to work fine and nothing is change beyond this. > > What I did so far: > git clone mda > git branch vsites_fragments > git add . > git commit -m 'adds vsites as base-vsite bonds in the ITP reader > (Molecule) using vsite parsers' > gith push??? (where and how ><) > > Now how would I send this for review in the MDA project? I am sorry > for being a noob, but it has been a long long time since I added a > contribution to some one else's library. > > Cheers, > > Bart > > PS > Changed pieces of code in ITPParser.Molecule an example is given for > parse_virutal_sites3: > > > self.parsers = { > 'atoms': self.parse_atoms, > 'bonds': self.parse_bonds, > 'angles': self.parse_angles, > 'dihedrals': self.parse_dihedrals, > 'constraints': self.parse_constraints, > 'virtual_sites1' : self.parse_virtual_sites1, #BMHB > 'virtual_sites2' : self.parse_virtual_sites2, #BMHB > 'virtual_sites3' : self.parse_virtual_sites3, #BMHB > 'virtual_sites4' : self.parse_virtual_sites4, #BMHB > 'virtual_sitesn' : self.parse_virtual_sitesn, #BMHB > 'settles': self.parse_settles > } > > def parse_virtual_sites3(self, line): # BMHB > # [ virtual_sites[n1234] ] are formats of constructing > # a particle geometrically based on the base particles > # that define it, with some underlying values to > # specify the distance of the vsite with respect to the > # base particles. In GMX there are multiple forms of > # vsites, but we are going to parse them all in a similar > # manner. We add a bonds between the vsite and all of its > # base particles. Example: vsite3 has first three elements as > # the the base particles the 4th element is the vsite. Thus > # 3 bonds are created 1-4, 2-4, 3-4. This is the same for > # all vsites with a different number so vsite2 adds 2 bonds. > n_funct = 3 > elements = line.split() > base_indices = elements[:n_funct] > vsite_index = elements[n_funct] > funct_value = elements[n_funct+1] > # Add every element vsite pair as a bond. > for base_index in base_indices: > temp_line = f'{base_index} {vsite_index} {funct_value} 0.0 0.0' > self.add_param(temp_line, self.bonds, n_funct=2, > funct_values=(1, 2, 3, 4,)) > > > > Op ma 9 okt 2023 om 15:44 schreef Bart Bruininks @.>: > > > > Ok I'll give it a go. > > > > Op ma 9 okt. 2023 15:21 schreef Oliver Beckstein @.>: > >> > >> @BartBruininks contributions are very welcome! Start with a PR and then it's a back-and-forth on the PR with multiple revisions until it has the right shape. > >> > >> — > >> Reply to this email directly, view it on GitHub, or unsubscribe. > >> You are receiving this because you were mentioned.Message ID: @.***>

My macbook Pro has the same problem

"ImportError: MDAnalysis not installed properly. This can happen if your C extensions have not been built.".

GPT4 says

"The error you're encountering is due to the architecture incompatibility of the MDAnalysis library with your ARM-based Mac (Apple Silicon). "

Indeed, I use Apple M1 Max Chip.

@orbeckst
Copy link
Member

@kimdn I think GPT4 is relying on outdated trainings data... have a look at https://www.mdanalysis.org/2022/06/27/applem1packages/. MDAnalysis builds just fine on Apple silicon and we have conda packages for osx-arm64.

Furthermore, please open a new issue as installation problems are not relevant to the discussion here. Having separate issues really helps us to focus the discussion. However, for installation problems I recommend you ask in our Discussions forum or in our Discord, where people are happy to engage in helpful discussions. Thank you!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

7 participants