Adding Critical Points (Vectorized) and Finding Bond Paths

[Ali-Tehrani]

Description

• Added my own approach of finding critical points that is significantly faster than the current method.
• Added an option to use the logarithm of the density and (hence gradient and hessian) when finding critical points using the Newton-Ralphson equation. This should increase speed of convergence since it becomes quadratic near the centers
• Added finding bond paths (using vectorization approach)
• Added ode.py with Runge-Kutta 4(5) with adaptive step-size and added bond-path finding.
• Added tests for critical points and bond paths which weren't there before.

The critical points algorithm is as follows: 1. Points that are close to center and that have small gradients are included to be initial guesses. Points whose density values is small are not. 2. While-loop is done until each point from the initial guess converges to a critical point using Newton-Ralphson (Armijo backtracing is not used). 3. Points whose density values are too small are removed from the iteration process. 4. When the while-loop is finished. Points that are identical are merged together by first rounding to some decimal places and then a KDTree is used to group the points into small balls and the points with the smallest gradient are picked out of each ball.

• The poincare-Hopf equation is checked, if not then a warning is raised.
• If any of the gradients are smaller than 1e-5, then a warning is raised.

The API follows the previous critical point finder:

mol = Molecule.from_file(file_path)
centers = mol.coordinates

# Generate points
pts = MolecularGrid.from_molecule(mol).points

# Construct Topology Objectt
result = Topology(
    mol.compute_density,
    mol.compute_gradient,
    mol.compute_hessian,
    pts
 )
result.find_critical_points_vectorized(centers, verbose=False)
result.find_bond_paths(max_ss=1e-3)

Examples

It takes four seconds (with ChemtoolsCUDA) to find critical points of dipeptide ALA_ALA (23 atoms) and 20 seconds for PHE_TRP (53 atoms and def2-SPVD basis-set) It took 5 seconds to find the bond paths and note the bcp and rcp found between pi-teeing (T-shaped) side-chains :

Here is 2014 GPU implementation of other people code (with basis-set 6-31++G(2p,2d) ). Note that our method seems to be significantly faster than theirs)

Here is a bond-path version

Type of Changes

Please remove the lines that don't represent the type of your PR.

:sparkles: New Feature