docs: def nuclear_attraction(alpha1,lmn1,A,alpha2,lmn2,B,C):
      """
      Full form of the nuclear attraction integral
      >>> isclose(nuclear_attraction(1,(0,0,0),array((0,0,0),'d'),1,(0,0,0),array((0,0,0),'d'),array((0,0,0),'d')),-3.141593)
      True
      """
      l1,m1,n1 = lmn1
      l2,m2,n2 = lmn2
      gamma = alpha1+alpha2
  
      P = gaussian_product_center(alpha1,A,alpha2,B)
      rab2 = norm2(A-B)
      rcp2 = norm2(C-P)
  
      dPA = P-A
      dPB = P-B
      dPC = P-C
  
      Ax = A_array(l1,l2,dPA[0],dPB[0],dPC[0],gamma)
      Ay = A_array(m1,m2,dPA[1],dPB[1],dPC[1],gamma)
      Az = A_array(n1,n2,dPA[2],dPB[2],dPC[2],gamma)
  
      total = 0.
      for I in range(l1+l2+1):
          for J in range(m1+m2+1):
              for K in range(n1+n2+1):
                  total += Ax[I]*Ay[J]*Az[K]*Fgamma(I+J+K,rcp2*gamma)
                  
      val= -2*pi/gamma*exp(-alpha1*alpha2*rab2/gamma)*total
      return val