Generated system is sometime unstable, potential issue with the estimation of protein and cosolvent radii?

[Cucco23]

Hi,

I tried to run the simulation with the aa.json cosolvents.
My impression is that in certain cases, if there are too many molecules, there is a problem in defining their position or allocating the space.
The system is generated but the simulation does not run and crashes already during minimization.
This is the error I got.

Traceback (most recent call last):
File "/home/ema/miniforge3/envs/cosolvkit/bin/create_cosolvent_system", line 33, in
sys.exit(load_entry_point('cosolvkit', 'console_scripts', 'create_cosolvent_system')())
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/media/ema/Memory_1/cosolvkit-master/cosolvkit/cli/create_cosolvent_system.py", line 168, in main
run_simulation(
File "/media/ema/Memory_1/cosolvkit-master/cosolvkit/simulation.py", line 131, in run_simulation
simulation.step(int(warming_steps / nT))
File "/home/ema/miniforge3/envs/cosolvkit/lib/python3.11/site-packages/openmm/app/simulation.py", line 147, in step
self._simulate(endStep=self.currentStep+steps)
File "/home/ema/miniforge3/envs/cosolvkit/lib/python3.11/site-packages/openmm/app/simulation.py", line 212, in _simulate
self.integrator.step(10) # Only take 10 steps at a time, to give Python more chances to respond to a control-c.
^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/ema/miniforge3/envs/cosolvkit/lib/python3.11/site-packages/openmm/openmm.py", line 15287, in step
return _openmm.LangevinMiddleIntegrator_step(self, steps)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
openmm.OpenMMException: Particle coordinate is NaN. For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#nan

Now the interesting part is that I followed this idea that may be is a problem of organization of the system
In cosolvent_system.py:
self.protein_radius = 3.5 * openmmunit.angstrom
self.cosolvents_radius = 2.5 * openmmunit.angstrom
I changed the parameter 3.5 to 4 and 2.5 to 3.
These parameters seems to control the radius of the protein and of the cosolvents. How these two numbers where chosen since proteins and cosolvents change?
How the simulation will be affected? With this modification the simulation is now running...

Thanks for the help

Ema

[nbruciaferri]

Hi, the system is populated with cosolvent molecules randomly, both in terms of position and orientations. We chose 3.5 and 2.5 because those are reasonable cutoff for the most general and usual cases, where 3.5 is the distance of the cosolvent from every receptor's atom and 2.5 is the distance between the atoms of different cosovlent molecules.

That being said, yes you are right that some systems are more challenging than others and increasing the padding around the receptor and between the cosolvent molecules should result in a more "relaxed" system that should not suffer the MD explosion as much.

I believe that there's no perfect recipe that can be applied to every system/cosolvent, hence the high flexibility of the API and the config file to allow the exploration of different parameters that better fit each special case.
It could be worth adding the receptor and cosolvents padding in the config file so that the user can adjust it.

@diogomart, @mallanos and @mattholc what are your thoughts on this?

[mattholc]

From the traceback it looks like the NaN's occurred during warming. It may be that there are nearby molecules that have bad initial coords and the combination is leading to instabilities. That could be a reason why increased padding would help, the alternative could be that it's actually just luck, and a different random seed at the lower padding would also run. Its hard to say without inspecting the systems just how much crowding is playing a role, and how much its a pair of bad particles.

Out of curiosity, can you provide more system details? Ideally it would be good to know all the cosolvent identities, target concentrations, and actual number placed. If those details are sensitive, total molality (mass of all placed cosolvent molecules/mass of all placed waters) would be a good proxy. I just want to understand how crowded the system was.

Regarding exposing these parameters in the config. I would go with modifiable from config, but with defaults so users don't have to specify them. But they do seem to be a useful debugging tool, so it makes sense to expose them.