Add support for HMR

[xiki-tempula]

This is something that could be done later on. It is quite often (in our lab) to do hydrogen mass repartition, where you move the mass of the carbon or oxygen to nearby hydrogens, such that one could use a timestep of 4 fs. I think this could be a function that one could use later on.

[selimsami]

Interesting, do you have some literature on this? Does this work better than simply using LINCS to constrain the hydrogens? (based on the timestep, it seems so...)

[xiki-tempula]

So using LINCS to constrain the hydrogens would only permit a time step of 2 fs but HMR would permit a time step of 4 fs so it is quite a significant acceleration.

The original paper is purposed by Adrian E. Roitberg https://pubs.acs.org/doi/10.1021/ct5010406

This kind of thing is commonly used in free energy calculations (based on the 400 citations that this paper has).

[selimsami]

Interesting. Would need to scale the mass-weighted Hessian for this, but other than that, it should be relatively simple? Or am I missing something?

[xiki-tempula]

@selimsami You do have a point that the mass-weight Hessian would be different. I think the "common" strategy is that they just scale the mass and leave the bond, angle and dihedral unchanged, so the potential energy is the same. Because the goal of this method is to slow down the oscillation of the C-H bond so that one could use a larger time step, if we keep the real Hessian, the oscillation won't get slowed down.

My suggestion is that the current Hessian routine should be the same and only apply the mass-scaling when eventually writing the topology. So the mass-weighted Hessian is fitted to reproduce the Hessian with the real mass not the scaled mass.

Like when one does

[ff]
HMR = False

Everything is the same.

When

[ff]
HMR = 4.032

mol_qforce_HMR.itp will be generated together with mol_qforce.itp. mol_qforce.itp is the current file and mol_qforce_HMR.itp is the same topology but only with mass rescaled. The HMR = 4.032 means the mass of the hydrogen will be scaled by 4 times.

The tricky thing is that the gromacs implementation scale the mass of the hydrogen by a factor of 4 but the amber implementation is to scale by a factor of 3. In practice, this makes no difference but would be good to allow the user to choose to scale it by a factor of 3 or a factor of 4.

The number 4.032 is taken to mimic the interface of https://parmed.github.io/ParmEd/html/parmed.html#hmassrepartition and https://ambermd.org/tutorials/basic/tutorial12/index.php. We could do 3 or 4 if you felt it is more clear.

[selimsami]

Hmm if it is just about changing the mass numbers on the output FF file, I am less excited about this :-D

You are then generating a high accuracy FF with Q-Force that matches to QM vibrational frequencies, then proceed to make them all wrong.

Maybe I am still missing something though. I'll have a more detailed look later.

The tricky thing is that the gromacs implementation scale the mass of the hydrogen by a factor of 4 but the amber implementation is to scale by a factor of 3. In practice, this makes no difference but would be good to allow the user to choose to scale it by a factor of 3 or a factor of 4.

Letting the user decide the value is not an issue. This could be done easily with an optional input like you suggested.

[xiki-tempula]

@selimsami This is just a suggestion and I do agree that HMR will screw the kinetics up. I could do the mass rescaling easily with parmed so other users should be able to do this by themselves if they want.

[selimsami]

Well, if there is a community who is interested in HMR, then it might be worth to have it just as an optional input, no harm done and maybe more users for Q-Force :)