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=1, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=0, 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=4, molecule=Molecule(name=CH6N, charge=1) with 0 conformers), Atom(index=12, 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=2, molecule=Molecule(name=C7H14N2O2, charge=0) with 0 conformers), Atom(index=2, molecule=Molecule(name=CH6N, charge=1) 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=14, 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=0, molecule=Molecule(name=C7H14N2O2, charge=0) 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],
....: )
....:
In [23]: print(len(mol_constraints.charge_sum_constraints))
2
And the sp3 equivalences are added:
In [24]: print(len(mol_constraints.charge_equivalence_constraints))
4