Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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'])
Add Comment
Please, Sign In to add comment