I have tried to reproduce the problem by running millions of 39.1776 GeV proton collisions on Molybdenum target nucleus, with different versions of Geant4 (10.6.p02, 10.6.p03, 10.7), but I don't see any crash/problem.
Therefore, I cannot do anything, sorry about that!

Since the error persists in version 10.7.2, and eventually we want to launch a new production campaign, I come back to this.
In my case, it is not an issue of running million of proton events, a few is enough. It is also not related to the incoming proton (400 GeV in our case), it happens with particles produced when traveling through the thick target (1m). See below examples for anti-proton and for pion. I suspect that it has to do with using a thick target and maybe the energy loss associated with it. 

===> G4GeneralPhaseSpaceDecay::Pmx energy in cms < mass1+mass2
Call for FTFP
Target element Mo  Z= 42  A= 98
Unrecoverable error in the method ApplyYourself of anti_protonInelastic
TrackID= 5  ParentID= 0  anti_proton
Ekin(GeV)= 5.24208;  direction= (0.122321,-0.0544709,0.990995)
Position(mm)= (38.369,-16.5703,-3601.58);  material molybdenummisis
PhysicalVolume  <Target_9>
 ApplyYourself failed

*** Fatal Exception *** core dump ***
G4Track (0x4a98be20) - track ID = 5, parent ID = 0
 Particle type : anti_proton - creator process : not available
 Kinetic energy : 5.24208 GeV - Momentum direction : (0.122321,-0.0544709,0.990995)
 Step length : 1.70838 cm  - total energy deposit : 22.3057 MeV
 Pre-step point : (36.2973,-15.6412,-3618.51) - Physical volume : Target_9 (molybdenummisis)
 - defined by : Transportation - step status : 1
 Post-step point : (38.369,-16.5703,-3601.58) - Physical volume : Target_9 (molybdenummisis)
 - defined by : anti_protonInelastic - step status : 4
 *** Note: Step information might not be properly updated.

or

===> G4GeneralPhaseSpaceDecay::Pmx energy in cms < mass1+mass2
Call for FTFP
Target element Mo  Z= 42  A= 98
Unrecoverable error in the method ApplyYourself of pi-Inelastic
TrackID= 162299  ParentID= 4  pi-
Ekin(GeV)= 12.4731;  direction= (-0.0652905,-0.00295326,0.997862)
Position(mm)= (-3.88987,0.30987,-3590.59);  material molybdenummisis
PhysicalVolume  <Target_9>
 ApplyYourself failed

*** Fatal Exception *** core dump ***
G4Track (0x496ba428) - track ID = 162299, parent ID = 4
 Particle type : pi- - creator process : pi-Inelastic, creator model : Undefined
 Kinetic energy : 12.4731 GeV - Momentum direction : (-0.0652905,-0.00295326,0.997862)
 Step length : 1.4463 mm  - total energy deposit : 1.72405 MeV
 Pre-step point : (-3.79564,0.313972,-3592.03) - Physical volume : Target_9 (molybdenummisis)
 - defined by : hIoni - step status : 4
 Post-step point : (-3.88987,0.30987,-3590.59) - Physical volume : Target_9 (molybdenummisis)
 - defined by : pi-Inelastic - step status : 4

What mitigates the problem is commenting the line
throw G4HadronicException(__FILE__, __LINE__, "G4GeneralPhaseSpaceDecay::Pmx energy in cms < mass1+mass2"); in G4GeneralPhaseSpaceDecay::Pmx

this results in warning messages see below, but the program continues and says it will re-sample.


 pmx called with (129.785 139.570 139.570)
-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : had012
      issued by : G4HadronicProcess:CheckResult()
Warning: Bad energy non-conservation detected, will re-sample the interaction
 Process / Model: pi+Inelastic / QGSP
 Primary: pi+ (211),  E= 106419, target nucleus (42, 95)
 E(initial - final) = -77429.8 MeV.

*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

pmx called with (79.020 139.570 139.570)

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : had012
      issued by : G4HadronicProcess:CheckResult()
Warning: Bad energy non-conservation detected, will re-sample the interaction
 Process / Model: pi-Inelastic / FTFP
 Primary: pi- (-211),  E= 14663.1, target nucleus (74, 186)
 E(initial - final) = -8220.53 MeV.

After adding print statements and switching some debug options on, I think I understand the source of the problem:

G4GeneralPhaseSpaceDecayrho0 94.111 148.928 2 pi+ pi-   0x7fffffff5740 
or G4GeneralPhaseSpaceDecayrho0 91.561 148.928 2 pi+ pi-   0x7fffffff5740 

A rho0 has
 Mass [GeV/c2] : 0.77526     Width : 0.1491

If a Breit-Wigner function is used to generate the actual mass of a rho0, it is mathematically possible that one gets ~90MeV.
Have not figured out yet where the actual mass is calculated.

The "mistake" happens in G4VPartonStringModel::Scatter. Example below, a rho0 is being created with a mass of 42MeV which is below the mass of two pions. This eventually leads to a G4HadronicException and crash of the program. 

-----------------------Parton-String model is runnung ------------
Projectile Name Mass pi+ 139.5701
           Momentum  (0,0,158328.36286877;158328.42438589)
Target nucleus   A Z 98 42

Initial baryon number 98
Initial charge        43
-------------- Parton-String model:  Generation of strings -------

String No 1 (-551.17306902763,-44.557426694119,22775.305075267;23102.862753861) 3836.921817811 Partons   2          2101
String No 2 (239.02138938293,567.10154664467,77946.026070812;78442.353452211) 8788.6915693395 Partons   1          -1
String No 3 (393.76378872586,-233.07760525519,52786.29392673;53818.868801724) 10481.814608487 Partons   1          1103
String No 4 (-289.82541158578,-232.96028145553,4226.3812081322;4551.9975988902) 1649.2770421951 Partons   2          -1
String No 5 (-361.00449229046,-23.739475167759,24.019661435367;1007.0987121343) 939.56536 Partons   1          2103
String No 6 (120.25992969809,-293.05765167554,53.610994838555;992.97655681092) 939.56536 Partons   2          1103
String No 7 (-22.057448683393,309.2673224664,-134.86624959013;997.33448924899) 938.272013 Partons   2          2101
0 : rho0 (-146.50698340786,171.68827771744,3654.4884735792;3661.6986289055) 42.540386837166 775.26
...

G4HadronicException:  G4GeneralPhaseSpaceDecay::Pmx energy in cms < mass1+mass2"

One more step closer to the origin of the problem:

G4ExcitedStringDecay::FragmentString:

G4cout<<"Resonance *** "<<aTrack+1<<" min mass "<< BrW.GetMinimumMass( TrackDefinition ) <<" "
                     <<(*generatedKineticTracks)[aTrack]->GetDefinition()->GetParticleName()<<" "

Resonance *** 4 min mass 274.5467 rho- 
Resonance *** 1 min mass 274.5467 rho+ 


Resonance *** 3 min mass 1.02199782 rho0
Resonance *** 2 min mass 1.02199782 omega


G4SampleResonance::GetMinimumMass 
This loops over all decay channel of a resonance and calculates the minimum invariant mass adding up the daughter masses. 
By default, a rho0 in GEANT4 only decays to pi+pi-. But this is not good for our purpose, rho0 can also decay to mu+mu- and e+e-. I had modified the decay list of a rho0 to include all known channels. Therefore, MinimumMass is as low as 2*electron mass. But a rho0 of such low mass can of course not decay to pi+pi-. But this is not checked by the procedure and causes eventually a G4Exception.

After your detailed information, I was finally able to reproduce the problem and studied it.

The problem was due to the fact that rare, light decay channels - such as e+ e- - can push down the mass threshold of resonance decays, which can cause problems in the case of wide resonance - such as rho0 - to find a kinematically acceptable decay within a few attempts.

A simple, practical solution that I would propose is to consider only the decay channels whose branching ratio is above a certain threshold (e.g. 10%) to define the mass range of resonances - but letting all decay channels to be available when sampling the decay channel.
In the case of future, weird new resonances with several decay channels, but without any of them with a branching ratio above the threshold, the channel with the highest branching ratio is considered to determine the minimum mass of the resonance.

Note that before Geant4 10.4, the hadronic resonances produced by Geant4 string models were all set to their nominal mass. Only starting with Geant4 10.4, Breit-Wigner smearing was introduced (as required by the PANDA experiment).
This explains why you did not get any problem with Geant4 10.3.

The fix is in attachment: G4SampleResonance.cc
This should be copied in the directory:
    geant4/source/processes/hadronic/util/src/

This file can be used for any Geant4 version starting with 10.4, including any patch of 10.6 and 10.7.

Let me know if it works for you. If this is the case, and you are fine with the solution, it will be included in the next patch of 10.7 (i.e. 10.7.p03, whose release date is not yet scheduled) and version 11.0 (scheduled for December).

Many thanks!

Created attachment 706 [details]
Patch

The practical solution is in principle OK for our purpose. Only question, what happens for resonances with many decay channels and none above 10%?
Another issue could come up if somebody starts to modify the decay table and enhances rare decays for dedicated studies.

If there is no channel with branching ratio above 10%, then only the channel with the biggest branching ratio - whatever its value - is used to determine the minimum mass (i.e. the value returned by the method GetMinimumMass).

The current solution works also when rare decays are artificially enhanced if both of the following conditions hold:
1. the enhanced rare decays have branching ratios less than 10%
2. the decay with the highest branching ratio is a "natural" decay, i.e. not a rare decay which has been artificially enhanced.
I hope that these two restrictions can be accepted, else we need to find another, more complicated solution, e.g. excluding channels with leptons for determining the minimum mass...

Let me know.

Looks good!