Your manual implementation of the binary cross entropy loss seems to work perfectly fine.
Here is a comparison using a dummy model:
model = nn.Sequential(
    nn.Linear(1, 1),
    nn.ReLU(),
    nn.Linear(1, 1)
)

x = torch.randn(10, 1)
target = torch.randint(0, 2, (10, 1)).float()
criterion = nn.BCEWithLogitsLoss()

output = model(x)

loss_pytorch = criterion(output, target)
print(loss_pytorch)
loss_pytorch.backward(retain_graph=True)
print([p.grad for p in model.parameters()])

model.zero_grad()

loss_manual = bce_loss(output, target)
print(loss_manual)
loss_manual.backward()
print([p.grad for p in model.parameters()])

If you run the code, you’ll see that both losses and the corresponding gradients are equal.
If you are summing both losses together, you might need to lower the learning rate or play around with some other hyperparameters.