pore radius from gromacs traj

[Nikhil Maroli]

Dear all,

i want to plot pore radius vs channel co-ordinate
and pore radius vs time ,i have done up to this,kindly tell me what to do next.

i wanted to get results in table as .dat file


import numpy as np
import os
import MDAnalysis as mda
from MDAnalysis.analysis.hole import HOLE, HOLEtraj
from MDAnalysis.tests.datafiles import PDB_HOLE
u = mda.Universe('ab.gro','abc.

xtc')
H = HOLEtraj(u)
H = HOLE(PDB_HOLE, executable="~/hole2/exe/hole")
H.run()
H.collect()

[Oliver Beckstein]

See below.

> On 28 Jul, 2016, at 06:30, Nikhil Maroli <scin...@gmail.com> wrote:
>
> Dear all,
>
> i want to plot pore radius vs channel co-ordinate
> and pore radius vs time ,i have done up to this,kindly tell me what to do next.
>
> i wanted to get results in table as .dat file

It would help if you told us what values you would like to have in the output file. Python is very flexible and in principle any output format can be written. If you are not overly familiar with Python then I understand that this can be a bit daunting at first so we're more than happy to help. But it does make our life a lot easier (and makes better use of our very limited time) if you are as specific as you can.

Now, given that I have to guess at what you want, I am just giving you a basic example and you can start modifying it to your own ends.

For a start, I don't know if you want to run HOLE on a trajectory or a single structure. I *assume* a trajectory because in the commands below you have your own "ab.gro" and "abc.xtc" for the trajectory analysis but the example PDB_HOLE for the single frame. Thus, I will ignore the single frame analysis – if this is what you want, say so clearly and paste the exact commands that you are using

TRAJECTORY ANALYSIS

I commented out lines that you don't need.

> import numpy as np
> # import os

> import MDAnalysis as mda

> # from MDAnalysis.analysis.hole import HOLE, HOLEtraj
# only need HOLEtraj for trajectory analysis
from MDAnalysis.analysis.hole import HOLEtraj

# don't need exampl file
> # from MDAnalysis.tests.datafiles import PDB_HOLE

> u = mda.Universe('ab.gro','abc.xtc')

# I assume you read the docs and you don't need any other options to run HOLE?
> H = HOLEtraj(u)


But you still need to actuall *run* the trajectory analysis:

H.run()

You can then write it to a file; here I separate each frame with '&' but you can adjust the code to whatever you need:

with open("profiles.xvg", "w") as xvg:
for frame, profile in H:
for _, zeta, radius in profile:
xvg.write("{0} {1}\n".format(zeta, radius))
xvg.write("&\n")

For more details, see the updated notebook at http://nbviewer.jupyter.org/gist/orbeckst/64c0bd5a037b3e434cc8ee6732030252#Extracting-data-to-a-file


NOTE: the rest is NOT needed and in fact, overwrites your trajectory analysis:

> H = HOLE(PDB_HOLE, executable="~/hole2/exe/hole")
> H.run()
> H.collect()

Hope that helps.

Oliver

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst * orbe...@gmail.com

[Nikhil Maroli]

Thanks for your detailed reply, Yeah , I'm new to python language and i am able to  generate pore radius for the structure and trying to perform other analysis which is available in mdanalysis .

thanks again. 

[Jeebu]

Hi Dr. Beckstein,

I am new to MDanalysis, this is a great software to analyze MD trajectory files.  I want to use HOLEtraj.

Continuing to your recent post on writing  profiles.xvg to a separate file for each frames you have mentioned to "One can use the above code to write the profiles in any format; for instance, it would be easy to write each frame's profile to a separate file by moving the with statement inside the outer for loop"

#

for frame, profile in H:

                 with open("profile.xvg", "w") as xvg:

                        for _, zeta, radius in profile:
                            xvg.write("{0} {1}\n".format(zeta, radius))
                xvg.write("&\n)

#

What else shall I change to write files profile1.xvg profile2.xvg .... for frames 1 2... separately.

Thank you.

[Oliver Beckstein]

Hello Jeebu,

Welcome to MDAnalysis. Glad to hear that it is useful to you.

On Feb 7, 2018, at 11:10 AM, Jeebu <gcjee...@gmail.com> wrote:

Hi Dr. Beckstein,

Just to be clear: most of the work is done by everyone else (the MDAnalysis core developers and many, many users who contribute code). You can just say “hi!” or “Hi MDAnalysis community”. Everyone on the list can answer (and in general you’d rather want anyone to answer than someone in particular, who might be busy… we all are).

I am new to MDanalysis, this is a great software to analyze MD trajectory files.  I want to use HOLEtraj.

Just fyi (and in case you haven’t seen these two notebooks on HOLEtraj:

HOLE with order parameters: https://gist.github.com/orbeckst/a0c856f58059112538591af108df6d59

HOLE hacking: https://gist.github.com/orbeckst/a0c856f58059112538591af108df6d59

They might be helpful… and we might actually put them up on the MDAnalysis Binder, come to think… https://mybinder.org/v2/gh/MDAnalysis/binder-notebook/master )

To your question:

Continuing to your recent post on writing  profiles.xvg to a separate file for each frames you have mentioned to "One can use the above code to write the profiles in any format; for instance, it would be easy to write each frame's profile to a separate file by moving the with statement inside the outer for loop"

#
for frame, profile in H:
                 with open("profile.xvg", "w") as xvg:
                        for _, zeta, radius in profile:
                            xvg.write("{0} {1}\n".format(zeta, radius))
                xvg.write("&\n)
#

What else shall I change to write files profile1.xvg profile2.xvg .... for frames 1 2... separately.

I think you have done it almost. You just have to use a different name for each file, something like

for frame, profile in H:

    filename = profile_{0:04d}.xvg".format(frame)   # produces profile_0000.xvg, profile_0001.xvg, ... 

    with open(filename, "w") as xvg:

         for _, zeta, radius in profile: 
              xvg.write("{0} {1}\n".format(zeta, radius)) 

You should get one profile_NNNN.xvg file per frame. You don’t need the “&” separator between data sets. (I haven’t tested the code, just wrote it down, so you might have to adjust things but it should be roughly what you need.)

Hope that helps.

Oliver

Thank you.

--

Oliver Beckstein

orbe...@gmail.com * orbe...@gmx.net

[Jeebu]

Thank you.

Your suggestions worked well.

filename = "profile_{0:04d}.xvg".format(frame)   ## added a close " in your suggestion.