import arcpy
from collections import defaultdict
from itertools import combinations

#Change the four lines below
arcpy.env.workspace = r'C:Default.gdb'
grid = 'fishnet_500m'
resultfield = 'meanpointdistance' #New field that will be added and populated with mean point distance
points = 'randompoints'

#Intersect grid with points
inter = 'grid_inter_points'
arcpy.Intersect_analysis(in_features=[grid,points], out_feature_class=inter, output_type='POLYGON')
fid_grid_field = 'FID_'+grid

#Group points together in a defaultdict(list) using fishnet id
allpoints = defaultdict(list)
with arcpy.da.SearchCursor(inter,[fid_grid_field,'SHAPE@']) as cursor:
    for row in cursor:
        allpoints[row[0]].append(row[1])

#Calculate distances between all point combinations and store in a dictionary
d = {}
for k,v in allpoints.items(): #d.iteritems in py2/ArcMap
    if len(v)>1:
        pointdistances = [c[0].distanceTo(c[1]) for c in combinations(v,2)]
        d[k] = sum(pointdistances)/len(pointdistances)

#Add new field and calculate as mean point distance
arcpy.AddField_management(in_table=grid, field_name=resultfield, field_type='DOUBLE')
with arcpy.da.UpdateCursor(grid,['OID@',resultfield]) as cursor:
    for row in cursor:
        if row[0] in d:
            row[1] = d[row[0]]
            cursor.updateRow(row)