The directories hold 3 similar sims
411: CDM
412: WDM 7keV
413: WDM 5.5keV
414: WDM 5.5keV low mass clusters
417: CDM low mass clusters

the 411,412,413 sims have star particles of 1 M_sun

simfile has all the particles, including dark. Quite big.
starfile has just the star particles
np.savez(simfile,time=time,rc3=rc3,rd3=rd3,rs3=rs3,vc3=vc3,vd3=vd3,vs3=vs3,ics=ics)
np.savez(starfile,time=time,rc3=rc3,rs3=rs3,vc3=vc3,vs3=vs3,ics=ics)

time: time in "Gyr" = 1.022 Gyr
rc3, vc3: positions (xyz) and velocities of cluster centers
rs3, vs3: pos and v of stars
ics: cluster of origin for every star
rd3, vd3: pos and vel of dark matter particles

Some example code that reads star positions and velocities and plots:
#
import numpy as np
import astropy as ast
import scipy as sc
import matplotlib.pyplot as plt

import numpy as np
import math
import astropy
import numpy.linalg

from tqdm import tqdm

dir = '411'
ifile = 1310
#simfile = dir + '/sim'+dir+'file' + f'{ifile:04}' + '.npz'
simfile = dir + '/stars'+dir+'file' + f'{ifile:04}' + '.npz'
npzfile= np.load(simfile)
print(npzfile.files)
time = npzfile['time']
rc3 = npzfile['rc3']
rs3 = npzfile['rs3']
vs3 = npzfile['vs3']
ics = npzfile['ics']
nc = len(rc3)
print(time,nc)


plt.rcParams["figure.figsize"] = (18,8)
fig, (ax1,ax2) = plt.subplots(1,2)
ax1.set_aspect(aspect=1)
rplot = 200
vplot = 300
ax1.axis([-rplot,rplot,-rplot,rplot])
ax1.set_xlabel("x [kpc]")
ax1.set_ylabel("y [kpc]")
ax2.set_xlim(0,rplot)
ax2.set_ylim(-vplot,vplot)
ax2.set_xlabel("x [kpc]")
ax2.set_ylabel("v_z [km/sec]")
for i in tqdm(range(nc)):
    pin = ics == i
    ax1.scatter(rs3[pin,0],rs3[pin,1],s=0.2,lw=0)
    strip = np.abs(rs3[:,1]) < rplot
    pin2 = pin & strip
    ax2.scatter(rs3[pin2,0],vs3[pin2,1],s=1,lw=0)


#np.savez(simfile,time=time,rc3=rc3,rs3=rs3,rd3=rd3,vc3=rd3,vs3=vs3,vd3=vd3,ics=ics)
dsm = np.linalg.norm(rs3,axis=1)
lon = np.arctan2(rs3[:,1],rs3[:,0])
lat = np.arcsin(rs3[:,2]/dsm)
hd = np.sqrt(1.+np.cos(lat)*np.cos(lon/2))
hx = 2.*math.sqrt(2)*np.cos(lat)*np.sin(lon/2)/hd
hy = math.sqrt(2)*np.sin(lat)/hd
dcm = np.linalg.norm(rc3,axis=1)

dmin = 10
dmax = 60
Lmin = 1000
plt.rcParams["figure.figsize"] = (18,9)
fig, (ax1) = plt.subplots()
ax1.axis([-math.pi,math.pi,-math.pi/2,math.pi/2])
ax1.set_xlabel('Hammer Azimuth',size='x-large')

for i in range(nc):
    dplot1 = dcm[i] > dmin
    dplot2 = dcm[i] < dmax
    if dplot1 & dplot2:
        pin = ics == i
        ax1.scatter(hx[pin],hy[pin],s=0.2,lw=0)