hbonds analysis with lammps

[yann claveau]

Hello,

My system is water in a nanopore (O is type 1, H is type 2). I would like to plot the number of hbonds with respect to the radial position (pore centered). I use lammps.
.

I am quite new in python and MDA so I would like

If I understood correctly :

u = MDAnalysis.Universe(lammps.data, traj.lammpsdump)

hbonds = HBA(universe=u, hydrogens_sel='type 1', acceptors_sel='type 2') hbonds.run()

will give me

results = [ [ <frame>, <donor index (0-based)>, <hydrogen index (0-based)>, <acceptor index (0-based)>, <distance>, <angle> ], ... ]

From here, I have to get the number of occurence of each donor and its xy coordinates. However, the "donor index" is not the same as the lammps' one, since it is 0-based, is it?
 How can I do that ?

Thanks !

[Oliver Beckstein]

Hi Yann,

Welcome to the MDAnalysis list.

The “donor index” and “acceptor index” are indices into the universe.atoms AtomGroup: If you have donor index d_ix then the corresponding atom is

a = u.atoms[d_ix]

Then you can get the position of this Atom (a.position) and check where it is in your nano pore. In this way you could create a spatially-resolved hydrogenbond density.

Note that if you collect ALL the indices in a frame in a list or array d_indices then

donors = u.atoms[d_indices]

is an AtomGroup and donors.positions are all the positions. This might be easier to handle than looping over atoms.

Oliver

--

Oliver Beckstein (he/his/him)

email: orbe...@mdanalysis.org
twitter: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

[yann claveau]

Thanks a lot !

I tested several post process packages and to me, MDA is the simplest to use, it's a very great tool !

[yann claveau]

Hello again,

Other question: The default angle cutoff in VMD plugin is 20° while it is 150° in MDA. Does the 0° in VMD is 180° in MDA ? I guess it is not but just to be sure...

I have a remark. For someone (as for me) who begins in Python, the doc for hbonds could be clearer. It took me a lot of time to understand how to access the outputs of hbonds. Indeed:

     from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA

     hbonds = HBA(universe=u, donors_sel='name O', hydrogens_sel='name H', acceptors_sel='name O', d_h_a_angle_cutoff=30 ,d_a_cutoff=3.0)
     results= hbonds.run()

     print(results)

gives

     <MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis object at 0x7fa4a60f6580>

I have to try a dir(results) to understand that results can be accessed through results.hbonds

Thanks,

Yann

[Oliver Beckstein]

Hi Yann,

On Jun 18, 2021, at 9:38 AM, yann claveau <yann.c...@gmail.com> wrote:

Hello again,

Other question: The default angle cutoff in VMD plugin is 20° while it is 150° in MDA. Does the 0° in VMD is 180° in MDA ? I guess it is not but just to be sure...

Yes, I am pretty sure that this is the case. You need to check how VMD defines the angle, though. MDAnalysis defines it as the angle with the apex at the hydrogen between donor-H…acceptor:

acceptor

   \                

    \   (angle)     

     \       

      hydrogen————donor

I have a remark. For someone (as for me) who begins in Python, the doc for hbonds could be clearer. It took me a lot of time to understand how to access the outputs of hbonds.

Sorry, there’s a better hbonds tutorial for the User Guide in the works (see PR 133 — you’re very welcome to add there any comments, it helps us a lot to learn what users really want to see).

Indeed:

     from MDAnalysis.analysis.hydrogenbonds.hbond_analysis importHydrogenBondAnalysis as HBA

     hbonds = HBA(universe=u, donors_sel='name O', hydrogens_sel='name H', acceptors_sel='name O', d_h_a_angle_cutoff=30 ,d_a_cutoff=3.0)
     results= hbonds.run()

There’s a misunderstanding here. hbdonds.run() returns itself, namely results == hbonds. This is a design decision that makes it easier to chain method calls if necessary or combine setup and run (see below). Most of the analysis tools work the same way:

1. set up the analysis class instance x

2. do x.run()

3. pull data from x.results as indicated by the docs for the analysis tool

Often 1 and 2 are combined:

# set up and run

x = XAnalysis(…).run()

# analyze “observable"

plt.plot(x.results.time, x.results.observable)

I see that you’re using 2.0.0b. Starting with MDAnalysis 2.0.0 we are trying to make clearer where the computed results are so in 2.x, you’ll find them in hbonds.results.hbonds — all standard analysis tools get the “results” attribute to hold computed data. Once you know this, it makes it easier to discover the output, even if the docs are less clear.

     print(results)

gives

     <MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis object at 0x7fa4a60f6580>

I have to try a dir(results) to understand that results can be accessed through results.hbonds

HydrogenBondAnalysis is a complicated analysis tool. We strive to make the docs complete and we appreciate it when you tell us what’s missing from https://docs.mdanalysis.org/2.0.0-dev0/documentation_pages/analysis/hydrogenbonds.html . As I said above, beyond the API docs, we’re also working on tutorial-style docs and we love any contributions — be it comments or otherwise.

Oliver

-- 
You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/691a7b17-6ef0-4776-ab50-261de4192aa1n%40googlegroups.com.

[yann claveau]

Hi Oliver,

Le vendredi 18 juin 2021 à 19:36:44 UTC+2, orbe...@mdanalysis.org a écrit :

Hi Yann,

On Jun 18, 2021, at 9:38 AM, yann claveau <yann.c...@gmail.com> wrote:

Hello again,

Other question: The default angle cutoff in VMD plugin is 20° while it is 150° in MDA. Does the 0° in VMD is 180° in MDA ? I guess it is not but just to be sure...

Yes, I am pretty sure that this is the case. You need to check how VMD defines the angle, though. MDAnalysis defines it as the angle with the apex at the hydrogen between donor-H…acceptor:

acceptor

   \                

    \   (angle)     

     \       

      hydrogen————donor

I have a remark. For someone (as for me) who begins in Python, the doc for hbonds could be clearer. It took me a lot of time to understand how to access the outputs of hbonds.

Sorry, there’s a better hbonds tutorial for the User Guide in the works (see PR 133 — you’re very welcome to add there any comments, it helps us a lot to learn what users really want to see).

Thanks, it is very clear !

 

Indeed:

     from MDAnalysis.analysis.hydrogenbonds.hbond_analysis importHydrogenBondAnalysis as HBA

     hbonds = HBA(universe=u, donors_sel='name O', hydrogens_sel='name H', acceptors_sel='name O', d_h_a_angle_cutoff=30 ,d_a_cutoff=3.0)
     results= hbonds.run()

There’s a misunderstanding here. hbdonds.run() returns itself, namely results == hbonds. This is a design decision that makes it easier to chain method calls if necessary or combine setup and run (see below). Most of the analysis tools work the same way:

1. set up the analysis class instance x

2. do x.run()

3. pull data from x.results as indicated by the docs for the analysis tool

Often 1 and 2 are combined:

# set up and run

x = XAnalysis(…).run()

# analyze “observable"

plt.plot(x.results.time, x.results.observable)

I see that you’re using 2.0.0b. Starting with MDAnalysis 2.0.0 we are trying to make clearer where the computed results are so in 2.x, you’ll find them in hbonds.results.hbonds — all standard analysis tools get the “results” attribute to hold computed data. Once you know this, it makes it easier to discover the output, even if the docs are less clear.

I am using v1.1.1. I tried v2 for MSD but the result was strange so I get back to the stable version.

hbonds.results gives 'HydrogenBondAnalysis' object has no attribute 'results'

I have to write h.hbonds to get the results.

 

     print(results)

gives

     <MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis object at 0x7fa4a60f6580>

I have to try a dir(results) to understand that results can be accessed through results.hbonds

HydrogenBondAnalysis is a complicated analysis tool. We strive to make the docs complete and we appreciate it when you tell us what’s missing from https://docs.mdanalysis.org/2.0.0-dev0/documentation_pages/analysis/hydrogenbonds.html . As I said above, beyond the API docs, we’re also working on tutorial-style docs and we love any contributions — be it comments or otherwise.

I will gladly help ! The tool is really great.

I just have another question / remark :

When I loop over frames to get the positions (for instance to calculate the radial position), it took one minutes, sometimes more (without the universe loading).

   Ow = u.select_atoms("type 1")
   R=[]
   for ts in u.trajectory:
       posxy = Ow.ts.positions[:,0:2]
       r=np.linalg.norm(posxy,axis=1)
       R=np.append(R,r)

While, if I loop on an array containing the same datas, it took 1s. I guess it is because of the use of classes such as atoms. I was wondering if it can be improved ? Or if I did something wrong.

Yann

[Oliver Beckstein]

Hi Yann,

Paul (and Lily) just finished the hbond tutorial notebooks and they are available at https://userguide.mdanalysis.org/2.0.0-dev0/examples/analysis/hydrogen_bonds/README.html

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/c09b151d-0ae1-4d2d-aa1d-9244be4b0d9an%40googlegroups.com.

[yann claveau]

Thanks, I 've seen it. The final version explains why I had some problem. The distance cutoof should be 3.5 A and I have to multiply by 2 the results. Everything is fine, now.