There are 2 kinds of NPT MD - with and without restriction on the shape of the computational box. If the software allows you, keep the shape constant. Apparently, you will get a metastable structure - since removing the shape restriction makes your crystal transform to another symmetry. You might get a new structure if your system is large eneough even at the restricted shape.  If you are absolutely sure that the ab Initio NPT MD should produce the same orthorhombic phase - check whether the results are converged on cut-off, k-points etc. It could be also that your volume is too far from the one you should get at these particular P and T and system moves into a new structure. May be the timestep is too large. May be a the mass of barostat and thermostat are too small or too large. There are quite a few ways to get wrong results.

If you run MD, you'll always lower symmetry because of (small or large) fluctuations. If the off dimensional elements are sufficiently small, and (after a long run) average to zero, they are simply just that - a fluctuation.

In my experience if the forcefield you are using is appropriate then MD will cause very little fluctuation in your overall crystal structure and symmetry. Unlike polymers/amorphous structures, crystal with its inherent symmetry is already in its most stable form and MD minimization and equilibration causes very minor changes to it. If you are still seeing large fluctuations it is most probably because of faulty simulation parameters. 

Dear Anatoly Belonoshko

Thank you very much for your suggestions. I am using the VASP simulation package which do not allow the NPT molecular dynamical simulation at constant shape.  I have checked the energy cuttoff, K-point and volume convergence and there is no problem in this. The time step of simulation is 1 fs which I think should be fine. I do not get much information regarding the mass of barostat and thermostat from available sources and I have taken it to be 

LANGEVIN_GAMMA = 5.0 

LANGEVIN_GAMMA_L = 5.0

PMASS= 20

But now I have checked that the off diagonal elements are small (0.1- 1%) when simulation is run for large time and final structure is obtained from the average over configurations during simulation time. 

However, I suspect that there could be some difficulty to study the phase transition ( like orthorhombic to monoclinic) using this technique. 

Thank you Kristof M. Bal & Aleena Alex for your reply. Yes I think the force field (pseudopotentials) are fine and the fluctuations are small and can be neglected as long as we are not dealing with the phase transitions occurring in the system.