Evidence for the existence of Li₂S₂ clusters in lithium–sulfur batteries: ab initio Raman spectroscopy simulation

Abstract

Using density functional theory calculations and ab initio molecular dynamics simulations we have studied the structures and the Raman spectra of Li₂S₄ clusters, which are believed to be the last polysulfide intermediates before the formation of Li₂S₂/Li₂S during the discharge process in Li–S batteries. Raman spectra have been obtained using a new technique to estimate polarizabilities using Wannier functions. We have observed clear evidence of Li₂S₄ → Li₂S₂ transition by studying systematic changes in the simulated Raman spectra of (Li₂S₄)_n, n = 1, 4, and 8 towards that of (Li₂S₂)₈. Furthermore, we have shown that the dominant Raman peak of the Li₂S₂ cluster at ∼440 cm⁻¹ arises from sulfur–sulfur stretching mode. This peak has been experimentally observed in the discharged state of Li–S batteries and has also been attributed to the formation of Li₂S₂. We have also demonstrated that the transition is mainly due to the strong electrostatic interactions between Li₂S₄ monomers, which results in energy lowering by arranging the local Li^+δ–S^−δ dipole moments in an anti-parallel fashion.

---

1 Introduction

Lithium–sulfur (Li–S) batteries are promising energy storage systems with a high theoretical energy density of 2500 W h kg⁻¹ (or 2800 W h L⁻¹).¹ During the discharge process in Li–S batteries, lithium ions react with sulfur on the cathode, leading to the formation of soluble Li–polysulfides, for example Li₂S₈, Li₂S₆, and Li₂S₄ clusters, followed by the formation of insoluble Li₂S₂, and eventually Li₂S crystals.^2,3 There have been many experimental studies, such as in situ Raman,^4–8 X-ray diffraction,^4,6,9 scanning electron,^4,6 transmission electron,^6,9 electron paramagnetic resonance,¹⁰ and ultraviolet-visible spectroscopies,¹¹ all aimed at understanding the processes occurring during discharge/charge cycles, and investigating the cathode structure. Li₂S₄ clusters are expected to be the last intermediates in the S₈ → Li₂S transition.^10,12,13 Nevertheless, it is not fully understood yet whether any other product, for example Li₂S₂, exists besides Li₂S crystals in the discharged state.^14–18 Li₂S₂, however, has not been clearly identified in experimental measurements on Li–S batteries during discharge.^19,20 It is partially due to the fact that it is experimentally difficult to isolate intermediate Li–polysulfide clusters and perform the measurements on them separately. However, theoretical simulations on structures, energetics, and vibrational properties enable one to study the structural transitions during the discharge process with molecular resolution, and can provide a complementary picture to experimental results. There are important studies based on density functional theory (DFT) calculations as well as many-body techniques aiming at understanding the structures of the Li–polysulfides and the mechanism of their transition to Li₂S₂/Li₂S (for the most recent studies, see ref. 17, 18 and 21–23). There have also been promising advances in simulation of the vibrational spectra of systems with general symmetries.^24–32 In particular, ab initio simulation of the Raman spectra of candidate structures combined with experimental Raman data can unambiguously identify the most probable atomic configurations. It has been shown recently that the electronic properties of a system introduce distinct features in the vibrational spectra.^29,32 Therefore, the electronic properties of the system need to be explicitly taken into account, and as such, classical simulations, usually based on interatomic parametric potentials, may not be the most suitable choice. Moreover, in order to have a realistic simulation, directly comparable to experimental measurements, finite-temperature effects need to be taken into account as well, a feature which is absent in perturbative methods at zero temperature.^33,34 All these requirements can be met by ab initio molecular dynamics (AIMD) simulations, in which forces are calculated “on-the-fly” using DFT.^26,35 To the best of our knowledge, ab initio simulation of Raman spectroscopy on individual Li–polysulfide clusters has not been performed so far.

In this article, using DFT calculations and AIMD simulations we first determine the atomic structures of (Li₂S₄)₁, (Li₂S₄)₄, and (Li₂S₄)₈ polysulfides. Furthermore, using a new method based on the Wannier function technique along with the second-generation Car–Parrinello method of Kühne et al.,^36,37 we obtain the polarizabilities of the structures as a function of time and thereby the Raman spectra of the clusters. We observe a clear tendency to form the Li₂S₂ structure as the size of the cluster grows. The crystallization is found to be due to the strong electrostatic interaction between Li₂S₄ monomers, which results in anti-parallel arrangement of the local Li^+δ–S^−δ dipole moments.

2 Computational details

The minimum-energy structures of (Li₂S₄)_n clusters have been obtained at the DFT level of theory using the mixed Gaussian and plane-wave code CP2K/QUICKSTEP^38,39 in conjunction with a very accurate TZV2PX Gaussian basis set,⁴⁰ Goedecker–Teter–Hutter pseudo potentials⁴¹ and the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functional,^42,43 plus semi-empirical correction for the long-range dispersion interactions (DFT-D3).⁴⁴ A real-space grid was represented by a plane-wave energy cutoff of 270 Ry, and the convergence criterion for the self-consistent field was set to 10⁻⁶. The clusters were modeled using supercells with unit cells of 16 × 16 × 16 ų, 21 × 21 × 21 ų, and 37 × 37 × 37 ų for (Li₂S₄)₁, (Li₂S₄)₄, and (Li₂S₄)₄, respectively, and therefore, the Brillouin zone was sampled only using the Γ-point. The minimum-energy structures were obtained using a BFGS optimizer.^45–49 Convergence criteria for the maximum geometry change, and the maximum force component were set to 1.6 × 10⁻³ Å, and 0.02 eV Å⁻¹, respectively. Vibrational normal mode analysis⁵⁰ was carried out with an atomic displacement of 5 × 10⁻³ Å.

The intensity of the isotropic Raman scattering, σ, as a function of the frequency of the incident beam, ν, is related to the mean polarizability, Ā = 1/3Tr[Â] through^31,51

where 〈⋯〉_cl denotes the statistical average in classical mechanics. Furthermore, by applying a periodic electric field using the Berry phase approach^52,53 one can obtain the polarizability tensor,from the total dipole moment of the system M. E_j denotes Cartesian component j of the applied periodic electric field, and M_i is the ith component of the total dipole moment. M_i can be estimated unambiguously using Maximally Localized Wannier Functions (MLWFs), which allow the total electronic density to be partitioned into individual fragments.^54–57 However, calculating the derivatives numerically based on density functional perturbation theory (DFPT)^33,34 or using higher-order finite difference (FD) methods^58,59 is computationally rather expensive. Alternatively, inspired by the fact that molecular polarizability changes linearly with the volume of a molecule,⁶⁰ one can assume that the total mean polarizability of the system can be expressed as a sum over the polarizabilities assigned to each Wannier function in the system,where β is a proportionality constant and S_i is the spread of the ith Wannier function. In a separate study (see ref. 61 for details) we have shown that in the new approach, which we hereafter refer to as the Wannier polarizability (WP) method, the best value for β is 0.30. In the FD approach, however, during an AIMD run with N_MD steps after equilibrium, one needs to perform 6N_MD extra single-point calculations to obtain the polarizability. The WP method, which has been used to simulate the Raman spectra in the present work, facilitates to routinely calculate finite-temperature spectra with only minimal extra computational cost. We have demonstrated that for typical test systems the simulated Raman spectra obtained using FD and WP methods are essentially the same in terms of relative intensity of the peaks and their lineshapes.⁶¹

All Raman spectra in this work were obtained by performing 20 ps AIMD simulations in the canonical ensemble to achieve equilibration, followed by 10 ps AIMD in the micro-canonical ensemble to remove thermostat effects and to sample polarizabilities. A time-step of 1 fs was used in all the simulations and the polarizabilities were sampled for every 5 fs. All AIMD simulations were also carried out using the CP2K/QUICKSTEP software package with the same simulation setup as mentioned before.

3 Results and discussion

We first start by calculating the isotropic Raman spectrum of the S₈ molecule in the gas phase. The spectrum is shown in Fig. 1. The S₈ molecule is known to have eleven fundamental vibrational modes, among which six peaks at 86, 152, 218, 248, 437, and 475 cm⁻¹ are Raman active. The two dominant peaks are observed at 218 and 475 cm⁻¹.⁶² Our simulated spectrum shows six peaks at around 63, 145, 218, 232, 427, and 470 cm⁻¹, two of which, namely 218 cm⁻¹ and 470 cm⁻¹ being the dominant ones. The Raman spectrum presented here using the WP method is in agreement with previous experimental^63–65 and theoretical¹⁷ studies, showing a mean absolute error of ∼10 cm⁻¹ with respect to the experimental data reported in ref. 62. It is worth mentioning here that the choice of basis set is crucial, especially for the case of low-frequency modes below 150 cm⁻¹. These frequencies are absent when the DZVP basis set is employed, while they emerge when more diffuse basis functions, such as TZV2PX, are used. Vibrational normal mode analysis at zero temperature reveals that the most dominant peak is due to doubly degenerate asymmetric bending mode of the three consecutive sulfur atoms in the S₈ molecule. However, the stretching modes occur in the frequency range of 350 to 470 cm⁻¹ with much lower intensity. Similar findings have also been reported experimentally.^63–66

Fig. 1 Atomic structure and isotropic Raman spectrum of S₈ molecules in the gas phase. Insets are zoomed views of weak Raman activities.

In order to find the structural transition of Li₂S₄ to Li₂S₂/Li₂S, we have examined a variety of possible arrangements of Li₂S₄ clusters with four and eight monomers at the DFT level. The charge population analysis on Li and S atoms of (Li₂S₄)_n structures using Mulliken population analysis (MPA)⁶⁷ reveals that the minimum-energy structures are those with an anti-parallel arrangement of local dipole moments of Li^+δ–S^−δ in Li₂S₄ molecules. The electrons are partially transferred from Li to S atoms, leaving the −δ charge on S and the +δ charge on Li, leading to the lowest electrostatic repulsion and highest electrostatic attraction between Li₂S₄ monomers (see Fig. 2). For a single Li₂S₄ molecule, the electron charges are transferred from Li₁ and Li₂ atoms to all S atoms (see Fig. 2). However, the major part (∼85%) of electron charges is transferred to the nearest S1 and S2 atoms and only small part to the next nearest atoms, namely S3 and S4. The averaged electron transfer (using MPA) is calculated to be 0.49 |e|. The calculated dipole moment for the isolated Li₂S₄ molecule is 2.63 Debye, which is comparable to that of a water molecule in bulk liquid water.⁶⁸ However, the residual total dipole moments of 0.11 Debye and 0.31 Debye for (Li₂S₄)₄ and (Li₂S₄)₈ clusters, respectively, show a large cancellation of dipole moments which would lead to a zero net polarization for (Li₂S₄)_n as n → ∞. These results are also in line with the energy gain of the clusters as a function of the number of Li₂S₄ monomers. The formation energy of (Li₂S₄)_n as a function of the cluster size n can be written as

where E^cluster_tot, , and E^Li_tot are the total energy of (Li₂S₄)_n clusters, a single S₈ molecule, and a Li atom, respectively. n_S and n_Li are the number of S and Li atoms in (Li₂S₄)_n clusters, respectively. As shown in Fig. 3, we find that the formation energy increases with n. The energy gain with increasing size of (Li₂S₄)_n clusters is comparable to the Coulomb interaction energy between Li^+δ and S^−δ point charges of interacting Li₂S₄ molecules, namely , while . As can be seen in Fig. 3, the formation energy may still increase with the size of clusters, which is due to the fact that the number of dangling bonds around the (Li₂S₄)₈ cluster is still large compared to the saturated S bonds. We therefore expect that the formation energy of (Li₂S₄)_n, n → ∞ is even lower than that for (Li₂S₄)₈. This clearly indicates a tendency towards forming an iconic crystal.

Fig. 2 Structures of (Li₂S₄)_n clusters. S and Li atoms are in yellow and purple, respectively. Initial and final configurations refer to structures before and after geometry optimization, respectively.

Fig. 3 Formation energy as a function of the cluster size, n. The structures are illustrated in Fig. 2.

The simulated Raman spectra are presented in Fig. 4. In comparison with the Raman spectrum of the S₈ molecule (Fig. 1), we find three new large peaks at around 330, 380, and 450 cm⁻¹, for the Li₂S₄ monomer (Fig. 4(a)). As the size of the (Li₂S₄)_n increases from n = 1 (Fig. 4(a)) to n = 8 (Fig. 4(c)), the peak at ∼200 cm⁻¹ loses its relative intensity with respect to the new peaks at 330, 380, and 450 cm⁻¹. As stated before, a S₈ molecule has a dominant vibrational mode at ∼200 cm⁻¹ which comes from the asymmetric bending of the three consecutive sulfur atoms in the ring. Therefore, the peak at around 200 cm⁻¹ in Fig. 4 should also arise from the collective vibrations of sulfur atoms in Li₂S_x clusters with x ≥ 3. Note that the stretching vibrations of the S₈ ring have frequencies above 350 cm⁻¹ with rather weak intensities. Moreover, in the case of (Li₂S₄)₂ (Fig. 4(b)), the peak at 330 cm⁻¹ becomes less notable, while the peaks at 380 and 450 cm⁻¹ turn more significant. When the (Li₂S₄)₈ cluster forms, we find a dominant peak at ∼440 cm⁻¹ and a smaller peak at ∼380 cm⁻¹ (Fig. 4(c)). Recently, these two peaks have been experimentally attributed to the formation of crystalline Li₂S₂/Li₂S structures: it has been shown by Raman spectroscopy studies that during the discharge process, when a Li–S battery is discharged to 1.7 V, two dominant peaks emerge for the cathode at ∼380 cm⁻¹ and ∼440 cm⁻¹.^5,16 The former peak is very close to the experimentally observed characteristic peak of the Li₂S crystal at 375 cm⁻¹.⁶⁹ Using normal mode analysis for a Li₂S₄ monomer at zero temperature we also observe Raman-active modes mainly in the frequency range of 350–500 cm⁻¹. This shows that the Raman activity between 350 and 500 cm⁻¹ in Fig. 4 mostly arises from lithium–sulfur interactions. Additionally, the decline in the relative intensity of the ∼200 cm⁻¹ peak with respect to the ones at 380 and 440 cm⁻¹ indicates a clear transition from Li₂S_x clusters with x ≥ 3 to structures with x < 3. The possible structures would be Li₂S₂ or Li₂S. We find that the peak at 440 cm⁻¹ is due to sulfur–sulfur stretching vibration (see insets in Fig. 4(a)–(c)). This can be directly seen by replacing Ā in eqn (1) with A_i, where i is the index of a Wannier function centered along a typical S–S bond in (Li₂S₄)_n structures. Here we have set i to be the index of the Wannier function along the S1–S3 (or S2–S4) bond shown in Fig. 2.

Fig. 4 Simulated Raman spectra of (a) (Li₂S₄)₁, (b) (Li₂S₄)₄, and (c) (Li₂S₄)₈. The structures are illustrated in Fig. 2. The insets show the contribution of S–S stretching vibration to the total Raman spectra.

We also simulated the Raman spectrum of the (Li₂S₂)₈ cluster which was in this case obtained after 60 ps of AIMD simulation in the canonical ensemble at 300 K, followed by a 20 ps production run in the micro-canonical ensemble. The Raman spectrum is shown in Fig. 5(a). The peak at ∼440 cm⁻¹ is also present in the (Li₂S₂)₈ cluster which, in turn, implies the existence of the S–S covalent bond. This peak has been experimentally attributed to Li₂S₂.^5,16Fig. 5(b) shows the S–S stretching contribution to the total spectrum along with the corresponding Wannier function used to calculate the spectrum. Additionally, by studying the S–S bond length during our long AIMD simulation (60 ps), we do not observe any S–S bond breaking which shows that (Li₂S₂)₈ is stable. Furthermore, there is no Raman activity observed at around 200 cm⁻¹ in Fig. 5(a), which is in line with our previous assignment of this vibration to more than two covalently bonded sulfur atoms. However, in all the insets of Fig. 4 there is a weak activity slightly above 200 cm⁻¹ which can be due to the coupling of collective sulfur vibrations with S–S stretching. Based on these observations, we therefore ascribe the peak at ∼440 cm⁻¹ in Fig. 4(b) and (c) to the formation of the Li₂S₂ structure in (Li₂S₄)₄ and (Li₂S₄)₈ clusters.

Fig. 5 Simulated Raman spectra of (a) (Li₂S₂)₈ clusters considering all Wannier functions and (b) partial Raman spectrum arising from the Wannier function centered along the S–S covalent bond. A typical covalent bond connecting two sulfur atoms and the Wannier function centers in the (Li₂S₂)₈ cluster are also depicted in yellow and blue, respectively.

4 Conclusions

In summary, we have performed DFT calculations and AIMD simulations to investigate the structure of Li₂S₄ and Li₂S₂ clusters. We have also simulated the Raman spectra of the most favorable structures, based on the Wannier function technique to study Li₂S₄ → Li₂S₂ structural transition in Li–S batteries during the discharge process. Our findings show that Li₂S₄ monomers bind strongly to each other due to the electrostatic interactions between Li^+δ–S^−δ dipoles of the interacting monomers and form larger clusters. In particular, we observe a dominant Raman peak at ∼440 cm⁻¹ in (Li₂S₄)_n clusters which we then assign to S–S stretching mode. The other peak at ∼200 cm⁻¹, which corresponds to the vibration of more than two covalently bonded sulfur atoms, loses its intensity as the cluster size grows. Simultaneously, the peak at 440 cm⁻¹ becomes more notable with the cluster size. Particularly, in the case of the (Li₂S₂)₈ cluster, we find a dominant peak at ∼440 cm⁻¹ with nearly no Raman activity at around 200 cm⁻¹. Therefore, we have been able to identify the structural transition to Li₂S₂ by following the trend in the Raman spectra. Moreover, we unambiguously ascribe the experimentally observed peak at around 440 cm⁻¹ to the existence of the Li₂S₂ structure. This work has implications in understanding the structure of the sulfur cathode in Li–S batteries. The simulation method presented in this work also enables one to routinely simulate the Raman spectrum “on-the-fly” based on ab initio molecular dynamics simulations, with minimal extra computational cost. Therefore, the method can be used as an efficient approach to study the structure of complex systems relevant to energy storage.

Acknowledgements

P. K. acknowledges support from the “Bundesministerium für Bildung und Forschung” (BMBF).

References

1. Y. Yang, G. Zheng and Y. Cui, Chem. Soc. Rev., 2013, 42, 3018–3032.
2. L.-X. Miao, W.-K. Wang, A.-B. Wang, K.-G. Yuan and Y.-S. Yang, J. Mater. Chem. A, 2013, 1, 11659–11664.
3. G. Zhou, D.-W. Wang, F. Li, P.-X. Hou, L. Yin, C. Liu, G. Q. M. Lu, I. R. Gentle and H.-M. Cheng, Energy Environ. Sci., 2012, 5, 8901–8906.
4. Z. W. Seh, H. Wang, P.-C. Hsu, Q. Zhang, W. Li, G. Zheng, H. Yao and Y. Cui, Energy Environ. Sci., 2014, 7, 672–676.
5. H. Yao, K. Yan, W. Li, G. Zheng, D. Kong, Z. W. Seh, V. K. Narasimhan, Z. Liang and Y. Cui, Energy Environ. Sci., 2014, 7, 3381–3390.
6. Z. W. Seh, J. H. Yu, W. Li, P.-C. Hsu, H. Wang, Y. Sun, H. Yao, Q. Zhang and Y. Cui, Nat. Commun., 2014, 5, 1–8.
7. M. A. Lowe, J. Gao and H. D. Abruña, RSC Adv., 2014, 4, 18347–18353.
8. K. H. Wujcik, J. Velasco-Velez, C. H. Wu, T. Pascal, A. A. Teran, M. A. Marcus, J. Cabana, J. Guo, D. Prendergast and M. Salmeron, et al. , J. Electrochem. Soc., 2014, 161, A1100–A1106.
9. Y. Yang, M. T. McDowell, A. Jackson, J. J. Cha, S. S. Hong and Y. Cui, Nano Lett., 2010, 10, 1486–1491.
10. Q. Wang, J. Zheng, E. Walter, H. Pan, D. Lv, P. Zuo, H. Chen, Z. D. Deng, B. Y. Liaw and X. Yu, et al. , J. Electrochem. Soc., 2015, 162, A474–A478.
11. C. Barchasz, F. Molton, C. Duboc, J.-C. Leprêtre, S. Patoux and F. Alloin, Anal. Chem., 2012, 84, 3973–3980.
12. M. Cuisinier, P.-E. Cabelguen, S. Evers, G. He, M. Kolbeck, A. Garsuch, T. Bolin, M. Balasubramanian and L. F. Nazar, J. Phys. Chem. Lett., 2013, 4, 3227–3232.
13. J. Song, S. J. Lee, Y. Kim, S.-S. Kim, K. T. Lee and N.-S. Choi, ECS Electrochem. Lett., 2014, 3, A26–A29.
14. X. Ji and L. F. Nazar, J. Mater. Chem., 2010, 20, 9821–9826.
15. A. Manthiram, Y. Fu and Y.-S. Su, Acc. Chem. Res., 2012, 46, 1125–1134.
16. J.-T. Yeon, J.-Y. Jang, J.-G. Han, J. Cho, K. T. Lee and N.-S. Choi, J. Electrochem. Soc., 2012, 159, A1308–A1314.
17. L. Wang, T. Zhang, S. Yang, F. Cheng, J. Liang and J. Chen, J. Energy Chem., 2013, 22, 72–77.
18. G. Yang, S. Shi, J. Yang and Y. Ma, J. Mater. Chem. A, 2015, 3, 8865–8869.
19. J. Nelson, S. Misra, Y. Yang, A. Jackson, Y. Liu, H. Wang, H. Dai, J. C. Andrews, Y. Cui and M. F. Toney, J. Am. Chem. Soc., 2012, 134, 6337–6343.
20. S. Waluś, C. Barchasz, J.-F. Colin, J.-F. Martin, E. Elkam, J.-C. Leprêtre and F. Alloin, Chem. Commun., 2013, 49, 7899–7901.
21. H. Park, H. S. Koh and D. J. Siegel, J. Phys. Chem. C, 2015, 119, 4675–4683.
22. B. Wang, S. M. Alhassan and S. T. Pantelides, Phys. Rev. Appl., 2014, 2, 034004.
23. R. S. Assary, L. A. Curtiss and J. S. Moore, J. Phys. Chem. C, 2014, 118, 11545–11558.
24. R. Ramírez, P. Kumar and D. Marx, J. Chem. Phys., 2004, 121, 3973.
25. R. Iftimie and M. E. Tuckerman, J. Chem. Phys., 2005, 122, 214508.
26. R. Iftimie, P. Minary and M. E. Tuckerman, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6654–6659.
27. G. Cicero, J. C. Grossman, E. Schwegler, F. Gygi and G. Galli, J. Am. Chem. Soc., 2008, 130, 1871–1878.
28. A. Morita and T. Ishiyama, Phys. Chem. Chem. Phys., 2008, 10, 5801–5816.
29. D. Donadio, G. Cicero, E. Schwegler, M. Sharma and G. Galli, J. Phys. Chem. B, 2009, 113, 4170–4175.
30. C. Zhang, D. Donadio and G. Galli, J. Phys. Chem. Lett., 2010, 1, 1398–1402.
31. T. Ishiyama, V. V. Sokolov and A. Morita, J. Chem. Phys., 2011, 134, 024510.
32. Q. Wan, L. Spanu, G. A. Galli and F. Gygi, J. Chem. Theory Comput., 2013, 9, 4124–4130.
33. P. Giannozzi and S. Baroni, J. Chem. Phys., 1994, 100, 8537–8539.
34. S. Baroni, S. De Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515.
35. D. Marx and J. Hutter, Ab initio molecular dynamics: basic theory and advanced methods, Cambridge University Press, 2009.
36. T. D. Kühne, M. Krack, F. R. Mohamed and M. Parrinello, Phys. Rev. Lett., 2007, 98, 066401.
37. T. D. Kühne, WIREs Comput. Mol. Sci., 2014, 4, 391.
38. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–488.
39. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
40. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
41. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703.
42. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
43. M. Krack, Theor. Chem. Acc., 2005, 114, 145–152.
44. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
45. C. G. Broyden, IMA J. Appl. Math., 1970, 6, 76–90.
46. C. G. Broyden, IMA J. Appl. Math., 1970, 6, 222–231.
47. R. Fletcher, Comput. J., 1970, 13, 317–322.
48. D. Goldfarb, Math. Comput., 1970, 24, 23–26.
49. D. F. Shanno, Math. Comput., 1970, 24, 647–656.
50. R. M. Martin, Electronic structure: basic theory and practical methods, Cambridge University Press, 2004.
51. P. Atkins and J. De Paula, Atkins' physical chemistry, Oxford University Press, 2010.
52. I. Souza, J. Íñiguez and D. Vanderbilt, Phys. Rev. Lett., 2002, 89, 117602.
53. P. Umari and A. Pasquarello, Phys. Rev. Lett., 2002, 89, 157602.
54. N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza and D. Vanderbilt, Rev. Mod. Phys., 2012, 84, 1419.
55. R. Resta, Phys. Rev. Lett., 1998, 80, 1800.
56. P. L. Silvestrelli, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 9703–9706.
57. G. Berghold, C. J. Mundy, A. H. Romero, J. Hutter and M. Parrinello, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 10040.
58. I. Vasiliev, S. Öğüt and J. R. Chelikowsky, Phys. Rev. Lett., 1997, 78, 4805.
59. A. Putrino and M. Parrinello, Phys. Rev. Lett., 2002, 88, 176401.
60. J. S. Murray, P. Lane, T. Brinck, K. Paulsen, M. E. Grice and P. Politzer, J. Phys. Chem., 1993, 97, 9369–9373.
61. P. Partovi-Azar, T. D. Kühne, 2015, arXiv preprint arXiv:1504.03639.
62. B. Meyer, Chem. Rev., 1976, 76, 367–388.
63. R. Steudel and H.-J. Mäusle, Angew. Chem., Int. Ed. Engl., 1977, 16, 112–113.
64. A. Talyzin, L.-E. Tergenius and U. Jansson, J. Cryst. Growth, 2000, 213, 63–69.
65. S. N. White, Chem. Geol., 2009, 259, 240–252.
66. P. D. Harvey and I. S. Butler, J. Raman Spectrosc., 1986, 17, 329–334.
67. R. S. Mulliken, J. Chem. Phys., 1955, 23, 1833–1840.
68. T. D. Kühne, M. Krack and M. Parrinello, J. Chem. Theory Comput., 2009, 5, 235–241.
69. Z. Lin, Z. Liu, N. J. Dudney and C. Liang, ACS Nano, 2013, 7, 2829–2833.

---