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

_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion

Reply via email to