code: def overlap(alpha1,lmn1,A,alpha2,lmn2,B):
      
      l1,m1,n1 = lmn1
      l2,m2,n2 = lmn2
      rab2 = norm2(A-B)
      gamma = alpha1+alpha2
      P = gaussian_product_center(alpha1,A,alpha2,B)
  
      pre = pow(pi/gamma,1.5)*exp(-alpha1*alpha2*rab2/gamma)
  
      wx = overlap1d(l1,l2,P[0]-A[0],P[0]-B[0],gamma)
      wy = overlap1d(m1,m2,P[1]-A[1],P[1]-B[1],gamma)
      wz = overlap1d(n1,n2,P[2]-A[2],P[2]-B[2],gamma)
      return pre*wx*wy*wz
docs: Full form of the overlap integral. Taken from THO eq. 2.12
      >>> isclose(overlap(1,(0,0,0),array((0,0,0),'d'),1,(0,0,0),array((0,0,0),'d')),1.968701)
      True