search

MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.32k stars 651 forks source link

Topology readers are not aware of virtual sites #1954

Open jbarnoud opened 6 years ago

jbarnoud commented 6 years ago

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 commented 5 years ago

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 commented 1 year ago

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 https://github.com/MDAnalysis/mdanalysis/issues/3168#issuecomment-1755942135 , @hejamu 's https://github.com/MDAnalysis/mdanalysis/issues/3168#issuecomment-1757282196 , and @mattwthompson 'https://github.com/MDAnalysis/mdanalysis/issues/3168#issuecomment-1763222187 )

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 commented 1 year ago

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 https://github.com/MDAnalysis/mdanalysis/issues/3168 (based @BartBruininks https://github.com/BartBruininks initial implementation of virtual sites as additional bonded interactions in PR #4318 https://github.com/MDAnalysis/mdanalysis/pull/4318 and on @jbarnoud https://github.com/jbarnoud 's #3168 (comment) https://github.com/MDAnalysis/mdanalysis/issues/3168#issuecomment-1755942135 , @hejamu https://github.com/hejamu 's #3168 (comment) https://github.com/MDAnalysis/mdanalysis/issues/3168#issuecomment-1757282196 , and @mattwthompson https://github.com/mattwthompson '#3168 (comment) https://github.com/MDAnalysis/mdanalysis/issues/3168#issuecomment-1763222187 )

  • Reading vsites and adding the dummy particles as bonded particles to the topology is straightforward and works, PR #4318 https://github.com/MDAnalysis/mdanalysis/pull/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 https://github.com/MDAnalysis/mdanalysis/issues/588 .)

In order to support the functionality required to solve #3168 https://github.com/MDAnalysis/mdanalysis/issues/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 https://github.com/MDAnalysis/mdanalysis/issues/1954#issuecomment-1763439788, or unsubscribe https://github.com/notifications/unsubscribe-auth/AALLQBYYKYSWUF3APTB5GCLX7QFGDAVCNFSM4FGTSNO2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCNZWGM2DGOJXHA4A . You are receiving this because you were mentioned.Message ID: @.***>

jbarnoud commented 1 year ago

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 commented 11 months ago

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

BartBruininks commented 8 months ago

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 https://github.com/MDAnalysis/mdanalysis/pull/4318 ?

— Reply to this email directly, view it on GitHub https://github.com/MDAnalysis/mdanalysis/issues/1954#issuecomment-1858329517, or unsubscribe https://github.com/notifications/unsubscribe-auth/AALLQB4QKM2BJ54T4NF74RDYJSKPFAVCNFSM4FGTSNO2U5DIOJSWCZC7NNSXTN2JONZXKZKDN5WW2ZLOOQ5TCOBVHAZTEOJVGE3Q . You are receiving this because you were mentioned.Message ID: @.***>