-
Notifications
You must be signed in to change notification settings - Fork 114
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
PDBFixer.addSolvent
often incorrectly neutralises systems with charged ligands as though they were neutral
#299
Comments
That's an interesting corner case I hadn't thought about. Lines 1358 to 1364 in 92044f3
Notice that Modeller lets you pass in a force field to use for identifying charges, but PDBFixer generates its own force field, which of course only knows about standard residues. Here are a few possibilities.
A complication with the third one is that it often doesn't know what the formal charges are. Registering a template is only required if you're going to add new atoms to a residue. If all you're doing is preserving a molecule from the input file, it will happily copy it over even if it has no idea what that molecule is. It doesn't need to know anything about it... except in this case, it kind of does because it affects how many ions get added. |
Thanks @peastman, that clarifies the issues really well for me! I might take a crack at option 3 in a fork - would you accept a PR if I can come up with a good solution? |
If you can come up with a good solution, sure, but let's discuss it first. What do you have in mind for figuring out the charges? If we can't figure out anything robust, I'd vote for the second one. It's easy to implement, and if you plan to simulate the system with OpenMM, you presumably have a force field that can parameterize it. |
I'd be happy to try using @Yoshanuikabundi's fork as a way to evaluate implementing Peter's option 3 via either taking formal charges from the CCD, or using the MDAnalysis chemical information guesser at scale. These would both depend on the user knowing their protonation state beforehand, so users could either use the already-implemented protonation functionality, or they could come with fully-protonated inputs. Unfortunately, Peter's options 1 and 2 would hit a circular dependency in SMIRNOFF force field land. We can't make a FF for a molecule without already knowing its formal charges. |
I think that's orthogonal to this issue? This is only about how to determine the total charge for purposes of adding counter ions. Parameterizing the system is much more complicated. And if you have a solution for parameterizing, then allowing a force field to be passed to |
Option 2 would be fine if we had an OpenMM force field that could charge, say, anything in the PDB... but unfortunately we don't. And as you point out, if we did have an appropriate force field we could just use Modeller to add solvent. The I'd propose adding a method to the Then I think @j-wags point is that OpenFF does have a solution for parametrizing, but (a) it isn't in the form of an OpenMM force field, and (b) it requires formal charges. I'm not familiar enough with the new code to know off the top of my head how much of a refactor the above is, but I'm happy to experiment in a fork! |
Feel free to give it a try and see if you can work out something robust. Doing it based on the templates will have two main issues you have to overcome. The first, as noted above, is that we often don't have a template. Ligands may not be in the CCD, or the name may be nonstandard so we don't know what it is. The other is that formal charges depend on protonation. For example, here's the CCD template for LYS. Formal charge is the 5th column.
It lists a +1 charge on NZ, and no charges anywhere else. But at high pH NZ can be neutral, and terminal residues can be charged. So you can't assume the charges listed in the template are correct for the current model. |
The
addSolvent
method uses a force field to assign charges, which it then uses to add ions that neutralise the net charge of the box. Unfortunately, the force field seems to assign formal charges of zero to many charged groups in small molecules; for example the citrate ion, which has a -3 formal charge:The
citrate.pdb
file generated above is available here: citrate.pdb.zipWith the current
main
branch (commit 92044f3), you can also see this behavior with a PDB file like 1PO0, which has 2 citrate ions (each with -3 formal charge).It would be great if the correct formal charges could be stored when the template is downloaded! Then they could be used when neutralising the system. This might also give an opportunity for users to manually specify the formal charges of
UNK
residues.Thanks!
(also the ACE and NME residues in the second example above don't seem to be getting hydrogens - let me know if you'd like me to open a second issue for that!)
The text was updated successfully, but these errors were encountered: