Gaussian Iterative Stockholder Analysis (GISA) method
[1]:
import numpy as np
from scipy.optimize import least_squares, minimize
from setup import prepare_argument_dict, prepare_grid_and_dens, print_results
from horton_part import GaussianISAWPart
Quadprog Solver
The problem of this solver is that the initials values are set to zeros for all pro-atom parameters and no threshold setting is availiable.
Using API
[2]:
mol, grid, rho = prepare_grid_and_dens("data/h2o.fchk")
def quadprog_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = "quadprog"
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
quadprog_solver()
The number of electrons: 10.000003764139395
Coordinates of the atoms
[[ 0. 0. 0.224]
[-0. 1.457 -0.896]
[-0. -1.457 -0.896]]
Atomic numbers of the atom
[8 1 1]
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-08
Maximum iterations : 1000
lmax : 3
Solver : quadprog
Grid type : 1
Iteration Change Entropy
1 2.11650e-01 2.11396e-01
2 2.49517e-02 1.12635e-01
3 1.32558e-02 8.50056e-02
4 8.55599e-03 6.90944e-02
5 5.91191e-03 5.95412e-02
6 4.24249e-03 5.35617e-02
7 3.12465e-03 4.96967e-02
8 2.34690e-03 4.71329e-02
9 1.79038e-03 4.53930e-02
10 1.38337e-03 4.41873e-02
11 1.08043e-03 4.33355e-02
12 8.51611e-04 4.27227e-02
13 6.76617e-04 4.22746e-02
14 5.41337e-04 4.19419e-02
15 4.35763e-04 4.16914e-02
16 3.52679e-04 4.15003e-02
17 2.86803e-04 4.13530e-02
18 2.34217e-04 4.12381e-02
19 1.91987e-04 4.11477e-02
20 1.57888e-04 4.10759e-02
21 1.30221e-04 4.10186e-02
22 1.07673e-04 4.09724e-02
23 8.92255e-05 4.09350e-02
24 7.40810e-05 4.09045e-02
25 6.16098e-05 4.08796e-02
26 5.13121e-05 4.08592e-02
27 4.27890e-05 4.08424e-02
28 3.57201e-05 4.08285e-02
29 2.98466e-05 4.08169e-02
30 2.49587e-05 4.08074e-02
31 2.08855e-05 4.07994e-02
32 1.74872e-05 4.07928e-02
33 1.46491e-05 4.07873e-02
34 1.22769e-05 4.07827e-02
35 1.02925e-05 4.07788e-02
36 8.63147e-06 4.07756e-02
37 7.24041e-06 4.07729e-02
38 6.07489e-06 4.07706e-02
39 5.09794e-06 4.07687e-02
40 4.27879e-06 4.07671e-02
41 3.59174e-06 4.07658e-02
42 3.01536e-06 4.07647e-02
43 2.53171e-06 4.07637e-02
44 2.12581e-06 4.07630e-02
45 1.78511e-06 4.07623e-02
46 1.49911e-06 4.07617e-02
47 1.25898e-06 4.07613e-02
48 1.05737e-06 4.07609e-02
49 8.88068e-07 4.07606e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.863 0.432 0.432]
cartesian multipoles:
[[-0.863 0. -0. 0.113 -5.313 -0. -0. -4.718 -0. -5.013 -0. -0. -0.039 0. -0. 0. -0. 0.322 -0. 0.155]
[ 0.432 0. 0.023 -0.001 -0.309 -0. -0. -0.282 0.002 -0.288 0. 0.013 0.003 0. 0. 0. 0.039 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.309 0. 0. -0.282 -0.002 -0.288 0. -0.013 0.003 0. -0. 0. -0.039 0.009 -0.015 0.02 ]]
Customized Implementation
[3]:
import quadprog
from horton_part import check_pro_atom_parameters
def _prepare_quantities(bs_funcs, rho, propars, weights, alphas):
nprim = len(propars)
G = (
2
/ np.pi**1.5
* (alphas[:, None] * alphas[None, :]) ** 1.5
/ (alphas[:, None] + alphas[None, :]) ** 1.5
)
G = (G + G.T) / 2
# a = np.zeros(nprim, float)
# for k in range(nprim):
# a[k] = 2 * np.einsum("i,i,i", weights, bs_funcs[k], rho)
a = 2 * np.einsum("i,ni,i->n", weights, bs_funcs, rho)
# Construct linear equality or inequality constraints
C = np.zeros([nprim, nprim + 1])
# First column : corresponds to the EQUALITY constraint sum_{k=1..Ka} c_(a,k)= N_a
C[:, 0] = np.ones(nprim)
# Other K_a columns : correspond to the INEQUALITY constraints c_(a,k) >=0
C[0:nprim, 1 : (nprim + 1)] = np.identity(nprim)
b = np.zeros(nprim + 1)
# First coefficient : corresponds to the EQUALITY constraint sum_{k=1..Ka} c_(a,k) = N_a
pop = np.einsum("i,i", weights, rho)
b[0] = pop
return G, a, C, b
def opt_propars_qp_quadprog(bs_funcs, rho, propars, points, weights, alphas, threshold):
r"""
Optimize parameters for proatom density functions using quadratic programming.
.. math::
N_{Ai} = \int \rho_A(r) \frac{\rho_{Ai}^0(r)}{\rho_A^0(r)} dr
G = \frac{1}{2} c^T S c - c^T b
S = 2 \int \zeta(\vec{r}) \zeta(\vec{r}) d\vec{r}
= \frac{2}{\pi \sqrt{\pi}} \frac{(\alpha_k \alpha_l)^{3/2}}{(\alpha_k + \alpha_l)^{3/2}}
b = 2 * \int \zeta(\vec{r}) \rho_a(\vec{r}) d\vec{r}
Parameters
----------
bs_funcs : 2D np.ndarray
Basis functions array with shape (M, N), where 'M' is the number of basis functions
and 'N' is the number of grid points.
rho : 1D np.ndarray
Spherically-averaged atomic density as a function of radial distance, with shape (N,).
propars : 1D np.ndarray
Pro-atom parameters with shape (M). 'M' is the number of basis functions.
points : 1D np.ndarray
Radial coordinates of grid points, with shape (N,).
weights : 1D np.ndarray
Weights for integration, including the angular part (4πr²), with shape (N,).
threshold : float
Convergence threshold for the iterative process.
alphas : 1D np.ndarray
The Gaussian exponential coefficients.
threshold : float
The convergence threshold of the optimization method.
Returns
-------
1D np.ndarray
Optimized proatom parameters.
Raises
------
RuntimeError
If the inner iteration does not converge.
"""
G, a, C, b = _prepare_quantities(bs_funcs, rho, propars, weights, alphas)
result = quadprog.solve_qp(G=G, a=a, C=C, b=b, meq=1)[0]
check_pro_atom_parameters(result, total_population=float(b[0]))
return result
Now we can call the customized solver through follows:
[4]:
def customized_quadprog_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = opt_propars_qp_quadprog
kwargs["inner_threshold"] = 1e-8
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
customized_quadprog_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-08
Maximum iterations : 1000
lmax : 3
Solver : <function opt_propars_qp_quadprog at 0x10880fe20>
Grid type : 1
Iteration Change Entropy
1 2.11650e-01 2.11396e-01
2 2.49517e-02 1.12635e-01
3 1.32558e-02 8.50056e-02
4 8.55599e-03 6.90944e-02
5 5.91191e-03 5.95412e-02
6 4.24249e-03 5.35617e-02
7 3.12465e-03 4.96967e-02
8 2.34690e-03 4.71329e-02
9 1.79038e-03 4.53930e-02
10 1.38337e-03 4.41873e-02
11 1.08043e-03 4.33355e-02
12 8.51611e-04 4.27227e-02
13 6.76617e-04 4.22746e-02
14 5.41337e-04 4.19419e-02
15 4.35763e-04 4.16914e-02
16 3.52679e-04 4.15003e-02
17 2.86803e-04 4.13530e-02
18 2.34217e-04 4.12381e-02
19 1.91987e-04 4.11477e-02
20 1.57888e-04 4.10759e-02
21 1.30221e-04 4.10186e-02
22 1.07673e-04 4.09724e-02
23 8.92255e-05 4.09350e-02
24 7.40810e-05 4.09045e-02
25 6.16098e-05 4.08796e-02
/Users/runner/work/horton-part/horton-part/src/horton_part/gisa.py:192: UserWarning: Customized solver is used, the argument `inner_threshold` is not used.
warnings.warn("Customized solver is used, the argument `inner_threshold` is not used.")
26 5.13121e-05 4.08592e-02
27 4.27890e-05 4.08424e-02
28 3.57201e-05 4.08285e-02
29 2.98466e-05 4.08169e-02
30 2.49587e-05 4.08074e-02
31 2.08855e-05 4.07994e-02
32 1.74872e-05 4.07928e-02
33 1.46491e-05 4.07873e-02
34 1.22769e-05 4.07827e-02
35 1.02925e-05 4.07788e-02
36 8.63147e-06 4.07756e-02
37 7.24041e-06 4.07729e-02
38 6.07489e-06 4.07706e-02
39 5.09794e-06 4.07687e-02
40 4.27879e-06 4.07671e-02
41 3.59174e-06 4.07658e-02
42 3.01536e-06 4.07647e-02
43 2.53171e-06 4.07637e-02
44 2.12581e-06 4.07630e-02
45 1.78511e-06 4.07623e-02
46 1.49911e-06 4.07617e-02
47 1.25898e-06 4.07613e-02
48 1.05737e-06 4.07609e-02
49 8.88068e-07 4.07606e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.863 0.432 0.432]
cartesian multipoles:
[[-0.863 0. -0. 0.113 -5.313 -0. -0. -4.718 -0. -5.013 -0. -0. -0.039 0. -0. 0. -0. 0.322 -0. 0.155]
[ 0.432 0. 0.023 -0.001 -0.309 -0. -0. -0.282 0.002 -0.288 0. 0.013 0.003 0. 0. 0. 0.039 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.309 0. 0. -0.282 -0.002 -0.288 0. -0.013 0.003 0. -0. 0. -0.039 0.009 -0.015 0.02 ]]
CVXOPT Solver
Another sovler implemented in CVXOPT can also be applied.
Using API
[5]:
def cvxopt_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = "cvxopt"
kwargs["solver_options"] = {"eps_rel": 1e-8}
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
opt_propars = part._cache["propars"]
print("Difference between total electrons:")
print(grid.integrate(rho) - np.sum(opt_propars))
cvxopt_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-08
Maximum iterations : 1000
lmax : 3
Solver : cvxopt
Grid type : 1
Iteration Change Entropy
1 2.11603e-01 2.11396e-01
2 2.49842e-02 1.12501e-01
3 1.32921e-02 8.48478e-02
4 8.57237e-03 6.89801e-02
5 5.91785e-03 5.94610e-02
6 4.24425e-03 5.35056e-02
7 3.12533e-03 4.96581e-02
8 2.34836e-03 4.71079e-02
9 1.79441e-03 4.53802e-02
10 1.37842e-03 4.41882e-02
11 1.07512e-03 4.33343e-02
12 8.48063e-04 4.27203e-02
13 6.74378e-04 4.22717e-02
14 5.39965e-04 4.19391e-02
15 4.34960e-04 4.16889e-02
16 3.52247e-04 4.14983e-02
17 2.86606e-04 4.13515e-02
18 2.34163e-04 4.12373e-02
19 1.92009e-04 4.11475e-02
20 1.57947e-04 4.10765e-02
21 1.30281e-04 4.10198e-02
22 1.07710e-04 4.09742e-02
23 8.92219e-05 4.09375e-02
24 7.40284e-05 4.09077e-02
25 6.15061e-05 4.08834e-02
26 5.11599e-05 4.08635e-02
27 4.25939e-05 4.08472e-02
28 3.54894e-05 4.08338e-02
29 2.95885e-05 4.08227e-02
30 2.46814e-05 4.08136e-02
31 2.05967e-05 4.08060e-02
32 1.71939e-05 4.07997e-02
33 1.43571e-05 4.07945e-02
34 1.19911e-05 4.07902e-02
35 1.00168e-05 4.07866e-02
36 8.36874e-06 4.07836e-02
37 6.99267e-06 4.07811e-02
38 5.84343e-06 4.07790e-02
39 4.88343e-06 4.07772e-02
40 4.08141e-06 4.07758e-02
41 3.41127e-06 4.07746e-02
42 2.85128e-06 4.07736e-02
43 2.38330e-06 4.07727e-02
44 1.99218e-06 4.07720e-02
45 1.66528e-06 4.07714e-02
46 1.39205e-06 4.07709e-02
47 1.16366e-06 4.07705e-02
48 9.72756e-07 4.07702e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.863 0.432 0.432]
cartesian multipoles:
[[-0.863 0. -0. 0.112 -5.312 -0. -0. -4.718 -0. -5.013 0. -0. -0.04 0. 0. 0. -0. 0.322 -0. 0.153]
[ 0.432 0. 0.023 -0.001 -0.309 -0. -0. -0.283 0.002 -0.288 0. 0.013 0.003 0. 0. 0. 0.04 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.309 0. 0. -0.283 -0.002 -0.288 0. -0.013 0.003 0. -0. 0. -0.04 0.009 -0.015 0.02 ]]
Difference between total electrons:
1.99248908909766e-05
Customized Implementation
[6]:
import cvxopt
def opt_propars_qp_cvxopt(bs_funcs, rho, propars, points, weights, alphas, threshold):
"""
Optimize pro-atom parameters using quadratic-programming implemented in the `CVXOPT` package.
Parameters
----------
bs_funcs : 2D np.ndarray
Basis functions array with shape (M, N), where 'M' is the number of basis functions
and 'N' is the number of grid points.
rho : 1D np.ndarray
Spherically-averaged atomic density as a function of radial distance, with shape (N,).
propars : 1D np.ndarray
Pro-atom parameters with shape (M). 'M' is the number of basis functions.
points : 1D np.ndarray
Radial coordinates of grid points, with shape (N,).
weights : 1D np.ndarray
Weights for integration, including the angular part (4πr²), with shape (N,).
threshold : float
Convergence threshold for the iterative process.
alphas : 1D np.ndarray
The Gaussian exponential coefficients.
threshold : float
The convergence threshold of the optimization method.
Returns
-------
1D np.ndarray
Optimized proatom parameters.
Raises
------
RuntimeError
If the inner iteration does not converge.
"""
nprim, npt = bs_funcs.shape
S = (
2
/ np.pi**1.5
* (alphas[:, None] * alphas[None, :]) ** 1.5
/ (alphas[:, None] + alphas[None, :]) ** 1.5
)
P = cvxopt.matrix(S)
vec_b = np.zeros((nprim, 1), float)
for k in range(nprim):
vec_b[k] = 2 * np.einsum("i,i", weights * bs_funcs[k], rho)
q = -cvxopt.matrix(vec_b)
# Linear inequality constraints
G = cvxopt.matrix(0.0, (nprim, nprim))
G[:: nprim + 1] = -1.0
h = cvxopt.matrix(0.0, (nprim, 1))
# Linear equality constraints
A = cvxopt.matrix(1.0, (1, nprim))
pop = np.einsum("i,i", weights, rho)
b = cvxopt.matrix(pop, (1, 1))
# initial_values = cvxopt.matrix(np.array([1.0] * nprim).reshape((nprim, 1)))
opt_CVX = cvxopt.solvers.qp(
P,
q,
G,
h,
A,
b,
# initvals=propars,
initvals=np.zeros_like(propars),
options={"feastol": threshold, "show_progress": 0},
# options={"show_progress": log.do_medium, "feastol": threshold},
)
new_propars = np.asarray(opt_CVX["x"]).flatten()
check_pro_atom_parameters(new_propars, total_population=float(pop))
return new_propars
[7]:
def customized_cvxopt_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = opt_propars_qp_cvxopt
kwargs["inner_threshold"] = 1e-8
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
opt_propars = part._cache["propars"]
print("Difference between total electrons:")
print(grid.integrate(rho) - np.sum(opt_propars))
customized_cvxopt_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-08
Maximum iterations : 1000
lmax : 3
Solver : <function opt_propars_qp_cvxopt at 0x1088948b0>
Grid type : 1
Iteration Change Entropy
1 2.11603e-01 2.11396e-01
2 2.49842e-02 1.12501e-01
3 1.32921e-02 8.48478e-02
4 8.57237e-03 6.89801e-02
5 5.91785e-03 5.94610e-02
6 4.24425e-03 5.35056e-02
7 3.12533e-03 4.96581e-02
8 2.34836e-03 4.71079e-02
9 1.79441e-03 4.53802e-02
10 1.37842e-03 4.41882e-02
11 1.07512e-03 4.33343e-02
12 8.48063e-04 4.27203e-02
13 6.74378e-04 4.22717e-02
14 5.39965e-04 4.19391e-02
15 4.34960e-04 4.16889e-02
16 3.52247e-04 4.14983e-02
17 2.86606e-04 4.13515e-02
18 2.34163e-04 4.12373e-02
19 1.92009e-04 4.11475e-02
20 1.57947e-04 4.10765e-02
21 1.30281e-04 4.10198e-02
22 1.07710e-04 4.09742e-02
23 8.92219e-05 4.09375e-02
24 7.40284e-05 4.09077e-02
25 6.15061e-05 4.08834e-02
26 5.11599e-05 4.08635e-02
27 4.25939e-05 4.08472e-02
28 3.54894e-05 4.08338e-02
29 2.95885e-05 4.08227e-02
30 2.46814e-05 4.08136e-02
31 2.05967e-05 4.08060e-02
32 1.71939e-05 4.07997e-02
33 1.43571e-05 4.07945e-02
34 1.19911e-05 4.07902e-02
35 1.00168e-05 4.07866e-02
36 8.36874e-06 4.07836e-02
37 6.99267e-06 4.07811e-02
38 5.84343e-06 4.07790e-02
39 4.88343e-06 4.07772e-02
40 4.08141e-06 4.07758e-02
41 3.41127e-06 4.07746e-02
42 2.85128e-06 4.07736e-02
43 2.38330e-06 4.07727e-02
44 1.99218e-06 4.07720e-02
45 1.66528e-06 4.07714e-02
46 1.39205e-06 4.07709e-02
47 1.16366e-06 4.07705e-02
48 9.72756e-07 4.07702e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.863 0.432 0.432]
cartesian multipoles:
[[-0.863 0. 0. 0.112 -5.312 -0. -0. -4.718 -0. -5.013 -0. -0. -0.04 0. -0. 0. -0. 0.322 -0. 0.153]
[ 0.432 0. 0.023 -0.001 -0.309 -0. -0. -0.283 0.002 -0.288 0. 0.013 0.003 0. 0. 0. 0.04 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.309 0. 0. -0.283 -0.002 -0.288 0. -0.013 0.003 0. -0. 0. -0.04 0.009 -0.015 0.02 ]]
Difference between total electrons:
1.9924890892752956e-05
ECOS Sovler
[8]:
def ecos_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = "ecos"
kwargs["solver_options"] = {"reltol": 1e-10}
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
opt_propars = part._cache["propars"]
print("Difference between total electrons:")
print(grid.integrate(rho) - np.sum(opt_propars))
ecos_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-08
Maximum iterations : 1000
lmax : 3
Solver : ecos
Grid type : 1
Iteration Change Entropy
1 2.11662e-01 2.11396e-01
2 2.49254e-02 1.12672e-01
3 1.32492e-02 8.50267e-02
4 8.53556e-03 6.91032e-02
5 5.93385e-03 5.95247e-02
6 4.24356e-03 5.35553e-02
7 3.11140e-03 4.96930e-02
8 2.33603e-03 4.71239e-02
9 1.82012e-03 4.53760e-02
10 1.39577e-03 4.41781e-02
11 1.08352e-03 4.33397e-02
12 8.45891e-04 4.27194e-02
13 6.72398e-04 4.22667e-02
14 5.38656e-04 4.19314e-02
15 4.34270e-04 4.16797e-02
16 3.52063e-04 4.14884e-02
17 2.86808e-04 4.13414e-02
18 2.34613e-04 4.12271e-02
19 1.92580e-04 4.11376e-02
20 1.57877e-04 4.10667e-02
21 1.36816e-04 4.10081e-02
22 1.13458e-04 4.09653e-02
23 9.30038e-05 4.09299e-02
24 7.62189e-05 4.09004e-02
25 6.26717e-05 4.08757e-02
26 5.17460e-05 4.08552e-02
27 4.28451e-05 4.08381e-02
28 3.55703e-05 4.08239e-02
29 2.96656e-05 4.08120e-02
30 2.47262e-05 4.08022e-02
31 2.06659e-05 4.07939e-02
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/qpsolvers/solvers/ecos_.py:126: UserWarning: warm-start values are ignored by this wrapper
warnings.warn("warm-start values are ignored by this wrapper")
32 1.73065e-05 4.07871e-02
33 1.45088e-05 4.07813e-02
34 1.21624e-05 4.07765e-02
35 1.02060e-05 4.07725e-02
36 8.54192e-06 4.07692e-02
37 7.18210e-06 4.07664e-02
38 6.06573e-06 4.07640e-02
39 5.04520e-06 4.07621e-02
40 4.27853e-06 4.07604e-02
41 3.55020e-06 4.07590e-02
42 3.03717e-06 4.07578e-02
43 2.51777e-06 4.07569e-02
44 2.12495e-06 4.07560e-02
45 1.79963e-06 4.07554e-02
46 1.49222e-06 4.07548e-02
47 1.26242e-06 4.07543e-02
48 1.07009e-06 4.07539e-02
49 8.93342e-07 4.07535e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.864 0.432 0.432]
cartesian multipoles:
[[-0.864 -0. 0. 0.113 -5.313 -0. -0. -4.719 -0. -5.014 0. -0. -0.039 -0. -0. 0. -0. 0.323 -0. 0.156]
[ 0.432 0. 0.023 -0.001 -0.309 -0. -0. -0.282 0.002 -0.288 0. 0.013 0.003 0. 0. 0. 0.039 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.309 0. 0. -0.282 -0.002 -0.288 0. -0.013 0.003 0. -0. 0. -0.039 0.009 -0.015 0.02 ]]
Difference between total electrons:
1.9957000491288568e-05
OSQP Solver
[9]:
def osqp_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = "osqp"
kwargs["solver_options"] = {"eps_rel": 1e-14}
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
opt_propars = part._cache["propars"]
print("Difference between total electrons:")
print(grid.integrate(rho) - np.sum(opt_propars))
osqp_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-08
Maximum iterations : 1000
lmax : 3
Solver : osqp
Grid type : 1
Iteration Change Entropy
1 2.11649e-01 2.11396e-01
2 2.49710e-02 1.12965e-01
3 1.32661e-02 8.52882e-02
4 8.56305e-03 6.93508e-02
5 5.91729e-03 5.97838e-02
6 4.24685e-03 5.37974e-02
7 3.12836e-03 4.99294e-02
8 2.35018e-03 4.73647e-02
9 1.79338e-03 4.56249e-02
10 1.38617e-03 4.44201e-02
11 1.08306e-03 4.35694e-02
12 8.54102e-04 4.29579e-02
13 6.78981e-04 4.25110e-02
14 5.43600e-04 4.21797e-02
15 4.37940e-04 4.19309e-02
16 3.54765e-04 4.17414e-02
17 2.88797e-04 4.15952e-02
18 2.36114e-04 4.14815e-02
19 1.93779e-04 4.13921e-02
20 1.59569e-04 4.13211e-02
21 1.31788e-04 4.12644e-02
22 1.09127e-04 4.12188e-02
23 9.05683e-05 4.11818e-02
24 7.53150e-05 4.11518e-02
25 6.27388e-05 4.11271e-02
26 5.23410e-05 4.11069e-02
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/qpsolvers/conversions/ensure_sparse_matrices.py:28: SparseConversionWarning: Converted matrix 'P' of your problem to scipy.sparse.csc_matrix to pass it to solver 'osqp'; for best performance, build your matrix as a csc_matrix directly.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/qpsolvers/conversions/ensure_sparse_matrices.py:28: SparseConversionWarning: Converted matrix 'G' of your problem to scipy.sparse.csc_matrix to pass it to solver 'osqp'; for best performance, build your matrix as a csc_matrix directly.
warnings.warn(
/Library/Frameworks/Python.framework/Versions/3.10/lib/python3.10/site-packages/qpsolvers/conversions/ensure_sparse_matrices.py:28: SparseConversionWarning: Converted matrix 'A' of your problem to scipy.sparse.csc_matrix to pass it to solver 'osqp'; for best performance, build your matrix as a csc_matrix directly.
warnings.warn(
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.94063E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.78186E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.67535E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.60093E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.54745E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.50816E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.47877E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.45646E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.44044E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.42800E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.41825E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.41053E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40618E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40489E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40397E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40331E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40326E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40336E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40348E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40360E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40372E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40384E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40395E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40405E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40413E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40421E-04
warnings.warn(info)
27 4.37231e-05 4.10903e-02
28 3.65650e-05 4.10765e-02
29 3.06083e-05 4.10651e-02
30 2.56433e-05 4.10556e-02
31 2.14990e-05 4.10477e-02
32 1.80354e-05 4.10411e-02
33 1.51377e-05 4.10356e-02
34 1.27113e-05 4.10310e-02
35 1.06778e-05 4.10271e-02
36 8.97254e-06 4.10239e-02
37 7.54170e-06 4.10212e-02
38 6.34051e-06 4.10189e-02
39 5.33170e-06 4.10170e-02
40 4.48416e-06 4.10154e-02
41 3.77188e-06 4.10141e-02
42 3.17314e-06 4.10129e-02
43 2.66971e-06 4.10120e-02
44 2.24634e-06 4.10112e-02
45 1.89026e-06 4.10105e-02
46 1.59072e-06 4.10099e-02
47 1.33872e-06 4.10095e-02
48 1.12669e-06 4.10091e-02
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40428E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40434E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40439E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40443E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40447E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40450E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40453E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40455E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40457E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40459E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40461E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40462E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40463E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40464E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40465E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40466E-04
warnings.warn(info)
/Users/runner/work/horton-part/horton-part/src/horton_part/utils.py:454: UserWarning: WARNING: The sum of pro-atom parameters is not equal to reference population.
WARNING: The difference is -1.40467E-04
warnings.warn(info)
49 9.48273e-07 4.10087e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.865 0.432 0.432]
cartesian multipoles:
[[-0.865 0. -0. 0.114 -5.316 -0. -0. -4.722 0. -5.017 -0. -0. -0.037 0. -0. 0. -0. 0.326 -0. 0.162]
[ 0.432 0. 0.023 -0.001 -0.308 -0. -0. -0.281 0.002 -0.287 0. 0.013 0.003 0. 0. 0. 0.038 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.308 0. 0. -0.281 -0.002 -0.287 0. -0.013 0.003 0. -0. 0. -0.038 0.009 -0.015 0.02 ]]
Difference between total electrons:
0.00030096872035656475
More supported sovlers can be found in `qpsovlers <https://qpsolvers.github.io/qpsolvers/supported-solvers.html>`__ package.
Customized SLSQP Solver
A equivalent method is using minimization with the ‘SLSQP’ sovler implemented in SciPy package.
[10]:
def opt_propars_qp_slsqp(bs_funcs, rho, propars, points, weights, alphas, threshold):
r"""
Optimize parameters for proatom density functions using quadratic programming.
.. math::
N_{Ai} = \int \rho_A(r) \frac{\rho_{Ai}^0(r)}{\rho_A^0(r)} dr
G = \frac{1}{2} c^T S c - c^T b
S = 2 \int \zeta(\vec{r}) \zeta(\vec{r}) d\vec{r}
= \frac{2}{\pi \sqrt{\pi}} \frac{(\alpha_k \alpha_l)^{3/2}}{(\alpha_k + \alpha_l)^{3/2}}
b = 2 * \int \zeta(\vec{r}) \rho_a(\vec{r}) d\vec{r}
Parameters
----------
bs_funcs : 2D np.ndarray
Basis functions array with shape (M, N), where 'M' is the number of basis functions
and 'N' is the number of grid points.
rho : 1D np.ndarray
Spherically-averaged atomic density as a function of radial distance, with shape (N,).
propars : 1D np.ndarray
Pro-atom parameters with shape (M). 'M' is the number of basis functions.
points : 1D np.ndarray
Radial coordinates of grid points, with shape (N,).
weights : 1D np.ndarray
Weights for integration, including the angular part (4πr²), with shape (N,).
threshold : float
Convergence threshold for the iterative process.
alphas : 1D np.ndarray
The Gaussian exponential coefficients.
threshold : float
The convergence threshold of the optimization method.
Returns
-------
1D np.ndarray
Optimized proatom parameters.
Raises
------
RuntimeError
If the inner iteration does not converge.
"""
G, a, C, b = _prepare_quantities(bs_funcs, rho, propars, weights, alphas)
meq = 1
def _obj_func(x):
return 0.5 * np.dot(x, G).dot(x) - np.dot(a, x)
constraints = []
if C is not None:
constraints = [
{
"type": "eq" if i < meq else "ineq",
"fun": lambda x, C=C, b=b, i=i: (np.dot(C.T, x) - b)[i],
}
for i in range(C.shape[1])
]
result = minimize(
_obj_func,
x0=np.zeros(len(G)),
method="SLSQP",
constraints=constraints,
tol=threshold,
options={"maxiter": 2000},
)
return result.x
[11]:
def customized_slsqp_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = opt_propars_qp_cvxopt
# the value is from the test file of quadprog package.
kwargs["inner_threshold"] = 1e-10
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
opt_propars = part._cache["propars"]
print("Difference between total electrons:")
print(grid.integrate(rho) - np.sum(opt_propars))
customized_slsqp_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-10
Maximum iterations : 1000
lmax : 3
Solver : <function opt_propars_qp_cvxopt at 0x1088948b0>
Grid type : 1
Iteration Change Entropy
1 2.11603e-01 2.11396e-01
2 2.49842e-02 1.12501e-01
3 1.32921e-02 8.48478e-02
4 8.57237e-03 6.89801e-02
5 5.91785e-03 5.94610e-02
6 4.24425e-03 5.35056e-02
7 3.12533e-03 4.96581e-02
8 2.34836e-03 4.71079e-02
9 1.79441e-03 4.53802e-02
10 1.37842e-03 4.41882e-02
11 1.07512e-03 4.33343e-02
12 8.48063e-04 4.27203e-02
13 6.74378e-04 4.22717e-02
14 5.39965e-04 4.19391e-02
15 4.34960e-04 4.16889e-02
16 3.52247e-04 4.14983e-02
17 2.86606e-04 4.13515e-02
18 2.34163e-04 4.12373e-02
19 1.92009e-04 4.11475e-02
20 1.57947e-04 4.10765e-02
21 1.30281e-04 4.10198e-02
22 1.07710e-04 4.09742e-02
23 8.92219e-05 4.09375e-02
24 7.40284e-05 4.09077e-02
25 6.15061e-05 4.08834e-02
26 5.11599e-05 4.08635e-02
27 4.25939e-05 4.08472e-02
28 3.54894e-05 4.08338e-02
29 2.95885e-05 4.08227e-02
30 2.46814e-05 4.08136e-02
31 2.05967e-05 4.08060e-02
32 1.71939e-05 4.07997e-02
33 1.43571e-05 4.07945e-02
34 1.19911e-05 4.07902e-02
35 1.00168e-05 4.07866e-02
36 8.36874e-06 4.07836e-02
37 6.99267e-06 4.07811e-02
38 5.84343e-06 4.07790e-02
39 4.88343e-06 4.07772e-02
40 4.08141e-06 4.07758e-02
41 3.41127e-06 4.07746e-02
42 2.85128e-06 4.07736e-02
43 2.38330e-06 4.07727e-02
44 1.99218e-06 4.07720e-02
45 1.66528e-06 4.07714e-02
46 1.39205e-06 4.07709e-02
47 1.16366e-06 4.07705e-02
48 9.72756e-07 4.07702e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.863 0.432 0.432]
cartesian multipoles:
[[-0.863 0. 0. 0.112 -5.312 -0. -0. -4.718 -0. -5.013 -0. -0. -0.04 0. -0. 0. -0. 0.322 -0. 0.153]
[ 0.432 0. 0.023 -0.001 -0.309 -0. -0. -0.283 0.002 -0.288 0. 0.013 0.003 0. 0. 0. 0.04 0.009 0.015 0.02 ]
[ 0.432 0. -0.023 -0.001 -0.309 0. 0. -0.283 -0.002 -0.288 0. -0.013 0.003 0. -0. 0. -0.04 0.009 -0.015 0.02 ]]
Difference between total electrons:
1.9924890892752956e-05
Customized Least-Square Sovler
Users can easily apply customized solvers. The first step involves defining a customized solver, for example, customized_solver. The corresponding arguments and their meanings are as follows:
[12]:
def opt_propars_lstsq(bs_funcs, rho, propars, points, weights, alphas, threshold):
r"""
Optimize pro-atom parameters using quadratic-programming implemented in the `CVXOPT` package.
.. math::
N_{Ai} = \int \rho_A(r) \frac{\rho_{Ai}^0(r)}{\rho_A^0(r)} dr
G = \frac{1}{2} c^T S c - c^T b
S = 2 \int \zeta(\vec{r}) \zeta(\vec{r}) d\vec{r}
= \frac{2}{\pi \sqrt{\pi}} \frac{(\alpha_k \alpha_l)^{3/2}}{(\alpha_k + \alpha_l)^{3/2}}
b = \int \zeta(\vec{r}) \rho_a(\vec{r}) d\vec{r}
Parameters
----------
bs_funcs : 2D np.ndarray
Basis functions array with shape (M, N), where 'M' is the number of basis functions
and 'N' is the number of grid points.
rho : 1D np.ndarray
Spherically-averaged atomic density as a function of radial distance, with shape (N,).
propars : 1D np.ndarray
Pro-atom parameters with shape (M). 'M' is the number of basis functions.
points : 1D np.ndarray
Radial coordinates of grid points, with shape (N,).
weights : 1D np.ndarray
Weights for integration, including the angular part (4πr²), with shape (N,).
threshold : float
Convergence threshold for the iterative process.
alphas : 1D np.ndarray
The Gaussian exponential coefficients.
threshold : float
The convergence threshold of the optimization method.
Returns
-------
1D np.ndarray
Optimized proatom parameters.
Raises
------
RuntimeError
If the inner iteration does not converge.
"""
def _obj_func(x):
pro = np.sum(x[:, None] * bs_funcs, axis=0)
return np.sqrt(weights) * np.abs(pro - rho)
result = least_squares(
_obj_func, x0=np.zeros_like(propars), bounds=(0, np.inf), ftol=threshold
)
return result.x
Next, we need to pass this customized solver to our partitioning method.
[13]:
def customized_lstsq_solver():
kwargs = prepare_argument_dict(mol, grid, rho)
kwargs["solver"] = opt_propars_lstsq
kwargs["inner_threshold"] = 1e-14
part = GaussianISAWPart(**kwargs)
part.do_all()
print_results(part)
opt_propars = part._cache["propars"]
print("Difference between total electrons:")
print(grid.integrate(rho) - np.sum(opt_propars))
customized_lstsq_solver()
================================================================================
Information of integral grids.
--------------------------------------------------------------------------------
Grid size of molecular grid: 18460
************************************ Atom 0 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18 18
18 18 18 18 18 18 26 26 26 26 38 38 38 38 38 50 50 50 50 50
86 86 110 110 110 110 170 194 194 434 590 590 434 434 434 302 302 302 194 194
170 110 110 110 110 110 110 110 110 110 110 110 86 50 50 18 18 18 18 18
************************************ Atom 1 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
************************************ Atom 2 ************************************
|-- Radial grid size: 120
|-- Angular grid sizes:
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6
6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 18 18 18 18 18
18 18 18 18 18 18 18 18 18 26 26 26 38 38 38 38 38 38 38 50
50 50 50 50 50 86 86 86 110 110 110 170 170 170 194 194 194 194 170 170
170 110 110 110 110 110 110 110 110 110 110 86 86 50 50 50 18 18 18 18
--------------------------------------------------------------------------------
================================================================================
Performing a density-based AIM analysis with a wavefunction as input.
Molecular grid : MolGrid
Using local grids : True
Scheme : Gaussian Iterative Stockholder Analysis (GISA)
Outer loop convergence threshold : 1.0e-06
Inner loop convergence threshold : 1.0e-14
Maximum iterations : 1000
lmax : 3
Solver : <function opt_propars_lstsq at 0x108897be0>
Grid type : 1
Iteration Change Entropy
1 2.11298e-01 2.11396e-01
2 2.46519e-02 1.80014e-01
3 1.28132e-02 1.26698e-01
4 8.42982e-03 9.80376e-02
5 5.91102e-03 8.25970e-02
6 4.27284e-03 7.36446e-02
7 3.15178e-03 6.83471e-02
8 2.36034e-03 6.52028e-02
9 1.78897e-03 6.33447e-02
10 1.36930e-03 6.22637e-02
11 1.05671e-03 6.16543e-02
12 8.21161e-04 6.13315e-02
13 6.41884e-04 6.11797e-02
14 5.04238e-04 6.11293e-02
15 3.97855e-04 6.11379e-02
16 3.15079e-04 6.11756e-02
17 2.50290e-04 6.12273e-02
18 1.99396e-04 6.12843e-02
19 1.59204e-04 6.13391e-02
20 1.27371e-04 6.13904e-02
21 1.02074e-04 6.14360e-02
22 8.18963e-05 6.14757e-02
23 6.58049e-05 6.15106e-02
24 5.29191e-05 6.15398e-02
25 4.26226e-05 6.15643e-02
26 3.42929e-05 6.15839e-02
27 2.76582e-05 6.16019e-02
28 2.23147e-05 6.16157e-02
29 1.80136e-05 6.16268e-02
30 1.44952e-05 6.16356e-02
31 1.17148e-05 6.16445e-02
32 9.45362e-06 6.16508e-02
33 7.63823e-06 6.16559e-02
34 6.19796e-06 6.16599e-02
35 4.97943e-06 6.16625e-02
36 4.02984e-06 6.16658e-02
37 3.24430e-06 6.16681e-02
38 2.63320e-06 6.16703e-02
39 2.13223e-06 6.16714e-02
40 1.70096e-06 6.16723e-02
41 1.41018e-06 6.16739e-02
42 1.11664e-06 6.16739e-02
43 9.19287e-07 6.16750e-02
Computing atomic populations.
Computing atomic charges.
Computing density decomposition for atom 0
Computing density decomposition for atom 1
Computing density decomposition for atom 2
Computing cartesian and pure AIM multipoles and radial AIM moments.
Storing proatom density spline for atom 0.
Storing proatom density spline for atom 1.
Storing proatom density spline for atom 2.
charges:
[-0.829 0.415 0.415]
cartesian multipoles:
[[-0.829 -0. 0. 0.08 -5.229 -0. -0. -4.611 -0. -4.916 -0. 0. -0.116 0. -0. -0. 0. 0.229 0. -0.061]
[ 0.415 0. 0.028 -0.003 -0.351 -0. -0. -0.313 0.001 -0.322 0. 0.024 -0.005 0. 0. 0. 0.066 0.004 0.025 0.001]
[ 0.415 0. -0.028 -0.003 -0.351 0. 0. -0.313 -0.001 -0.322 0. -0.024 -0.005 0. -0. 0. -0.066 0.004 -0.025 0.001]]
Difference between total electrons:
0.018706938820722385
The result of the least squares method differs slightly from the quadratic programming method due to the constraints. We can observe that the deviation in the total electrons obtained by the least squares method is larger than that obtained by the quadratic programming method.