After studying the article and theses, as well as article Ground-state properties and elementary excitations of quantum droplets in dipolar Bose-Einstein condensates from which the ansatz was taken, I came to the conclusion that the ansatz may not have been calibrated to find the magnitude of the energy. Therefore, we calibrate the ansatz using the data in Fig. 1:
ClearAll["Global`*"] Idip[X_?NumericQ, Y_?NumericQ, u_?NumericQ] := Re[Exp[-u^2/2] - 3 X*Y/(1 - X^2)^(3/2) NIntegrate[ v^2 Exp[-u^2 X^2 v^2/(2 (1 - X^2) (1 - v^2))]/(Sqrt[ 1 - v^2] Sqrt[1 - v^2 (1 - Y^2)/(1 - X^2)]), {v, 0, Sqrt[1 - X^2]}]]; f[X_?NumericQ, Y_?NumericQ] := Re[Idip[X, Y, 0]]; Iqf[u_?NumericQ] := Re[2 Exp[-5 u^2/8]/Sqrt[Pi] NIntegrate[ Exp[-l^2] Cosh[Sqrt[2/5] u*l]^(5/2), {l, 0, Infinity}]]; ettot[x_, y_, z_, u_, wy_] := h*w[wy]/Kb (1/4 (1/x^2 + 1/y^2 + 1/z^2) + 1/4 (wx[wy]^2 x^2 (1 + u^2/2) + wy^2 y^2 + wz[wy]^2 z^2) + as[wy]*n/(2 Sqrt[2 Pi] (x*y*z)) (1 + Exp[-u^2/2]) + add[wy]*n/(2 Sqrt[2 Pi] (x*y*z)) (-f[x/z, y/z] - Idip[x/z, y/z, u]) + 512*Sqrt[ 2] as[wy]^(5/ 2) n^(3/2)/(75 Sqrt[5] Pi^(7/4) (x*y*z)^(3/2)) (1 + 3/2 add[wy]^2/as[wy]^2) Iqf[u]); a0 = 5.29*10^(-11);(*Bohr radius*)h = 1.054*10^(-34);(*Reduced Planck's constant*)M = 163.9*1.66*10^(-27);(*Dy-164 mass in kg*)wxx = 2 Pi*70;(*Experimental value of x-frequency*)wzz = 2 Pi*1000;(*Experimental value of z-frequency*) w[wy_] := (wxx*wy*wzz)^(1/3); a[wy_] := Sqrt[h/(M*w[wy])]; wx[wy_] := wxx/w[wy];(*Dimensionless x-frenquency*) wz[wy_] := wzz/w[wy];(*Dimensionless w-frenquency*) add[wy_] := 131*a0/a[wy];(*Dimensionless Dipolar length*) as[wy_] := 70*a0/a[wy];(*Dimensionless contact length*)n = 10^4; Kb = 1.38*10^(-23)*10^-9(*nK*);(*Boltzmann constant*)ddata1 = Table[minsol1 = FindMinimum[{ettot[x, y, z, 0, wy/w[wy]], x > 0 && y > 0 && z > 0}, {{x, 1.01}, {y, 1.012}, {z, 0.14}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, wy/(2 Pi), minsol1[[1]]} /. Last[minsol1]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; ddata2 = Table[ minsol2 = FindMinimum[{ettot[x, y, z, d/x, wy/w[wy]], x > 0 && y > 0 && z > 0 && d > x}, {{x, 1.01}, {y, 1.012}, {z, 0.14}, {d, 1.02}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, d, wy/(2 Pi), minsol2[[1]]} /. Last[minsol2]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; Compose a new ansatz $E_{total}=k E_{tot}+b$ and find the constants $b,k$
et1 = Interpolation[ddata1[[All, {4, 5}]]]; et2 = Interpolation[ddata2[[All, {5, 6}]]]; sol = NSolve[{b1 + k1 et1[200] == b2 + k2 et2[200], b1 + k1 et1[160] == 0, b2 + k2 et2[144] == 0, b2 + k2 et2[220] == 20}, {b1, b2, k1, k2}] Now we optimize a new ansatz
{b2, k2} = First[{b2, k2} /. sol] ddata2 = Table[ minsol2 = FindMinimum[{b2 + k2 ettot[x, y, z, d/x, wy/w[wy]], x > 0 && y > 0 && z > 0 && d > x}, {{x, 1.01}, {y, 1.012}, {z, 0.14}, {d, 1.02}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, d, wy/(2 Pi), minsol2[[1]]} /. Last[minsol2]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; {b1, k1} = First[{b1, k1} /. sol] ddata1 = Table[ minsol1 = FindMinimum[{b1 + k1 ettot[x, y, z, 0, wy/w[wy]], x > 0 && y > 0 && z > 0}, {{x, 1.01}, {y, 1.012}, {z, 0.14}}, Method -> Automatic, PrecisionGoal -> Automatic, AccuracyGoal -> Automatic]; Re[{x, y, z, wy/(2 Pi), minsol1[[1]]} /. Last[minsol1]], {wy, 2 Pi*100, 2 Pi*250, 2 Pi*20}]; Finally, build the curves and plot the data. It looks similar to Fig. 1, but there are differences in curvature. Note that you can not do re-optimization, because the data of re-optimization accurately lie on the curves obtained in the first optimization 