Molecular Sciences Modeling Nmr Chemical Shifts: Crystal Potential Derived Point Charge (cppch) Model to Calculate Solid State Effects on 31 P Chemical Shifts Tensors

This paper presents a new method to calculate solid-state effects on NMR chemical shifts. Using full crystal potentials, this new method (CPPCh) eliminates the need to arbitrarily select the point charges that are included in the calculations of the NMR chemical shieldings to take into account intermolecular effects. By eliminating the arbitrary selection of the point charges, the method provides a mechanism to systematically improve the simulation of intermolecular effects on chemical shielding calculations. This new method has been applied to the calculation of the 31 P NMR chemical shifts tensors in P 4 S 3. The shielding calculations were done using the DFT approach with the BLYP gradient corrected exchange correlation functional. This method was selected to calculate the 31 P chemical shifts because it includes electron correlation effects at a reasonable cost.


Introduction
Solid state nuclear magnetic resonance is nowadays an important tool to elucidate the structure of organic, biological macromolecules and inorganic materials.This is a consequence of the increasing progress in the field of high-resolution solid state NMR spectroscopy (SSNMR).Due to the popularization of techniques such as cross-polarization magic angle spinning (CP/MAS), SSNMR is routinely available in many structural analysis laboratories [1].With the increasing use of SSNMR data to solve structural problems in the solid state there is a greater need to enhance its usefulness by combining the experimental NMR data with molecular modeling of the NMR chemical shifts [1].Recently, molecular modeling and simulation techniques have achieved notable improvements in the description of intermolecular interactions and bulk phenomena.For instance, the solid-state effects observed in both, 13 C and 15 N chemical shifts have been accurately described by several modeling methods [2,6].
Quantum mechanical methods commonly employed to compute NMR chemical shifts do not include crystalline effects in the calculations.Recently Mauri et al. [7,8] have introduced a new method to calculate NMR chemical shifts in infinite periodic systems, which of course includes all the intermolecular interactions.Unfortunately, this method has been implemented only using the LDA approximation [9,11], which reproduces chemical shifts with much less accuracy than other DFT (density functional theory) methods that use hybrid exchange correlation functionals [12].Therefore, when these functionals are used in the calculation of chemical shieldings the intermolecular effects must be taken into account by explicit modeling techniques.
Currently, there are two general schemes to model intermolecular effects: i) the cluster models and ii) the charge point models.A detailed discussion of both schemes is available in the introductory section of Ref. [13] .The cluster model explicitly includes in the calculation of the chemical shielding tensors the neighboring molecules or some of their significant fragments.The difficulties with this method arise from the selection of the significant fragments, which is somehow arbitrary, and the considerable increase in the computational resources needed when large fragments are added to the calculations.Charge models, which describe the intermolecular effects by a discrete distribution of point charges, are very attractive because they derive from a very simple idea: the crystalline electrostatic potential can be represented by a finite distribution of point charges, bond dipoles, and/or distributed multipoles.The great advantage of these models is that the cost of shielding calculations does not significantly increase when a finite charge distribution is added to the calculations.Unfortunately, current charge models suffer from one of the shortcomings of the cluster model: the selection of charges to be included in the calculations is also arbitrary.
This paper shows that using a distribution of point charges fitted to the electrostatic potential given by the full crystalline wavefunction, it is possible to eliminate the problems associated with the arbitrary selection of charges representing the crystalline environment.Moreover, the description of the potential can be systematically improved by increasing the number of charges included in the calculations.The crystalline electrostatic potential can be calculated using the full crystal wave-function like one available from the CRYSTAL98 package [14].
Phosphorous compounds are a constituent of many important systems such as minerals, glasses, zeolites, bio-organic materials, etc. [15].The high sensitivity of 31 P, which has natural abundance of 100 %, makes solid-state 31 P NMR spectroscopy a valuable tool to gather information about the structure, dynamics, and properties of these materials.Consequently, it is appropriate to choose a phosphorous compound to test the new charge point model reported here.
This contribution presents the results obtained using the new CPPCh method to calculate the 31 P NMR chemical shifts in tetra phosphorous trisulfide, P 4 S 3 .The results are compared with those obtained using the GRID point charge method [16] and with experimental values.

Crystal potential derived point charge (CPPCh) model
The CPPCh model calculates the intermolecular effects on the chemical shifts by performing the calculations on a molecule (the central molecule) surrounded by a distribution of point charges that represents the crystalline potential generated by the rest of the crystal lattice.This model differs from other methods of modeling intermolecular interactions by using a distribution of point charges chosen to reproduce the electrostatic potential of the full crystal.Other methods, used to calculate chemical shieldings, select the distribution of point charges to reproduce the potential generated by the electronic distribution of individual molecules [14].
The crystalline potential can be evaluated employing the CRYSTAL98 code [15], using the following methodology.Considering a region of space, C, in which the charge density is zero and in which the electrostatic potential is produced by a charge distribution, ρ(r), lying entirely outside C; the Green's theorem [17,18] states that the electrostatic potential inside C, V el (r), is independent of the charge distribution, for any distribution that generates the same potential on the boundary surface S of the volume C. Note that the electrostatic potentials generated outside C will be in general different for different charge distributions.This is not relevant to our problem because we are only interested in the influence of the crystal lattice on the nuclei for which the chemical shifts are calculated.All these nuclei are inside C due to the way in which C is constructed (see bellow).If we consider C as the cavity where the central molecule is located (see Fig. 1), after removing the central molecule the charge density in this region is zero, therefore the Laplace equation [18,19] is satisfied for the electrostatic potential V ¨9r)=0 (1) for r inside C; using Green's theorem it is possible to replace the electrostatic field generated by the infinite crystal lattice at the location of the central molecule by an equivalent field generated by a finite charge distribution outside C.Such finite charge distribution should generate the same potential over the boundary of C than the infinite crystal lattice.A Van der Waals molecular surface (of the central molecule) is a convenient way to create the boundary surface for the cavity C. In this work we employed the Gepol routine [19], improved by C.
Chipot, Gepol92 [20], to determine the surface S. To evaluate the lattice potential over this surface, S, we calculated both the crystalline electrostatic potential and the molecular electrostatic potential on points of the surface S.These points were selected using the Gepol92 routine.Both electrostatic potentials were calculated using the CRYSTAL98 program employing the same basis set and the same quantum mechanical method.The electrostatic potential that gives the boundary condition necessary to determine the finite distribution of charges is given by where V(r) crys the crystalline electrostatic potential corresponding to the infinite lattice and V(r) mol is the electrostatic potential generated by the electronic distribution of the central molecule.The desired distribution of charges can be obtained by minimazing the function f(q 1 , r 1 , q 2 , r 2 ,….. q nq , r nq ) in Eq. ( 3), with respect to the position and value of the charges.
f(q 1 , r 1 , q 2 , r 2 ,….. q nq , r nq ) If the position of the charges is defined a priori, see bellow, Eq. ( 3) becomes a linear least square problem on the charges.In this paper we have used two different approaches to set the position of the charges, these are identified as Distribution I and II.Distribution I corresponds to the situation in which the charges are located on a second surface, S 2 , enclosing S. The second surface, S 2 , is another Van der Waals surface which was built duplicating the atomic radii.This is a convenient choice because it avoids penetration problems.The number and position of point charges in S 2 were selected employing the Gepol92 routine and it can be increased as necessary, to improve the accuracy of the fit.The value of the charges was determined using Eq.(3).Distribution II was designed by assigning point charges to each atomic position belonging to first neighbor molecules outside of the region C.These charges are fitted to the crystalline electrostatic potential using Eq. ( 3) by a least squares procedure, but it is important to realize that they do not have any relation to the concept of atomic charges in the crystal.Second and higher order neighbors can be added to improve the fitting of the crystalline potential by this charge distribution.A similar approach was applied recently to fit the Madelung potentials generated by Ewald summations in ionic crystals [21].
The selected distribution of point charges was used thereafter to perform shielding calculations using the Gaussian98 [22] suite of programs.The shielding calculations were performed according to the DFT (Density Functional Theory) approach proposed by Cheesman et al. [13], using the BLYP exchange correlation functional [23,24] with a D95** basis set [25].This method makes use of an efficient implementation [26] of the gauge including atomic orbitals (GIAO) method [27,28] and it has been extensively used to calculate chemical shieldings [2].

Results and Discussion
The calculated principal components of the 31 P shielding in P 4 S 3 , are given in Tables 1 and 2. The results for the three basal 31 P nuclei are in Table 1, the corresponding results for the apical phosphorous nucleus are entered in Table 2.The first and second row in these tables correspond to the calculations for the isolated molecule employing the molecular geometries optimized with a 6-31G** basis set [29] and determined by neutron diffraction (ND) [30], respectively.The ND geometry was also employed in all the calculations for which a finite distribution of charges was included.The results obtained using the GRID model, with a configuration of 56 point charges, are given in the third row of the tables.Finally, the results obtained with the CPPCh model are presented in the fourth and fifth rows of the tables.The results in the fourth row correspond to calculations using the charges given by the Distribution I: 190 point charges distributed on the surface S 2 , which was described in the preceding section.The fifth row depicts the results corresponding to the calculations including the charges given by the Distribution II: 56 point charges located at the atomic positions of the 56 atoms (the same number employed in the calculations using the GRID model) corresponding to first neighbor molecules to the central one.The RMS between the electrostatic potential calculated using the ab initio method, V SCF , and those generated by the finite distributions of charges are 10 -3 au. for the GRID, 13×10 -3 au. for the CPPCh-I, and 4×10 -3 au. for the CPPCh-II models, respectively.As discussed above the number of charges included in the CPPCh models could be increased to lower the corresponding RMS, securing a better description of the crystalline electrostatic potential.Unfortunately, convergence problems were encountered when a large number of very small charges were included in the shielding calculations.The large difference between the optimized molecular geometry and the one corresponding to the ND determination, explains the large discrepancy (more than 100 ppm) observed between the chemical shift values calculated using these two geometries.The principal values of the chemical shift tensor of the apical phosphorous are better reproduced than those of the basal phosphorous.The calculated values show that the apical shift tensor has axial symmetry when the optimized geometry is used in the calculations.This is a consequence of the existence of a molecular C 3 axis normal to the base of the cage and in the direction of the apical phosphorous.The results calculated using the ND geometry, regardless of the inclusion of intermolecular effects, show a small difference between the δ 22 and δ 33 components.This difference, which is a consequence of the breaking of the C 3 symmetry by intermolecular interactions and packing, is much smaller than the one observed in the experimental data.All the basal phosphorous are equivalent in the experimental spectra [31,32].This is also the case when the optimized molecular geometry is used in the shielding calculations, but the basal phosphorous become non-equivalent when the ND geometry is used in the calculations.The difference among the chemical shift principal values calculated for different basal nuclei is less than 10 ppm.This value is of the order of the errors observed in the experimental measurements, which may explain the lack of resolution of these shifts in the experimental spectra.The calculated values for the basal phosphorous have been averaged for comparison with the experimental values.
The RMS between the experimental [32,33] and calculated 31 P chemical shifts are given in Table 3.The values given in the first column correspond the comparison of the calculated values with the experimental data from Ref. [32].The RMS values calculated with respect to the experimental data of Gibby et al. [33], which exhibit a similar behavior, are given in the second column of the table.
The results obtained with the different charge models show very similar trends and represent a significant improvement over the conventional calculations for an isolated molecule using the optimized geometry.The differences among the different models are quite modest compared with the large improvements obtained with respect to the bare calculation using the optimized geometry.For the compound studied here the results of the CPPCh method do not show a significant improvement over those obtained using the GRID charge method, this may be due to the large experimental errors in the 31 P chemical shifts of P 4 S 3 .There are formal advantages in using the CPPCh method over the previous point charge methods These include its ability to be extended to ionic systems for which the previous methods do not work [33] and its capability to systematically improve the description of the crystalline potential by increasing the number of charges in the distribution.
The CPPCh method must be checked not only in other phosphorous compounds but also in other nuclei.The case of P 4 S 3 has been employed as an example to present the model, but a more comprehensive study is underway

Figure 1 .
Figure 1.a) crystal fragment, b) the cavity in the fragment represents the C region.

Table 1 .
31P NMR chemical shifts for the basal nuclei of P 4 S 3 .Chemical shifts in ppm with respect to phosphoric acid.a

Table 2 .
31P NMR chemical shifts for the apical nucleus of P 4 S 3 .Chemical shifts in ppm with respect to phosphoric acid.a

Table 3 .
Comparison between the 31 P experimental and calculated chemical shifts principal values in P 4 S 3 .All values in ppm.