RESP Charges#1552
Conversation
proteneer
left a comment
proteneer
left a comment
There was a problem hiding this comment.
Two major things to discuss offline
- openff toolkit/recharge dependency - I think this is mostly from _elf generation.
- Auto3D/AIMNET dependency (which adds torch, torchvision, torchaudio).
| from Auto3D.ASE.geometry import opt_geometry | ||
| from jax import jit, vmap | ||
| from numpy.typing import NDArray | ||
| from openff.recharge.aromaticity import AromaticityModel |
There was a problem hiding this comment.
is this overwritten by from timemachine.ff.handlers.bcc_aromaticity import AromaticityModel down below?
| from numpy.typing import NDArray | ||
| from openff.recharge.aromaticity import AromaticityModel | ||
| from openff.recharge.utilities.molecule import extract_conformers | ||
| from openff.toolkit import unit, RDKitToolkitWrapper |
There was a problem hiding this comment.
I think unit here is just openmm.unit
| molecule: Molecule = Molecule.from_rdkit(rdmol) | ||
| # Apply ELF conformer selection to get diverse, representative conformers | ||
| # This helps ensure charges are computed from a representative ensemble | ||
| molecule.apply_elf_conformer_selection( |
There was a problem hiding this comment.
ugh - I think this is the biggest dependency on the openff
| cache_key = nonbonded.RESP | ||
| am1h = nonbonded.RESPHandler(smirks, params, props) | ||
| mol = Chem.AddHs(Chem.MolFromSmiles("C1CNCOC1F")) | ||
| AllChem.EmbedMolecule(mol) |
There was a problem hiding this comment.
If RESP to be independent of the initial conformer, we should explicitly be calling mol.RemoveAllConformers() beforehand as opposed to EmbedMolecule()
| # assert vjp_fn(charges_adjoints) == None | ||
|
|
||
|
|
||
| def test_resp_parameterization(): |
There was a problem hiding this comment.
We should add a test for determinism (either in the presence or absence of a seed) since the underlying algorithm uses a stochastic conformer generator. i.e. we want to make sure that calling RESP multiple items returns identical charges on fresh non-cached molecules.
| from openff.toolkit import unit, RDKitToolkitWrapper | ||
| from openff.toolkit.topology import Molecule | ||
| from openff.toolkit.utils import AntechamberNotFoundError | ||
| from openff.units import Quantity |
There was a problem hiding this comment.
again, I think this is just what's in OpenMM
| mol.build() | ||
|
|
||
| # Perform Restricted Hartree-Fock calculation on GPU | ||
| mf = scf.RHF(mol).to_gpu() |
There was a problem hiding this comment.
Is there a meaningful speed difference of cpu vs gpu? We've never tested downstream timemachine code robustly in multi-CUDA context settings, and generally assumes that tm is the only library within the process that creates/initializes the underlying cuda context.
| partial_charges = np.mean(am1_partial_charges, axis=0) | ||
|
|
||
| # Ensure total charge equals the formal charge of the molecule | ||
| # Small numerical errors can accumulate, so we redistribute any discrepancy |
There was a problem hiding this comment.
How small are we talking about here? Is the source of error numerical or algorithmic?
Numerical: using float32 vs float64 etc.
Algorithm: use of restraints, constraints, poor convergence, etc.
There was a problem hiding this comment.
And would it be better to do this on a per-individual-conformer basis first and then average?
is avg(spread(charges)) == spread(avg(charges)) i.e. do the operators spread and avg commute?
| params.pruneRmsThresh = 1.0 | ||
| params.clearConfs = True | ||
|
|
||
| AllChem.EmbedMultipleConfs(mol, 800, params) |
There was a problem hiding this comment.
It is just using the same number of conformers the existing oechem code uses. I don't know the origin of that choice.
2 participants
Uses pyscf to compute RESP charges as an alternative to OpenEye's AM1BCC