I write the dimensions in the comments. Given:
z = torch.randn(7,5)   # i, d use torch.stack([list of z_i], 0) if you don't know how  to get this otherwise.
mu = torch.randn(6,5)  # j, d
nu = 1.2

you do
# I don't use norm. Norm is more memory-efficient, but possibly less numerically stable in backward
q_raw = (1 + ((z.unsqueeze(1)-mu.unsqueeze(0))**2).sum(2) / nu)**(-(nu+1) / 2)  # i, j
q_sum = q_raw.sum(1, keepdim=True) # i, 1 --> will be broadcast
q = q_raw / q_sum # i, j

p_raw = q**2 / q.sum(0, keepdim=True) # i, j
p_sum = p_raw.sum(1, keepdim=True)    # 1, j --> will be broadcast
p = p_raw / p_sum

kl_div = (p * (p.clamp(min=1e-7).log()  - q.clamp(min=1e-7).log())).sum()

Kindly send your invoicing address by PM. 
Best regards
Thomas