code: def V(a,b,C):
      
      if b.contracted:
          return sum(cb*V(pb,a,C) for (cb,pb) in b)
      elif a.contracted:
          return sum(ca*V(b,pa,C) for (ca,pa) in a)
      return a.norm*b.norm*nuclear_attraction(a.exponent,a.powers,a.origin,
                                              b.exponent,b.powers,b.origin,C)
docs: Simple interface to the nuclear attraction function.
      >>> from pyquante2 import pgbf,cgbf
      >>> s = pgbf(1)
      >>> isclose(V(s,s,(0,0,0)),-1.595769)
      True
  
      >>> sc = cgbf(exps=[1],coefs=[1])
      >>> isclose(V(sc,sc,(0,0,0)),-1.595769)
      True
  
      >>> sc = cgbf(exps=[1],coefs=[1])
      >>> isclose(V(sc,s,(0,0,0)),-1.595769)
      True
  
      >>> isclose(V(s,sc,(0,0,0)),-1.595769)
      True