The cluster variation method (CVM) coupled with the constrained simulated annealing (CSA) algorithm for global free energy minimization is proposed as an efficient approach for the study of alloy phase diagram. The chemical order–disorder transitions in Fe–Co alloy system between 300 and 1100K are investigated using the proposed CSA–CVM approach using the irregular tetrahedron approximation in a bcc lattice. The calculations have been carried out for a hypothetical paramagnetic Fe–Co alloy, i.e. with chemical interactions only, as well as for Fe–Co alloy with both magnetic and chemical interactions acting simultaneously, i.e., considering the quaternary Fe ↑ Fe ↓ Co ↑ Co ↓ system. The calculated A2/B2 boundary is found to match well with that determined experimentally. The incorporation of magnetic interaction is also found to have a large influence on the order–disorder transitions. The main advantages of the CSA–CVM approach is that it works efficiently with no initial value dependence on optimization parameters, and in contrast to some of the previously used minimization techniques, it is always found to converge to the global minimum. Another new aspect of the present investigation is that the phase diagram is obtained here by plotting long range order versus composition for the most ordered phase (B2 in this study) and not the commonly used grand potential versus chemical potential curves for various phases, thereby saving a lot of labor. Precise predictions for the variation of heat capacity with temperature have also been made.