Modern approaches to nonlinear regression which are available in R

Question

I would like to fit a complex nonlinear regression model: basically, I have a complex computer code which has an input vector $\mathbf{x}$, a vector of calibration parameters $\boldsymbol{\theta}$ and produces an output $s=f(\mathbf{x},\boldsymbol{\theta})$. I also have iid observations $y$ for this output, which I model as contaminated by Gaussian, zero mean noise:

$$y = f(\mathbf{x},\boldsymbol{\theta})+\epsilon$$

Let $D=\{\mathbf{x}_i,y_i\}_{i=1}^N$ be the data set of my observations. I would like to find the vector $\boldsymbol{\theta}$ which maximizes the likelihood of observing $D$. Because of the assumption I made on the noise, maximizing the likelihood is equivalent to minimizing the residual sum of squares. This means which I can find $\boldsymbol{\theta}$ by means of Nonlinear Least Squares.

I need to use R for this project. In R, nls is available to implement NLS, which I think implements the Newton algorithm, but I've had convergence issues using it. Also, the interface is (for me) a bit awkward.

Are there other methods for Nonlinear Least Square regression, which are readily available in R, and have better convergence properties?

EDIT I already received great suggestions (nsl2, nlstools and nlmrt) which may hopefully solve my issues. Anyway, I'm providing some extra details on my case. $f(\mathbf{x},\boldsymbol{\theta})$ is a black box, i.e., it's the output of a computer code for which I don't have the source. Thus I can't provide you with an analytical, explicit expression for $f(\mathbf{x},\boldsymbol{\theta})$. I can tell you that the code computes damping and stiffness for a particular kind of turbomachinery seal, by numerically integrating the 1D isothermal compressible Navier-Stokes equations in an annular geometry. The code has been tested and validated multiple times in the course of years, so it's reliable.

Parametrization: to respond to @Roland's comment, the code is an executable with an input file, which contains many input values. Among these values, with the support of seal experts I chose some values to be predictors (vector $\mathbf{x}$) and some others to be calibration parameters (vector $\boldsymbol{\theta}$). Of course, since I don't have access to the code source, I cannot be sure that my parameters will be identifiable. But I think there must be a way to at least get an idea: I cannot open the code, but I can run it multiple times, for arbitrary values of $\mathbf{x}$ and $\boldsymbol{\theta}$ (as long as the code converges :) Can this knowledge be used to address @Roland's points? I.e.,

1. whether the data support the problem (I don't understand this point)
2. if I need to reparametrize (this is clear and I think it's related to identifiability)
3. if I have a partial linear problem (what's this?)

EDIT2 sadly I don't have access to gradients. It's a '70 code, and writing adjoint codes was not at all common in the engineering community at that time.