Molecules
PsiRESP is built on MolSSI’s QC stack. The psiresp.molecule.Molecule
,
psiresp.conformer.Conformer
, and psiresp.orientation.Orientation
classes each wrap a qcelemental.models.Molecule
with the qcmol
attribute.
This has several advantages; molecules can be easily hashed, for example, and they interact natively with the engine-agnostic QCEngine and QCFractal packages. A QCElemental molecule can be created from, and corresponds directly to, the commonly-used XYZ format.
There is a hierarchy to these molecules in the PsiRESP API.
A psiresp.molecule.Molecule
contains
psiresp.conformer.Conformer
s, which contains
psiresp.orientation.Orientation
s.
Conformers
RESP methods are known to be very conformation-dependent, so including more conformers increases the likelihood of getting a charge profile that’s more suited to the conformers explored in simulation.
Adding conformers manually
You can define your own conformers with QCElemental molecules:
In [1]: import qcelemental as qcel
In [2]: import psiresp as sip
In [3]: dmso_spec = """\
...: 10
...: ! from R.E.D. examples/2-Dimethylsulfoxide
...: C 3.7787218489 0.7099460365 -8.4358800149
...: H 3.7199547803 1.7353551063 -8.0912740044
...: H 3.3727951491 0.0324988805 -7.6949604622
...: H 3.2200309413 0.6050316563 -9.3564160613
...: S 5.4843687520 0.2699083657 -8.7873057009
...: O 5.4949906162 -1.1820495711 -9.0993845212
...: C 6.1255577314 0.4439602615 -7.1184894575
...: H 7.1686363815 0.1575024332 -7.1398743610
...: H 6.0389920737 1.4725205702 -6.7894907284
...: H 5.5889517256 -0.2186737389 -6.4509246884
...: """
...:
In [4]: qcdmso = qcel.models.Molecule.from_data(dmso_spec, dtype="xyz")
In [5]: dmso = sip.Molecule(qcmol=qcdmso)
# No conformers are generated automatically
In [6]: assert dmso.n_conformers == 0
In [7]: dmso.add_conformer(qcmol=qcdmso)
In [8]: assert dmso.n_conformers == 1
In [9]: print(dmso.conformers)
[Conformer(qcmol=Molecule(name='C2H6OS', formula='C2H6OS', hash='d8fbedf'), orientations=[], is_optimized=False)]
Or directly with coordinates:
In [10]: import numpy as np
In [11]: coordinates = np.array([
....: [0.00000000, 0.00000000, 0.00000000],
....: [-0.37443200, -1.01033397, -0.11267908],
....: [-0.36445856, 0.63526149, -0.79767520],
....: [-0.32471908, 0.40230888, 0.95038637],
....: [1.79620837, 0.00000000, -0.00000000],
....: [2.22305520, 1.42249416, 0.00000000],
....: [2.03592873, -0.61092095, -1.67202714],
....: [3.10077808, -0.62557805, -1.86283992],
....: [1.63738124, -1.61401863, -1.76489819],
....: [1.55810229, 0.05835618, -2.37659975],
....: ])
....:
In [12]: dmso.add_conformer_with_coordinates(coordinates, units="angstrom")
In [13]: assert dmso.n_conformers == 2
In [14]: print(dmso.conformers)
[Conformer(qcmol=Molecule(name='C2H6OS', formula='C2H6OS', hash='d8fbedf'), orientations=[], is_optimized=False), Conformer(qcmol=Molecule(name='C2H6OS', formula='C2H6OS', hash='e02f359'), orientations=[], is_optimized=False)]
Automatically generating conformers
However, automatically generating conformers is probably easiest
and likely to get better results. The conformers generated depend
on the psiresp.conformer.ConformerGenerationOptions
passed to a psiresp.molecule.Molecule
.
The process of generating and selecting conformers is as follows:
Use RDKit to generate
n_conformer_pool
initial conformers at leastrms_tolerance
angstrom apart in RMSDKeep only the conformers within a certain energy window in kcal/mol. This means only those conformers within
energy_window
kcal/mol of the lowest energy conformer are considered for the next step.Select a set with, at most,
n_max_conformers
maximally diverse conformers from the remaining pool. Diversity is calculated by heavy atom RMSD.
It is recommmended to geometry optimize these conformers before
generating Orientations from them. psiresp.job.Job.run()
will
do this automatically, providing
psiresp.molecule.Molecule.optimize_geometry = True.
Orientations
It is also recommended to include multiple orientations
for each conformer in the RESP calculation.
The orientations are controlled by the
psiresp.molecule.Molecule.reorientations
,
psiresp.molecule.Molecule.rotations
, and
psiresp.molecule.Molecule.translations
attributes, as well as
psiresp.molecule.Molecule.keep_original_orientation
.
psiresp.molecule.Molecule.reorientations
and
psiresp.molecule.Molecule.rotations
are lists of atom indices, whereas
psiresp.molecule.Molecule.translations
is a translation vector.
Reorientations
Three atom indices must be specified. The first atom becomes the new origin; the second defines the x-axis from the origin; and the third defines the xy plane.
Rotations
Three atom indices must be specified. The first two atoms define a vector parallel to the x-axis, while the third defines a plane parallel to the xy-plane.
Translations
Three floats must be given, as the translation in the x, y, and z axes.
Automatically generating transformations
As with Conformers, Orientation specifications can be automatically generated with
psiresp.molecule.Molecule.generate_transformations()
.
Note
This method does not generate the Orientations themselves, but rather
fills the reorientations
,
rotations
, and
translations
lists. This means that
you can, and should, generate the transformations before generating
conformers.
If given a desired number of reorientations or rotations, combinations of atoms will be generated to reorient the molecule around. The method first combines heavy atoms, before including hydrogens.
If given a desired number of translations, random translation vectors will be generated between -5 to 5 angstrom on each axis.