Hirshfeld atom refinement based on projector augmented wave densities with periodic boundary conditions

Deriving the atomic form factors from Hirshfeld-partitioned periodic projector augmented wave calculations shows a great benefit for H-atom bond lengths and a smaller benefit for H-atom atomic displacement parameters.

In a previous application of PAW-DFT by Wall (Wall, 2016) the author postulated, that the multiplication by a phase factor "is not appropriate for translations by fractional grid points on a fixed rectilinear grid such as is used here." Instead, he proposed a different factor that can be found in the same paper in equation 3. His derivation was demonstrated with the one-dimensional case, culminating in the one-dimensional equation 15: In this equation = Δ = + 1 , with 0 being the integer component and 1 the fractional remainder in the range 0 ≤ 1 ≤ 1, h is our reciprocal index and N is the number of lattice points in our one-dimensional partition.
We however found the application of this work to lead to worse results in application. We want to give two additional reasons, why we finally opted to use a traditional phase factor with FFT in our study.

S1.1. The factor for Δx and -Δx do not offset each other for small Δx
Let us now assume, we have only a small Δ ≤ 1/ . In this case u0 is zero and u1 is NΔx and the equation simplifies to: However, if we enter − Δ into the equation with | − Δ | ≤ 1/ , u0 is -1 and 1 = 1 − Δ and this then can be simplified to: The product Δ (ℎ) ⋅ −Δ (ℎ) should be one so that the operation is reversible, but this is not the case. For N=10, Δ = 0.05 and h = 1 the product becomes: GPAW does offer the possibility of using different descriptions for the electron density. The description as a linear combination of atomic orbitals does offer the possibility to compare the expansion of the atomic densities on the rectangular grid with an expansion onto a spherical grid. We can therefore directly compare the difference in results from the two expansions.
The expansion was written in Python using the grid implementation of HORTON (Verstraelen et al., 2015) with the integration of HORTON used for discrete Fourier Transform on the spherical grid.
The dataset we used for comparison was the L-Alanine at 23K. As the LCAO calculation in GPAW is only available for calculations on the GGA level, the RPBE functional was used for comparison.
Grid-spacing for wavefunction calculations was 0.12 Å. Expansion on the rectangular grid was done with four-fold interpolation resulting in a grid of 0.03 Å. Atomic grids of the "insane" preset of HORTON were used for spherical expansion.

Figure S1
Differences in wR2(F 2 ) and agreement to neutron values for different expansions on grids for the calculation of atomic form factors. Distribution of parameters are displayed as box-whisker plots. All calculations were done using the RPBE functional.
The resulting quality indicators of the different refinement can be found in Figure S1. If we compare the two LCAO refinements, the wR2(F 2 ) is barely lower in the rectangular grid Fourier transform, the agreement to neutron data increases for the spherical grid Fourier transform.
However, using a finite difference real space description of the density shows a significantly larger benefit. This means from a practical point of view, even if there was a more beneficial alternative to the phase, its effect should be small. If we would use densities described as a linear combination of atomic orbitals with a spherical grid Fourier transform, we would discard the more beneficial density description, as well as the possibility to calculate densities with meta-GGA functionals.
As a consequence of all the afore mentioned points, this work uses a rectangular grid description as well as the traditional phase expression.

S2. Validation of Refinement against SHELXL
We used a completely new implementation within PYTHON for the refinement. In order to be sure that the implementation is actually correct, we compared IAM calculations of the 23 K L-Alanine dataset with the Python script to IAM calculations in SHELXL (Sheldrick, 2015). The resulting atomic parameters can be found in Tables S2 and S3. For all atoms the resulting parameters and estimated standard deviations are identical.

S3.1. PAW-HAR in GPAW
If not noted differently all calculations used the finite difference (FD) mode in GPAW. The convergence criterion for the density was tightened to 10 -7 . If applied, Monkhorst-Pack grids (Monkhorst & Pack, 1976) were shifted to include the Γ point. If the Monkhorst-Pack grid is listed as Γ, the calculation only included a reciprocal space calculation at the origin of the reciprocal grid. In deviation to the default, full symmetry was considered. The Hirshfeld Atom Refinement cycle was employed until the difference of positions was under 10 -6 Å. The Hirshfeld partitioning and atomic form factor calculation was calculated on a twice interpolated grid from the given wavefunction grid.
This means, that the grid spacing for the FFT was one fourth of the given value. Additional parameters for the individual calculations are listed in Table S3.

S3.2. Non-periodical refinement in NoSpherA2
For A23K, HMa-8HQ and Xy the fragment was the asymmetric unit itself. For Urea, the molecule was completed by applying the mirror plane that lies in the C,O axis but not in the complete molecule.
For the HMa-Mg dataset the inversion centre on the Magnesium ion was applied to complete the hexaqua magnesium unit. In order to keep the overall inversion symmetry, the operation was applied to the maleate ion as well.

Figure S2
Depiction of the employed hexaqua magnesium maleate fragment, at the example of the Tonto calculation. The marked Mg1 atom is located on the inversion centre that was applied to the atoms of the asymmetric unit.

S3.2.1. Tonto cluster charge calculations
All structures were refined with the B3LYP functional and the def2-TZVPP basis set. Integration accuracy was set to high. Clusters were always completed. Structures were pre-refined using 4 Å of 7 cluster charges and a DIIS convergence criterion of 10 -2 . The final refinement for all structures except A23K for was done with 8 Å of cluster charges and completed clusters. The DIIS convergence criterion was 10 -4 . Due to convergence problems A23K used 10 -5 as the convergence criterion, however full convergence could not be reached within 20 cycles. Additionally, as the overall performance of 4 Å of cluster charges was superior to the larger cluster, the smaller radius was used for the final comparison to other A23K refinements.

S3.2.2. Orca isolated molecule calculations
Structures were refined with the noted functional and the def2-TZVPP basis set. Integration accuracy was set to high. The SCF Threshold was set to "VeryTightSCF" the SCF Strategy was set to NormalConv.

S5. Quality indicators for Hirshfeld Atom Refinement structures
In addition to the quality indicators listed in the paper, the following indicators were evaluated for the evaluation of the agreement of the hydrogen atom description:

wRMSD(Δr):
This quality indicator is the weighted root mean square difference between the X-H bond distances derived from neutron diffraction and X-ray diffraction. The weighting employed is the combined estimated standard deviations of the X-ray and neutron derived values leading to the formula:

wRMSD(ΔUij):
This quality indicator is the weighted root mean square difference between the atomic displacement parameters derived from neutron diffraction and X-ray diffraction. The weighting employed is combined estimated standard deviations of the X-ray and neutron derived values leading to the formula: ‹ΔV/Vn›: To determine whether the deviations between the anisotropic descriptions are random or systematic, we determined the difference in between calculated volumes of the probability function from neutron and Hirshfeld Atom Refinement divided by the respective volume of the neutron probability function. If X and n are the Cartesian displacement matrices from X-ray and neutron refinement respectively, the relative difference in volume can be calculated as: Additionally aggregated DRK(Adam Stash, 2007) and Henn-Meindl-plots / residual density plots (Meindl & Henn, 2008) have been produced for all structures and functionals investigated.

Table S5
Aggregated quality indicators the agreement in X-ray intensities and hydrogen atom description from for all density sources for the non-periodical calculations of the L-Alanine at 23 K

Figure S7
DRK plots of all non-periodic calculations in this work for the L-Alanine at 23 K (A23K) dataset.

Figure S8
Henn-Meindl plots of all non-periodic calculations in this work for the L-Alanine at 23 K (A23K) dataset.              Figure S25 DRK plots of all periodic PAW calculations in this work for the xylitol at 122 K (Xy) dataset.

Figure S26
Henn-Meindl plots of all periodic PAW calculations in this work for the xylitol at 122 K (Xy) dataset.  Figure S27 DRK plots of all non-periodic calculations in this work for the xylitol at 122 K (Xy) dataset.

Figure S28
Henn-Meindl plots of all non-periodic calculations in this work for the xylitol at 122 K (Xy) dataset.     Figure S32 DRK plots of all periodic PAW calculations in this work for the Urea at 123 K dataset.

Figure S33
Henn-Meindl plots of all periodic PAW calculations in this work for the Urea at 123 K dataset.  Figure S34 DRK plots of all non-periodic calculations in this work for the Urea at 123 K dataset.

Figure S35
Henn-Meindl plots of all non-periodic calculations in this work for the Urea at 123 K dataset.

S6. Gram-Charlier refinement of urea.
In Section 3.3.4 of the publication, we postulated that the improvement in the density description of Urea reveals the difference electron density pattern indicating anharmonic vibration, which we then claimed to have successfully described with refinement of Gram-Charlier parameter refinement up to the fourth order. Here we want to describe our approach to check the correct refinement within our script and to then validate the Gram-Charlier description and finally want to give the resulting quality indicators.

S6.1. Validation against Olex Refine
Olex (Dolomanov et al., 2009) provides the possibility to read in pre-written .tsc files (Midgley et al.,20.11.2019) for use in olex.refine. At the same time, we can also refine Gram-Charlier parameters up to the fourth order. In order to validate the performance of our own script we first did a refinement until convergence. From the final positions we then calculated atomic form factors with GPAW and the SCAN functional used for refinement in our script and wrote a .tsc file. This was then used for one refinement in Olex. The resulting Gram-Charlier parameters can be found in Table S15 for the Third order Gram-Charlier parameter and in Table S16 and Table S17 for the fourth order. As third order Gram-Charlier refinement of carbon did not yield values, that were significantly different from zero, no Gram-Charlier parameters for carbon were refined in the calculation.

S6.2. Significance, positive probability density and Kuhs' rule
Additionally, we can see that for both refined atoms there are third and fourth order Gram-Charlier parameters with a value of more than three estimated standard deviations. Therefore, the first criterion for the Gram-Charlier refinement is met.
For checking Kuhs' rule (Kuhs, 1992)and excluding the possibility of negative probability density from the Gram-Charlier refinement, we wrote a xd parameter file and subsequently used XDPDF from the XD2016 Suite (Volkov et al., 2016) for validation. All probability density is positive. With a given sin( ) value of 1.405 Å -1 we are below the required value for the refinement of fourth order Gram-Charlier parameters for the carbon atom (1.43 Å -1 ) and well above the values for the fourth order Gram-Charlier refinement of nitrogen (1.12 Å -1 ) and oxygen (1.28 Å -1 ). Subsequently, we tried refining third-order Gram-Charlier parameters for carbon and third and fourth-order Gram-Charlier parameters for nitrogen and oxygen. Both oxygen and nitrogen have anharmonic vibration parameters of third and fourth order, which are significant (|c|/σ(c) > 3). However, no cijk-parameter of the carbon atom was significant. Therefore Gram-Charlier correction was only refined for the nitrogen and oxygen atoms.