ListContourPlot3D
has a VERY nasty default option; it joins the points creating a surface, not discrete 3D points and this feature cannot be removed. In order to deal with this we have to come up with a radical solution like the following:
First of all, let's read the data file
data = ReadList["data.out", Number, RecordLists -> True];
I suppose that the fourth column corresponds to the value of an index so, let's find out the minimum and the maximum values of that index in our sample data file
datat = Table[data[[i, 4]], {i, 1, Length[data]}];
min = Min[datat];
max = Max[datat];
Print["min = ", min]
Print["max = ", max]
Now let's paint the 3D points with colors according to the values of the index
valrange = {min, max};
data[[All, 4]] = Rescale[data[[All, 4]], valrange];
colfunc[x_, cf_] := ColorData[cf][1 - x[[4]]];
and then plot them
S0 = Graphics3D[{PointSize[0.004], Point[#[[1 ;; 3]], VertexColors -> colfunc[#, "Rainbow"]] & /@ data}]
If we like we can also create a custom-made colorbar
Clear[colorbar]
colorbar[{min_, max_}, colorFunction_: Automatic, divs_: 150] :=
DensityPlot[y, {x, 0, 0.1}, {y, min, max}, AspectRatio -> 10,
PlotRangePadding -> 0, PlotPoints -> {2, divs}, MaxRecursion -> 0,
Frame -> True, FrameLabel -> {{None, "index"}, {None, None}},
LabelStyle -> Directive[FontFamily -> "Helvetica", 17],
FrameTicks -> {{None, All}, {None, None}},
FrameTicksStyle -> Directive[FontFamily -> "Helvetica", 15, Plain],
ColorFunction -> colorFunction]
and combine the colorbar with the previous 3D plot
With[{opts = {ImageSize -> {Automatic, 600}}, cf = "Rainbow"},
Row[{Show[{S0}, Axes -> True,
AxesStyle -> Directive[FontSize -> 18, FontFamily -> "Helvetica"],
AxesLabel -> {"x", "y", "z"}, BoxRatios -> {1, 1, 1}, opts],
Show[colorbar[valrange, ColorData[cf][1 - #] &],
ImagePadding -> {{20, 60}, {55, 25}}, opts]}]]

ListInterpolation
$\endgroup$