Everything depends on your needs and a lot of cofactors, including hardware you can use.

I believe more accurate would the question about choosing force field for your system(s) - a MD engine can be usually changed, but you should keep on mind that some FFs work better with some kind of systems - at least in case of limited computational powers.

 Thanks Dr. Rafal for the valuable infromation

Actually, I'm working with gromacs on simulating cytochrome p450-Drug interaction, but I'm facing troubles with finding a forcefield containing the parameters of heme-oxo Cpd I of cytochrome p450, so I think about switching to a software other than gromacs 

It's rather forcefiled issue than gromacs itself (which is a MD engine, NOT ff).

Anyway, I'm pretty sure that the heme-oxo should be somehow done - maybe not in standard gromos topologies, but simply don't stop googling (I'm not a gromos expert as I use mainly charmm ff).

Thanks a lot Dr. Rafal for your cooperation

As mentioned before, Gromacs can do this (as can Amber, lammps, namd, etc...)

I've used Gromacs for simulating lipids, helical proteins, and enzymes with metal centre. All of this comes down to picking the right forcefield for your needs, and in some cases adding or modifying existing forcefields.

For examples, I had to recalculate all of the partial charges, bond distances and bond angles (including respective force constants) for the amberffsb99 forcefield within Gromacs for the simulation of a matrix metalloproteinase with two zinc binding domains. 

Thanks Dr. Anthony,

Unfortunately, I'm suffering finding the heme-oxo cpd I (cytochrome p450) forcefield parameters compatible with gromacs 

If you find something compatible with Amber then it is also compatible with gromacs, it will be simply a case of converting the force field parameters into something which is readable by what ever MD engine you are using. 

Remember, gromacs is only a for-loop which integrates newtonian equation of motion. The same goes for Amber. 

Thanks again Dr. Anthony for your kind reply and cooperation

How should I proceed? My system contains protein and a single mg ion. The protein structure is modeled and so is Mg ion. I want to simulate the system in water and see the movements in Mg ion and to check if it is stable in particular region. I done a literature search which suggested use of QM/MM for which we need to define a bond and dummy atoms but in my case I am not sure to which atom I should define/restrain the ion as my protein and ion is modeled. Please suggest me

If any one of you can provide me some basic reference that would be extremly helpful.

I know this post is old, but if someone is searching for modeling Cpd I for heme systems, this can be done easily in AMBER by following this tutorial: http://ambermd.org/tutorials/advanced/tutorial20/mcpbpy_heme.htm