How to use the chempy.cpv function in chempy

def centroid(selection='all', center=0, quiet=1):

    model = cmd.get_model(selection)
    nAtom = len(model.atom)

    centroid = cpv.get_null()

    for a in model.atom:
        centroid = cpv.add(centroid, a.coord)
    centroid = cpv.scale(centroid, 1. / nAtom)

    if not int(quiet):
        print ' centroid: [%8.3f,%8.3f,%8.3f]' % tuple(centroid)

    if int(center):
        cmd.alter_state(1, selection, "(x,y,z)=sub((x,y,z), centroid)",
                        space={'centroid': centroid, 'sub': cpv.sub})

    return centroid

raise Exception("Unexpected error!")

    if (nstates==1):
        if not quiet: print("Creating one state object!")

    if startframe > endframe:
        startframe, endframe = endframe, startframe
        if not quiet: print("Inverted start and end frames!")


    ########## BEGIN OF FUNCTIONAL SCRIPT ##########

    #normalize and get orthogonal vector

    # define vectors from points
    xyz2 = cpv.sub(xyz2, xyz1)
    xyz3 = cpv.sub(xyz3, xyz1)

    #NB! cpv.get_system2 outputs normalized vectors [x,y,z]
    xyz4 = cpv.get_system2(xyz2,xyz3)
    xyz2 = xyz4[0]
    xyz3 = xyz4[1]
    for x in range(0,3):
        for z in range(0,3):
            if x==z:
                continue
            if xyz4[x]==xyz4[z]:
                raise Exception("Illegal vector settings!")
    xyz4 = cpv.negate(xyz4[2]) #overwrites original

    # transform origin to corner
    if mode==0:

if method in ['bahar1996', '1', 1]:
        modelAll = cmd.get_model('(%s) and resn %s' % (selection, '+'.join(sidechaincenteratoms)))
        for at in modelAll.atom:
            if at.name in sidechaincenteratoms[at.resn]:
                atmap.setdefault((at.segi, at.chain, at.resn, at.resi), []).append(at)
    elif method in ['centroid', '2', 2]:
        modelAll = cmd.get_model('(%s) and polymer and not (hydro or name C+N+O)' % selection)
        for at in modelAll.atom:
            atmap.setdefault((at.segi, at.chain, at.resn, at.resi), []).append(at)
    else:
        print('Error: unknown method:', method)
        raise CmdException

    model = models.Indexed()
    for centeratoms in atmap.values():
        center = cpv.get_null()
        for at in centeratoms:
            center = cpv.add(center, at.coord)
        center = cpv.scale(center, 1./len(centeratoms))
        atom = Atom()
        atom.coord = center
        atom.index = model.nAtom + 1
        atom.name = name
        for key in ['resn','chain','resi','resi_number','hetatm','ss','segi']:
            atom.__dict__[key] = at.__dict__[key]
        model.add_atom(atom)
    model.update_index()
    if object in cmd.get_object_list():
        cmd.delete(object)
    cmd.load_model(model, object)
    return model

regions_iter = qdelaunay((a.coord for a in model.atom), 3, len(model.atom),
            qdelaunay_exe=qdelaunay_exe)
    edges = set(tuple(sorted([region[i-1], region[i]]))
            for region in regions_iter for i in range(len(region)))

    edgelist=[]
    r = []

    minco = 9999
    maxco = 0

    for edge in edges:
        ii, jj = edge
        a = model.atom[ii]
        b = model.atom[jj]
        co = cpv.distance(a.coord, b.coord)
        if cutoff > 0.0 and co > cutoff:
            continue
        if as_cgo:
            minco=min(co,minco)
            maxco=max(co,maxco)
            edgelist.append(a.coord + b.coord + [co])
        else:
            bnd = Bond()
            bnd.index = [ii, jj]
            model.add_bond(bnd)
        r.append((a,b,co))

    if not as_cgo:
        cmd.load_model(model, name, 1)
        return r

obj.extend([
        radius,
        0.8, 0.8, 0.8,
    ])
    obj.extend(color_list)
    if symmetric:
        start = cpv.sub(center, cpv.scale(direction, scale))
        obj.append(cgo.SAUSAGE)
        obj.extend(center)
        obj.extend(start)
        obj.extend([
            radius,
            0.8, 0.8, 0.8,
        ])
        obj.extend(color2_list)
    coneend = cpv.add(end, cpv.scale(direction, 4.0*radius))
    obj.append(cgo.CONE)
    obj.extend(end)
    obj.extend(coneend)
    obj.extend([
        radius * 1.75,
        0.0,
    ])
    obj.extend(color_list * 2)
    obj.extend([
        1.0, 1.0, # Caps
    ])
    cmd.load_cgo(obj, cmd.get_unused_name('oriVec'), zoom=0)

if 'C' in x and 'O' in x:
            vec_list.append(cpv.sub(x['O'], x['C']))
            count += 1
    if count == 0:
        print('warning: count == 0')
        raise CmdException
    vec = _vec_sum(vec_list)
    if count > 2 and sigma_cutoff > 0:
        angle_list = [cpv.get_angle(vec, x) for x in vec_list]
        angle_mu, angle_sigma = _mean_and_std(angle_list)
        vec_list = [vec_list[i] for i in range(len(vec_list))
                    if abs(angle_list[i] - angle_mu) < angle_sigma * sigma_cutoff]
        if not quiet:
            print('Dropping %d outlier(s)' % (len(angle_list) - len(vec_list)))
        vec = _vec_sum(vec_list)
    vec = cpv.normalize(vec)
    return _common_orientation(selection, vec, visualize, quiet)

def cgo_tetrahedron(x, y, z, r):
    vertices = [cpv.add((x, y, z), cpv.scale(v, r)) for v in [
        [0., 1., 0.],
        [0.0, -0.3338068592337708, 0.9426414910921784],
        [0.8163514779470693, -0.3338068592337708, -0.471320745546089],
        [-0.816351477947069, -0.3338068592337708, -0.4713207455460897]
    ]]
    return [
        cgo.BEGIN, cgo.TRIANGLES,
        cgo.NORMAL, 0.8165448970931916, 0.33317549135767066, 0.4714324161421696,
        cgo.VERTEX] + vertices[0] + [
        cgo.VERTEX] + vertices[1] + [
        cgo.VERTEX] + vertices[2] + [
        cgo.NORMAL, 0., 0.3331754913576707, -0.9428648322843389,
        cgo.VERTEX] + vertices[0] + [
        cgo.VERTEX] + vertices[2] + [
        cgo.VERTEX] + vertices[3] + [
        cgo.NORMAL, -0.8165448970931919, 0.3331754913576705, 0.4714324161421693,