Charge constraints

The charge constraints used to fit charges to the ESPs can have a substantial effect on the resulting charges. You can both add manual charge constraints yourself (described in the section below, “Custom charge constraints”) and set general options.

The key options are:

  • symmetric_methyls (default: True)

  • symmetric_methylenes (default: True)

  • symmetric_atoms_are_equivalent (default: False)

  • split_conformers (default: False)

  • constrain_methyl_hydrogens_between_conformers (default: False)

The first three options deal with symmetry. They constrain methyl hydrogens, methylene hydrogens, and symmetric atoms to have the same charge, respectively. Symmetric atoms are determined from the graph representation only, rather than from 3D coordinates. It is generally a good idea to have them all as True, although in some cases where the 3D coordinates differentiate symmetric groups, you may want to keep symmetric_atoms_are_equivalent=False.

split_conformers affects how the constraint matrices are constructed. When split_conformers=False (the default), conformers are treated separately, and equivalence constraints ensure that atoms are given the same charge across multiple molecules. When split_conformers=True, all conformers are merged into one matrix so that they are essentially averaged, without equivalence constraints needed.

When conformers are split, constrain_methyl_hydrogens_between_conformers affects whether equivalence constraints are set between conformers for methyl and methylene hydrogen atoms, for the first stage of a two-stage fit. All atoms are always given equivalence constraints between conformers during the last stage of the fit. Therefore, this option is only relevant if split_conformers=True. Using split_conformers=True, constrain_methyl_hydrogens_between_conformers=False is the approach first proposed by Cornell et al., 1993 in one of the original RESP papers, where the charges methyl/ene hydrogens are free to vary in the first stage.

However, using split_conformers=False generates charges more in line with other packages, and for this reason that is the recommended default. It also saves considerably on memory because the matrices are smaller.

Custom charge constraints

PsiRESP offers two forms of constraints: a ChargeSumConstraint and a ChargeEquivalenceConstraint.

A ChargeSumConstraint constrains one or a group of atoms to a single charge. For example, this is very helpful if you are parametrizing a residue in a large molecule, and you need the cap atoms to sum to 0.

A ChargeEquivalenceConstraint constrains every atom specified to the same charge. For example, you could manually constrain all hydrogens around a single carbon to the same charge. For this reason, a ChargeEquivalenceConstraint must contain at least two atoms.

All of these are contained in ChargeConstraintOptions. ChargeConstraintOptions are passed to a Job. As ChargeConstraintOptions may contain charge constraints for molecules that are not considered in the Job, the job will create MoleculeChargeConstraints to work on. Users should typically not interact with MoleculeChargeConstraints themselves, but let the Job handle it. MoleculeChargeConstraints contains helper methods such as adding symmetry constraints for sp3 hydrogens.

Creating a constraint requires specifying the atoms that the constraint applies to. This can be done in multiple ways. For example, intra-molecular constraints can be created from the molecule and atom indices:

In [1]: import psiresp

In [2]: nme2ala2 = psiresp.Molecule.from_smiles("CC(=O)NC(C)(C)C(NC)=O")

In [3]: constraints = psiresp.ChargeConstraintOptions()

In [4]: constraints.add_charge_sum_constraint_for_molecule(nme2ala2,
   ...:                                                    charge=0,
   ...:                                                    indices=[0, 1])
   ...: 

In [5]: print(constraints.charge_sum_constraints)
[ChargeSumConstraint(atoms={Atom(index=0, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=1, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers)}, charge=0.0)]

These indices can be obtained through SMARTS matching:

In [6]: nme_smiles = "CC(=O)NC(C)(C)C([N:1]([H:2])[C:3]([H:4])([H:5])([H:6]))=O"

In [7]: nme_indices = nme2ala2.get_smarts_matches(nme_smiles)

In [8]: print(nme_indices)
[(8, 21, 9, 22, 23, 24)]

In [9]: constraints.add_charge_equivalence_constraint_for_molecule(nme2ala2,
   ...:                                                            indices=nme_indices[0])
   ...: 

Alternatively, you can pass a list of atoms. This is especially useful for inter-molecular constraints, that involve multiple molecules:

In [10]: methylammonium = psiresp.Molecule.from_smiles("C[NH3+]")

In [11]: methyl_atoms = methylammonium.get_atoms_from_smarts("C([H])([H])([H])")

In [12]: ace_atoms = nme2ala2.get_atoms_from_smarts("C([H])([H])([H])C(=O)N([H])")

In [13]: constraint_atoms = methyl_atoms[0] + ace_atoms[0]

In [14]: constraints.add_charge_sum_constraint(charge=0, atoms=constraint_atoms)

In [15]: constraints.charge_sum_constraints[-1]
Out[15]: ChargeSumConstraint(atoms={Atom(index=14, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=11, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=1, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=0, molecule=Molecule(name=CH6N, charge=1) with 0 conformers), Atom(index=3, molecule=Molecule(name=CH6N, charge=1) with 0 conformers), Atom(index=0, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=13, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=3, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=2, molecule=Molecule(name=CH6N, charge=1) with 0 conformers), Atom(index=12, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=2, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=4, molecule=Molecule(name=CH6N, charge=1) with 0 conformers)}, charge=0.0)

You can also indirectly add constraints with the symmetric_methylenes and symmetric_methyls terms. These add a ChargeEquivalenceConstraint for the appropriate hydrogens.

Note

For now, detecting sp3 carbons requires accurate chemical perception. For reliable symmetry detection, it is highly advisable to create Molecules from SMILES, RDKit molecules, or QCElemental molecules with the connectivity specified.

While the actual constraints are not created in ChargeConstraintOptions, they are specified in MoleculeChargeConstraints. MoleculeChargeConstraints are created by a Job; users should not typically create their own or interact with it. They contain methods for detecting and merging redundant constraints. For example, we create constraint options where a constraint for nme2ala2 is added twice, and a constraint is added that includes atoms from both nme2ala2 and methylammonium:

In [16]: constraints = psiresp.ChargeConstraintOptions(symmetric_methyls=True,
   ....:                                               symmetric_methylenes=True)
   ....: 

# add this constraint twice
In [17]: constraints.add_charge_sum_constraint_for_molecule(nme2ala2,
   ....:                                                    indices=nme_indices[0])
   ....: 

In [18]: constraints.add_charge_sum_constraint_for_molecule(nme2ala2,
   ....:                                                    indices=nme_indices[0])
   ....: 

# add constraint with both nme2ala2 and methylammonium
In [19]: constraints.add_charge_sum_constraint(charge=0, atoms=constraint_atoms)

In [20]: print(len(constraints.charge_sum_constraints))
3

In [21]: print(len(constraints.charge_equivalence_constraints))
0

When we create MoleculeChargeConstraints with only the nme2ala2 molecule, the redundant constraint is removed:

In [22]: mol_constraints = psiresp.charge.MoleculeChargeConstraints.from_charge_constraints(
   ....:                     constraints,
   ....:                     molecules=[nme2ala2],
   ....:                     )
   ....: 
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [22], in <cell line: 1>()
----> 1 mol_constraints = psiresp.charge.MoleculeChargeConstraints.from_charge_constraints(
      2                     constraints,
      3                     molecules=[nme2ala2],
      4                     )

File ~/checkouts/readthedocs.org/user_builds/psiresp/conda/latest/lib/python3.8/site-packages/psiresp-0.4.2+6.g9518fa4.dirty-py3.8.egg/psiresp/charge.py:363, in MoleculeChargeConstraints.from_charge_constraints(cls, charge_constraints, molecules)
    351 eqvs = [constr.copy(deep=True)
    352         for constr in charge_constraints.charge_equivalence_constraints
    353         if constr.molecule_set & molecule_set]
    355 dct = {
    356     k: getattr(charge_constraints, k)
    357     for k in [
   (...)
    360     ]
    361 }
--> 363 constraints = cls(
    364     charge_sum_constraints=sums,
    365     charge_equivalence_constraints=eqvs,
    366     molecules=molecules,
    367     **dct
    368 )
    370 accepted = []
    371 if charge_constraints.symmetric_methyls:

File ~/checkouts/readthedocs.org/user_builds/psiresp/conda/latest/lib/python3.8/site-packages/psiresp-0.4.2+6.g9518fa4.dirty-py3.8.egg/psiresp/base.py:42, in Model.__init__(__pydantic_self__, **data)
     40 __pydantic_self__.__pre_init__(**data)  # lgtm[py/init-calls-subclass]
     41 super().__init__(**data)
---> 42 __pydantic_self__.__post_init__(**data)  # lgtm[py/init-calls-subclass]

File ~/checkouts/readthedocs.org/user_builds/psiresp/conda/latest/lib/python3.8/site-packages/psiresp-0.4.2+6.g9518fa4.dirty-py3.8.egg/psiresp/charge.py:308, in MoleculeChargeConstraints.__post_init__(self, **kwargs)
    306 self.clean_charge_sum_constraints()
    307 self.clean_charge_equivalence_constraints()
--> 308 self._generate_molecule_increments()

File ~/checkouts/readthedocs.org/user_builds/psiresp/conda/latest/lib/python3.8/site-packages/psiresp-0.4.2+6.g9518fa4.dirty-py3.8.egg/psiresp/charge.py:324, in MoleculeChargeConstraints._generate_molecule_increments(self)
    322     self._molecule_increments[hash(mol)] = inc_list
    323 self._n_total_atoms = increment
--> 324 self._generate_edges()

File ~/checkouts/readthedocs.org/user_builds/psiresp/conda/latest/lib/python3.8/site-packages/psiresp-0.4.2+6.g9518fa4.dirty-py3.8.egg/psiresp/charge.py:328, in MoleculeChargeConstraints._generate_edges(self)
    326 def _generate_edges(self):
    327     increments = list(self._molecule_increments.values())
--> 328     self._n_molecule_atoms = np.array([x[0] for x in increments] + [self._n_total_atoms])
    329     self._edges = []
    330     for mol in self.molecules:

File ~/checkouts/readthedocs.org/user_builds/psiresp/conda/latest/lib/python3.8/site-packages/psiresp-0.4.2+6.g9518fa4.dirty-py3.8.egg/psiresp/charge.py:328, in <listcomp>(.0)
    326 def _generate_edges(self):
    327     increments = list(self._molecule_increments.values())
--> 328     self._n_molecule_atoms = np.array([x[0] for x in increments] + [self._n_total_atoms])
    329     self._edges = []
    330     for mol in self.molecules:

IndexError: list index out of range

In [23]: print(len(mol_constraints.charge_sum_constraints))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [23], in <cell line: 1>()
----> 1 print(len(mol_constraints.charge_sum_constraints))

NameError: name 'mol_constraints' is not defined

And the sp3 equivalences are added:

In [24]: print(len(mol_constraints.charge_equivalence_constraints))
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [24], in <cell line: 1>()
----> 1 print(len(mol_constraints.charge_equivalence_constraints))

NameError: name 'mol_constraints' is not defined