Confusion about the unit of CMAPTorsionForce in force field file

Hi!

I find that the parameters of <CMAPTorsionForce> in force field files have the unit of kcal. That seems inconsistent with the default setting of openmm units described in this page.

Here is my experiment.

I have RNA CMAP parameters from an modified amber force filed in chamber format, which has the unit in kcal. I calculated the energy of an RNA molecule by tLEaP and it give me the result of -9.31, again in kcal/mol.

'cmap': -9.313757515318306
# this has unit kcal/mol

Then I move the numbers into <CMAPTorsionForce> in .xml file. Basically, this was done by two methods, Parmed and my self-implement. At this time, the unit in forcefield remains kcal. Then I used the following codes:

# label force groups
for i in range(system.getNumForces()):
    force = system.getForce(i)
    print(i, force)
    force.setForceGroup(i)
integrator = mm.LangevinIntegrator(300, 1/unit.picosecond,  0.0005*unit.picoseconds)
context = mm.Context(system, integrator)
context.setPositions(pdb.positions)
for i in range(system.getNumForces()):
    energy = context.getState(getEnergy=True, groups={i}).getPotentialEnergy()
    print(energy)

This gives me:

<openmm.openmm.CMAPTorsionForce; proxy of <Swig Object of type 'OpenMM::CMAPTorsionForce *' at 0x13694df00> >
-38.94414138793945 kJ/mol

which is exactly the same result as -9.313 kcal/mol.

I also used the inferface of Parmed to calculate the energy, which gave me the same result.

parmed.openmm.energy_decomposition(context)

So does this indicate that the unit of CMAP paramaters in .xml file is 'kcal'? All the other parameters are in kJ, so this totally confused me.

Since I can barely find any instruction about implement of CMAP in openmm, I hope to get some help here.

Thanks in advance!

openmm / openmm

Confusion about the unit of CMAPTorsionForce in force field file #3893