[BUG] Coordinates from GROMACS vacuum simulations can cause overflow as input to AMBER

[lohedges]

For GROMACS we implement vacuum simulations by using pseudo-PBC conditions, i.e. running calculations in a near infinite box, as described here. (This approach works with all versions of GROMACS that we support.) However, a result of this is that absolute molecular coordinates can end up being really large, so can't be used as input to other engines, e.g. AMBER, since they overflow the formatting restrictions of the input file. For example:

In [1]: import BioSimSpace as BSS

INFO:numexpr.utils:Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
WARNING:root:Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.

In [2]: BSS.setVerbose(True)

In [3]: m = BSS.Parameters.openff_unconstrained_2_0_0("CC").getMolecule()

In [4]: p = BSS.Protocol.Minimisation()

In [5]: process = BSS.Process.Gromacs(m.toSystem(), p)
/home/lester/Code/openbiosim/biosimspace/python/BioSimSpace/Process/_gromacs.py:294: UserWarning: No simulation box found. Assuming gas phase simulation.
  _warnings.warn("No simulation box found. Assuming gas phase simulation.")
/home/lester/Code/openbiosim/biosimspace/python/BioSimSpace/_Config/_config.py:96: UserWarning: No simulation box found. Assuming gas phase simulation.
  _warnings.warn("No simulation box found. Assuming gas phase simulation.")

In [6]: process.start()
Out[6]: BioSimSpace.Process.Gromacs(<BioSimSpace.System: nMolecules=1>, BioSimSpace.Protocol.Minimisation(steps=10000, restraint=None, force_constant=10.0 kcal_per_mol/angstrom**2), exe='/home/lester/.conda/envs/openbiosim/bin.AVX2_256/gmx', name='gromacs', work_dir='/tmp/tmpyuhwqevj', seed=None)

In [7]: process.wait()

In [8]: s = process.getSystem()

In [9]: s[0].coordinates()
Out[9]:
[(9998.2400 A, 0.0200 A, 9999.0200 A),
 (9999.7600 A, -0.0200 A, 9998.9800 A),
 (9997.8600 A, -0.8300 A, 9999.6000 A),
 (9997.8900 A, 0.9500 A, 9999.4900 A),
 (9997.8300 A, -0.0400 A, 9998.0100 A),
 (1.0000e+04 A, -0.9200 A, 9998.4600 A),
 (1.0000e+04 A, 0.8600 A, 9998.4500 A),
 (1.0000e+04 A, -0.0300 A, 9999.9900 A)]

In [10]: BSS.IO.saveMolecules("test", s, "rst7")
╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ /home/lester/Code/openbiosim/biosimspace/python/BioSimSpace/IO/_io.py:730 in saveMolecules       │
│                                                                                                  │
│    727 │   │   │   │   _os.rename(file, new_file)                                                │
│    728 │   │   │   │   file = [new_file]                                                         │
│    729 │   │   │   else:                                                                         │
│ ❱  730 │   │   │   │   file = _SireIO.MoleculeParser.save(                                       │
│    731 │   │   │   │   │   system._sire_object, filebase, _property_map                          │
│    732 │   │   │   │   )                                                                         │
│    733                                                                                           │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
OSError: SireError::io_error: Cannot write the (perhaps some of the ) files as the following errors occurred:
Failed to write the file '/home/lester/Downloads/vac_issue/test.rst7' using the parser for fileformat 'RST7'. Errors reported by the parser
are:
Errors converting the system to a Amber Rst7 format...
Could not write the float at index 15, value '10000.1' into the specified format AmberFormat( 6 x float[width = 12, precision = 7] ).
Could not write the float at index 18, value '10000.1' into the specified format AmberFormat( 6 x float[width = 12, precision = 7] ).
Could not write the float at index 21, value '10000.2' into the specified format AmberFormat( 6 x float[width = 12, precision = 7] ). (call
sire.error.get_last_error_details() for more info)

The above exception was the direct cause of the following exception:

╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <module>:1                                                                                    │
│                                                                                                  │
│ /home/lester/Code/openbiosim/biosimspace/python/BioSimSpace/IO/_io.py:742 in saveMolecules       │
│                                                                                                  │
│    739 │   │   except Exception as e:                                                            │
│    740 │   │   │   msg = "Failed to save system to format: '%s'" % format                        │
│    741 │   │   │   if _isVerbose():                                                              │
│ ❱  742 │   │   │   │   raise IOError(msg) from e                                                 │
│    743 │   │   │   else:                                                                         │
│    744 │   │   │   │   raise IOError(msg) from None                                              │
│    745                                                                                           │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
OSError: Failed to save system to format: 'RST7'

I'm confused why this is happening since the molecule should be placed at the centre of the box, but perhaps this isn't working correctly. For a minimisation it should be easy enough to translate it back to some reasonable place afterwards, but maybe not for other protocols where a user might want a trajectory.

[lohedges]

After a single step the coordinates are completely sensible. Perhaps the vacuum box origin is at (0, 0, 0) so that we need to translate the molecule to the middle, then back again at the end. This is what I need to do for solvation too.

In [1]: import BioSimSpace as BSS

INFO:numexpr.utils:Note: NumExpr detected 20 cores but "NUMEXPR_MAX_THREADS" not set, so enforcing safe limit of 8.
INFO:numexpr.utils:NumExpr defaulting to 8 threads.
WARNING:root:Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.

In [2]: m = BSS.Parameters.openff_unconstrained_2_0_0("CC").getMolecule()

In [3]: p = BSS.Protocol.Minimisation(steps=1)

In [4]: process = BSS.Process.Gromacs(m.toSystem(), p)
/home/lester/Code/openbiosim/biosimspace/python/BioSimSpace/Process/_gromacs.py:294: UserWarning: No simulation box found. Assuming gas phase simulation.
  _warnings.warn("No simulation box found. Assuming gas phase simulation.")
/home/lester/Code/openbiosim/biosimspace/python/BioSimSpace/_Config/_config.py:96: UserWarning: No simulation box found. Assuming gas phase simulation.
  _warnings.warn("No simulation box found. Assuming gas phase simulation.")

In [5]: process.start()
Out[5]: BioSimSpace.Process.Gromacs(<BioSimSpace.System: nMolecules=1>, BioSimSpace.Protocol.Minimisation(steps=1, restraint=None, force_constant=10.0 kcal_per_mol/angstrom**2), exe='/home/lester/.conda/envs/openbiosim/bin.AVX2_256/gmx', name='gromacs', work_dir='/tmp/tmpoiefdyqf', seed=None)

In [6]: process.wait()

In [7]: s = process.getSystem()

In [8]: s[0].coordinates()
Out[8]:
[(-0.7500 A, 0.0400 A, 0.0100 A),
 (0.7500 A, 0.0000e+00 A, 0.0000e+00 A),
 (-1.1300 A, -0.6400 A, 0.8100 A),
 (-1.1800 A, 1.0300 A, 0.2000 A),
 (-1.2000 A, -0.3300 A, -0.9400 A),
 (1.0800 A, -0.7400 A, -0.7700 A),
 (1.2300 A, 0.9600 A, -0.2700 A),
 (1.2000 A, -0.3200 A, 0.9500 A)]

[lohedges]

No, shouldn't be needed since you can see that negative coordinates are present above. I guess it's just wandered off a long way! Looks like it goes bonkers after about 50 steps.

[jmichel80]

annoying...maybe the hydration free energy tests can all use OpenMM as engine for equilibration for the time being?

[lohedges]

Yes, I'd suggest using OpenMM for vacuum minimisation/equilibration for now.