How is Poier and Jensen's Bond Capacity Model Parameterized and Optimized?

Question

I am working on adding a charge polarization model into my own research and have been exploring a few approaches. One of the most attractive options is the Bond Capacity (BC herein) model of Paolo and Jensen (DOI: 10.1021/acs.jctc.8b01215). Alas, there is no example implementation available and though it has been added to Tinker-HP (DOI: 10.1021/acs.jctc.9b00721) the version that might contain it has yet to see a release.

I have been prototyping my own implementation, but I have ran into some roadblocks when trying to reproduce results in the paper. The primary issue is a lack of information regarding the model's parameterization and handling of integrals. I understand the questions are narrow, but I'm taking a shot in the dark and looking for any valuable insights.

I will briefly highlight the purpose and key equations to save you from reading the paper. I should note, the BC model was designed to describe molecular polarization and polarizabilities using fluctuating charges.

Now for the equations. I am going to presume we are working with a neutral system of $N$ atoms acting under the influence of an external electric field $F$. In the BC model, the electrostatic potential for atom $i$, $V_i$, is the sum of the intrinsic potential of atomic electronegativity ($\chi_i$), the cumulative external potential of neighboring atoms ($\phi_i$), and the potential due to the external field ($\psi_i$). The amount of charge that is transferred between a pair of atoms $p_{ij}$ is taken to be proportional to a pair capacity $\xi_{ij}$ and the difference in electrostatic potential for the atoms: $$p_{ij} = \xi_{ij}(V_i - V_j) \\ p_{ij} = \xi_{ij}\left[(\chi_i+\phi_i+\psi_i)-(\chi_j+\phi_j+\psi_j)\right]$$

This pair capacity parameter is geometry dependent and is defined as: $$\xi_{ij}(R_{ij}) = \xi_{ij}^{0} g(R_{ij}) = \xi_{ij}^{0} \left(\frac{1-S_{ij}}{R_{ij}}\right)$$ Where $\xi_{ij}^{0}$ is pair capacity at equilibrium geometry and $g(R_{ij})$ is an attenuation function that is one at equilibrium and decays to zero as inter atomic distance grows. This is where my first two questions come in. In the paper, they choose the overlap between Slater orbitals ($S_{ij}$) as part of the attenuation function.

1. By chance, does anyone have any idea if they are using STOs with angular momenta greater than $l=0$ to evaluate the overlap of atoms?

2. In similar models, only one s Slater type orbital is used per atom. I presume the same in my implementation and evaluate the overlap using equations (5) and (6) from (DOI:10.1002/qua.560210612). Though, this form is not analytically derived. Does anyone have a reference that gives an exact form for the two-centered overlap of S-type slater orbitals?

Okay moving forward, the total charge on any atom $Q_i$ due to the aforementioned electrostatic potential components can be determined via, $$Q_i = \sum_j^N -\xi_{ij}(V_i - V_j)$$ thus the charge on every atom in the system can summarized in the matrix form, $$\textbf{Q}=-\textbf{CV}$$

Where $\textbf{C}$ is symmetric ($\xi_{ij}=\xi_{ji}$) and contains negative pair capacities of atoms on the off-diagonal and the diagonal is a sum of these pair capacities in the corresponding column or row, $\xi_{ii} = \sum_{k\neq j}^N \xi_{ik}$. $\textbf{V}$ is a vector of the atomic electrostatic potentials defined above.

The components of $\textbf{V}$ should be defined briefly. To determine the polarization potential of atom $i$ due to its neighbors, $\phi_i$, you must evaluate a Coulomb interaction, ($J_{ik}$), with a similar form to $\xi_{ij}(R_{ij})$: $$\phi_{i}=\sum_{k=1}^{N} J_{ik} Q_{k}=\sum_{k \neq i}^{N} \frac{1-S_{i k}}{R_{ik}} Q_{k}+J_{ii} Q_{i} \\ \boldsymbol{\phi = JQ}$$ Where again the overlap between Slater orbitals ($S_{ik}$) attenuation function is applied. Meanwhile, the intrinsic potential of atomic electronegativity, $\boldsymbol{\chi}$, is taken to be a vector of free parameters. The potential due to the external electric field has it's usual definition $$\psi_i = \sum_{\sigma=x,y,z}R_{i \sigma}F_{\sigma} \\ \boldsymbol{\psi = R^TF}$$ where $\textbf{R}$ is a (N,3) matrix of atomic positions and $\textbf{F}$ is the electric field vector.

Combining the above equations and performing some rearrangements allows one to define the 3x3 molecular polarizability tensor as, $$\boldsymbol{\alpha = R^T\left( I + CJ \right)^{-1}CR}$$ where all $\textbf{I}$ is the identity matrix and all other terms have been defined above. This brings me to my final question:

3. The authors are using an, "in house," code to perform the fitting of their parameters to reproduce ab initio polarizabilities. I have no way of knowing what routines they are using. I am trying to determine what STO exponential constants $\zeta$ and electronegativities $\chi_i$ were used to attain their results. I am currently using scipy's least_squares to optimize these parameters via minimizing the residuals between the components of my polarizability matrix and theirs (iteratively of course). I am not a statistician nor a fitting expert. Since this is just a system of equations, does anyone know of a better way of fitting the parameters to ab initio data in python or C++?

I apologize for the long winded post, but these details have been vexing for a while now and I'm getting tired of perusing Google for answers. If you made it this far, thanks for reading.

Answer 1

The specific questions:

1. Only $l = 0$ overlaps are (currently) considered, in the spirit of atoms being represented by (only) a spherical charge distribution (point charges).

2. The overlap/attenuation function for the bond capacities is very (!) preliminary. It should be designed to reflect the actual physics, and this is work in progress. However, any reasonable function that models how the bonding decays as a function of distance should at least give a fair model. The same holds for the Coulomb screening function. In the paper, it was taken to be the same as the screening function for the bond capacity, but one could (and should) considered these to have different functional forms. We opted for the simple solution of using the overlap calculated using the s-functions of the STO-3G basis set. This involves only overlap of s-type Gaussian functions which is trivial to calculate.

3. We used a Newton-Raphson optimizer with numerical gradients, a piece of code that had been developed for another purpose, and could easily be modified to do the job. We assigned the bond capacities from calculated polarizabilities, and then fitted the electronegativity parameters to give the desired atomic charges. One also need to assign the atomic Coulomb self-interaction (aka hardness) parameters to give the desired polarization response.

As inferred from the above, the parameterization of the BC model is currently very rudimentary, and there is certainly room for refinements, but it does have some fundamental features that look attractive.