Regarding Angular distribution function (ADF's) in MDAnalysis

[cy16f0...@nitk.edu.in]

Hello,
I have a .gro file containing glycine+water (1 glycine,100 water molecules)mixture. Now i need to calculate the Angular distribution function of the water around the first hydration shell of N-terminal (NH3+) and the distance which is 0.35 nm (the value is taken from the RDF 1st minima).
The .gro file is:-
1GLY N 1 1.273 1.239 1.343
1GLY H1 2 1.239 1.214 1.436
1GLY H2 3 1.373 1.205 1.339
1GLY H3 4 1.223 1.196 1.266
1GLY CA 5 1.296 1.388 1.327
1GLY HA1 6 1.243 1.437 1.409
1GLY HA2 7 1.255 1.416 1.230
1GLY C 8 1.448 1.416 1.332
1GLY OT1 9 1.485 1.536 1.329
1GLY OT2 10 1.512 1.306 1.340
2SOL OW 1035 0.225 0.275 0.996
2SOL HW1 1036 0.260 0.258 1.088
2SOL HW2 1037 0.137 0.230 0.984
3SOL OW 1038 1.755 0.607 0.231
3SOL HW1 1039 1.743 0.594 0.132
3SOL HW2 1040 1.725 0.526 0.280
4SOL OW 1041 0.850 0.798 1.823
4SOL HW1 1042 0.846 0.874 1.888
4SOL HW2 1043 0.872 0.834 1.732
5SOL OW 1044 0.685 1.012 0.665
5SOL HW1 1045 0.754 0.996 0.735
5SOL HW2 1046 0.612 1.069 0.703
6SOL OW 1047 1.603 0.447 0.737
6SOL HW1 1048 1.529 0.493 0.687
6SOL HW2 1049 1.654 0.515 0.790
7SOL OW 1050 0.231 1.713 0.483
7SOL HW1 1051 0.265 1.790 0.537
7SOL HW2 1052 0.275 1.713 0.393
8SOL OW 1053 1.127 1.341 1.690
8SOL HW1 1054 1.174 1.341 1.778
8SOL HW2 1055 1.079 1.254 1.679
.
.
.

And the code that i have tried is -
import MDAnalysis
from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD
from matplotlib import pyplot as plt

u = MDAnalysis.Universe('glywat.gro', 'nptmd.trr')
selection = "byres name OW and sphzone 0.35 (protein and (resid 42 or resid 26) )"
bins = 30
AD_analysis = AD(u,selection,bins)
AD_analysis.run()
#now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector ,
#the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns
#are P(cos(theta)) vs cos(theta) for dipole vector
#for bin in range(bins):
print("{AD_analysisOH} {AD_analysisHH}
{AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], #AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin]))

#and if we want to graph our results
plt.figure(1,figsize=(18, 6))

#AD OH
plt.subplot(131)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for OH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[0][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[0][:-1]])

#AD HH
plt.subplot(132)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for HH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[1][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[1][:-1]])

#AD dip
plt.subplot(133)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for dipole')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[2][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[2][:-1]])

plt.show()

But i am getting error as:- "print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin]))

SyntaxError: keyword can't be an expression".

Any suggestions are highly appreciated

[Oliver Beckstein]

I'm not 100% sure if the following fixes your problem but try as the first line of your script

from __future__ import print_function

*Maybe* that was the issue. If not I would run the script manually line by line in ipython and make sure every line executes. Once you get to one that doesn't then work on this line. The reason why I'm suggesting this is because sometimes one forgets a parenthesis somewhere and then Python tries to include whole lines as expressions in previous lines. 

--

Oliver Beckstein

email: orbe...@gmail.com

[Naveen Michaud-Agrawal]

Hi,

Your issue is within the format call in the print function:

print("{AD_analysisOH} {AD_analysisHH}
{AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], #AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin]))

You have something that looks like: "string".format(a.b = something, ...). The issue is that keyword arguments to a function have to be a name, not an attribute. Try this:

print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysisOH=AD_analysis.graph[0][bin], AD_analysisHH=AD_analysis.graph[1][bin], AD_analysisHH=AD_analysis.graph[2][bin]))

--

-----------------------------------
Naveen Michaud-Agrawal

[cy16f0...@nitk.edu.in]

Hello,

I treid with the following command also but now i am getting error as:-

File "<ipython-input-78-a8fc0bc674a2>", line 15

    print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysisOH=AD_analysis.graph[0][bin], AD_analysisHH=AD_analysis.graph[1][bin], AD_analysisHH=AD_analysis.graph[2][bin]))

        ^
IndentationError: expected an indented block

1] I dont find any source here...

2] Or should i define the vectors/dipole ..?? or is this by default taking the correct requirements..??

Any suggestions..??

[Oliver Beckstein]

Not an MDAnalysis error. This is a Python error (or more precisely: a programming error). The error message tells you what happened. 

Have you worked through an introductory Python tutorial? If yes: work through it again until you understand why this is an error in Python. 

Quick answer: Remove leading white space. 

--

Oliver Beckstein

email: orbe...@gmail.com

--
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 post to this group, send email to mdnalysis-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mdnalysis-discussion.
For more options, visit https://groups.google.com/d/optout.

[cy16f0...@nitk.edu.in]

I have tried with the following modifications...

import MDAnalysis

from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD

from matplotlib import pyplot as plt

u = MDAnalysis.Universe('glyammw.gro', 'nptmd.trr')

selection = "byres name OW and sphzone 0.35 (resid 42 or resid 26)"

bins = 100

AD_analysis = AD(u,selection,bins)

AD_analysis.run()

#now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector ,

#the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns

#are P(cos(theta)) vs cos(theta) for dipole vector

for bin in range(bins):

      

    print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysisOH=AD_analysis.graph[0][bin],AD_analysisHH=AD_analysis.graph[1][bin], AD_analysisDip=AD_analysis.graph[2][bin]))

#and if we want to graph our results

plt.figure(1,figsize=(18, 6))

#AD OH

plt.subplot(131)

plt.xlabel('cos theta')

plt.ylabel('P(cos theta)')

plt.title('PDF cos theta for OH')

plt.plot([float(column.split()[0]) for column in AD_analysis.graph[0][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[0][:-1]])

#AD HH

plt.subplot(132)

plt.xlabel('cos theta')

plt.ylabel('P(cos theta)')

plt.title('PDF cos theta for HH')

plt.plot([float(column.split()[0]) for column in AD_analysis.graph[1][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[1][:-1]])

#AD dip

plt.subplot(133)

plt.xlabel('cos theta')

plt.ylabel('P(cos theta)')

plt.title('PDF cos theta for dipole')

plt.plot([float(column.split()[0]) for column in AD_analysis.graph[2][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[2][:-1]])

plt.show()

And the images of the graphs are attached.

What are the inferences that i can make from this.. or any changes that i need to make..??

Thank you.

[Oliver Beckstein]

On Dec 5, 2017, at 5:33 AM, cy16f0...@nitk.edu.in wrote:

What are the inferences that i can make from this.. or any changes that i need to make..??

That is YOUR job as a researcher.

If you run analysis then we assume you know why you’re doing this and you should also have some sense of what to expect. People on this mailing list are not going to do your research for you. We help with MDAnalysis and related technical issues. Read the papers that used this method, discuss with your PI or with colleagues that are familiar with the research. Present at conferences and discuss with other researchers there.

--

Oliver Beckstein, DPhil * oliver.b...@asu.edu
http://becksteinlab.physics.asu.edu/

Assistant Professor of Physics
Arizona State University
Center for Biological Physics and Department of Physics
Tempe, AZ 85287-1504
USA

Office: PSF 348
Phone: +1 (480) 727-9765
FAX: +1 (480) 965-4669

Department of Physics: https://physics.asu.edu/content/oliver-beckstein
Center for Biological Physics: https://cbp.asu.edu/content/oliver-beckstein

[cy16f0...@nitk.edu.in]

Yes, 

But as in the literature's the ADF's/orientation profile/cos theta distribution functions shows as in the following images...

So what are the modifications that i have to do in the code..??

Do i have to define any other things..?

Suppose i want to calculate the orientation profile ie., cos theta v/s Pcos theta of NH3+ of glycine around 1st solvation shell (taken from rdf and it is 0.35 nm) w.r.t water (ie., water molecules surrounded around NH3+ of glycine throughout the simulation).

[cy16f0...@nitk.edu.in]

1] In the manual, it is given as " vector ˆn parallel to a chosen axis (z is the default value)". So should i define the 'n' vector..?? or how is it choosing the z-axis vector.

2] should i consider the shell (sphzone 0.35 ) or bulk is also good enough..??

Since i am getting the same type of graph i am indeed asking this questions...

Any suggestions, please.

[cy16f0...@nitk.edu.in]

And in the code if i give as:

 

u = MDAnalysis.Universe('nptmd.tpr', 'nptmd.xtc')

selection = "name OW and resname GLY"

bins = 30

AD_analysis = AD(u,selection,bins)

AD_analysis.run()

for bin in range(bins):

      

    print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysisOH=AD_analysis.graph[0][bin],AD_analysisHH=AD_analysis.graph[1][bin], AD_analysisDip=AD_analysis.graph[2][bin]))

#and if we want to graph our results

plt.figure(1,figsize=(18, 6))...

it gives the output as:

0.0 nan 0.0 nan 0.0 nan
0.01 nan 0.01 nan 0.01 nan
0.02 nan 0.02 nan 0.02 nan
0.03 nan 0.03 nan 0.03 nan
0.04 nan 0.04 nan 0.04 nan
0.05 nan 0.05 nan 0.05 nan
0.06 nan 0.06 nan 0.06 nan
0.07 nan 0.07 nan 0.07 nan
0.08 nan 0.08 nan 0.08 nan
0.09 nan 0.09 nan 0.09 nan
0.1 nan 0.1 nan 0.1 nan
0.11 nan 0.11 nan 0.11 nan
0.12 nan 0.12 nan 0.12 nan
0.13 nan 0.13 nan 0.13 nan
0.14 nan 0.14 nan 0.14 nan

So what does this nan means..?? 

Thank you.