The latest (?) paper on the subject is this, but see also some well-thought criticism here.

Assuming that you use the AMBER force field, preparing a simulation using HMR is a piece of cake : the implementation has been completely automated through the program 'parmed'. It goes like this :

* Prepare the system as you would normally do (AMBERTOOLS).

* Once you have the prmtop file, type :

# parmed original.prmtop 

                                  /^^^/           /]
                                 /   ]           / ]
                         _______/    ]___       /  ]
                        /                \_____/   /
                      _/   [@]  \ \                \
                     /..         | |                ]
                      VVVvvv\    | |         _/\    ]
          P A R M E D       |               /    \  ]
                     AAA^^^^^              /       \]
                      \_________\   \_____/
                                \    \
                                 \____\

ParmEd: a Parameter file Editor


Loaded Amber topology file original.prmtop
Reading input from STDIN...
> HMassRepartition
Repartitioning hydrogen masses to 3.024 daltons. Not changing water hydrogen masses.
> outparm new.prmtop
Outputting Amber topology file new.prmtop
> quit
Done!

That's it. The 'new.prmtop' file is ready for your HMR simulation with a documented 4fs timestep. This procedure will work fine for a pure water simulation, but if you are using a mixed solvent (for example TFE/water), you are alone : I do not have a clue whether you should repartition everything (and what the effect of increased viscosity would be). Better safe than sorry : don't do HMR with mixed solvents. There, I said it.

In the NAMD configuration file don't forget to change a few things in timesteps so that your DCD (and .temps file if you do adaptive tempering) make sense with a 4fs timestep. Here is a relevant diff (new parameters first) :

24c23
< adaptTempOutFreq        250
---
> adaptTempOutFreq        400
37c36
< dcdFreq                 250
---
> dcdFreq                 400
44,46c43,45
< outputEnergies          250
< outputTiming            1500
< xstFreq                 500
---
> outputEnergies          400
> outputTiming            1600
> xstFreq                 400
51c50
< timestep                4.0
---
> timestep                2.5

The steps above have been selected with long peptide folding runs in mind : it simplifies following the simulation (1,000,000 frames is 1 μs) but without becoming too coarse (1 frame is 1 picosecond).

This is not trivial, especially for folding simulations. To answer this you have to demonstrate that the systems (2.5fs vs 4.0fs) sample the same space with similar frequencies and with similar sampling efficiency. This is difficult to demonstrate. Having said that, because the gains are significant (something like 400 vs 640 nanoseconds per day for 5,000 atom system), we have performed a couple of simulations on established systems and compare the results as follows.

See this page for a more extended discussion and other possible solutions.

The αLa peptide is this one : http://dx.doi.org/10.1016/j.bpj.2011.08.044. The sampling (2fs vs 4fs) with respect to Good-Turing statistics appears to be comparable :

The secondary structure preferences are more difficult to analyze with any confidence due to limited sampling :

Extending the 4fs simulation to a total length of 3.5 μs, gave the following secondary structure analysis :

The peptide sequence is the one of the PDB entry 2n9n. This is a highly flexible peptide. The sampling (2.5fs vs 4fs) with respect to Good-Turing statistics appears to be comparable :

The secondary structure preferences also look comparable :