Applying scipy and weave to simple calculations on large sets

scipyweave_atoms.py

#First the standard python list style code:

#Calculate the distance between 2 atoms:
def dist(p1,p2):
    return sum([(a-b)*(a-b) for a,b in zip(p1,p2)])**0.5

#Calculate the angle between 3 atoms:
def ang(p1,p2,p3):
    P12=dist(p1,p2)
    P13=dist(p1,p3)
    P232=dist(p2,p3)**2
    x=(P12*P12 + P13*P13 - P232)/(2*P12*P13)
    if fabs(fabs(x) - 1.0) <= 1e-9:
        return 0.0
    return acos(x)


#Now in scipy with numpy.array and weave:

#Calculate the distance between 2 atoms:
distcode = """
double c=0.0;
for(int i=0;i<3;i++)
  c = c + (a[i]-b[i])*(a[i]-b[i]);
return_val=sqrt(c);
"""
def dist(a,b):
    return weave.inline(distcode,['a','b'])

#Calculate the angle between 3 atoms:
angcode = """
double ang,x,dab=0.0,dac=0.0,dbc=0.0;
for(int i=0;i<3;i++)
  dab = dab + (a[i]-b[i])*(a[i]-b[i]);
for(int i=0;i<3;i++)
  dac = dac + (a[i]-c[i])*(a[i]-c[i]);
for(int i=0;i<3;i++)
  dbc = dbc + (b[i]-c[i])*(b[i]-c[i]);
x=(dab + dac - dbc)/(2.0*sqrt(dab)*sqrt(dac));
if(int(fabs(x*1e9)) == int(1e-9))
  return_val=0.0;
else
  return_val=acos(x);
"""
def ang(a,b,c):
    return weave.inline(angcode,['a','b','c'])