wk_minimize_trajectory.py fails: kcal not accepted unit by openmm

[blakemertz]

The minimize step for running waterkit fails because openmm will not read in units of energy in kcal/mol:

python /media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py -p protein_system.prmtop -t protein_system.nc -o protein_min.nc
Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.
Traceback (most recent call last):
  File "/media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py", line 145, in <module>
    main()
  File "/media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py", line 141, in main
    m.minimize_trajectory(prmtop_filename, traj_filename, output_filename)
  File "/media/bak11/binaries/git/waterkit/scripts/wk_minimize_trajectory.py", line 110, in minimize_trajectory
    simulation.minimizeEnergy(maxIterations=self._n_steps, tolerance=tolerance)
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/app/simulation.py", line 143, in minimizeEnergy
    mm.LocalEnergyMinimizer.minimize(self.context, tolerance, maxIterations, reporter)
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/openmm.py", line 4422, in minimize
    tolerance = tolerance.value_in_unit(unit.kilojoules_per_mole/unit.nanometer)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/unit/quantity.py", line 626, in value_in_unit
    val = self.in_units_of(unit)
          ^^^^^^^^^^^^^^^^^^^^^^
  File "/media/bak11/binaries/miniconda3/envs/waterkit/lib/python3.11/site-packages/openmm/unit/quantity.py", line 662, in in_units_of
    raise TypeError('Unit "%s" is not compatible with Unit "%s".' % (self.unit, other_unit))
TypeError: Unit "kilocalorie/mole" is not compatible with Unit "kilojoule/(nanometer*mole)".

Looks like the culprit is in the import calls at the beginning of wk_miniimize_trajectory.py (simtk is deprecated in openmm 7/8 btw):

     17 from simtk.unit import Quantity, picoseconds, kilocalories_per_mole, angstroms, kelvin
     18 from simtk.openmm import CustomExternalForce, LangevinIntegrator, Platform
     19 from simtk.openmm.app import AmberPrmtopFile, HBonds, CutoffNonPeriodic, Simulation

and here is the code snippet from wk_minimize_trajectory.py that is causing it to fail:

     63     def minimize_trajectory(self, prmtop_filename, traj_filename, output_filename):
     64         nonbondedMethod = CutoffNonPeriodic
     65         nonbondedCutoff = 9 * angstroms
     66         rigidWater = True
     67         constraints = HBonds
     68         dt = 0.002 * picoseconds
     69         temperature = 300 * kelvin
     70         friction = 1.0 / picoseconds
     71         tolerance = 1 * kilocalories_per_mole
     72         K = self._restraint * kilocalories_per_mole / angstroms**2

Should be a simple enough fix -- import kilojoules per mol instead of kcal (http://docs.openmm.org/latest/userguide/theory/01_introduction.html) and change the minimize_trajectory function for tolerance and K to use kJ/mol instead of kcal/mol. However, are there any other steps in the waterkit workflow that would be affected that utilize kcal/mol? wk_minimize_trajectory.py is the only script that explicitly references kcal/mol.

EDIT: on second thought, it must be something in how tolerance is being defined as a variable, since I don't get an error in python when I import unit from simtk.openmm. Is it something similar to this issue? https://github.com/openmm/openmm/issues/4069

[jeeberhardt]

(3 months later...)

The issue is now fixed!

[blakemertz]

@jeeberhardt thank you -- how was the issue fixed? I dont see any updates to the source code reflected in github.

[jeeberhardt]

Here the modifications: https://github.com/forlilab/waterkit/commit/fa61e50d15ee6496fad6a9e6802bbb7b5b03ba50#diff-59cd211306dec5adc28683b6041057cacfa65a2ae0bf7bc2c5cd8e85e6de270dR70-R75

[blakemertz]

Thanks for the clarification and fixing this issue. Much appreciated.