Support removal of Hydrogens when using gaff forcefield #60
Conversation
Codecov Report
@@ Coverage Diff @@
## master #60 +/- ##
==========================================
- Coverage 98.95% 98.91% -0.04%
==========================================
Files 6 6
Lines 191 184 -7
==========================================
- Hits 189 182 -7
Misses 2 2
|
planckton/init.py
Outdated
| @@ -182,7 +177,7 @@ def pack(self, box_expand_factor=5): | |||
| typed_system : ParmEd structure | |||
| ParmEd structure of filled box | |||
| """ | |||
| if self.remove_hydrogen_atoms: | |||
| if self.remove_hydrogen_atoms and self.ff == FORCEFIELD["gaff-custom"]: | |||
| self._remove_hydrogen() | |||
There was a problem hiding this comment.
Does this work on a typed system? It may be better to type the system, then check if we want to remove hydrogens, then do it the gaff way or the non gaff way instead of doing it in two different parts of the code
There was a problem hiding this comment.
Well, the existing function to remove hydrogens does so from the mbuild compound before the forcefield is applied by foyer. This only works if someone is using one of the pre-existing *_typed.mol2 compounds along with the gaff-custom forcefield.
The other scenario is if someone is generating a compound from SMILES, or a file that hasn't been typed using the old ambertools method. These compounds can't have hydrogens removed from the mbuild compound before being typed by foyer since they are typed using SMARTS matching.
We could check if removing the hydrogens from the parmed structure works for the *_typed.mol2 compounds, and get rid of or change the _remove_hydrogens to only use this method.
typed_system.strip([a.atomic_number == 1 for a in typed_system.atoms])
Ultimately, this is just a temporary solution that allows someone to do UA sims if they are simulating systems outside of the existing group we have (e.g. the REUs this summer), but I think the goal is to move away from having the pre-typed mol2 files and the custom gaff forcefield #59. Once that happens then there would only be one way of
There was a problem hiding this comment.
It looks like typed_system.strip([a.atomic_number == 1 for a in typed_system.atoms]) would still work on one of the typed mol2 files after playing around real quick in a jupyter notebook. I'll try it out in the PR and see how it works.
There was a problem hiding this comment.
Ok, I tested it out in planckton and it looks like that works. Keep the hydrogens on until after the forcefield is applied by foyer and just remove them from the resulting parmed structure.
There was a problem hiding this comment.
Perfect that should be much more simple!
There was a problem hiding this comment.
This is so much cleaner--thanks @chrisjonesBSU! We might notice a slowdown with large systems (removing H from one compound vs. from whole system) but I'm not sure. Anyway, using the element number is much better than maintaining a list of all the H types, and this gives us much more flexibility! :)
This PR should allow users to run united atom simulations while using the gaff force field. Since we still have the custom-gaff force field to deal with, this basically results in 2 different methods of hydrogen being removed, but this really only adds one extra if statement. This can be cleaned up once #59 is resolved.
This solves #45