Topology readers are not aware of virtual sites

[jbarnoud]

Virtual sites are atoms the coordinates of which are built frame by frame from the coordinates of other atoms. MDAnalysis reads the coordinates of these atoms as for any other ones. Yet, it does not read the connection between the virtual sites and the atoms they are built from. Therefore, the virtual sites appear as not connected to the rest of the molecule which breaks the expected behaviour of ag.fragments and make_whole.

A solution would be to have a VirtualSites topology attribute that connects the relevant atoms. With the relationship described, the appropriate edges can be added to the graph used by ag.fragments and make_whole can translate the virtual sites with the building atoms. The difficulty is that virtual sites can be built from any number of atoms, yet from my tests, TopologyGroups store the indices as an array of ints. I do not know if it is possible to make a TopologyGroup with an array of lists (dtype=object) instead.

Note that virtual sites are especially used by gromacs, and I have the TPR reader in mind, but they seem to exist in other simulation engines.

See #1949 and #588.

[orbeckst]

The issue that TIP4P is not recognized as a fragment (details in #1949) makes the "make whole" transformation PR #2038 not correctly wrap typical OPLS-AA simulation systems (such as our TPR/XTC test trajectory).

[orbeckst]

Let me summarize the discussion on #3168 (based @BartBruininks initial implementation of virtual sites as additional bonded interactions in PR #4318 and on @jbarnoud 's #3168 (comment) , @hejamu 's #3168 (comment) , and @mattwthompson '#3168 (comment) )

• Reading vsites and adding the dummy particles as bonded particles to the topology is straightforward and works, PR adds vsites as bonds between base particles and vsite #4318.
• However, there's broad consensus that it's important to not confuse dummies/vsites with real atoms and chemical bonds because
  1. it makes the actual chemistry confusing,
  2. implementations of vsites differ between codes so a flexible/customizable solution is needed

How to implement vsites? — essentially what @jbarnoud said at the top of this issue:

vsites should be tracked with a separate TopologyAttr, i.e., add the virtual site to a new VirtualSites topology attribute (to be defined in mda.core.topologyattrs) that would look very similar to the Bonds topology attribute. (This would resolve #588 .)

In order to support the functionality required to solve #3168 (unwrapping), the function that builds the fragments (which is a method of the Bonds topology attribute defined in mda.core.topologyattrs) is changed to account for the virtual sites in addition to the bonds.

[BartBruininks]

Sounds good to me. This leaves me with some questions about the vbonds attribute. If the vsites are in a separate attribute. What is the format of this. Vsites come in many flavours, meaning putting them in a single matrix will be hard. I would be in favour of making the vsites attribute to be like the bonds, always linking 2 atoms. So a virtualsites4 (a,b,c,d,e) in gromacs would add 4 vbonds ((a,e), (b,e), (c,e) (d,e)). And the virtualsitesn adds n vbonds. Then all virtual sites are obtainable by np.unique(vbonds[:,1]). Or should there also be a vatoms attribute? Cheers Op zo 15 okt. 2023 18:27 schreef Oliver Beckstein ***@***.***

…

: Let me summarize the discussion on #3168 <#3168> (based @BartBruininks <https://github.com/BartBruininks> initial implementation of virtual sites as additional bonded interactions in PR #4318 <#4318> and on @jbarnoud <https://github.com/jbarnoud> 's #3168 (comment) <#3168 (comment)> , @hejamu <https://github.com/hejamu> 's #3168 (comment) <#3168 (comment)> , and @mattwthompson <https://github.com/mattwthompson> '#3168 (comment) <#3168 (comment)> ) - Reading vsites and adding the dummy particles as bonded particles to the topology is straightforward and works, PR #4318 <#4318>. - However, there's broad consensus that it's important to not confuse dummies/vsites with real atoms and chemical bonds because 1. it makes the actual chemistry confusing, 2. implementations of vsites differ between codes so a flexible/customizable solution is needed *How to implement vsites?* — essentially what @jbarnoud <https://github.com/jbarnoud> said at the top of this issue: vsites should be tracked with a separate TopologyAttr, i.e., add the virtual site to a new *VirtualSites* topology attribute (to be defined in mda.core.topologyattrs) that would look very similar to the Bonds topology attribute. (This would resolve #588 <#588> .) In order to support the functionality required to solve #3168 <#3168> (unwrapping), the function that builds the fragments (which is a method of the Bonds topology attribute defined in mda.core.topologyattrs) is changed to account for the virtual sites in addition to the bonds. — Reply to this email directly, view it on GitHub <#1954 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AALLQBYYKYSWUF3APTB5GCLX7QFGDAVCNFSM4FGTSNO2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCNZWGM2DGOJXHA4A> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[jbarnoud]

As you said, virtual sites do not all involve the same number of particles, and therefore cannot be stored all in one array if we store them one site per row. The way I understand it, you came to the same conclusion as I did: store vsites as pairs of particles with the first particle being the site and the second being one of the particles from which it is built.

[orbeckst]

@BartBruininks is this sufficient information to make progress on PR #4318 ?

[BartBruininks]

I think it is pretty clear what is intended, however, I do not know if I have time and the understanding to add the additional vsites attributes and make sure it is handled correctly by all the other methods in MDA (almost everything that can make use of a bond should be altered to also take the vsites into account). This was the reason I was just adding them to the bonds. In my opinion adding them to the bonds and then having a list of the vsites is much more trivial, then the user can always use the vsiste selection to prune the bonds taken into account, where maybe default behavior is to remove all the bonds which have a link to vsite? Op vr 15 dec 2023 om 19:39 schreef Oliver Beckstein < ***@***.***>:

…

@BartBruininks <https://github.com/BartBruininks> is this sufficient information to make progress on PR #4318 <#4318> ? — Reply to this email directly, view it on GitHub <#1954 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AALLQB4QKM2BJ54T4NF74RDYJSKPFAVCNFSM4FGTSNO2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCOBVHAZTEOJVGE3Q> . You are receiving this because you were mentioned.Message ID: ***@***.***>