Hello Meryem,
your expressions are indeed somewhat complicated, so I choose a much simpler example which may, at least I hope so, give you a hint.
First of all you should check the "mixed derivatives" of your gradient, that means check the "integrability conditions". That will show if your function f exists.
Now for the example.
Define
grad[f_] := {D[f, x], D[f, y], D[f, z]}
Let's say you have a gradient
gradf = {(2 x)/z, (2 y)/z, -((x^2 + y^2)/z^2)};
Check the mixed derivatives
D[gradf[[1]], y] == D[gradf[[2]], x]
D[gradf[[1]], z] == D[gradf[[3]], x]
D[gradf[[2]], z] == D[gradf[[3]], y]
Each answer is True, so there is indeed a function giving the gradient in question.
Now chose a straight line from { x0, y0, z0 } to { x1, y1, z1 } and the derivative of the line element to calculate
Integrate [ gradf . dr ]
along this line. Set
rule = {x -> x0 (1 - t) + t x1, y -> y0 (1 - t) + t y1, z -> z0 (1 - t) + t z1};
xp = D[{x, y, z} /. rule, t];
df = (gradf /. rule).xp // FullSimplify;
Switch off the output of conditions in Integrate and do the integration, change the variables back to x, y, z and get a function yielding the given gradient
SetOptions[Integrate, GenerateConditions -> False];
ff = Integrate[df, {t, 0, 1}] /. {x1 -> x, y1 -> y, z1 -> z} // FullSimplify
grad[ff]
Of course you may as well integrate along the edges of the cuboid:
fff = Integrate[(gradf[[1]] /. {x -> u, y -> y0, z -> z0}), {u, x0, x1}] +
Integrate[(gradf[[2]] /. {x -> x1, y -> u, z -> z0}), {u, y0, y1}] +
Integrate[(gradf[[3]] /. {x -> x1, y -> y1, z -> u}), {u, z0, z1}]
fff = (fff // FullSimplify) /. {x1 -> x, y1 -> y, z1 -> z}
grad[fff]
gradf
I hope you got an idea how to calculate your function f.