openforcefield / openff-nagl

OpenFF NAGL
https://docs.openforcefield.org/projects/nagl/en/latest/?badge=latest
MIT License
13 stars 1 forks source link

Weird effect when charging cations alone or together: #136

Closed mrshirts closed 1 month ago

mrshirts commented 1 month ago

Expected behavior

A ionic compound can be written in smiles with a . connecting them (like [Na+].[Cl-]). I would EXPECT, because there are no actual bonds, the charges generated when passing an organic cation and anion together would be same as when they are processed apart. But that's not the case, interestingly!

Actual behavior

See attached Python notebook for an example. When charged as an ionic compound (a single smiles with a . joining the cation and anion) , the charges are +/- 0.81, rather than +1/-1, and all of the charges are slightly off. When processed separately, the formal charges are zero.

I hypothesize this is because the total formal charge OVERALL is 0, rather than being +1 and -1, so the constraints are different, and it does something weird.

The REASON for trying to do it this way is that this is to make an ionic liquid system, and for packing molecules, it would be better to pack them TOGETHER, so the anions and cations don't have to find each other.

Note that there's also another issue, in that the two molecules are output to GROMACS/LAMMPS as a single molecule. This is not necessarily horrible (file should still run - there are no bonds), though obviously not great for analysis since the cation and anion are not separate. Though I haven't actually tried to run the file, it could be that GROMACS gags eventually as some intramolecular interactions become really long. I can try that if people are interested.

What is the priority? Probably low - knowing this exists, it can be worked around - just have to work harder in equilibration. But I thought it was interesting that something you COULD do without any warnings gets easily messed up.

Code to reproduce the behavior

Attached notebook (remove the .txt): IL_tests_nagl_charging.ipynb.txt

Current environment

Everything installed via conda recently, conda list attached: conda.list.txt

lilyminium commented 1 month ago

Ah, yup, you're exactly right -- this arises because if the OpenFF toolkit considers something a single Molecule, NAGL accepts it as such. The individual partial charges are solved for with the total charge of the "molecule" as the constraint, so they will be different if the anion and cation are combined to have a total formal charge of 0, vs being split.

Happy to put in a fix that splits up unbound fragments. I'll just check first with the OpenFF toolkit though to see if the Molecule spec should be accepting these as single Molecules, or if the warning should occur upstream instead.

Thanks for raising this issue!

j-wags commented 1 month ago

cc https://github.com/openforcefield/openff-toolkit/issues/1587

j-wags commented 1 month ago

I hate to push the patch down to NAGL since this is inherently a Toolkit problem, but there are a few important pieces of functionality that rely on this behavior/loophole, and NAGL is so close to release that I don't want to delay it any more than absolutely necessary. So I think implementing a solution here in NAGL is unfortunately the right way to go.

mrshirts commented 1 month ago

Yeah, I don't want to dictate what the right way to solve this is (disallow it, or break up the molecule), but just wanted to point out that the current functionality is probably not something that should happen longer term. But it's also probably a lower priority!

mrshirts commented 1 month ago

Note that there's also another issue, in that the two molecules are output to GROMACS/LAMMPS as a single molecule. This is not necessarily horrible (file should still run - there are no bonds), though obviously not great for analysis since the cation and anion are not separate. Though I haven't actually tried to run the file, it could be that GROMACS gags eventually as some intramolecular interactions become really long. I can try that if people are interested.

This is the toolkit problem that persists even if the NAGL problem is fixed. I think there are workarounds for now, so can be prioritized lower. Something to address eventually.

lilyminium commented 1 month ago

@mrshirts I think I have a fix running in #137, if this is currently affecting your workflows I'll put out a release tomorrow that you could run with -- if not I'll target next week to see if I can bundle some other nagl fixes along with.

mrshirts commented 1 month ago

No, does not affect current workflows; we're just testing functionality for writing future workflows. This can wait.