Using python to create a curve and attach its endpoints with hooks to two spheres
Based on the answer linked by Seppo Silaste I have come up with the following example code which you can modify to fit your application:

import bpy

def hookCurve(o1, o2, scn):
    curve = bpy.data.curves.new("link", 'CURVE')
    spline = curve.splines.new('BEZIER')

    spline.bezier_points.add(1)
    p0 = spline.bezier_points[0]
    p1 = spline.bezier_points[1]
    p0.co = o1.location
    p0.handle_right_type = 'VECTOR'
    p1.co = o2.location
    p1.handle_left_type = 'VECTOR'

    obj = bpy.data.objects.new("link", curve)

    m0 = obj.modifiers.new("alpha", 'HOOK')
    m0.object = o1
    m1 = obj.modifiers.new("beta", 'HOOK')
    m1.object = o2

    scn.objects.link(obj)
    scn.objects.active = obj

    # using anything in bpy.ops is a giant pain in the butt
    bpy.ops.object.mode_set(mode='EDIT')

    # the mode_set() invalidated the pointers, so get fresh ones
    p0 = curve.splines[0].bezier_points[0]
    p1 = curve.splines[0].bezier_points[1]

    p0.select_control_point=True
    bpy.ops.object.hook_assign(modifier="alpha")

    p0.select_control_point = False
    p1.select_control_point = True
    bpy.ops.object.hook_assign(modifier="beta")

    return obj

o1 = bpy.data.objects['atom 1']
o2 = bpy.data.objects['atom 2']

obj = hookCurve(o1,o2, bpy.context.scene)