Python program to calculate Kirkwood−Buff Integrals

[DHN]

I want to write a python program to calculate the Kirkwood−Buff Integrals.
The equation is as:
G=4π∫₀ ^∞ r ² [g(r)-1] dr
I have the g(r) values.

Any suggestions are highly appreciated.

Thank you.

 

[jedishrfu]

You will have to use a numerical integration scheme:

https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html

https://plot.ly/python/numerical-integration/

where you first calculate the values for the G function and then use custom code to integrate it or use the scipy package integrate() function or some variation thereof.

Please be aware that we can't write this program for you but can only provide hints and help along the way if you show us what you have and where you are stuck.

 

[DHN]

I didn't mean to write the code for me @ jedishrfu.
I have the don the following code:
A=open ('abc.dat','r')
C=open('def.dat','r')
D=open('ghi','w+')
f=4*pi
for row, col in izip (C, A):
cols=col.strip().split()
if(float(row)!=0):
x2 = ((float(cols[0])**2)*(float(row)-1))
print(integrate.quad(x2, 0, 4))
D.write(cols[0] + '\t' + (float(f)*(integrate.quad(x2,0,np.inf))))
D.close()
C.close()
A.close()

I get the following error as:
D.write(cols[0] + '\t' + (float(f)*(integrate.quad(x2,0,np.inf))))
File "/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.py", line 316, in quad points)
File "/usr/lib/python2.7/dist-packages/scipy/integrate/quadpack.py", line 383, in _quad return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
quadpack.error: quad: first argument is not callable
[Finished in 0.1s with exit code 1]

Any suggestions are highly appreciated.

Thank you.

 

[Ibix]

You can use code tags to allow indentation in the forum - critical in python. Put

 before the code and

after and you get:

A=open ('abc.dat','r')
C=open('def.dat','r')
D=open('ghi','w+')
f=4*pi
for row, col in izip (C, A):
    cols=col.strip().split()
    if(float(row)!=0):   
    x2 = ((float(cols[0])**2)*(float(row)-1))
    print(integrate.quad(x2, 0, 4))
    D.write(cols[0] + '\t' + (float(f)*(integrate.quad(x2,0,np.inf))))
D.close()
C.close()
A.close()

...which won't run because the indentation is invalid. My guess is that the two lines after the if should be further indented because the print of quad ought to fail too and it didn't, so presumably it didn't execute.

Did you read the first link @jedishrfu provided, which tells you what kind of thing the first argument to quad must be? Given the inputs you have, is quad an appropriate approach? Also, what output are you expecting from this? You will get a huge list of numbers, but $\int_0^\infty 4\pi r^2(g(r)-1)dr$ is a definite integral so should be a single number.

Finally, naming things in meaningless ways makes everyone's lives harder. Not correcting your code, just renaming your variables and files, I'd have written

rfile=open ('r.dat','r')
gfile=open('g.dat','r')
outfile=open('integral','w+')
f=4*pi
for g,r in izip (gfile, rfile):
    rs=r.strip().split()
    if(float(g)!=0):   
    x2 = ((float(rs[0])**2)*(float(g)-1))
    print(integrate.quad(x2, 0, 4))
    outfile.write(rs[0] + '\t' + (float(f)*(integrate.quad(x2,0,np.inf))))
outfile.close()
rfile.close()
gfile.close()

That way I can see at a glance where you are reading your r and g values from, and can more easily relate the code to the maths.

 

[DHN]

@lbix, I apologize for not making the code simpler as you have written.
Actually, i need to get a graph as:
https://imageshack.com/a/img922/2125/rC7rQL.png
The graph is r v/s G_αβ (in the image it is N_αβ). So from the values of r and g it should do the integration and give the values in which i plot the same r(x-axis) and v/s G_αβ on y-axis. It is not a single number.

 

[PeterDonis]

My guess is that the two lines after the if should be further indented

Assuming you only want that code to execute if the if statement's condition is true, yes. As the code stands now, you have an if statement with nothing after the condition, which is a syntax error; but the interpreter won't see anything wrong syntactically with the indentation of the rest of the statements (logically they probably should be indented, but syntactically they would be correct as they are).

 

[Ibix]

Assuming you only want that code to execute if the if statement's condition is true, yes.

I'm simply observing that neither call to integrate.quad should work, so if the error comes from the second one, as reported, the first must not execute for some reason. Thus those two lines must be in a conditional branch that the later code isn't (edit: although I'd need to check the scoping rules to determine if x2 even exists outside the if branch if it's defined there). I agree that I don't know why the program's set up that way.

@DHN - my problem is that the maths you've posted doesn't produce the output you want. And the code you've posted doesn't work. So you are asking me to guess at what you are trying to do, and guess at how you are trying to do it. This isn't going to get far.

 

[Ibix]

@PeterDonis - to be clear, the indentation in the code in my #4 is present in DHN's #3, but suppressed by the forum. It's not my interpretation - I just hit reply and added code tags.