Partial charges from hits

[matteoferla]

This idea is interesting, but challenging and potentially of little effect.

Hypothesis: Some compounds (with large/delocalised) orbitals are polarisable, ie. shuffle partial charges around: copying these partial charges from the hits may improve energy calculations.

Rosetta does not have a Drude particle and the partial charges are defined in the topology definition of the residue (residuetype, defined in the params file), so are permanent. Rosetta virtual atoms do not have charge and were they charged they would not affect LJ by definition.
Whereas nothing can be done for soft atoms like bromine, for heterocycles RDKit can give the resonance forms, each with a different set of partial charges. I figured out how to switch in PyRosetta the partial charges of a residue type, so one could score each resonance form of a hit and find the best one. Then find the resonance form in the followup that best matches these.

There are two Qs: technical feasibility and actual utility.

Technical

Aim: get the gasteiger chargers of the atoms in the various resonance forms and find the combination which scores the highest in PyRosetta.

RDKit can give resonant forms:

from rdkit import Chem
from rdkit.Chem import AllChem

original = Chem.MolFromSmiles('c12c(cc[nH]2)cccc1')

flags = Chem.rdchem.ResonanceFlags.ALLOW_CHARGE_SEPARATION + \
         Chem.rdchem.ResonanceFlags.ALLOW_INCOMPLETE_OCTETS +\
        Chem.rdchem.ResonanceFlags.UNCONSTRAINED_ANIONS + \
        Chem.rdchem.ResonanceFlags.UNCONSTRAINED_CATIONS

for i, mol in enumerate(AllChem.ResonanceMolSupplier(original, flags=flags)):
    AllChem.ComputeGasteigerCharges(mol)
    from rdkit.Chem.Draw import SimilarityMaps
    AllChem.ComputeGasteigerCharges(mol)
    contribs = [mol.GetAtomWithIdx(i).GetDoubleProp('_GasteigerCharge') for i in range(mol.GetNumAtoms())]
    SimilarityMaps.GetSimilarityMapFromWeights(mol, contribs, colorMap='jet', contourLines=10)
    #.savefig(f"resonance_form_{i}.png", bbox_inches='tight')

As far as I can tell, in PyRosetta one has to generate a new residuetype and place it anew (while copying over the coordinates and crosslinks).
Even though one can mutate residuetypes in PyRosetta, the mutable type is converted to a standard one in the residue type set —the ResidueFactory does not accept mutable residues:

from fragmenstein import Igor
import pyrosetta
Igor.init_pyrosetta()
from rdkit_to_params import Params
# make a pose with that an indole
params = Params.from_smiles('c12c(cc[nH]2)cccc1')
pose = pyrosetta.Pose()
rts = pose.conformation().modifiable_residue_type_set_for_conf(pyrosetta.rosetta.core.chemical.FULL_ATOM_t)
buffer = pyrosetta.rosetta.std.stringbuf(params.dumps())
stream = pyrosetta.rosetta.std.istream(buffer)
new = pyrosetta.rosetta.core.chemical.read_topology_file(stream, params.NAME, rts)
# rts.add_base_residue_type(new)
print(f'old charge: {params.ATOM[1].partial}')
# make a mutable type and change it
mrt = pyrosetta.rosetta.core.chemical.MutableResidueType(new)
sac = pyrosetta.rosetta.core.chemical.SetAtomicCharge(params.ATOM[0].name, 0.8)
sac.apply(mrt)
# add it to the residue type set —here the issue lies
rts.add_base_residue_type(mrt)
pose.conformation().reset_residue_type_set_for_conf(rts)
# add new residue. `rts.name_map( params.NAME )` returns a `pyrosetta.rosetta.core.chemical.ResidueType`
lig: pyrosetta.rosetta.core.conformation.Residue = pyrosetta.rosetta.core.conformation.ResidueFactory.create_residue( rts.name_map( params.NAME ) )
print(lig.atom_name(1), lig.atomic_charge(1))
# is it a mutable type set?
print(type(mrt), type(lig.type())) # <class 'pyrosetta.rosetta.core.chemical.MutableResidueType'> <class 'pyrosetta.rosetta.core.chemical.ResidueType'>

Therefore, the only option is switching out the residue for every partial charge set.

Conclusion

Given the hassle, it would need to be proven first that this matters in the first place: given that the atoms are constrained anyway, the benefits may be minor. A test case could be PLP which is has a heterocycle nitrogen interaction and several structures (covalently bound inhibitors, e.g. 2-chloroalanine) in the PDB.