Methods to assign charge during or after atomtyping about foyer HOT 8 OPEN

Excellent. That sounds great. LMK if it would be helpful to have a call with us to discuss strategies for this. Otherwise, please just tag me as you get going as I'll be an interested spectator and potentially interested in using as well.

from foyer.

So,

antechamber is installable in conda:

conda install -c omnia ambermini

And isn't too hard to run. Here, I run it on a propylene carbonate molecule and store the charges in an output mol2 file.

$ time antechamber -i prcarb.pdb -fi pdb -o prcarb.mol2 -fo mol2 -c bcc -nc 0
Total number of electrons: 54; net charge: 0

Running: //anaconda/bin/../bin/sqm -O -i sqm.in -o sqm.out


real	0m0.297s
user	0m0.264s
sys	0m0.025s

And we can obviously grab these charges in parmed:

>>> import parmed as pmd
>>> pc = pmd.load_file('prcarb.mol2', structure=True)
>>> [atom.charge for atom in pc.atoms]
[0.1121, 0.0797, 0.0797, -0.4069, 0.8487, -0.5525, -0.4059, 0.0964, 0.0857, -0.1111, 0.0587, 0.0587, 0.0587]

So it's easy enough to generate these partial charges and store them. But given AM1's inability to scale past small systems, this would need to be run on single prototype molecules, which, as far as I know, isn't really how atomtyping works right now.

One idea I have (say I have an existing gaff.xml file with correct information and partial charges of 0 in <NonBondedForce>):

1. Atomtype the system as normal so the rest of the topology has bonded and LJ terms
2. Pass a list of residue names to a function that slices out a single residue of each type
3. Run an antechamber wrapper on this slice
4. Grab the generated partial charges and map them back onto each copy of this residue

This works in my head and should work for a box of small molecules (or a more complex system in which only its small molecules are passed to this function). This certainly wouldn't work on big slabs and might not work for polymers unless each monomer is accounted for as a single residue.

Any thoughts? Worth pursuing this approach?

from foyer.

The basic procedure that one needs to apply here is conceptually simple but hard to implement in a general way. It's most established for proteins and nucleic acids (and is where the charges in those force fields typically come from) so you may want to see what's been done there. Basically, you "slice" (to use your term) a piece of the system out, cap it in a way so that the caps don't affect its electronic structure much, compute charges, remove the caps, distribute residual charge from the caps back across the slice so that it still has the proper formal charge, then link it back into the larger chain.

It's conceptually simple. The main problems, as I understand it, relate to finding a general way to pick where to slice so that you end up with electronically relatively decoupled regions so that your charges aren't (a) artificial and (b) totally dependent on where you picked to slice.

If you have an OpenEye license, I think OpenEye is working on a general approach for this it might be possible to try out.

I'm also happy to discuss briefly by phone if you guys want.

from foyer.

Thanks for your comment, looks like I wasn't too far off the mark. Could you share some references that demonstrate this? I looked around a bit but I'm not familiar with the biophysical literature and I only found its use in small molecules.

I'm not working with any polymers/biopolymers at the moment, so I was thinking about doing a simpler implementation that was limited to small molecules. Obviously we would want to support arbitrarily large molecules in the futures but this might be useful as a first step (and it would be more useful for my research in the immediate future, so I can start using GAFF and forcefields that rely on setting partial charges).

Thoughts guys? Want to wait until @ctk3b is back to chat or is this worth implementing now?

from foyer.

from foyer.

Thanks for the note on ATB - I hadn't realized they discussed these issues on there so I'll have to look.

Let me know if I can be helpful at all.

Regarding having files versus having charging procedures -- I have always felt knowing the origin of the files is even more important than having the files themselves. The worst case scenario is, as you describe it, papers which just say "We got charges from ATB" without saying which/how, but it's still less than ideal to say, "We got charges from ATB and here are our files," as one still doesn't know where they came from. It's sort of like buying a "diamond ring" on a street corner. Who knows what it really is? I'd much rather have, "Here are my files, and here's the algorithm I used to generate them."

If you do implement something along these lines, can you try and do it in a way which is not Foyer-specific? We'll ultimately need something for blocking too, and it would be nice to be able to reuse rather than rewrite.

from foyer.

from foyer.

Alright, to sum up some of these discussions, I'll be taking a stab at building another package outside of foyer that calls a QM method to assign charge. I'll need to figure out a handful of details along the way (blocking being the biggest hurdle) and I'll tag people in issues elsewhere if I need help.

I agree that documenting the underly algorithms, not necessarily just spitting out files, is the best approach for reproducibility and openness.

If you do implement something along these lines, can you try and do it in a way which is not Foyer-specific?

My vision is for it to take and return parmed structures, which fits right into the foyer workflow but should also be generally useful to people who have other workflows.

from foyer.