2010/7/2 Charles R Harris <charlesr.harris@gmail.com>

On Fri, Jul 2, 2010 at 12:15 PM, Nicolas Bigaouette <nbigaouette@gmail.com

...

wrote:

...

Hi all,

I don't really know where to ask, so here it is.

I was able to vectorize the normalization calculation in quantum mechanics: <phi|phi>. Basically it's a volume integral of a scalar field. Using:

...

norm = 0.0 for i in numpy.arange(len(dx)-1): for j in numpy.arange(len(dy)-1): for k in numpy.arange(len(dz)-1): norm += psi[k,j,i]**2 * dx[i] * dy[j] * dz[k]

if dead slow. I replaced that with:

...

norm = (psi**2 * dx*dy[:,numpy.newaxis]*dz[:,numpy.newaxis,numpy.newaxis]).sum()

which is almost instantanious.

I want to do the same for the calculation of the kinetic energy: <phi|p^2|phi>/2m. There is a laplacian in the volume integral which complicates things:

...

K = 0.0 for i in numpy.arange(len(dx)-1): for j in numpy.arange(len(dy)-1): for k in numpy.arange(len(dz)-1): K += -0.5 * m * phi[k,j,i] * ( (phi[k,j,i-1] - 2.0*phi[k,j,i] + phi[k,j,i+1]) / dx[i]**2 + (phi[k,j-1,i] - 2.0*phi[k,j,i] + phi[k,j+1,i]) / dy[j]**2 + (phi[k-1,j,i] - 2.0*phi[k,j,i] + phi[k+1,j,i]) / dz[k]**2 )

My question is, how would I vectorize such loops? I don't know how I would manage the "numpy.newaxis" code-foo with neighbours dependency... Any idea?

Thanx!

Wouldn't the Laplacian look better in momentum space? That is, do a 3d Fourier transform, multiply by the appropriate weight, and do the integral in the transform space.

Chuck

Thanx all for suggestions. The best approach was the convolution. I couldn't do a single convolution since the cell sizes are the same in x, y and z. I had to separate the laplacian into a sum over the three dimensions and convolve each individually. There is also terms appearing because of the non-uniform grid. Not that I work in atomic units, so m=1. I'm pasting the code I've finally written:

def Energy_Kinetic(phi, J1x, J1y, J1z, J2x, J2y, J2z):

kernel_x1 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], [[0,0,0], [1,-2,1], [0,0,0]], [[0,0,0], [0,0,0], [0,0,0]]]) kernel_y1 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], [[0,1,0], [0,-2,0], [0,1,0]], [[0,0,0], [0,0,0], [0,0,0]]]) kernel_z1 = numpy.array([[[0,0,0], [0,1,0], [0,0,0]], [[0,0,0], [0,-2,0], [0,0,0]], [[0,0,0], [0,1,0], [0,0,0]]])

kernel_x2 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], [[0,0,0], [1,0,-1], [0,0,0]], [[0,0,0], [0,0,0], [0,0,0]]]) kernel_y2 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], [[0,1,0], [0,0,0], [0,-1,0]], [[0,0,0], [0,0,0], [0,0,0]]]) kernel_z2 = numpy.array([[[0,0,0], [0,1,0], [0,0,0]], [[0,0,0], [0,0,0], [0,0,0]], [[0,0,0], [0,-1,0], [0,0,0]]])

convolution = ( \ scipy.signal.convolve(phi, kernel_x1, mode='same') / J1x**2

+ scipy.signal.convolve(phi, kernel_y1, mode='same') /

J1y[:,numpy.newaxis]**2 + scipy.signal.convolve(phi, kernel_z1, mode='same') / J1z[:,numpy.newaxis,numpy.newaxis]**2

+ 0.5*(J2x / J1x**3)*scipy.signal.convolve(phi, kernel_x2,

mode='same') \ + 0.5*(J2y[:,numpy.newaxis] / J1y[:,numpy.newaxis]**3)*scipy.signal.convolve(phi, kernel_y2, mode='same') \ + 0.5*(J2z[:,numpy.newaxis,numpy.newaxis] / J1z[:,numpy.newaxis,numpy.newaxis]**3)*scipy.signal.convolve(phi,kernel_z2, mode='same') \ )

K = -0.5 * (numpy.conjugate(phi) * convolution * J1x * J1y[:,numpy.newaxis] * J1z[:,numpy.newaxis,numpy.newaxis]).real.sum()

return K #