There is a way to add custom split normal use Python API?
Custom Split Normals
Elaborating on this answer and here
Test script. Run in object mode.  Switch to edit mode to see results.
import bpy

context = bpy.context
ob = context.object
me = ob.data
# visualize in edit mode    
me.show_normal_vertex = me.show_normal_loop = True
me.use_auto_smooth = True
# normal custom verts on each axis
me.normals_split_custom_set([(0, 0, 0) for l in me.loops])


end up with custom split normals in each of the corner directions. Blue are the vert normals
# all custom split normals pointing up.
normals = [(0, 0, 1) for v in me.vertices]
# make csn's all face up.  
me.normals_split_custom_set_from_vertices(normals)


all custom split normals facing up (0, 0, 1)
I am not aware if custom split normals can be accessed from bmesh.