Dear yt users,

I am running a binary star simulation with enzo, and in my box there is a star made up of a gas envelope and a particle as core, and a secondary star represented by a single particle (no gas).

I would like to plot the equipotential surfaces of the system on a slice plot, hence, since the default Grav_Potential field does not include the particles potential, I created a custom made field as follows:

  def _Total_Grav_Potential(field,data): 
#G in cgs
grav_const_cgs = 6.674e-8
#data for gas
gas_grav_potential_cgs = data['Grav_Potential']  * (length_unit1/time_unit1)**2.0 #in units of velocity^2
gas_position_x_cgs = data['x'] * length_unit1
gas_position_y_cgs = data['y'] * length_unit1
gas_position_z_cgs = data['z'] * length_unit1
#initialize the array for the total grav potential to the gas grav potential
total_grav_potential_cgs = gas_grav_potential_cgs
#loops on the particles
for i in range(len(data['ParticleMassMsun'])):
#data of the current particle
current_particle_mass_cgs= data['ParticleMassMsun'][i] * mass_unit
current_particle_position_x_cgs = data['particle_position_x'][i] * length_unit1
current_particle_position_y_cgs = data['particle_position_y'][i] * length_unit1
current_particle_position_z_cgs = data['particle_position_z'][i] * length_unit1
#computes the array of grav potential for the current particle
current_particle_grav_potential_cgs = -(grav_const_cgs * current_particle_mass_cgs)/(((gas_position_x_cgs - current_particle_position_x_cgs)**2.0 + (gas_position_y_cgs - current_particle_position_y_cgs)**2.0 + (gas_position_z_cgs - current_particle_position_z_cgs)**2.0)**0.5)
#adds the current particle potential to the total potential
total_grav_potential_cgs = total_grav_potential_cgs + current_particle_grav_potential_cgs
return total_grav_potential_cgs
#adds the field to the available yt fields
add_field("Total_Grav_Potential", function=_Total_Grav_Potential, take_log=True ,units=r"\rm{erg}/\rm{g}")

where the length unit and time unit variables have the function to convert my quantities in cgs.
If I then call my "Total_Grav_Potential" field, I get reasonable values, but when I try to inspect the shape of the equipotential surfaces on the xy-plane by plotting the following:

sp = SlicePlot(pf, 'z', "Density", width = 1.0)
sp.annotate_velocity(factor=16, normalize=True)
sp.annotate_particles(1.0, p_size=50.0, marker='o', col='black')
sp.annotate_text((0.7,1.05), "time = "+str(pf.current_time)+"yr")"density_slice_z_"+str(index)+".png")

nothing happens, and the plot does not show any additional contour; while if I do the same with the default "Grav_Potential" field, the gas equipotential contours get plotted without problems. 
I also tried to play with the "annotate_contour" options "ncont", "factor" and "clim" to refine my plot but the result is always the same.
Additionally, if I do not set the "clim" option the following error shows up:

ValueError: zero-size array to minimum.reduce without identity

I would like to ask if my custom definition of the Total_Grav_Potential field is correct from a yt coding point of view and if I should modify it (or if there is any improvement I can do).
Also, I would like to ask if you have any idea of what could be the problem with the plotting.

Thanks for the help,