Inconsistent results between GROMACS/2016.4 and GROMACS/2021.3

GROMACS version: 2016.4 and 2021.3
GROMACS modification: No

Dear GROMACS users and developers,
I am simulating a crystal structure using both GROMACS/2016.4 and GROMACS/2021.3, adopting the GAFF2 force field, at 100 K and 1 atm.

The results that I obtain using two different versions of GROMACS are remarkably different, more than one could expect. I attach below a plot of the potential energy as a function of the simulation time.
Results are related to a 10 ns simulation using Parrinello-Rahman barostat and anisotropic scaling, with tau-p = 5. I have previously equilibrated the system for 10 ns using Berendsen barostat (with the same tau-p) and it worked fine for both GROMACS versions.

The purple line is referred to 2016 version while the green one concerns the 2021 version.
In more detail, the 2016 version seems much more stable, also for what regards volume fluctuations (which I do not report here but I can share if useful).
For the sake of completeness, I observed the same issue simulating an amorphous drug structure at room temperature using OPLS force field and isotropic scaling with Parrinello-Rahman barostat (1 atm). GROMACS/2021 led to an unacceptable energy drift, which I did not observe with GROMACS/2016.

This is the input file for the simulation:

integrator = md
nsteps = 10000000
dt = 0.001

nstxout = 10000
nstvout = 0
nstenergy = 10000
nstlog = 10000

continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4

cutoff_scheme = Verlet
ns_type = grid
nstlist = 10
rcoulomb = 1.0
rvdw = 1.0
DispCorr = EnerPres

coulombtype = PME
pme_order = 4
fourierspacing = 0.16

tcoupl = V-rescale
tau_t = 1.0
tc-grps = System
ref_t = 100.0

pcoupl = Parrinello-Rahman
pcoupltype = anisotropic
tau_p = 5.0
ref_p = 1.0 1.0 1.0 0.0 0.0 0.0
compressibility = 4.5e-5 4.5e-5 4.5e-5 0.0 0.0 0.0
refcoord_scaling = com

pbc = xyz

Have you observed similar issues?
Thanks and best regards,
Tommaso