Here is some code to calculate the energy of the star particles in the galactic potential. This potential does not include the potential of the star cluster.  The velocities are in units of km/sec, distances kpc, masses 1x10^10 M_sun.  Newton's G is 43018 in these units and the time units is kpc/(km/sec) ~1.022 Gyr. 

    KE = 0.5*(vs3[:,0]**2 + vs3[:,1]**2 + vs3[:,2]**2)
    KEc = 0.5*(vc3[:,0]**2 + vc3[:,1]**2 + vc3[:,2]**2)

    rad2 = (rs3[:,0]**2 + rs3[:,1]**2 + rs3[:,2]**2)
    rad = np.sqrt(rad2)
    rad2c = (rc3[:,0]**2 + rc3[:,1]**2 + rc3[:,2]**2)
    radc = np.sqrt(rad2c)

    Grav = 43018.7
    NFW_M                     =92.2541
    NFW_Rs                    =20
    NFW_c                     =13.839
    NFW_Tb                    =0.97
    NFW_Tc                    =0.83
    MN_M                      =6.8
    MN_a                      =3.0
    MN_b                      =0.28
    MWbulge_M                 =0.5
    MWbulge_a                 =0.25

    P_plum = -Grav*MWbulge_M/np.sqrt(rad2 + MWbulge_a**2)
    ra = np.sqrt(rs3[:,0]**2 + (rs3[:,1]/NFW_Tb)**2 + (rs3[:,2]/NFW_Tc)**2)
    gc = 1/(math.log(1+NFW_c) - NFW_c/(1+NFW_c))
    P_NFW = -Grav*gc*NFW_M*np.log(1.0+ ra/NFW_Rs)/ra
    MHaloT = 0.3
    MHaloN = 1.037
    P_NFW *=  MHaloN*2/math.pi * math.atan((MHaloT*time)**2)
    Rc2 = (rs3[:,0]**2 + rs3[:,1]**2)
    P_MN = -Grav*MN_M/np.sqrt(Rc2 + (np.sqrt(rs3[:,2]**2 + MN_b**2) + MN_a)**2)
    E = KE + P_plum + P_MN + P_NFW

    P_plumc = -Grav*MWbulge_M/np.sqrt(rad2c + MWbulge_a**2)
    rac = np.sqrt(rc3[:,0]**2 + (rc3[:,1]/NFW_Tb)**2 + (rc3[:,2]/NFW_Tc)**2)
    #gc = 1/(math.log(1+NFW_c) - NFW_c/(1+NFW_c))
    P_NFWc = -Grav*gc*NFW_M*np.log(1.0+ rac/NFW_Rs)/rac
    MHaloT = 0.3
    MHaloN = 1.037
    P_NFWc *=  MHaloN*2/math.pi * math.atan((MHaloT*time)**2)
    Rc2c = (rc3[:,0]**2 + rc3[:,1]**2)
    P_MNc = -Grav*MN_M/np.sqrt(Rc2c + (np.sqrt(rc3[:,2]**2 + MN_b**2) + MN_a)**2)
    Ec = KEc + P_plumc + P_MNc + P_NFWc
    E -= Ec