The complete and original calculation scheme beyond the finite difference approximation, for the atomic (and orbital) Fukui function (FF) indices is proposed. The method explores an expansion for derivatives of LCAO coefficients, ∂C/∂N=CU. The separation scheme for the U matrix has been elaborated at the ab initio level. Nucleophilic and electrophilic FF indices, as well as atomic softness, have been derived from the standard result of SCF HF ab initio calculations. The indices reproduce two effects; the change in orbital occupancy and the relaxation of the electronic system. The molecular hardness (softness) provided by this scheme explicitly includes these two effects.

In this paper, we verify the usefulness of the alchemical derivatives in the prediction of chemical properties. We concentrate on the stability of the transmutation products, where the term "transmutation" means the change of the nuclear charge at an atomic site at constant number of electrons. As illustrative transmutations showing the potential of the method in exploring chemical space, we present some examples of increasing complexity starting with the deprotonation, continuing with the transmutation of the nitrogen molecule, and ending with the substitution of isoelectronic B-N units for C-C units and N units for C-H units in carbocyclic systems. The basis set influence on the qualitative and quantitative accuracies of the alchemical predictions was investigated. The alchemical deprotonation energy (from the second order Taylor expansion) correlates well with the vertical deprotonation energy and can be used as a preliminary indicator for the experimental deprotonation energy. The results of calculations for the BN derivatives of benzene and pyrene show that this method has great potential for efficient and accurate scanning of chemical space.

We present an analytical approach to treat higher order derivatives of Hartree-Fock (HF) and Kohn-Sham (KS) density functional theory energy in the Born-Oppenheimer approximation with respect to the nuclear charge distribution (so-called alchemical derivatives). Modified coupled perturbed self-consistent field theory is used to calculate molecular systems response to the applied perturbation. Working equations for the second and the third derivatives of HF/KS energy are derived. Similarly, analytical forms of the first and second derivatives of orbital energies are reported. The second derivative of Kohn-Sham energy and up to the third derivative of Hartree-Fock energy with respect to the nuclear charge distribution were calculated. Some issues of practical calculations, in particular the dependence of the basis set and Becke weighting functions on the perturbation, are considered. For selected series of isoelectronic molecules values of available alchemical derivatives were computed and Taylor series expansion was used to predict energies of the "surrounding" molecules. Predicted values of energies are in unexpectedly good agreement with the ones computed using HF/KS methods. Presented method allows one to predict orbital energies with the error less than 1% or even smaller for valence orbitals.

With the idea of using alchemical derivatives to explore in an efficient, computer- and cost-effective way Chemical Space was launched several years ago. In the context of Conceptual DFT response functions, these energies vs nuclear charge derivatives permit the estimatation of the energy of transmutants of a given starting or reference molecule showing different nuclear compositions. After an explorative study on small and planar molecules ( Balawender et al. J. Chem. Theory Comput. 2013 , 9 , 5327 ) by the present authors of this paper, the present study fully exploits the computational advantages of the alchemical derivatives in larger three-dimensional systems. Starting from a single reference calculation on C, the complete BN substitution pattern, from single substituted CBN via the belt (C(BN) and the ball C(BN) structures to the fully substituted (BN), is explored. Successive and simultaneous substitution strategies are followed and compared, indicating that both techniques yield identical results up to 13 substitutions but that for higher substitutions the simultaneous approach needs to be taken. Due to the cost-efficiency of the algorithm this path can indeed be followed as opposed to earlier work in the literature where for each step a full SCF calculation was at stake leading to prohibitively large computational demands for adopting the simultaneous approach. Previously formulated rules governing the substitution pattern by Kar and co-workers are scrutinized in this context and reformulated giving chemical insight in the gradual substitution process and the relative energies of the isomers. In its present form the method offers an interesting venue to study BN substitution patterns in higher fullerenes and graphene and in general paves the way for more efficient exploration of the Chemical Space.

A calculation scheme of the nuclear Fukui function via a coupled perturbed Hartree–Fock approach is proposed avoiding the finite difference approach in DFT-based descriptors. Nucleophilic and electrophilic nuclear Fukui functions are compared with the numerical approximation for the nuclear Fukui function (FF) as the negative derivative of the chemical potential with respect to the atomic coordinates and as the derivative of the Helman–Feynman force with respect to the total number of electrons. The results for a set of diatomic molecules are shown. Analytical and numerical techniques do show a high correlation. Overall, values from both numerical methods are larger than those from the analytical one. The analytical results can be interpreted in terms of the character of the orbital involved during ionization or adding of electrons; the change in the equilibrium bond length upon ionization, which is positive for bonding orbitals and negative for antibinding orbitals is connected with the negative or positive values of the left-hand-side nuclear FF, respectively. The nucleophilic nuclear FF is positive for all cases except CO indicating a systematical increase of the bond length after addition of an electron.

The properties of the derivative of the total binding function (the virial of the forces) with respect to the number of electrons and its decomposition at local and atomic level have been analyzed. At local level the binding function is expressed by the Berlin function fv(r) and the electronic Fukui function f(r). The atomic analog is expressed in terms of the nuclear Fukui function (FF) and the nuclear position vectors. A relationship between the local maps of fv(r)f(r), the nuclear FF vectors, and the Jahn–Teller distortion direction is discussed. It is predicted that upon ionization the symmetry of the nearest local stationary point for BH3 is C2v, for AH4 molecules (CH4 and SiH4) D2d, and for C3H6 C2v. For the benzene anion a D2h symmetry is predicted.

The effective fragment potential (EFP) model has been used to study the effect of adding increasing numbers of the water molecules on several DFT-based reactivity descriptors of NH 3 . The HOMO-LUMO gap and electrophilic hardness are seen to increase with addition of water molecules. The importance on the wave function relaxation in the solvent effect on ammonia's properties is shown when analyzing the relaxation part in the electrophilic hardness and condensed Fukui function for the nitrogen atom. An increase in the atomic softness for the nitrogen atom with decreasing the global softness is observed. The saturation point for solvatation of ammonia was located around a cluster with 16 molecules of water. Atomic properties such as the Mulliken population, condensed Fukui function, and atomic softness for nitrogen and electrophilic global properties such as the hardness and its components for dilute solutions are predicted faithfully.

The effect of solvent on the electronegativity, hardness, and condensed Fukui function, and atomic softness for a set of diatomic and small polyatomic molecules and ions has been studied using the effective fragment potential (EFP) model. The binding function was used for monitoring the solvation of the molecule. We do not observe a decrease in the HOMO-LUMO gap in the solvent. All anions show a significant change in the chemical potential. Both HOMO and LUMO energy levels decrease in the solvent phase as compared to the gas phase. For the major part of the acids, the increase in the LUMO orbital energy is larger than in the HOMO orbital energy. For the group of salts, we observe an increase in the LUMO energy level and a similar decrease in the HOMO energy level, resulting in a small change in the chemical potential. The importance of the change in the wave function upon solvation was shown through an analysis of the relaxation part in the hardness and condensed Fukui function. Very close values found for the same ions in molecules such as LiH, LiF, NaH, NaF, and LiF indicate that in these cases very good separated ion pairs are present.