[Numpy-discussion] huge array calculation speed
Dan Lussier
dtlussier@gmail....
Thu Jul 10 11:38:43 CDT 2008
Hello,
I am relatively new to numpy and am having trouble with the speed of
a specific array based calculation that I'm trying to do.
What I'm trying to do is to calculate the total total potential
energy and coordination number of each atom within a relatively large
simulation. Each atom is at a position (x,y,z) given by a row in a
large array (approximately 1e6 by 3) and presently I have no
information about its nearest neighbours so each its position must be
checked against all others before cutting the list down prior to
calculating the energy.
My current numpy code is below and in it I have tried to push as much
of the work for this computation into compiled C (ideally) of numpy.
However it is still taking a very long time at approx 1 second per
row. At this speed even some simple speed ups like weave won't
really help.
Are there any other numpy ways that it would be possible to calculate
this, or should I start looking at going to C/C++?
I am also left wondering if there is something wrong with my
installation of numpy (i.e. not making use of lapack/ATLAS). Other
than that if it is relevant - I am running 32 bit x86 Centos 5.1
linux on a dual Xeon 3.0 GHz Dell tower with 4 GB of memory.
Any suggestions of things to try would be great.
Dan
========================================================================
=======================
def epairANDcoord(xyz,cutoff=4.0):
box_length = 160.0
cut2 = cutoff**2
ljTotal = numpy.zeros(xyz.shape[0])
coord = numpy.zeros(xyz.shape[0])
i = 0
for r0 in xyz:
r2 = r0-xyz
# application of minimum image convention so each atom interacts
with nearest periodic image w/i box_length/2
r2 -= box_length*numpy.rint(r2/box_length)
r2 = numpy.power(r2,2).sum(axis=1)
r2 = numpy.extract(r2<cut2, r2)
coord[i] = r2.shape[0]-1
r2 = numpy.extract(r2 != 0, r2) #
avoid problems of inf when considering own atom.
r2 = numpy.divide(1.0,r2)
ljTotal[i] = numpy.multiply(2,(numpy.power(r2,6)-
numpy.power(r2,3))).sum(axis=0)
# ljTotal[i] = ljPerAtom.sum(axis=0)
i += 1
del r2
# if i>25: break # used to cut down the length of
calculation as needed
return ljTotal,coord
More information about the Numpy-discussion
mailing list