How can I try to debug this discrepancy between Map and ParallelMap?

It is hard for me to present the actual set of functions that is prompting this question, because it is part of a very large code base. But to summarise, I have a setup like the following:

ClearAll[f, g1, g2];

f[mol_Molecule] := (* Transform mol and return a new molecule; f in pure *);

g1[lst : {__Molecule}] := Map[f, lst];

g2[lst : {__Molecule}] := ParallelMap[f, lst];


The argument lst usually has about 1,000 molecules, and the parallel version runs about 12 times faster on my computer. But while the function g1 always returns the correct result, the function g2 sometimes does not.

I understand that this is hardly enough detail. But I am quite baffled, and looking for pointers on what could cause such an error, and what I should try to inspect. The documentation clearly states that ParallelMap "...automatically distributes different applications of f among different kernels and processors," so I assumed that ParallelMap could be used as a drop-in replacement for Map.

flinty has brought up a very relevant point that I forgot to mention in the question. The function f is a pure function and does not mutate any variable.
• Is your f a pure function? If it has some state associated with it then ParallelMap might not give a result consistent with Map. Consider the following f[x_] := (y++; Return[x + y]); y = 0; Map[f, {1, 2, 3, 4}] then try: y = 0; ParallelMap[f, {1, 2, 3, 4}] and you'll probably get a different result. – flinty Nov 8 '20 at 16:28
• There could still be global state in a built-in function or package function somewhere. We'd need a minimal failing example f to test it. – flinty Nov 8 '20 at 16:31