How to calculate gradients for the intermediate outputs using symbol?

Is it possible to compute gradients of intermediate values of a computation?
For example:-

a = mx.sym.var(name = 'a', shape = (1, 1), dtype = 'float64')
b = mx.sym.var(name = 'b', shape = (1, 1), dtype = 'float64')

c = mx.sym.broadcast_add(a, b, name = 'c')
d = mx.sym.make_loss(mx.sym.broadcast_add(c, b, name = 'd'))

bind = d.simple_bind(ctx = mx.cpu(1))
bind.forward(a = mx.nd.ones((1, 1)), b = mx.nd.ones((1, 1)), grad_req = {'a': 'write', 'b':'write', 'c':'write'})
outputs = bind.outputs[0]
# print(outputs)
bind.backward()
print(bind.grad_dict)
#prints 
  # {'a': 
  # [[1.]]
  # <NDArray 1x1 @cpu(1)>, 
  # 'b': 
  # [[2.]]
  # <NDArray 1x1 @cpu(1)>}

Here bind.grad_dict doesn’t provide gradients of c even though I’ve written 'c':'write' in grad_req .
So how can I get the gradients of “d with respect to c” as well?