To successfully investigate this problem theoretically in some of its complexity one needs to use methods that are capable of modelling the reaction step as a quantum chemical process. However, the large size of the molecule precludes the use of ab-initio and even semi-empirical methods. The molecular dynamics methods that have shown some degree of success in simulating the properties of biological macromolecules in solution, have force fields that are parametrised far from the region of bond breakage and formation. As a result, their use for studying enzyamtic reactions are theoretically unjustified. One way out is to devise hybrid methods: those which perform quantum calculations for the reaction region and classical dynamics for the rest of the molecule and connect the two in an acceptable manner. Such methods are in the process of implementation in various molecular dynamics packages (and in particular, in CHARMM). So far I have used this method to characterise the transition state and calculate the reaction free energy profile of a small model compound, an acyclic pentacoordinated oxyphosphorane species, in gas phase, in a sphere of water molecule and in the presence and absence of Mg2+ ions. The results more or less corroborate ab initio results in that the solvated dianionic species is seen to have a stable intermediate state while the gas phase results are inconclusive yet.