Fitting a numerical integral via NonLinearModelFit to magnetic data

Question

I'm new to Mathematica and I'm currently trying to fit
$$m_T (H,T) = N_T \int\limits_0^{\infty} \frac{x k_\text{B} T}{\mu_0 H} \mathcal{L}(x) \text{pdf}(D_\text{mag}) \text{d}D_\text{mag}$$
with $x=\mu_0 M_S \pi D_\text{mag}^3 H/(6k_\text{B}T)$ and $\mathcal{L}(x)=\coth(x)-1/x$ (Langevin function)
to experimental data (this is a fit for the total magnetic moment of superparamagnetic nanoparticles).

So I tried the following:

kB = 1.38*10^-23; (*Boltzmann constant*)

pdf[d_?NumericQ, d0_?NumericQ, w_?NumericQ] = 
  1/(w*Sqrt[2 Pi])*1/d*E^(-1/(2*w^2)*(Log[d/d0])^2);

mt[b_?NumericQ, Nt_?NumericQ, Ms_?NumericQ, d0_?NumericQ, w_?NumericQ] := 
  Nt*NIntegrate[(Ms*Pi*d^3/6)*(Coth[b*Ms*Pi*d^3/(6*kB*300)] - 
    1/(b*Ms*Pi*d^3/(6*kB*300)))*pdf[d, d0, w], {d, 0, 1}];

As I think starting values are quite important here, I tried to find reasonable values and ended with

The blue dots are the experimental data and the red line is 1000*mt[0.0001*b, 0.35*10^14, 20150, 2.0246171576454905*10^-8, 0.3].

Note 1: The 1000*... and 0.0001*b are due to conversion of cgs units <-> SI units.

Note 2: I know that some of the nanoparticles exhibit ferromagnetic behaviour, thus the hysteresis, but for now I would be happy to just fit the superparamagnetic total magnetic moment to these points.

So I tried to fit mt to the data:

nlmB1 = 
  NonlinearModelFit[data, 
    {mt[b, Nt, Ms, d0, w], 
      {10^(15) > Nt > 10^(13) && 20200 > Ms > 20000 && 5*10^(-8) > d0 > 1*10^(-8) && 0.5 > w > 0.29}}, 
    {{Nt, 0.35*10^14}, {Ms, 20150}, {d0, 2.0246171576454905*10^-8}, {w, 0.3}}, 
    b]

But I get an error saying that There are no points that satisfy the constraints {}. So I removed most of the constraints and tried

nlmB1 = 
  NonlinearModelFit[data, 
    {mt[b, Nt, Ms, d0, w], {10^(15) > Nt && Nt > 10^(13)}}, 
    {{Nt, 0.35*10^14}, {Ms, 20150}, {d0, 2.0246171576454905*10^-8}, {w, 0.3}}, 
    b]

This time the fitting procedure takes forever and mathematica tells me

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in d near {d} = {2.07884*10^-8}. NIntegrate obtained 8.417808128264057*^-20 and 1.8814309233927094*^-16 for the integral and error estimates.
General::stop: "Further output of NIntegrate will be suppressed during this calculation."

What do I need to do to get Mathematica really fit this? I have also tried several Methods and I've changed the AccuracyGoal without getting a successful fit.

---

The data is included via

dataInput = Import["Measurement.dat", "Table", "FieldSeparators" -> ","];
data = dataInput[[All, {2, 4}]];

(If you want to try, here is the data: https://www.dropbox.com/s/z55iktrt0poki5y/Measurement.dat )

Answer 1

Following @Vitaliy's comment try this formulation:

kB = QuantityMagnitude@
    UnitConvert[Quantity["BoltzmannConstant"] , "Joule/Kelvin"]

pdf[dp_?NumericQ, d0p_?NumericQ, w_?NumericQ] = 
    kB^(-1/3)/(w Sqrt[2 Pi]) 1/dp *E^(-1/(2  w^2) (Log[dp /d0p])^2);
mt[b_?NumericQ, Nt_?NumericQ, Ms_?NumericQ, d0p_?NumericQ, w_?NumericQ] := 
    Nt kB^(4/3) Ms Pi /6  NIntegrate[  
        dp^3 (Coth[#] - 1/(#) &@(b Ms Pi dp^3/(6 300))) pdf[dp  , d0p, w],
          {dp, 0, Infinity}]

Note I've replaced d,d0 with dp,d0p which are normalized by kB^(1/3) (You should verify!!) Here is your fit plot, with the new d0p argument .844 ~ 2E-8/kB^(1/3)

Plot[ 1000*mt[0.0001*b, 0.35*10^14, 20150, .844 , 0.3] , {b, -10000, 10000}]

this could use some more work , but the essential issue here is your integrand was changing very rapidly for very small values of d. (The fact that setting the Infinite integration limit to 1 worked is a clue to that issue..). I didn't try the model fit but it should work much better. (note the integration limit now has to be infinity)

another thing..

that Coth[x]-1/x may be causing errors for small x (even though the limit behavior is well behaved and is x/3 + O(x^3). Try this:

Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), True}}]&

in place of:

Coth[#] - 1/(#)&