Added SolidAngle (#41)

* Added SolidAngle

- Added SolidAngle class
- Added solid angle documentation
- Added solid angle tests
- Updated changelog
- Upped version number
- Added symlinks for shared test data

CHANGELOG.md

@@ -6,6 +6,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.7.0] - 2022-06-22
+
+### Added
+- Added `SolidAngle` for solid angle calculations.
+
 ## [0.6.0] - 2022-03-28
 
 ### Added

README.md

@@ -54,6 +54,7 @@ For further information, see the separate [documentation](https://kjelljorner.gi
 * Exact ligand cone angle
 * Local force constant
 * Pyramidalization
+* Solid angle
 * Solvent accessible surface area
 * Sterimol parameters
 * XTB electronic descriptors

docs/index.rst

@@ -62,6 +62,7 @@ Features
 * Exact ligand cone angle
 * Local force constant
 * Pyramidalization
+* Solid angle
 * Solvent accessible surface area
 * Sterimol parameters
 * XTB electronic descriptors
@@ -109,6 +110,7 @@ both academic and commercial use.
    local_force
    pyramidalization
    sasa
+   solid_angle
    sterimol
    xtb
 

docs/refs.bib

@@ -1,9 +1,35 @@
+@article{bilbrey_solid_2013,
+  title     = {Exact ligand solid angles},
+  author    = {Bilbrey, Jenna A and Kazez, Arianna H and Locklin, Jason and Allen, Wesley D},
+  year      = {2013},
+  month     = dec,
+  journal   = {Journal of Chemical Theory and Computation},
+  volume    = {9},
+  number    = {12},
+  pages     = {5734--5744},
+  publisher = {{American Chemical Society}},
+  doi       = {10.1021/ct400426e}
+}
+
+@article{guzei_2006,
+  title   = {An Improved Method for the Computation of Ligand Steric Effects Based on Solid Angles},
+  author  = {Guzei, Ilia A. and Wendt, Mark},
+  year    = {2006},
+  journal = {Dalton Transactions},
+  number  = {33},
+  pages   = {3991},
+  issn    = {1477-9226, 1477-9234},
+  doi     = {10.1039/b605102b}
+}
+
 @article{durand_2019,
-  title   = {Computational ligand descriptors for catalyst design},
-  author  = {Durand, Derek J and Fey, Natalie},
+  title   = {Computational Ligand Descriptors for Catalyst Design},
+  author  = {Durand, Derek J. and Fey, Natalie},
   year    = {2019},
   journal = {Chemical Reviews},
-  issn    = {0009-2665},
+  volume  = {119},
+  number  = {11},
+  pages   = {6561--6594},
   doi     = {10.1021/acs.chemrev.8b00588}
 }
 

docs/solid_angle.rst

@@ -0,0 +1,73 @@
+###########
+Solid angle
+###########
+
+Solid angle and solid cone angles are implemented with a numerical recipe as
+described by Guzei and Wendt :footcite:`guzei_2006`.
+
+******
+Module
+******
+
+The SolidAngle class is provided to calculate and store the solid angles.
+
+.. code-block:: python
+  :caption: Example
+
+  >>> from morfeus import SolidAngle, read_xyz
+  >>> elements, coordinates = read_xyz("phosphines/PdPMe3.xyz")
+  >>> solid_angle = SolidAngle(elements, coordinates, 1)
+  >>> print(solid_angle.cone_angle)
+  112.53614024551955
+  >>> print(solid_angle.solid_angle)
+  2.794082404625142
+  >>> print(solid_angle.G)
+  22.234601305109027
+  >>> solid_angle.print_report()
+  Solid angle (sr): 2.794
+  Cone angle (°): 112.536
+  G: 22.235
+
+The Bondi vdW radii are used in reference :footcite:`bilbrey_solid_2013`, but
+radii from the CRC Handbook is the default here. It can be changed with
+``radii_type=<str>``. Custom radii can passed with ``radii=<list>``. The units
+for the solid cone angle are in degrees, steradians for the solid angle and the
+G parameter is in percent.
+
+For more detailed information, use ``help(SolidAngle)`` or see the API:
+:py:class:`SolidAngle <morfeus.solid_angle.SolidAngle>`.
+
+*******************
+Command line script
+*******************
+
+The command line script provides access to the basic functionality through the
+terminal.
+
+.. code-block:: console
+  :caption: Example
+
+  $ morfeus solid_angle PdPMe3.xyz - 1 - print_report
+  Solid angle (sr): 2.794
+  Cone angle (°): 112.536
+  G: 22.235
+
+**********
+Background
+**********
+
+The ligand solid angle, Ω is the solid angle of the complete shadow cast by a
+ligand when hypothetically illuminated from the metal center.
+:footcite:`bilbrey_solid_2013` From Ω, a solid "cone" angle Θ can also be
+calculated, which is analogous to Tolman's cone angle. The G parameter is a
+measure of how much of the metal coordination sphere is shielded by the
+ligands. :footcite:`guzei_2006`
+
+ᴍᴏʀғᴇᴜs uses a numerical method, in a similar way to the Solid-G__ program, but
+with different radii. The results have been benchmarked against the exact solid
+angle program by Allen and co-workers, :footcite:`bilbrey_solid_2013` and agree
+within numerical accuracy.
+
+.. __: https://xray.chem.wisc.edu/solid-g/
+
+.. footbibliography::

morfeus/__init__.py

@@ -1,6 +1,7 @@
 """Calculates steric descriptors for molecules.
 
 Modules:
+    bite_angle: Bite angle code
     buried_volume: Buried volume code
     calculators: Internal calculators
     cone_angle: Cone angle code
@@ -14,6 +15,7 @@
     pyramidalization: Pyramidalization code
     qc: Interface to quantum-chemical programs
     sasa: Solvent accessible surface area code
+    solid_angle: Solid angle code
     sterimol: Sterimol code
     typing: Typing code for arrays
     utils: Helper functions.
@@ -33,6 +35,7 @@
 from morfeus.local_force import LocalForce
 from morfeus.pyramidalization import Pyramidalization
 from morfeus.sasa import SASA
+from morfeus.solid_angle import SolidAngle
 from morfeus.sterimol import Sterimol
 from morfeus.visible_volume import VisibleVolume
 from morfeus.xtb import XTB
@@ -48,6 +51,7 @@
     "LocalForce",
     "Pyramidalization",
     "SASA",
+    "SolidAngle",
     "Sterimol",
     "VisibleVolume",
     "XTB",

morfeus/__main__.py

@@ -9,6 +9,7 @@
 import morfeus.local_force
 import morfeus.pyramidalization
 import morfeus.sasa
+import morfeus.solid_angle
 import morfeus.sterimol
 import morfeus.visible_volume
 import morfeus.xtb
@@ -26,6 +27,7 @@ def main() -> None:
             "local_force": morfeus.local_force.cli,
             "pyramidalization": morfeus.pyramidalization.cli,
             "sasa": morfeus.sasa.cli,
+            "solid_angle": morfeus.solid_angle.cli,
             "sterimol": morfeus.sterimol.cli,
             "visible_volume": morfeus.visible_volume.cli,
             "xtb": morfeus.xtb.cli,

morfeus/bite_angle.py

@@ -4,12 +4,12 @@
 
 from collections.abc import Iterable
 import functools
-from typing import Any
+from typing import Any, cast
 
 import numpy as np
 
 from morfeus.io import read_geometry
-from morfeus.typing import Array2DFloat, ArrayLike1D, ArrayLike2D
+from morfeus.typing import Array1DFloat, Array2DFloat, ArrayLike1D, ArrayLike2D
 
 
 class BiteAngle:
@@ -58,18 +58,18 @@ def __init__(
         v_2_norm = v_2 / np.linalg.norm(v_2)
 
         # Calculate angle between vectors
-        # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
         angle_rad = np.arctan2(
             np.linalg.norm(np.cross(v_1_norm, v_2_norm)), np.dot(v_1_norm, v_2_norm)
         )
-        angle = np.rad2deg(angle_rad)  # type: ignore
+        angle = np.rad2deg(angle_rad)
 
         # Check if angle should be inverted
         if ref_atoms is not None:
             ref_vector = (
                 np.mean(coordinates[[i - 1 for i in ref_atoms]], axis=0)
                 - coordinates[metal_index - 1]
             )
+        ref_vector = cast(Array1DFloat, ref_vector)
         inverted = False
         if ref_vector is not None:
             ref_vector /= np.linalg.norm(ref_vector)

morfeus/buried_volume.py

@@ -156,9 +156,8 @@ def __init__(
 
             v_1 = z_point - center
             v_2 = xz_point - center
-            # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
             v_3: Array1DFloat = np.cross(v_2, v_1)
-            real: Array2DFloat = np.vstack([v_1, v_3])  # type: ignore
+            real: Array2DFloat = np.vstack([v_1, v_3])
             real /= np.linalg.norm(real, axis=1).reshape(-1, 1)
             ref_1 = np.array([0.0, 0.0, -1.0])
             ref_2 = np.array([0.0, 1.0, 0.0])

morfeus/cone_angle.py

@@ -68,7 +68,7 @@ def __init__(  # noqa: C901
         radii_type: str = "crc",
         method: str = "libconeangle",
     ) -> None:
-        # Convert elements to atomic numbers if the are symbols
+        # Convert elements to atomic numbers if they are symbols
         elements = convert_elements(elements, output="numbers")
         coordinates: Array2DFloat = np.array(coordinates)
 

morfeus/geometry.py

@@ -600,9 +600,8 @@ def get_b_vector(self, coordinates: ArrayLike2D) -> Array1DFloat:
         phi = np.arccos(dot_product) * sgn
         g: Array1DFloat
         if abs(phi) > np.pi - 1e-6:
-            # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
             g = np.cross(w, a_1)
-            g = g / np.linalg.norm(g)  # type: ignore
+            g = g / np.linalg.norm(g)
             A = np.dot(v_1, ew) / np.linalg.norm(w)
             B = np.dot(v_2, ew) / np.linalg.norm(w)
             grad = [
@@ -612,9 +611,8 @@ def get_b_vector(self, coordinates: ArrayLike2D) -> Array1DFloat:
                 g / (np.linalg.norm(g) * np.linalg.norm(a_2)),
             ]
         elif abs(phi) < 1e-6:
-            # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
             g = np.cross(w, a_1)
-            g = g / np.linalg.norm(g)  # type: ignore
+            g = g / np.linalg.norm(g)
             A = np.dot(v_1, ew) / np.linalg.norm(w)
             B = np.dot(v_2, ew) / np.linalg.norm(w)
             grad = [
@@ -822,16 +820,14 @@ def rotate_coordinates(
 
     # Get imaginary dimensions and handle case of antiparallel vectors
     if real < 1e-6:
-        # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
         w: Array1DFloat = np.cross(vector, np.array([0, 0, 1]))
-        if np.linalg.norm(w) < 1e-6:  # type: ignore
+        if np.linalg.norm(w) < 1e-6:
             w = np.cross(vector, np.array([1, 0, 0]))
     else:
         w = np.cross(vector, axis)
 
     # Form quaternion and normalize
-    # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
-    q = np.concatenate((w, real))  # type: ignore
+    q = np.concatenate((w, real))
     q = q / np.linalg.norm(q)
 
     # Rotate atomic coordinates

morfeus/pyramidalization.py

@@ -5,7 +5,7 @@
 from collections.abc import Iterable, Sequence
 import functools
 import itertools
-from typing import Any
+from typing import Any, cast
 
 import numpy as np
 import scipy.spatial
@@ -127,9 +127,8 @@ def __init__(  # noqa: C901
         thetas: list[float] = []
         for v_1, v_2, v_3 in itertools.permutations([a, b, c], 3):
             # Calculate cos_alpha
-            # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
             normal: Array1DFloat = np.cross(v_1, v_2)
-            normal /= np.linalg.norm(normal)  # type: ignore
+            normal /= np.linalg.norm(normal)
             cos_alpha = np.dot(v_3, normal)
 
             # Test if normal vector is colinear with v_3
@@ -154,9 +153,9 @@ def __init__(  # noqa: C901
 
         # Calculate P
         v_1, v_2 = vectors[0]
-        # TODO: Remove type ignores when https://github.com/numpy/numpy/pull/21216 is released
         sin_theta = np.linalg.norm(np.cross(v_1, v_2))
-        P = sin_theta * cos_alphas[0]  # type: ignore
+        sin_theta = cast(float, sin_theta)
+        P = sin_theta * cos_alphas[0]
 
         # Correct P if pyramid is "acute" on average
         if np.mean(alphas) < 0: