triple cross prouct of tensor

Question

Im trying to compute a triple cross product of vectors a,b, and c in real space and integrate over the entire space. The result is a term in the hamiltonian for an electronic system so there are indices for the band number too (you can think of a and b as gradient of wavefunctions in real space and of c as a vector field): $$ h_{nm} = \int \boldsymbol{a}_n(x)\times\boldsymbol{c}(x)\times\boldsymbol{b}_m(x)\,d x % $$ where $n$ and $m$ are indices from i.e. 0 to 50. I can write the triple product as $$ \boldsymbol{A}\times\boldsymbol{B}\times\boldsymbol{C} = \boldsymbol{B}(\boldsymbol{A}\cdot\boldsymbol{C})-\boldsymbol{C}(\boldsymbol{A}\cdot\boldsymbol{B}) % $$ which should be faster.

So for the first term in python I'm doing

import numpy as np

N = 50   #  band index
v = 3    # direction x,y,z
g = 30000 # real space grid points

a = np.ones((N,v,g), dtype = 'complex')
b = np.ones((N,v,g), dtype = 'complex')
c = np.ones((v,g), dtype = 'complex')
H_vnn = np.zeros((N,v,N), dtype = 'complex')

for i in range(g):
         H_vnn += np.array([c[v,i]*(np.dot(a[:,:,i],b[:,:,i].T)) for v in range(3)])

which, when doesnt run into memory error, takes forever. Im pretty sure this is not the most efficient way of doing it, but the only alternative I can think of is looping over the g index which seems to be even more expansive. Thank to everybody and particular to those who can help me

Answer 1

Try this

import numpy as np

N = 50   #  band index
v = 3    # direction x,y,z
g = 30000 # real space grid points

a = np.ones((g,N,v), dtype = 'complex')
b = np.ones((g,N,v), dtype = 'complex')
c = np.ones((g,v), dtype = 'complex')
H_vnn = np.zeros((N,v,N), dtype = 'complex')

for i in range(g):
    adotb  = np.dot(a[i,:,:],b[i,:,:].T)
    H_vnn += np.array([c[i,v]*adotb) for v in range(3)])