This is an awful bug in ListContourPlot3D that got into version 13.1 and should have been caught before release. This has been fixed in development builds, and the fix will be included in the next Mathematica release.

Please try this function as a workaround in the meantime:

CubePlot::usage = "CubePlot[file] gives a list of surface plots for the given cube file. Any option for MoleculePlot3D or ContourPlot3D can be given.";
CubePlot[cubeFile_, opts___ ? OptionQ] := Module[
    {mol, data, range, molplot},
    mol = Import[cubeFile, "Molecule"];
    data = mol["MetaInformation"]["VolumetricData", "Data"];
    range = Part[mol @ "MetaInformation", "VolumetricData", "DataRange"];
    If[!MatchQ[data, {__NumericArray}],
       Return[Missing @ "NotAvailable", Module]
    ];
    molplot = MoleculePlot3D[mol, FilterRules[{opts}, Options @ MoleculePlot3D]];
    Map[cubePlot[molplot, range, #, opts]&, data]
];
cubePlot[molplot_, range_, data_, opts___] := Module[
    {interpFunc, plotVars, surfaces},
    interpFunc = ListInterpolation[Transpose[Normal @ data, {3, 2, 1}], range];
    plotVars = MapThread[Prepend, {range, {x, y, z}}];
    surfaces = ContourPlot3D[interpFunc[x, y, z],
       Evaluate[Sequence @@ plotVars],
       Evaluate @ FilterRules[{opts}, Options @ ContourPlot3D],
       Contours -> {-0.03, 0.03},
       ContourStyle -> Opacity[0.8],
       Mesh -> None
    ];
    Show[molplot, surfaces]
]