MD analysis

72 views
Skip to first unread message

Viswanath Vittaladevaram

unread,
Dec 21, 2021, 9:02:53 PM12/21/21
to MDnalysis discussion

I am seeing this error, need help on this!

warnings.warn("Failed to guess the mass for the following atom types: {}".format(atom_type)) Traceback (most recent call last): File "traj_residue-z.py", line 48, in protein_z=protein.centroid()[2] IndexError: index 2 is out of bounds for axis 0 with size 0

Irfan Alibay

unread,
Dec 21, 2021, 9:04:33 PM12/21/21
to MDnalysis discussion
Hi,

Welcome to MDAnalysis!

If possible could you provide more information on this please? Could you post the exact snippet of code that you are attempting to run?

Many thanks,

Irfan

Viswanath Vittaladevaram

unread,
Dec 21, 2021, 9:07:17 PM12/21/21
to MDnalysis discussion
#!/opt/local/bin/python

## TRAJ_RESIDUE_Z
##
## Program to determine the z separation between residue and surface
##
## Use terminal heavy atom for each residue

import sys, numpy, argparse
import MDAnalysis

p = argparse.ArgumentParser()
p.add_argument("-g","--grofile")
p.add_argument("-x","--xtcfile")
p.add_argument("-o1","--outfile1")
p.add_argument("-o2","--outfile2")
p.add_argument("--nres",type=int)
args=p.parse_args()


groFile=args.grofile
xtcFile=args.xtcfile
outFile1=args.outfile1
outFile2=args.outfile2
nprotein_res=args.nres

terminal_atom={'ALA':'CB','ARG':'CZ','ASN':'CG','ASP':'CG','CYS':'SG',
               'GLN':'CD','GLU':'CD','GLY':'CA','HIS':'CE1','ILE':'CD',
               'LEU':'CG','LYS':'NZ','MET':'CE','PHE':'CZ','PRO':'CG',
               'SER':'OG','THR':'CG2','TRP':'CH2','TYR':'OH','VAL':'CG1'}

u=MDAnalysis.Universe(groFile,xtcFile)
ouf1=open(outFile1,'w')
ouf2=open(outFile2,'w')

protein=u.select_atoms('resid 1-%d' % (nprotein_res))

nres=len(protein.residues)
for res in protein.residues:
    print ("##  %6d  %4s" % (res.resid,res.resname),file=ouf2)


z_ave=numpy.zeros((nres))
z_sqr=numpy.zeros((nres))
norm=0.0
for ts in u.trajectory:

    protein_z=protein.centroid()[2]
    print ("%12.3f %8.3f " % (ts.time,protein_z))
    print ("%12.3f %8.3f " % (ts.time,protein_z),file=ouf1)
    for ires,res in enumerate(protein.residues):
        res_z=res.atoms.centroid()[2]
        print ("%12.3f %6d %8.3f " % (ts.time,res.resid,res_z),file=ouf2)

Oliver Beckstein

unread,
Dec 21, 2021, 9:22:15 PM12/21/21
to mdnalysis-...@googlegroups.com
Hello,

Thanks for coming from StackOverflow to the mailinglist. Hopefully we can help you here. 

Does your selection contain any atoms?

What does 

print(protein)

show?

Also look at the output of centroid():

print(protein.centroid())

Perhaps that makes clearer what happens. 

Oliver

Am 12/21/21 um 19:07 schrieb Viswanath Vittaladevaram <viswa....@gmail.com>:

#!/opt/local/bin/python
--
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/9cfe96de-27f3-495f-b227-278e1d79e22fn%40googlegroups.com.

Viswanath Vittaladevaram

unread,
Dec 21, 2021, 9:33:09 PM12/21/21
to MDnalysis discussion
I am trying to run this 

python traj_residue-z.py -g fnIII-9_ps20_Nchain6_T298_nw.gro -x fnIII-9_ps20_Nchain6_T298_nw.xtc -o1 fnIII-9_ps20_Nchain6_run1.protein_z -o2 fnIII-9_ps20_Nchain6_run1.res_z --nres 89

Oliver Beckstein

unread,
Dec 21, 2021, 9:48:00 PM12/21/21
to mdnalysis-discussion
Hi,


On Dec 21, 2021, at 7:33 PM, Viswanath Vittaladevaram <viswa....@gmail.com> wrote:

I am trying to run this 

Did you write the script or did you get the script from someone else?



python traj_residue-z.py -g fnIII-9_ps20_Nchain6_T298_nw.gro -x fnIII-9_ps20_Nchain6_T298_nw.xtc -o1 fnIII-9_ps20_Nchain6_run1.protein_z -o2 fnIII-9_ps20_Nchain6_run1.res_z --nres 89

That does not really help me, I am sorry.. I don’t know what your data looks like. The only way we can help is if you provide some diagnostics, which means that you need to add code to your script to print out information that we need to be able to understand what’s happening.

Change your script by adding the lines in red:
----------

print("DEBUG: protein = “, protein)


nres=len(protein.residues)
for res in protein.residues:
    print ("##  %6d  %4s" % (res.resid,res.resname),file=ouf2)


z_ave=numpy.zeros((nres))
z_sqr=numpy.zeros((nres))
norm=0.0
for ts in u.trajectory:
    print("DEBUG: protein.centroid() = “, protein.centroid())

    protein_z=protein.centroid()[2]
    print ("%12.3f %8.3f " % (ts.time,protein_z))
    print ("%12.3f %8.3f " % (ts.time,protein_z),file=ouf1)
    for ires,res in enumerate(protein.residues):
        print("DEBUG: res = ", res)
        print(“DEBUG: res.centroid() = “, res.centroid())
        res_z=res.atoms.centroid()[2]
        print ("%12.3f %6d %8.3f " % (ts.time,res.resid,res_z),file=ouf2)


————

Then run it and copy and paste all output.

It is very likely that some of your selections return empty atomgroups and then centroid() returns an empty array that cannot be indexed. The debug print calls will show you what happens. You might need to add additional code to skip any selections that do not exist in your system of interest.

Oliver

p.s.: This is how I tested what centroid() does. 

In [1]: import MDAnalysis as mda; from MDAnalysis.tests import datafiles as data

In [2]: u = mda.Universe(data.TPR, data.XTC)

In [3]: u.select_atoms("resname foo”)   # residue foo does not exist
Out[3]: <AtomGroup with 0 atoms>

In [4]: ag = u.select_atoms("resname foo")

In [5]: ag.centroid()
Out[5]: array([], shape=(0, 3), dtype=float32)

In [6]: glutamates = u.select_atoms("resname GLU")

In [7]: glutamates.centroid()
Out[7]: array([61.15441038, 46.87385399, 26.82585312])

And this shows how you can check if an AtomGroup ag is empty (or check len(ag) > 0):

In [8]: if ag:
   ...:     print("AtomGroup contains atoms")
   ...: else:
   ...:     print("AtomGroup is empty")
   ...:
AtomGroup is empty
--
Oliver Beckstein (he/his/him)


MDAnalysis – a NumFOCUS fiscally sponsored project




Viswanath Vittaladevaram

unread,
Dec 21, 2021, 10:11:09 PM12/21/21
to MDnalysis discussion
Hi,
I tried as you said and actually it was not written by me. I was just doing analysis using code. 

As you highlighted the red ones, I tried and I see this!

File "traj_residue-z.py", line 38

    print("DEBUG: protein = “, protein)
                                       ^
SyntaxError: EOL while scanning string literal

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 6:19:02 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,
can I share the entire requirement file with code and data here!

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 8:02:13 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

Please find word document attached showing set of instructions that needs to be done. 
Fibronectin-Brush Analysis (1).docx

Irfan Alibay

unread,
Dec 22, 2021, 9:12:32 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

The string literal scanning error is due to improper quotes being used for this line.

Unfortunately emails aren't geared towards writing code, so Oliver's debug prints got mishandled by formatting.

Can you change all the print statements so that they have " quotes? This should hopefully fix this specific issue.

Please do report back with the debug information.

Best regards,

Irfan

--
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.

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 9:19:00 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

Can you please show me one example where to input commas

Irfan Alibay

unread,
Dec 22, 2021, 9:31:48 AM12/22/21
to mdnalysis-...@googlegroups.com
So in the code sections that Oliver provided in red in his previous email you should make sure you have the right commas in.

For example:

print("DEBUG: protein = ", protein) 

instead of:

print("DEBUG: protein = “, protein)

Best regards,

Irfan
 


Viswanath Vittaladevaram

unread,
Dec 22, 2021, 9:38:55 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

Thanks for your reply. But my concern is not with debug thing.

Here when I am trying to execute through MD analysis I am getting this error.
Traceback (most recent call last):
  File "traj_residue-z.py", line 53, in <module>

    protein_z=protein.centroid()[2]
IndexError: index 2 is out of bounds for axis 0 with size 0
You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/J8oJ0M9Rjb4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKVkbJgOAiVPg6naVuu-RcCOt4k8rWR9SpiH0bhGiQEuE%2BbM6A%40mail.gmail.com.

Irfan Alibay

unread,
Dec 22, 2021, 9:43:31 AM12/22/21
to mdnalysis-...@googlegroups.com
Yes I understood this, however Oliver's debug prints will help us understand what is going on with your code. Purely from the look of that one error all I can tell is that you AtomGroup "protein" is probably empty.

Oliver's debug prints were to verify this.

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 9:44:46 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,
I understand that but could you please elaborate how to fix this issue?

Irfan Alibay

unread,
Dec 22, 2021, 9:58:37 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

Here's how I generally approach debugging for something like this:

1) Our first thought is that the AtomGroup "protein" is likely empty. So first we need to be sure that this is indeed the case. Hence Oliver's debug prints.

2) Once we can tell that the issue is because the AtomGroup "protein" is empty, then the next step is to understand why that is the case. For the moment, the only thing I can think of given the information I have is that "protein=u.select_atoms('resid 1-%d' % (nprotein_res))" isn't doing what you intend it to so.

3) If it is that line, then it could be a number of different things. Maybe you don't have any residues in the range of 1 to 89. Maybe something is breaking argparse that I can't seen. etc...

Because debugging often ends up being a funnel based approach, the best thing we can do here is work our way down these steps to narrow down to the issue in the most efficient manner. In that interest, if you can confirm step 1 (i.e. that the AtomGroup is indeed empty using Oliver's debug prints, or otherwise), then we can progress towards a solution.

Best regards,

Irfan

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 10:02:51 AM12/22/21
to mdnalysis-...@googlegroups.com
Okay. I confirm and can you show how it works in both ways!

Irfan Alibay

unread,
Dec 22, 2021, 10:19:55 AM12/22/21
to mdnalysis-...@googlegroups.com
>   can you show how it works in both ways

I'm not sure what you mean by this sorry.

In any case if you confirm that the AtomGroup is indeed empty, then can you check your input GRO file to see if there are any atoms in the residue range 1-89?
You can either do this manually by opening the file with your favourite file reader (vim, etc...) or with MDAnalysis by doing the following:

"""
import MDAnalysis as mda

u = mda.Universe("fnIII-9_ps20_Nchain6_T298_nw.gro", "fnIII-9_ps20_Nchain6_T298_nw.xtc")

print(u.residues.resids[:10])
"""

If you can let us know what the first ~ 10 resid numbers are then that should help us work out if there's anything there.

Best regards,

Irfan

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 10:30:44 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,
1327LEU      N    1   2.013   3.349   8.848  0.4933 -0.2510  0.2982
 1327LEU     H1    2   1.953   3.277   8.893  0.0174  0.1791  0.3637
 1327LEU     H2    3   1.960   3.377   8.762  0.6275 -0.5669  0.1094
 1327LEU     H3    4   2.092   3.292   8.812  0.2711 -0.8049  0.6862
 1327LEU     CA    5   2.058   3.459   8.921  0.3604  0.7395  0.5225
 1327LEU     HA    6   1.975   3.524   8.944 -0.1585 -0.6378  2.9136
 1327LEU     CB    7   2.127   3.409   9.049  0.3145  0.2983 -0.1429
 1327LEU    HB1    8   2.062   3.345   9.112  0.8284 -0.2938 -0.2122
 1327LEU    HB2    9   2.217   3.352   9.018 -0.4983  0.4091 -2.9675
 1327LEU     CG   10   2.178   3.517   9.140  0.5360 -0.2089 -0.1033

Irfan Alibay

unread,
Dec 22, 2021, 10:42:37 AM12/22/21
to mdnalysis-...@googlegroups.com
Ok so your residue numbering doesn't start until 1327, so it's normal that if you selected for residues numbered (by the GRO file numbering) 1 to 89 you wouldn't get anything. Hence the empty AtomGroup.

Your option here is to either change the code to select the first nres residues (if that is specifically what you want to do), which you could do by changing the line:

protein=u.select_atoms('resid 1-%d' % (nprotein_res))

to

protein = u.residues[:protein_res]

or by changing your code to parse the residue ranges more accurately (this would require changing the argparse arguments accordingly).

Best regards,

Irfan


Viswanath Vittaladevaram

unread,
Dec 22, 2021, 11:15:21 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

As you said I did changes but instead I did this change as mentioned below and it is working

protein = u.residues[:nprotein_res]

Thank you very much.


Viswanath Vittaladevaram

unread,
Dec 22, 2021, 11:30:42 AM12/22/21
to mdnalysis-...@googlegroups.com
Further continuation to MD Analysis I tried to run next code and it shows error like this below.

traj_brush_prof.py: error: unrecognized arguments: -- brush_res_stop 244 --nbin 1400 --deltaz 0.1 --tskip 100000

Irfan Alibay

unread,
Dec 22, 2021, 11:36:20 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

Here again we would require information on the code contained within the `traj_brush_prof.py` script for us to understand what is going on.

That being said - "unrecognized arguments" errors is clearly an issue with argparse rather than the MDAnalysis library. It might be more fruitful to ask for support from the person who originally wrote this script?

Best regards,

Irfan

Viswanath Vittaladevaram

unread,
Dec 22, 2021, 11:41:19 AM12/22/21
to mdnalysis-...@googlegroups.com
Hi,

Okay. But I would like to have your view on the code once , find below
#!/opt/local/bin/python

"""
TRAJ_BRUSH_PROF

Program to calculate density profiles for
peptoid brushes

Uses C-alpha number density as hopefully reasonable
proxy
"""

import numpy, sys, argparse
import MDAnalysis

def anint(x):
    xlo=numpy.floor(x)
    xhi=xlo+1.0
    if xhi-x<x-xlo:
        return xhi
    else:
        return xlo

## get command line arguments

p = argparse.ArgumentParser()
p.add_argument("-g","--grofile")
p.add_argument("-x","--xtcfile")
p.add_argument("-o","--outfile")
p.add_argument("--brush_res_start",type=int)
p.add_argument("--brush_res_stop",type=int)
p.add_argument("--nbin",type=int)
p.add_argument("--deltaz",type=float)
p.add_argument("--tskip",type=float)
args=p.parse_args()


groFile=args.grofile
xtcFile=args.xtcfile
outFile=args.outfile


brush_res_start=args.brush_res_start
brush_res_stop=args.brush_res_stop

nbin=args.nbin
deltaz=args.deltaz
tskip=args.tskip

## initialise stuff
density_prof=numpy.zeros((nbin+1))
zmax=nbin*deltaz
norm=0.0

## open input files
u=MDAnalysis.Universe(grofile,xtcfile)

## make selection
ca_sel=u.select_atoms("resid %d:%d and name CA" % (brush_res_start,brush_res_stop))

for ts in u.trajectory:
    print (ts.time)
    if ts.time<tskip:
        continue

    norm+=1.0
    area=ts.dimensions[0]*ts.dimensions[1]

    ## loop over CA atoms
    for ca in ca_sel:
        zz=ca.position[2]
        if zz>zmax:
            continue

        ibin=int(anint(zz/deltaz))
        density_prof[ibin]+=1.0

dvol=area*deltaz
density_prof/=norm
density_prof/=dvol

ouf=open(outfile,'w')
for i in range(nbin):
    zz=(i+0.5)*deltaz
    print (zz,density_prof[i],file=ouf)

Irfan Alibay

unread,
Dec 22, 2021, 11:46:28 AM12/22/21
to mdnalysis-...@googlegroups.com
I can't 100% tell if that's the case or this is is just bad pasting into an email but:

"-- brush_res_stop 244 --nbin 1400 --deltaz 0.1 --tskip 100000"

There's a space between `--` and `brush_res_stop` which shouldn't be there?

Oliver Beckstein

unread,
Dec 22, 2021, 12:24:00 PM12/22/21
to MDnalysis discussion
I am glad it is working now.

I added an answer to you S/O post https://stackoverflow.com/a/70452793/ with links to this thread, just so that the S/O post does not look unanswered.

Best,
Oliver

Viswanath Vittaladevaram

unread,
Dec 23, 2021, 4:33:32 AM12/23/21
to mdnalysis-...@googlegroups.com
HI,

I am trying to run this(python3 traj_orientation-group.py -g fnIII-9_ps20_Nchain6_T298_nw.gro -x fnIII-9_ps20_Nchain6_T298_nw.xtc -o fnIII-9_ps20_Nchain6_run1.phrsn-orientation --protein_res_start 1 --protein_res_stop 89 –group 51 52 53 54 55) for MD analysis and I got this following below:

usage: traj_orientation-group.py [-h] [-g GROFILE] [-x XTCFILE] [-o OUTFILE] [--protein_res_start PROTEIN_RES_START] [--protein_res_stop PROTEIN_RES_STOP] [--group GROUP [GROUP ...]]
traj_orientation-group.py: error: unrecognized arguments: –group 51 52 53 54 55

Can you help on this!

--
You received this message because you are subscribed to a topic in the Google Groups "MDnalysis discussion" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mdnalysis-discussion/J8oJ0M9Rjb4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mdnalysis-discus...@googlegroups.com.

Irfan Alibay

unread,
Dec 23, 2021, 4:41:29 AM12/23/21
to MDnalysis discussion
Hi,

Similar to my previous comment, it's hard to tell based on email formatting but are you possibly making a typo in your arguments?

I.e. from your email it seems that you are using "–group" instead of "--group" (i.e. using a single dash rather than two as is customary for multi-character flags in argparse)

- Irfan

Viswanath Vittaladevaram

unread,
Dec 23, 2021, 4:59:42 AM12/23/21
to mdnalysis-...@googlegroups.com
Hi, thanks for your message. 

I am finding this error too and attached code to this email.

File "traj_orientation-group.py", line 49, in <module>
    norm=numpy.sqrt(numpy.dot(v,v))
  File "<__array_function__ internals>", line 5, in dot
ValueError: shapes (0,3) and (0,3) not aligned: 3 (dim 1) != 0 (dim 0)

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/95d90fa4-f879-4d90-ad52-4cf7295aa062n%40googlegroups.com.
print.pdf

Irfan Alibay

unread,
Dec 23, 2021, 5:04:20 AM12/23/21
to MDnalysis discussion
Hi,

Following previous discussions, are you sure that the "protein" AtomGroup is not empty?

In your previous .gro file, the starting residue was 1327. If it's still the same system, setting the starting and stopping residues at 1 and 89 in this script will give you an empty selection?

Viswanath Vittaladevaram

unread,
Dec 23, 2021, 5:06:01 AM12/23/21
to mdnalysis-...@googlegroups.com
Hi,

Yes you are right. Starting residue was 1327 and I think I need to change it then

Viswanath Vittaladevaram

unread,
Dec 23, 2021, 8:05:34 AM12/23/21
to mdnalysis-...@googlegroups.com
Hi,

It worked finally.Thank you.
Reply all
Reply to author
Forward
0 new messages