QTAIM Surface Algorithm

[Ali-Tehrani]

Description

Thursday, April 27, 2023: Mostly complete, needs refactoring, some features to be added, in progress

QTAIM is a theory that partitions the electron density based on the gradient path of the density. Each region is generally associated with an atom of the molecule, called the atomic basin of the atom. These basins all contain a boundary that borders each other. The boundary contains two parts, the interatomic surface where the boundary borders another atomic basin and the outeratomic surface which is the intersection of the isosurface of the electron density (at some value) and the atomic basin. The following is an example of a complicated hydrogen atomic basin that borders two atoms where the red points are the OAS points, black points are the IAS points, purple points are "refined" points, and green are the nuclei.

This pull-requests contains two algorithms (grid-based method of Henklemen [1] and Yu-trinkle [2] and a variant of Keith's PROMEGA algorithm [3] ) for finding atomic basins. The former is included in the file yu_trinkle.py and the latter in qtaim.py. This pull-requests requires heavy refactoring/cleaning and a couple of days of discussions before it should be merged. The rest of the pull request would be talking about the PROMEGA algorithm.

Requirements

1. Find the interatomic surface (IAS) and the outer atomic surface (OAS) of each atom.
2. Identify which point on the interatomic surface borders which atom.
3. Accurate IAS/OAS surfaces based on any accuracy.
4. Identify non-nuclear attractors.
5. Convergent and accurate integration over atomic basins.
6. Fast and reasonable speed.

General Algorithm

1. Over each atom, construct an angular grid, i.e. points on the sphere centered over each nucleus. The higher the degree, the more points are on the sphere.
2. Construct radial grids over each atom. This alongside the angular grids is used to construct rays centered over each nucleus.
3. Using a very small angular grid and dense radial grid, find the beta-sphere/capture-sphere for each atom: the sphere where all points inside it converge to that maximum. This reduces the need to backtrace all the way to the maximum along the gradient path. During this time, find any non-nuclear attractors.
4. Using a user-provided angular grid and coarse radial grid, classify each ray as either crossing another atomic basin (i.e. it crosses the IAS) or going towards infinity (crosses the OAS). Solving for which atomic basin is done by the Runge-Kutta of order 4 (5) with adaptive step size.
5. The rays that cross IAS are further refined until they satisfy some boundary error. This refinement is done by moving each point alongside the ray based on which atomic basin it converges to.
6. The rays that cross OAS are solved using in-house grid-based code for the isosurface point. This is preferred because scipy.optimize does not have a vectorized root-solver. Originally, the Newton-Ralpson was used since the logarithm of the electron density should be quadratic at large distances, where the Newton-Ralphson excels. However, this may not give the first root, but the second or third root.
7. Points on the IAS whose density values are smaller than the isosurface value, are re-assigned as OAS points.

The following tricks are used to improve performance:

• Everything is vectorized over all atoms. This requires the code to know how much available memory is there and how many points it can safely compute. This is done using in-house GPU code. There is a variant of an algorithm in qtaim_depreciated.py that doesn't require vectorization with memory management. However, the algorithm and its parameters are not optimized nor handle edge cases.
• This requires grouping all points of each ray of each atom together and having a way of tracking which points correspond to which atom. This is redone multiple times in different fashions and would benefit from being refactored into a single class that can work for all scenarios.
• If the ODE solver takes a step that decreases density values, then the ODE step is repeated with a smaller step size. This idea is obtained from the paper [4].
• Use of beta-sphere/capture-sphere/atomic trust sphere.
• Written my own Runge-Kutta and Newton-Ralphson solver, since SciPy lacks vectorized/independent ODE/root-solver.
• Can specify which atoms to find the surfaces for.
• Symmetric spherical t-design is used for the angular grid instead of Lebedev-Laikov grids. It can handle higher degrees and is evenly spaced out on the sphere.
• Refinement points (the purpose points in the pictures) are implemented that add points between the IAS and OAS. This is done by doing a convex hull on the IAS/OAS, then the distance between the vertices of the IAS and OAS is computed. The smallest distances are obtained, and for each point of the IAS, the two closest points on the OAS are obtained, where 9 points are constructed based on the triangle formed. These nine points are then solved to be either on IAS or OAS.

With P100 GPU, it takes anywhere from 5 to 15 minutes to solve for all surfaces for simple organic molecules up to 25 atoms with cc-pvtz basis set, high angular degrees of 107 for carbon, oxygen and nitrogen and 135 for hydrogen atoms and boundary error of 1e-5. This roughly corresponds to about 6K to 9K points on each surface.

Further Improvements Needed

• Newton-Ralphson solver is problematic where the ray intersects the isosurface of the electron density at more than one point close together or when the initial guess is far away. This motivated me to do a grid-search method for solving for the isosurface value.
• Multiple intersections of the IAS are not implemented but are found in step 4. It gives the NonImplementedError if the user wants them.
• Non-nuclear attractors are only found in step 3, and not in later points of the algorithm like in step 4, where it is more likely to find them. For now, the user needs to find the NNA before running this algorithm.
• The ODE solver can't handle the case where the points converge to a bond-critical or ring-critical point. Although it is unlikely to occur, it can happen and it just returns an error.
• A variant of an idea mentioned in [3] is to solve it for a smaller angular grid. Here it might be worthwhile to have an adaptive beta-sphere/capture-sphere. This could be done by solving for the surface for a smaller angular grid, then interpolated linearly or scipy.interpolate.LSQSphereBivariateSpline) or RectSphereBivariateSpline. Additionally, worthwhile to add padding to improve the trust.
• (Oct 23, 2023) There is a scenario where the beta-spheres (Neon in silica molecules sion) is set to such a large radius, that it's density values are less than the isosurface value. This causes problems in _classify_rays_as_ias_or_oas, where there are no pooints to classify the ray as, i.e. numb_rays_to_atom is zero at some places. This can be fixed by automatically classify the point as oas if the numb_rays_to_atom is zero because that's where it occurs. The other fix which is easier is to stop the beta-sphere at a smaller radius if any of the points on the radial shell has density less than the isosurface value, but this decreases computational performance, as it isn't the optimal beta-sphere.

API

Each calculation returns a Surface object from the file surface.py. The following is an example for methane:

from chemtools.topology.qtaim import qtaim_surface
from chemtools.wrappers import Molecule

mol = Molecule.from_file(METHANE FCHK HERE)
surface = qtaim_surface(
  angular=[70, 40, 40, 40, 40],                     # Angular degrees of each atom in methane.
  centers=mol.coordinates,                         # Atomic coordinates of methane
  dens_func=mol.compute_density,
  grad_func=mol.compute_gradient,
  iso_val=0.001,                                         # Isosurface value for the OAS
  bnd_err=1e-4,   # IAS Boundary error, controls the accuracy of the surface 
  iso_err = 1e-6,     # Controls the accuracy of the OAS
  beta_spheres=None,      # If None, finds beta-spheres automatically
  beta_sphere_deg=30,       #  The angular grid for finding beta-spheres, should be smaller than `angular`
  optimize_centers=True,   # if true, then optimizes the centers to the actual local maxima.
  hess_func=mol.compute_hessian,    # optional, used to find non-nuclear attractors.
  maximas_to_do=None         # Which atoms to find the surface for.  If None, then finds the surface for all.
)

# Atomic grid for integration
atomic_grid_carbon = surface.get_atom_grid_over_basin(0)
atomic_grid_hydrogen = surface.get_atom_grid_over_basin(1)
print("Density of carbon ", atomic_grid_carbon.integrate(Density values here ))

# Generate points on the surface, and include points from other atomic basins as part of the surface.
surface_carbon = surface.generate_pts_on_surface(0, include_other_surfaces=True)
ias_pts_carbon = surface.get_ias_pts_of_basin(0)
oas_pts_carbon = surface.get_oas_pts_of_basin(0)

# Get the surface of carbon that intersects the second hydrogen:
surface_pts = ias_pts_carbon[surface.basins_ias[0] == 3] 

# Save the surface object information
surface.save("methane_example")

# Two surfaces of the same molecule but with different angular degrees can be added together. The surfaces in `surface_new` replace the corresponding surface in `surface_old` 
#   Useful for avoiding to re-do surfaces that are considered "converged"
surface = surface_old + surface_new

Tests

The four tests in test_qtaim.py are performed over H2O, CH4, NH3 over a small angular grid size, the time for each test is around a minute except for the last one:

• Sum of the atomic density values is equal to the molecular charge up to two decimal places
• The integral of the Laplacian over the atomic basin (divided by four) is less than three decimal places
• The points on the OAS has the correct isosurface value
• Points inside and outside of the IAS correctly converge to the right atomic basin.

Checks

Please remove the checks that are not applicable to your PR, and replace [ ] with [x] to mark the
checks that are already completed. You can finish the other relevant checks after opening the PR.

• Each commit addresses one thing.
• Each commit message is informative and is less than 50 characters.
• Commits that can be grouped together are squashed.
• Each unit of code added is tested.
• Each unit of code added is documented.
• Commits are rebased onto master.

Type of Changes

✨ New Feature

References

[1] - Henkelman, Graeme, Andri Arnaldsson, and Hannes Jónsson. "A fast and robust algorithm for Bader decomposition of charge density." Computational Materials Science 36.3 (2006): 354-360.

[2] - Yu, Min, and Dallas R. Trinkle. "Accurate and efficient algorithm for Bader charge integration." The Journal of chemical physics 134.6 (2011): 064111.

[3] - T.A. Keith, Chapter 4 of Ph.D. Thesis: Molecules in Magnetic Fields, 176-213 (1993)

[4] - Rodríguez, Juan I., et al. "An efficient grid‐based scheme to compute QTAIM atomic properties without explicit calculation of zero‐flux surfaces." Journal of computational chemistry 30.7 (2009): 1082-1092.