MD analysis

[Viswanath Vittaladevaram]

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]

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]

#!/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]

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]

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]

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)

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

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

[Viswanath Vittaladevaram]

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]

Hi,

can I share the entire requirement file with code and data here!

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/f6b9a0c8-2d19-4c7a-8527-2b14e823fb2bn%40googlegroups.com.

[Viswanath Vittaladevaram]

Hi,

Please find word document attached showing set of instructions that needs to be done. 

[Irfan Alibay]

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.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/7896c558-12ca-42b1-ab60-c568d03b53f5n%40googlegroups.com.

[Viswanath Vittaladevaram]

Hi,

Can you please show me one example where to input commas

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKVkbJhMYRhomSRH7PLi7%2BaD36Fbio%2Bij2-%2B6q05wvBbC7sfyg%40mail.gmail.com.

[Irfan Alibay]

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

 

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZgnL%3DoFzAzXqU_fTNmcSwMCV%3DOb9srW0UkmadpaE52exw%40mail.gmail.com.

[Viswanath Vittaladevaram]

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]

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.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZjzcKQYqe9eOF4FYwz2q83xPQxzktoQSGb6Kvk79K7_uw%40mail.gmail.com.

[Viswanath Vittaladevaram]

Hi,

I understand that but could you please elaborate how to fix this issue?

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAJyMocv_kygySi5d4R9rWNpD0hBP3dKFf1Qkb_jdi_ncxwZHLA%40mail.gmail.com.

[Irfan Alibay]

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

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZhzVM%2BexfL_0RsfcFY0YgNw6UH7RjTCs_%2BzFuAU_zUWVw%40mail.gmail.com.

[Viswanath Vittaladevaram]

Okay. I confirm and can you show how it works in both ways!

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKVkbJgjb_ZwBx0JyPB%2B0F6_fMRv1iJXo55XwNoAhOSEFfhNkg%40mail.gmail.com.

[Irfan Alibay]

>   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

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZixrgXT4jo4mxp3Ti_eQ4_Xz%2BX98fJXwejC_%3Dyz7RU5gw%40mail.gmail.com.

[Viswanath Vittaladevaram]

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

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKVkbJhsn64%3D%3DWT-O7dzpYoopQKj%2BdSZjvY5Bx2eRKA24SwabA%40mail.gmail.com.

[Irfan Alibay]

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

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZg7qpuDGsBop-J9Bs1RsoudxZ879sixa6-BPfTktcv3rA%40mail.gmail.com.

[Viswanath Vittaladevaram]

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.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAJyMocsH3y%2Bs9YkEPU7G95NfcA1ReQKsOUVDCNkWV4gAxhgZ6A%40mail.gmail.com.

[Viswanath Vittaladevaram]

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]

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

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZi0U7T_mFbM6gik%3DTWYrj2QOUeHfr8CJMym4NJLCOSVCA%40mail.gmail.com.

[Viswanath Vittaladevaram]

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)

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAKVkbJihHRmJmwcrR-UT_qMBmpoaObF3qiBSBavyj%3DY8dVb%2B5w%40mail.gmail.com.

[Irfan Alibay]

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?

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/CAEQtsZiLio-7AuxdwRp50H-hkJyajgVtW7-W%3DDGvK3nZXdGXjw%40mail.gmail.com.

[Oliver Beckstein]

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]

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.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/dea2d4f3-2bc2-4279-bacc-960ea8e3f3f0n%40googlegroups.com.

[Irfan Alibay]

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]

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.

[Irfan Alibay]

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]

Hi,

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

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/a0202052-5e08-4831-afcb-4f45ea9712cbn%40googlegroups.com.

[Viswanath Vittaladevaram]

Hi,

It worked finally.Thank you.