Some information that may be helpful:
- Berendsen barostat will not correctly sample the isothermal-isobaric (NPT) ensemble, but it does lead to fast relaxation of the box vectors to equilibrium.
- Parrinello-Rahman (P-R) does correctly sample the NPT ensemble. However, it does lead to large oscillations if far away from equilibrium as Michele stated.
Because of these facts, generally the Berendsen barostat is favored over P-R for equilibration purposes (which is why GROMACS warns you if P-R is used for equilibration). However, in my own experience, sometimes P-R can be used for NPT equilibration as long as the initial box size is reasonable relative to the equilibration density and gen-vel = yes is used to make sure that the acceleration of the atoms is not ridiculously large. In these circumstances (and after an energy minimization run), my simulations tend to be stable enough not to crash during equilibration even when using P-R.
However, if you want to simulate in NVE, then I would suggest an Energy Minimization (steepest descent and then possibly conjugate gradient afterwards) → NPT (tcoupl=v-rescale, pcoupl=parrinello-rahman) → NVT (tcoupl=v-rescale, pcoupl=no) → NVE (tcoupl=no, pcoupl=no) workflow. If you experience crashing when using P-R for NPT equilibration, then add another NPT equilibration run with Berendsen barostat before that. It is not recommended to go straight from NPT to NVE because too many changes are taking place simultaneously.