Investigating the first hyperpolarizability of liquid carbon tetrachloride

Abstract

Sequential QMMM calculations have been carried out to investigate the first hyperpolarizability (β​_HRS,Liq) of liquid CCl₄. First, Monte Carlo simulations are performed to generate statistically uncorrelated snapshots representing the liquid structure. Then, the first hyperpolarizability of selected snapshots are evaluated using ab initio calculations. In these calculations the solvent effects are described either exclusively by point charges or a few neighboring CCl₄ molecules are also explicitly considered. In particular, it has been observed that considering small numbers of CCl₄ molecules, embedded in point charges, enables monitoring the emergence of the dipolar contribution to β_HRS,Liq and the increase of the depolarization ratio, confirming experimental results and substantiating that the dipolar contributions originate from intermolecular interactions between the CCl₄ molecules. Additional calculations using semi-empirical Hamiltonians for the QM part were performed on systems containing explicitly one or two solvation shells, and further confirmed the emergence of the dipolar contribution in liquid CCl₄.

---

1 Introduction

Carbon tetrachloride, CCl₄, has been used as a reference for hyper-Rayleigh scattering (HRS) first hyperpolarizability (β) measurements, both in gas and liquid phases.^1–6 In parallel to measurements, quantum chemical calculations have been performed on this benchmark compound, addressing various issues, including the importance of electron correlation to describe the wave-function and its responses, the vibrational contributions to these responses and their frequency dispersions, as well as the solvent effects.^6–9 For gas phase CCl₄ and an incident light wavelength λ = 1064 nm, a theoretical β_HRS value of 14.3 a.u. has been estimated,⁹ in good agreement with the experimental value of 12.8 ± 1.1 a.u. due to Shelton.⁵ The former was obtained by adding to the electronic dynamic coupled cluster singles and doubles – quadratic response function (CCSD-QRF) value (14.9 a.u.)⁹ the ZPVA (−1.1 a.u.) and pure vibrational (0.5 a.u.) contributions of ref. 7. Thus, the vibrational contributions to this second harmonic generation phenomenon are negligible.⁷ These calculations also demonstrated that the β_HRS ratio between the Hartree–Fock (HF) and CCSD with a perturbative estimate of the triples [CCSD(T)] methods amounts to 0.73.⁶ Moreover, frequency dispersion, evaluated by considering the β_HRS (λ = 1064 nm)/β_HRS (λ = ∞) ratio at the CCSD-QRF level, was shown to attain a value of 1.08.⁹

For liquid CCl₄, a β_Liq/β_Gas ratio of 1.43 has been estimated at the CCSD(T) level from calculations using the implicit polarizable continuum model (PCM).⁶ Experimentally, this ratio amounts to 1.57 ± 0.04,¹ evidencing an underestimation of β_HRS,Liq by the theoretical approximations. This comparison was however performed by considering the octupolar part (β_J=3) of β_HRS, i.e. by assuming that each CCl₄ molecule behaves as an independent light scatterer, giving rise to the so-called incoherent HRS signal. However, already 20 years ago, Kaatz and Shelton¹ have shown that the experimental HRS signal contains both a coherent (β_coh) and an incoherent (β_incoh) part, with β_coh/β_incoh ≈ 2/3. This coherent response originates from interactions between the CCl₄ molecules and attributes to liquid CCl₄ both dipolar (β_J=1) and octupolar (β_J=3) HRS responses. To our knowledge, the interplay between the β_J=1 and β_J=3 components of liquid CCl₄ has not yet been addressed from a theoretical chemistry viewpoint, which has motivated this investigation. On the basis of recent studies by the authors on the first hyperpolarizability of solutions,^10–12 a sequential quantum mechanics molecular mechanics (S-QMMM) method^13–17 was adopted to describe and analyze β_HRS of liquid CCl₄. In these recent works, the S-QMMM scheme has been employed to evaluate the effects of concentration on the β_HRS response of nitrobenzene in benzene solution^10,11 and to compare implicit and explicit solvation models for evaluating β_HRS of a biphotochromic nonlinear optical (NLO) switch.¹² The rest of this paper is organized as follows: the next section summarizes the key theoretical and computational aspects, the results are then presented and discussed in Section III, before conclusions and outlook are drawn in Section IV.

2 Computational procedures

The Monte Carlo (MC) simulations of liquid CCl₄ were performed following the standard procedure.¹⁸ A cubic box using periodic boundary conditions with the image method and the Metropolis sampling technique¹⁹ were applied. In the MC simulations the CCl₄ molecule was kept rigid at the geometry optimized at the B3LYP/6-311G(d) level in presence of solvent (CCl₄) described by the polarizable continuum model (PCM).^20,21 A total of 1000 CCl₄ molecules was considered in the NVT simulations carried out at 298 K. The volume was determined by the experimental density (1.584 g cm⁻³).²² The intermolecular interactions were described by the combination of Lennard-Jones (LJ) potentials and Coulomb terms. The OPLS²³ set of parameters (σ and ε) for the LJ potential were used as well as the atomic charges. The initial configuration, i.e. the atomic Cartesian coordinates for the molecules, was generated randomly. A new configuration (or snapshot) is generated after enabling all molecules in the box to move (translation and rotation around a chosen axis). Then, a first stage of 1.2 × 10⁸ MC steps was performed for thermalization. It was followed by an equilibrium simulation of 3.0 × 10⁸ MC steps. This part of the sequential QMMM method gives access to the configurations of the liquid and provides information about the structure and organization of CCl₄ molecules surrounding one CCl₄ molecule selected as reference. The radial distribution function g(r) was calculated to determine the first solvation shells as well as the number of molecules within each shell.

The solvent effects on the CCl₄ β_HRS response were evaluated using three different approaches. In the first, the reference CCl₄ molecule is embedded in an electrostatic potential represented by point charges (PC). In the second, it interacts explicitly with a few nearest-neighbor CCl₄ molecules and these explicit molecules are embedded in the PC electrostatic potential of the remaining CCl₄ molecules. In the third one, the full set of molecules constituting the first and second solvation shells around the reference molecule were considered explicitly, without PC embedding.

From the MC simulations described above a set of statistically uncorrelated snapshots representing the liquid structure are selected. The number of MC steps needed to uncorrelate these structures was determined by calculating the auto-correlation function of the energy.^13,14,17 The selected snapshots are used in the ab initio – or semi-emiprical – calculations to obtain the β responses.

When the solvent is represented by point charges only, we used the averaged solvent electrostatic configuration (ASEC) approach.²⁴ For the calculations where solvent molecules are considered explicitly, the results are presented under the form of simple averaged values with their standard deviations.

The static and dynamic first hyperpolarizabilities were evaluated at different levels of theory. Calculations using HF and DFT levels were performed with the aug-cc-pVQZ and d-aug-cc-pVTZ atomic basis sets. At these levels, the dynamic first hyperpolarizabilities were calculated using time-dependent schemes (TDHF and TDDFT).²⁵ The LC-BLYP and M05-2X exchange–correlation (XC) functionals have been adopted because they have been shown to be the most efficient, among a collection of about 20 XC functionals, to provide well-balanced dipolar and octupolar contributions to the first hyperpolarizability of reference molecules for nonlinear optics, including CCl₄.²⁶ The static first hyperpolarizability was also calculated by employing the high-level electron correlated CCSD(T) method within the finite field (FF) procedure.²⁷ The accuracy in the calculation of the numerical derivatives was improved and controlled using the Romberg procedure,^28,29 which allows to reduce the effects of higher-order contaminations in the calculation of the hyperpolarizabilities. Field amplitudes from ±0.0004 to ±0.0064 a.u. in a geometric progression ratio of 2 were used. To achieve high accuracy on the third-order energy derivatives, the convergence on the SCF (HF) energy was lowered to 10⁻¹¹ a.u. whereas to 10⁻⁹ a.u. in the CCSD(T) iterative procedure. In order to study the larger clusters, containing up to 61 explicit CCl₄ molecules (vide infra), semi-empirical Hamiltonians were employed in combination with the TDHF scheme. The most recent PM6 (ref. 30) and PM7 (ref. 31) parameterizations were selected.

The quantities of interest in this work are the first hyperpolarizability β_HRS and the depolarization ratio DR³²

where the decomposition of β_HRS in its two orientationally-averaged components, 〈β_ZZZ²〉 and 〈β_ZXX²〉, is evidenced. β_HRS is associated with hyper-Rayleigh scattering (HRS) experiments where a non-polarized/natural incident light beam propagates along the Y direction and the intensity of the vertically-polarized (along the Z axis) signal scattered at 90° (along the X axis) is detected. The β_HRS response can also be written in terms of its multipolar components, i.e., dipolar β_J=1 and octupolar β_J=3,

Detailed expressions for these quantities are presented in ref. 6. All reported β values are given in atomic units (1 a.u. of β = 3.6213 × 10⁻⁴² m⁴ V⁻¹ = 3.2064 × 10⁻⁵³ C³ m³ J⁻² = 8.639 × 10⁻³³ esu) and expressed within the T convention. The classical Monte Carlo simulations were performed with the DICE program,³³ the ab initio quantum chemistry calculations with Gaussian09 program,³⁴ the semi-empirical PM6 and PM7 calculations with MOPAC,³⁵ and the Romberg iterative scheme was carried out with a locally-developed program.²⁸

3 Results

3.1 Gas phase and electrostatic embedding

As a first step, Table 1 presents the static and dynamic (1064 nm) first hyperpolarizabilities of the isolated molecule, calculated using the LC-BLYP and M05-2X functionals in comparison with reference values obtained at the Hartree–Fock level as well as at coupled clusters levels of approximation, i.e. using CCSD quadratic response function (CCSD-QRF)⁹ and FF/CCSD(T).⁶ In agreement with ref. 26, the LC-BLYP functional appears to perform slightly better than M05-2X with an overestimation of the static CCSD(T) β_HRS value by 12% and close agreement with CCSD-QRF for the frequency dispersion factor, F(ω) = β (λ = 1064 nm)/β (λ = ∞).

Table 1 Static (λ = ∞) and dynamic (λ = 1064 nm) first hyperpolarizability β_HRS (a.u.) of the isolated CCl₄ molecule as determined at different levels of approximation. The last column gives F(ω), the dynamic/static ratios. At all levels the DR is equal to 1.50

	βHRS	βHRS	F(ω)
	λ = ∞	λ = 1064 nm
LC-BLYP/d-aug-cc-pVTZ	10.71	11.41	1.065
M05-2X/d-aug-cc-pVTZ	10.81	12.25	1.133
HF/d-aug-cc-pVTZ^6	6.96	7.03	1.010
HF/d-aug-cc-pVDZ^9	6.33	6.32	0.998
CCSD(T)/d-aug-cc-pVTZ^6	9.52	—	—
CCSD-QRF/d-aug-cc-pVDZ^9	13.84	14.89	1.076

---

Then, we present the β_HRS and DR of solvated CCl₄ as described by the electrostatic embedding scheme. The discrete structure of the liquid is analyzed using the radial distribution function g(r) (Fig. 1). Spherical integration of g(r) provides the average number of molecules within each solvation shell. In particular, there is an average of 12 CCl₄ molecules in the first shell, in agreement with experiment^36,37 (neutron and X-ray diffraction data) as well as with other theoretical simulations.^38–41 The second shell is composed of 48 molecules.

Fig. 1 Radial distribution function g(r) between the centers of mass of the CCl₄ molecules.

The selected snapshots are composed of 500 molecules forming a sphere of 26 Å of radius around the reference molecule. The QM values of the NLO properties are obtained by averaging over 100 overlapping statistically-uncorrelated snapshots (ASEC scheme). Table 2 reports the β_HRS, β_J=1, and β_J=3 values of the PC-embedded CCl₄ molecule evaluated at the DFT level using the LC-BLYP and M05-2X functionals. Six different integration grids were employed to sample the convergence of the β values as a function of the grid mesh. They highlight the independence of the first hyperpolarizability as a function of the integration grid for LC-BLYP while a small dependence is noticeable with the M05-2X functional when the (75,302) and (96,32,64) grids are used. On the basis of these results, a (99,590) grid was adopted for all DFT calculations.

Table 2 Static first hyperpolarizability β_HRS (a.u.) and β_J (a.u.) of one CCl₄ molecule embedded in the electrostatic field of 499 molecules described by point charges. The values were obtained using DFT with the LC-BLYP and M05-2X functionals and the d-aug-cc-pVTZ basis set for different grid meshes

Int grid	βHRS	βLiq/βGas	βJ=1	βJ=3
LC-BLYP
(75,302)	10.89	1.016	0.67	35.28
(99,590)	10.88	1.016	0.67	35.24
(99,974)	10.88	1.016	0.67	35.24
(120,974)	10.88	1.016	0.67	35.24
(150,974)	10.88	1.016	0.67	35.24
(96,32,64)	10.88	1.016	0.68	35.24
M05-2X
(75,302)	10.81	0.985	0.81	35.00
(99,590)	10.66	0.986	0.84	34.52
(99,974)	10.66	0.985	0.84	34.51
(120,974)	10.66	0.986	0.84	34.50
(150,974)	10.66	0.985	0.84	35.51
(96,32,64)	10.76	1.016	0.80	34.85

---

The HRS β_Liq/β_Gas ratios, also given in Table 2, demonstrate that the solvent effects as described by PC are very small, both on β_HRS and DR. All the DR values are equal to 1.50, which is consistent with the octupolar character of CCl₄. Nevertheless, owing to the non perfect isotropy of the point charge surrounding, the dipolar β_J=1 component is not exactly zero like for an isolated molecule. Still, the β_J=1 component is negligible with respect to the octupolar component β_J=3.

Table 3 lists the HRS quantities determined at the DFT level with these two functionals, as well as at the HF and CCSD(T) levels. The HF approximation underestimates β_HRS by 34% and 33% in comparison to DFT calculations using the LC-BLYP and M05-2X functionals, respectively. The reference CCSD(T) provides a value of 9.91 a.u. for β_HRS, about 7–9% smaller than DFT but larger than the HF value by 40%. On the other hand, the β_HRS,Liq/β_HRS,Gas ratios are close to one, whatever the employed method. This strong underestimation of the solvent effects in comparison to the experimental value¹ indicates that a simple electrostatic embbeding does not properly describe the impact of intermolecular interactions on the NLO properties of CCl₄.

Table 3 S-QMMM static and dynamic (1064 nm) first hyperpolarizability (a.u.) and depolarization ratio of a single solvated CCl₄ molecule as determined at different levels of approximation with the d-aug-cc-pVTZ basis set. The solvent is modeled by the PC of 499 CCl₄ molecules within the ASEC scheme. The β_HRS,Liq/β_HRS,Gas ratios are given in parentheses in column 3

	λ (nm)	βHRS	βJ=1	βJ=3	DR
HF	∞	7.19 (1.018)	0.57	23.27	1.50
	1064	7.29 (1.021)	0.63	23.62	1.50
LC-BLYP	∞	10.88 (1.016)	0.67	35.24	1.50
	1064	11.60 (1.017)	0.75	37.57	1.50
M05-2X	∞	10.66 (0.986)	0.84	34.52	1.50
	1064	12.29 (1.003)	0.95	39.76	1.50
CCSD(T)	∞	9.91 (1.041)	0.50	32.11	1.50

---

3.2 Explicit quantum embedding

The second approach to describe the solvent effects considers explicitly a few CCl₄ molecules interacting with the reference CCl₄. The calculations of the β responses were performed on the same structures generated in the MC simulations, but the QM region contains now two or five CCl₄ molecules, i.e., the reference molecule plus one or four nearest-neighbors within the first solvation shell. All interact with the electrostatic embedding of the remaining molecules. Owing to the substantial computational requirements (100 snapshots), the calculations were performed only at the HF and LC-BLYP levels. Table 4 presents the results for a pair of CCl₄ molecules embedded in point charges whereas Fig. 2 highlights the convergence of β_HRS and DR as a function of the number of snapshots. As quantified by the β_HRS (two molecules)/[2 × β_HRS (one molecule)] ratios, considering the explicit solvation has a little impact on the first hyperpolarizability per CCl₄ molecule. On the other hand, the β_J=1 dipolar contribution is no more negligible so that the depolarization ratio gets substantially larger than 1.50.

Table 4 S-QMMM dynamic (λ = 1064 nm) HRS quantities for a pair of solvated CCl₄ molecules as calculated at the Hartree–Fock and DFT (LC-BLYP) levels with the d-aug-cc-pVTZ basis set. Two molecules are considered explicitly in the QM calculations and are embedded in an electrostatic field of point charges. The quantities in parentheses are the ratios between the β_HRS of the CCl₄ molecule pairs and twice the value for a single PC-embedded CCl₄ molecule. All values are averages over 100 QM calculations and the uncertainties are the standard deviations

	HF	LC-BLYP
〈βHRS〉	14.93 ± 6.89 (1.024)	23.27 ± 10.45 (1.003)
〈βJ=1〉	16.99 ± 13.41	25.54 ± 20.30
〈βJ=3〉	39.28 ± 14.23	62.18 ± 21.99
〈DR〉	2.24 ± 0.61	2.17 ± 0.57

---

Fig. 2 Statistical convergence of β_HRS (top) and depolarization ratio (bottom) calculated at the HF/d-aug-cc-pVTZ level as a function of the number of configurations for a pair of CCl₄ molecules solvated by point charges embedding. The horizontal line indicates the averaged value reported in Table 4. The uncertainty corresponds to the statistical error.

Then, five CCl₄ molecules were considered in the QM region, again by considering 100 uncorrelated snapshots. However, this could only be performed using smaller basis sets, aug-cc-pVTZ and aug-cc-pVQZ. Indeed, doubly-augmented basis sets contain so many diffuse functions for 5 CCl₄ molecules that the SCF procedures has much difficulties to converge. Note that with the aug-cc-pVTZ and aug-cc-pVQZ basis sets, each of the CCl₄ molecule is represented by 246 and 416 atomic orbitals, respectively. The results given in Table 5 for one, two, and five molecules in the QM regions confirm the small variation of the β_HRS per molecule (without considering the aug-cc-pVTZ result for a single CCl₄ molecule) with respect to the size of the QM region, demonstrating saturation with respect to the number of molecules. Then, between 2 and 5 explicit molecules, DR is leveling off. DR is also smaller with aug-cc-pVQZ than with aug-cc-pVTZ. The percentage of dipolar contribution, estimated using the following expression

demonstrates that the dipolar or coherent response is substantial. Note that although the aug-cc-pVTZ basis set is not suitable to describe the first hyperpolarizability of a single CCl₄ molecule, its adequacy with respect to aug-cc-pVQZ is substantially improved when considering two or five CCl₄ molecules.

Table 5 S-QMMM dynamic (λ = 1064 nm) HRS quantities for one, two, and five solvated CCl₄ molecules as calculated at the Hartree–Fock level with the aug-cc-pVTZ (TZ) and aug-cc-pVQZ (QZ) basis sets. These molecules are considered explicitly in the QM calculations and are embedded in an electrostatic field of point charges. Quantities in parentheses are the values of β_HRS per molecule. See the caption of Table 4 for more details

	1 + PC		2 + PC		5 + PC
	TZ	QZ	TZ	QZ	TZ	QZ
〈βHRS〉	1.34	5.44	11.78 ± 6.84 (5.89)	14.01 ± 6.90 (7.00)	26.40 ± 9.19 (5.28)	28.38 ± 9.24 (5.68)
〈βJ=1〉	0.50	0.58	16.49 ± 13.26	16.93 ± 13.51	37.62 ± 18.30	37.91 ± 18.57
〈βJ=3〉	4.28	17.61	27.34 ± 12.46	35.79 ± 13.61	60.65 ± 21.01	68.69 ± 21.81
〈DR〉	1.56	1.50	2.73 ± 0.87	2.34 ± 0.66	3.03 ± 1.02	2.78 ± 0.92
% (dipolar)	3.05	0.24	43.5	32.5	45.1	39.7

---

Since ab initio calculations on larger CCl₄ clusters are out of the reach of actual computational resources, the last set of QM calculations were performed using PM6 and PM7 semi-empirical Hamiltonians. This allowed considering explicitly all molecules constituting the first (12 molecules), as well as the first and second (60 molecules) solvation shells. However, owing to the poor performance of these parameterizations to reproduce the isolated molecule β responses, the analysis is restricted to qualitative aspects, focusing on the description of the emergence of the dipolar contribution in the liquid phase. The results are given in Table 6 and were obtained after averaging over the same 100 uncorrelated snapshots without PC embedding. Consistent with ab initio approaches, semi-empirical calculations also predict the appearance of a nonzero dipolar character of the NLO response in the liquid phase of CCl₄. This contribution is however much smaller than that predicted using the higher levels of approximation described above, due to the inherent limitations of semi-empirical schemes. In particular, the use of a minimal Slater-type basis set is not sufficient to describe the full electronic response of a molecule subject to an external electric field. Still, it is noteworthy that these computationally cheap methods remain able to partly catch the emergence of the dipolar component induced by intermolecular interactions.

Table 6 S-QMMM static HRS quantities for CCl₄ with one and two solvation shells as calculated at the semi-empirical PM6 and PM7 levels. These 13 (one shell) or 61 (two shells) molecules are considered explicitly in the QM calculations. β values are given per CCl₄ molecule. All values are averages over 100 QM calculations

	1		1 + 1 solvation shell		1 + 2 solvation shells
	PM6	PM7	PM6	PM7	PM6	PM7
〈βHRS〉	32.85	43.65	9.36	8.56	4.72	4.17
〈βJ=1〉	0.02	0.02	6.19	4.00	3.26	2.09
〈βJ=3〉	106.46	141.44	28.34	26.92	14.19	13.05
〈DR〉	1.50	1.50	1.81	1.64	1.85	1.66
% (dipolar)	0.0	0.0	9.7	4.9	10.6	5.6

---

4 Conclusions and outlook

Sequential QMMM calculations have been carried out to investigate the first hyperpolarizability of liquid CCl₄. The method consists first in performing Monte Carlo simulations to generate uncorrelated snapshots to represent the liquid structure and then in calculating at the QM level the first hyper-polarizability for a selection of these snapshots. In these QM calculations, the solvent (surrounding molecules) is either exclusively described by point charges or by considering explicitly a few neighboring CCl₄ molecules. In all cases, the calculated β_HRS,Liq/β_HRS,Gas ratios remain close to one, much smaller than the experimental results,¹ as well as previous predictions obtained using the polarizable continuum model.⁶ On the other hand, considering a small number of CCl₄ molecules, embedded in point charges, enables to monitor the emergence of the dipolar contribution to β_HRS,Liq and the increase of the depolarization ratio. These calculations confirm the experimental data¹ and substantiate for the first time that the dipolar contributions originate from intermolecular interactions between the CCl₄ molecules. To get closer to the experimental dipolar/octupolar ratio as well as to fully reproduce the solvent effects in particular the β_HRS,Liq/β_HRS,Gas ratio, further investigations should consider much larger numbers of CCl₄ molecules, although this appears computationally challenging owing to the fact that highly extended basis sets are required to properly describe β_HRS of a single CCl₄ molecule. A first step towards that direction has been made by employing PM6 and PM7 semi-empirical calculations to describe the second-order NLO responses of systems containing explicitly one and two solvation shells. Though the absolute β values differ from the ab initio ones, these simulations also predict the appearance of a dipolar contribution to the HRS signal in liquid CCl₄. Another direction of investigation would be to combine the QMMM approach with electrostatic interaction schemes as done for other types of systems.^42,43 Finally, further refinements would also imply going beyond the rigid molecule approximation, by coupling the QM calculations to atomistic MD simulations allowing geometrical fluctuations of the CCl₄ molecules.

Acknowledgements

M. H. C. acknowledges the CNPq (Brazil) and the University of São Paulo for his postdoctoral grant as well as the Belgian Government (IUAP No. P7/5 “Functional Supramolecular Systems”) for financial support. This work is supported by funds from the Belgian Government (IUAP No. P7/5) and the Francqui Foundation. The calculations were performed on the computers of the Consortium des Équipements de Calcul Intensif, including those of the Technological Platform of High-Performance Computing, for which we gratefully acknowledge the financial support of the FNRS-FRFC (Conventions No. 2.4.617.07.F and 2.5020.11) and of the University of Namur as well as on the Mésocentre de Calcul Intensif Aquitain (MCIA) of the University of Bordeaux, financed by the Conseil Régional d'Aquitaine and the French Ministry of Research and Technology.

References

1. P. Kaatz and D. P. Shelton, Mol. Phys., 1996, 88, 683.
2. P. Kaatz, E. Donley and D. P. Shelton, J. Chem. Phys., 1998, 108, 849.
3. R. D. Pyatt and D. P. Shelton, J. Chem. Phys., 2001, 114, 9938–9946.
4. J. Campo, F. Desmet, W. Wenseleers and E. Goovaerts, Opt. Express, 2009, 17, 4587–4604.
5. D. P. Shelton, J. Chem. Phys., 2012, 137, 044312.
6. F. Castet, E. Bogdan, A. Plaquet, L. Ducasse, B. Champagne and V. Rodriguez, J. Chem. Phys., 2012, 136, 024506.
7. D. M. Bishop, F. L. Gu and S. M. Cybulsky, J. Chem. Phys., 1998, 109, 8407.
8. M. S. A. El-Kader and Y. N. Kalugina, Chem. Phys. Lett., 2015, 639, 93–98.
9. M. de Wergifosse, F. Castet and B. Champagne, J. Chem. Phys., 2015, 142, 194102.
10. M. Hidalgo Cardenuto and B. Champagne, J. Chem. Phys., 2014, 141, 234104.
11. M. Hidalgo Cardenuto and B. Champagne, Phys. Chem. Chem. Phys., 2015, 17, 23634.
12. J. Quertinmont, B. Champagne, F. Castet and M. Hidalgo Cardenuto, J. Phys. Chem. A, 2015, 119, 5496.
13. K. Coutinho, R. Rivelino, H. C. Georg, and S. Canuto, in Solvation Effects in Molecules and Biomolecules: Computational Methods and Applications, Springer, 2008, p. 159.
14. K. Coutinho and S. Canuto, Adv. Quantum Chem., 1997, 28, 89.
15. S. Canuto and K. Coutinho, Int. J. Quantum Chem., 2000, 77, 192.
16. K. Coutinho and S. Canuto, J. Chem. Phys., 2000, 113, 9132.
17. K. Coutinho, S. Canuto and M. Zerner, J. Chem. Phys., 2000, 112, 9874.
18. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendon Press, 1987.
19. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. Phys., 1953, 21, 1087.
20. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999.
21. J. Tomasi, Theor Chem Acc, 2004, 112, 184.
22. B. G. Rao and U. C. Singh, J. Am. Chem. Soc., 1991, 113, 4381.
23. E. M. Duffy, D. L. Severance and W. L. Jorgencen, J. Am. Chem. Soc., 1992, 114, 7535.
24. K. Coutinho, H. C. Georg, T. L. Fonseca, V. Ludwig and S. Canuto, Chem. Phys. Lett., 2007, 437, 148.
25. T. Helgaker, S. Coriani, P. Jorgensen, K. Kristensen, J. Olsen and K. Ruud, Chem. Rev., 2012, 112, 543.
26. F. Castet and B. Champagne, J. Chem. Theory Comput., 2012, 8, 2044.
27. H. D. Cohen and C. C. J. Roothaan, J. Chem. Phys., 1965, 43, S34.
28. M. de Wergifosse, V. Liégeois and B. Champagne, Int. J. Quantum Chem., 2014, 114, 900.
29. A. A. K. Mohammed, P. A. Limacher and B. Champagne, J. Comput. Chem., 2013, 34, 1497.
30. J. J. P. Stewart, J. Mol. Model., 2007, 13, 1173.
31. J. J. P. Stewart, J. Mol. Model., 2013, 19, 1.
32. T. Verbiest, K. Clays, and V. Rodriguez, Second-Order Nonlinear Optical Characterizations Techniques: An Introduction, CRC Press, New York, 2009.
33. K. Coutinho and S. Canuto, DICE (version 2.9): A general Monte Carlo program for liquid simulation, University of São Paulo, 2009.
34. J. Frisch, et al., Gaussian 09, Revision D01, Gaussian, Inc., Wallingford, CT, 2009.
35. J. J. P. Stewart, Stewart Computational Chemistry, MOPAC, Colorado Springs, CO, USA, 2016.
36. A. H. Narten, J. Chem. Phys., 1976, 65, 573.
37. R. W. Gruebel and G. T. Clayton, J. Chem. Phys., 1967, 46, 639.
38. T.-M. Chang, K. A. Peterson and L. X. Dang, J. Chem. Phys., 1995, 103, 7502.
39. J. C. Soetens, G. Jansen and C. Millot, Mol. Phys., 1999, 96, 1003–1012.
40. R. Rey, L. C. Prado, E. Llanta, K. Ando, D. O. López, J. L. Tamarit and M. Barrio, J. Chem. Phys., 2000, 112, 7505.
41. R. Rey, J. Chem. Phys., 2007, 126, 164506.
42. T. Seidler, K. Stadnicka and B. Champagne, J. Chem. Theory Comput., 2014, 10, 2114.
43. I. Harczuk, O. Vahtras and H. Ågren, Phys. Chem. Chem. Phys., 2016, 18, 8710.

---