First hyperpolarizability of water in bulk liquid phase: long-range electrostatic effects included via the second hyperpolarizability

Abstract

The molecular first hyperpolarizability β contributes to second-order optical non-linear signals collected from molecular liquids. For the Second Harmonic Generation (SHG) response, the first hyperpolarizability β(2ω, ω, ω) often depends on the molecular electrostatic environment. This is especially true for water, due to its large second hyperpolarizability γ(2ω, ω, ω,0). In this study we compute the electronic γ(2ω, ω, ω,0) and β(2ω, ω, ω) for water molecules in liquid water using QM/MM calculations. The average value of γ(2ω, ω, ω,0) is smaller than the one for the gaz phase, and its standard deviation among the molecules is relatively small. In addition, we demonstrate that the average bulk second hyperpolarizability 〈γ(2ω, ω, ω,0)〉 can be used to describe the electrostatic effects of the distant neighborhood on the first hyperpolarizability β(2ω, ω, ω). In comparison with more complex schemes to take into account long-range effects, the approximation is simple, and does not require any modifications of the QM/MM implementation. The long-range correction can be added explicitly, using an average value of γ for water in the condensed phase. It can also be easily added implicitly in QM/MM calculations through an additional embedding electric field, without the explicit calculation of γ.

---

1 Introduction

Non-Linear Optical (NLO) technique(s) are increasingly used to probe structural properties of matter. In Second Harmonic Generation (SHG) or Hyper-Rayleigh Scattering (HRS), two photons with the same fundamental frequency interact with a nonlinear material that generates a new photon with twice the energy of the initial photons.¹ SHG-based technologies have been developed to investigate simple and complex liquids,^2,3 biomimetic systems,^4–6 or even biological materials.^7,8 One interesting property of non-resonant SHG is the sensitivity of the response to the electrostatic environment.⁹ Historically, the prototypical application is the electric-field-induced second harmonic generation (EFISHG) of molecules in a gas phase.¹⁰ During such experiments, a macroscopic electrostatic field E^DC is applied and the SHG response of the system is described using second and third order susceptibility tensor: χ⁽²⁾(2ω, ω, ω) and χ⁽³⁾(2ω, ω, ω,0), respectively. For a mesoscopic unitary volume, the induced total dipole moment at the second harmonic frequency P^2ω expressed as a sum of two terms:

P^2ω ∝ χ⁽²⁾(2ω, ω, ω): E^ωE^ω + χ⁽³⁾(2ω, ω, ω,0): E^ωE^ωE^DC, (1)

where E^ω is the exciting electromagnetic field at the fundamental frequency, and E^DC is the external electrostatic field that potentially modifies (i) the orientation of dipoles, and (ii) the electronic-dependent molecular intrinsic SHG response. More recently, EFISHG has received renewed interest as a non-intrusive tool for probing electric fields in diluted media with spatial and temporal definition, using laser pulses.¹¹

Eqn (1) was also applied to condensed phases,^12,13 and the EFISHG is an established technique to determine the first and the second hyperpolarizabilities of compounds in solutions.¹⁴ Even more, studies on liquid/solid or liquid/air interfaces have reported the evolution of the Surface-SHG (S-SHG) response when the surface charge is modulated. Frameworks based on eqn (1) for the fluid nearby the surface have been commonly used in the S-SHG community to extract a surface potential or an effective surface charge.^15–19 But the interpretation of the different terms in the S-SHG intensity generated by aqueous solutions is still the subject of many recent works^20–23 requiring theoretical calculations at the molecular level.²⁴

To provide an interpretation at the microscopic scale, eqn (1) can also be rewritten at the molecular level:

where μ^2ω is the induced dipole moment of the molecule at the second harmonic frequency, e^ω the exciting electromagnetic field in the molecular frame, β⁽²⁾ (2ω, ω, ω) and γ⁽³⁾ (2ω, ω, ω,0) are respectively the first and second molecular hyperpolarizabilities, noted β and γ later on. The second hyperpolarizability γ describes the variations of the first hyperpolarizability value under a static external electrostatic field e^DC, present in the molecular frame. This approach permits to link the EFISHG macroscopic measurements for the gas phase to quantum calculations of the molecular first hyperpolarizability β.^25–30

The application of field expansions such as eqn (2) for condensed phases is very useful for the interpretation of S-SHG of aqueous solutions, and the values of β and γ for water in liquid water become parameters of the equation. To obtain values of β and γ for water in liquid water, one can use Quantum Mechanical/Molecular Mechanic (QM/MM) calculations.^{9,24,26,31,32} In such approaches, the electrostatic effect of the environment is included within the QM calculations of a molecular hyperpolarizability β^env for a given electrostatic environment, so that

where β^env depends on the electric field exerted by the environment e^env (T convention³³). Using the values of β^env, one can also calculate γ^env from its definition as , where β_ijk (e_l) is the ijk-component of the first hyperpolarizability of a molecule on which an extra electrostatic field e_l is applied along the l molecular direction. In the present work, we first provide the β^env and γ^env tensors of individual water molecules in the liquid phase using QM/MM approaches at an optical wavelength of 800 nm (typically used for experiments^34,35). Relatively to previous works,^9,26,32 we detail all the components of both tensors, and discuss the distribution of γ among water molecules.

In a second step, possessing a database of {β^env,γ^env} for individual molecules in liquid water, we can additionally scrutinize the validity of the field perturbation approach represented by eqn (2). Similar field expansions for the calculation of solvent effects have been tested for observables such as nuclear magnetic shieldings,³⁶ or linear and non linear optical properties.³⁷ However, the precision of the field expansion for the non-linear optical properties were limited,³⁷ and more terms are needed in the perturbation expression.³⁸ Concerning the first hyperpolarizability of water, Liang et al.⁹ have compared β^env of water within the liquid phase obtained by QM/MM calculations to the linear field expansion of eqn (2), where e^DC was the electrostatic field generated by the surrounding water molecules, noted e^env, and β and γ were the gas-phase values: a quantitative difference emerged. It is therefore expected that the strong solvation of water molecules within liquid water cannot be simplified to eqn (2). Nevertheless, in this work, we show that the perturbation field Taylor expansion can be useful to predict long-range electrostatic effects on the hyperpolarizability of water molecules in their liquid phase.

In the following, Section 2 describes the numerical details. In Section 3, the individual values of γ^env and β^env are first reported. Then, we show that the electrostatic field generated by the environment, namely e^env, is strongly heterogeneous in space. Based on this result, we propose to separate the effects of the electric field generated by the environment into a short- and a long-range part. Using this separation, we show that, for a water molecule in the bulk phase, the effect of the long-range electrostatic environment on β^env can be included using a linear correction proportional to γ.

2 Method

The calculations are based on a sequential QM/MM approach in two steps. First, we use a classical Molecular Dynamics (MD) simulation of bulk water to obtain typical structures of the liquid. Then, we compute the first and second hyperpolarizabilities of individual water molecules at the Density Functional Theory (DFT) level, within an electrostatic embedding framework. To investigate fluctuations due to the changing environment, statistics are performed over numerous configurations of the MD simulation.

2.1 Molecular dynamics

LAMMPS,³⁹ V.11.08.2017 is used to perform the MD simulation along with the rigid TIP4P/2005 water force field.⁴⁰ 15625 rigid TIP4P/2005 water molecules are placed in a simulation box (approx 7.8 × 7.8 × 7.8 nm³) to form a 3D-periodic bulk system. We have used an isothermal, isobaric ensemble (NPT) with Nose–Hoover thermostat at 300 K (τ_T = 0.4 ps) and Nose–Hoover barostat at 1 atm (τ_P = 2 ps).⁴¹ 1 ns of equilibration is performed before the 1 ns production run, both with a time step of 2 fs. Both electrostatic and Lennard-Jones intermolecular interactions are computed using the long-range Particle-Particle-Particle-Mesh (PPPM) formalism.^42,43 Neighbor lists are updated every time step within a radius of 10 Å. This simulation leads to a density of about 0.996 kg L⁻¹. Such a large system was necessary to investigate environmental effects up to 40 Å.

2.2 First hyperpolarizability in the liquid phase

To compute the first hyperpolarizability of water molecules within the liquid phase, an explicit environment composed of point charges is used: the Polarizable Embedding model at the zero^th-order (PE0),^44,45 implemented in the DALTON software,⁴⁶ release 2018.2 package. As in our previous recent work,³² the QM calculations were carried out on individual water molecules. Point charges represent the surrounding water molecules, and the same electrostatic description is used as for the MD (TIP4P/2005 model). The MD trajectories are used directly in the QM calculations without further optimization, in contra. We define the parameter R_c as the maximal distance upon which neighbors are included in the PE formalism. The electrostatic field generated by this environment is spatially-heterogeneous in the vicinity of the target molecule and creates a new potential for the target molecule electronic Hamiltonian.

The QM/MM approach could be refined by optimizing geometries,⁴⁷ or by adding dipole, quadrupole or polarizability of the environment,^45,48–51 but we shall see in the results that eqn (2) becomes relevant only for long-range effects, where the point-charges electrostatic field is dominant. We have therefore chosen to keep the environment description as simple and robust as possible. To optimize parts of the QM/MM routines, the home-made software FROG is used.‡

The whole procedure is composed of the following steps: (1) perform MD simulation using LAMMPS. (2) Read the MD trajectories using FROG to build the electrostatic environment for each water molecule defined into a sphere of radius R_c centered on the mass center of the target molecule. If a neighbor possesses at least one atom in the sphere, the whole molecule is considered in the environment. The output files for DALTON are created for 300 molecules for 8 selected MD configurations. (3) For each of the 2400 target molecules, calculate the electrostatic field, and its gradients, generated by the environment using FROG. (4) For each of the 2400 target molecules, calculate the value of β at the fundamental wavelength of 800 nm using DFT with the functional CAM-B3LYP⁵² and the basis d-aug-cc-pVTZ,⁵³ using the quadratic response formalism as implemented in DALTON.⁵⁴ (5) Read all DALTON output files and analyze the data using FROG.

The largest inaccuracy of our approach is due to the QM method and the corresponding functional. For the first hyperpolarizability of water molecules, the values in vacuum obtained with our method (DFT/CAM-B3LYP) are typically overestimated by about 10%³² in comparison to the golden standard for QM/MM calculation of small molecules (CCSD^55–57). Both calculation schemes use the TIP4P/2005 rigid force field, the MM size R_c = 15 Å and the basis is d-aug-pVTZ. The basis set d-aug-cc-pVTZ⁵³ used is widespread across the community^26,55–57 and presents a good balance between accuracy and cost. In a previous work,³² we have observed a smooth convergence of the first hyperpolarizability relative to the basis-set size for water in bulk water, and similar absolute values of the nuclear and electronic contributions to the electrostatic and polarization energies. Therefore, we estimate that spill-out effects are not dominant,⁵⁰ even if we use the large d-aug-cc-pVTZ basis set. We have not taken into account the nuclear degree of freedom (vibrational effects) since Beaujean et al.⁵⁸ have shown that this effect is small for the second hyperpolarizability of water, especially at small wavelengths.

2.3 Second hyperpolarizability in the liquid phase

In this work, we are interested in the evolution of the first hyperpolarizability under a static and homogeneous electrostatic field, which is reflected by γ(2ω, ω, ω, 0) noted here γ. For each individual molecule, the second hyperpolarizability γ can be computed from the first hyperpolarizability β using the Finite Field (FF) framework. For that, we use the PE formalism described above to compute β when an extra homogeneous electrostatic field is applied, and then compute γ according towhere β_ijk (e_l) is the first hyperpolarizability of a molecule on which an extra static electrostatic field e_l is applied along the l molecular direction. Numerically, the γ_ijkl is obtained for each molecule as the slope of a linear fit performed on seven calculations with e_l from 0 to 1.5 × 10⁻³ a.u. Noticeably, during the γ calculation, the target molecule is embedded in a total electrostatic field composed of a spatially-heterogeneous contribution from its PE environment plus the spatially-homogeneous e_l. More details regarding the FF procedure can be found in the ESI,† Section S1.1. The first and second hyperpolarizabilities are given in atomic units, our conventions are defined by eqn (2) and (4). They are expressed in the molecular frame {a, b, c} with the axis c along the dipole moment pointing from the Oxygen towards the Hydrogens, and the atoms within the {a, c} plane, see insert in Table 1. Noticeably, one could alternatively have used the response scheme to the third order, as implemented in Dalton, to obtain the second hyperpolarizability. For a given QM/MM conformation, the numerical cost of the response scheme is about ten times one of the finite field (FF) scheme, for similar results (see ESI,† Table S1). The comparison of the FF and the response schemes provides an estimation of the numerical error due to the FF procedure: about 5 to 20 a.u.

Table 1 Water second hyperpolarizability γ⁽³⁾ (2ω, ω, ω,0) components, either in the vacuum (γ^vac) or in its liquid phase (γ^env) for an exciting wavelength of 800 nm, in atomic units. For the bulk phase, the average value 〈γ^env〉 and the standard deviation σ [γ^env] are reported. The tensors components are expressed in the molecular frame {a, b, c}, see inserted figure. Comparing with response calculations, the error on γ due to the finite field differentiation is estimated at about 20 a.u. (see ESI, Table S1)

ijkl	Vacuum	Liquid
	γ ^vac	〈γ^env〉	σ[γ^env]
aaaa	1110	740	140
bbbb	4090	2940	650
cccc	2000	1400	230
aabb	930	640	200
bbaa	1060	710	150
baab	1010	680	220
abba	940	660	140
aacc	590	370	60
ccaa	610	380	80
caac	610	380	60
acca	590	370	80
bbcc	1120	770	170
ccbb	1010	720	180
cbbc	1010	730	150
bccb	1060	750	200

---

Frequency dispersion can be important for γ: for example, Beaujean et.al.⁵⁸ have shown a variation of about 25% between an excitation wavelength of 800 nm and at the infinite wavelength limit for γ(3ω, ω, ω, ω). Here, even if we are using a Finite Field approach, the frequency dispersion of β related to the fundamental frequency (here 800 nm) is correctly reproduced since we use the first hyperpolarizability obtained using the frequency dependent response scheme.

2.4 Embedding environment size: R_c

One of the objectives of this study is to understand the impact of the electrostatic environment on the first hyperpolarizability, distinguishing the effects of close and far neighbors. Therefore, different values of R_c are used for β calculations. In the following, the notation β^PE (R_c) denotes the hyperpolarizability obtained by a QM/MM calculation at a given R_c and e(R_c) the electric field generated by this PE environment. In this notation, R_c = 0 corresponds to a gas phase calculation. According to previous works,^32,57,59 a large environment (typically between 10 and 20 Å) is necessary to obtain a good convergence of the β values. Here, for each molecule, a reference calculation is performed using the PE formalism with an extremely large radius of R_f = 40 Å, and we consider that β^env = β_PE (40 Å). Similarly, we consider that the electrostatic field e^env = e (40 Å). Note that such a large radius still does not represent the complete MD liquid phase, described using periodic boundary conditions.

On the contrary, the variation of γ with the size of the environment is not investigated in details here. Therefore, the PE radius used to compute the second hyperpolarizability γ is set to 10 Å. Yet, we have verified that the γ does not evolve much using larger radius in agreement with Osted et al.⁵⁹ (data not shown) and have noted it γ^env.

2.5 Statistical averaging

The first and second hyperpolarizability in the liquid phase are obtained from 2400 configurations of the PE environment. These configurations are extracted from 8 MD snapshots, with 300 water molecules randomly selected in each frame. These frames are separated by 100 ps to ensure time decorrelation.³² We have verified that this number of configurations is sufficient to reach convergence; indeed, the numerical error due to the configuration number is about 1% (see ESI,† Fig. S2 and S3). This is smaller than the error due to the DFT Finite field calculation relative to CCSD calculations (see ESI,† Tables S2 and S3). The C_2v symmetry of the first and second hyperpolarizability obtained by QM/MM confirms this good convergence (see the Section 3 and ESI,† Section S2.1).

3 Results

In the following, Section 3.1 first presents the second hyperpolarizability tensor γ calculated using QM/MM approaches for single water molecules within the liquid bulk phase, and the standard deviation of the calculated values. Then, Section 3.2.1 presents some difficulties encountered when applying eqn (2) for condensed phases. We attribute them to the properties of the electrostatic field generated by the environment, e^env, described in 3.2.2. Based on these results, and on eqn (2), in Section 3.2.3, we propose a correction to β that permits to include long-range interactions. Finally, we show that this correction enhances the precision of standard QM/MM calculations.

3.1 QM/MM results: second hyperpolarizability of water in the bulk phase

Table 1 presents the values of selected γ components calculated either in the vacuum or in liquid phase – the other non-zero components of γ are presented in ESI,† Table S4.

For the vacuum phase, the molecular C_2v symmetry is fulfilled: the only non-vanishing components are the γ_iiii and γ_iijj for i and j in {a, b, c}, plus their permuted terms¹ (see ESI,† Table S4). In the liquid phase, this symmetry is not valid for individual values, but is recovered for the average 〈γ^env〉. Our results agree with the ones by Osted et al.⁵⁹ within about 15%. We impute the discrepancies to both the difference of methods (DFT vs. CCSD, different MD models) and the difference of excitation wavelengths (800 vs. 1080 nm). More detailed comparisons with literature are provided in ESI,† Section S1.1.

Noticeably, when transferred into the laboratory frame, the average second hyperpolarizability becomes centro-symmetric (see ESI,† Table S5), indicating that our sampling of molecular orientation in the liquid phase is sufficient. Interestingly, the average value of all C_2v-authorized γ^env components are positive. Moreover, in the liquid phase, the γ^env components are about 30% smaller than the ones in the vacuum – see ESI,† Fig. S4.

Beyond the average values, the environment induces fluctuations of γ^env, and of β^env. A typical and relevant case is illustrated by Fig. 1, displaying the join probabilities of {β_ccc,γ_cccc} obtained for the 2400 water molecules within their environment. Projections on the axes display the distributions of β_ccc and γ_cccc separately. As already discussed in the literature,^32,57,59 the water first hyperpolarizability fluctuations in the liquid phase are strong, with standard deviations similar to the absolute values of the averages – for instance 〈β^env_ccc〉 = 4.1 a.u. and σ [β^env_ccc] = 3.1 a.u. For individual molecules, some β components, null in average, can be larger than the one with a net average (see ESI,† Fig. S4).

Fig. 1 Join population of {β^env_ccc,γ^env_cccc} obtained by QM/MM approaches for bulk liquid water. The independent ones are presented on the top and right for the β^env_ccc and γ^env_cccc respectively. The values are normalized with respect to their average ones:〈 β^env_ccc〉 = 4.1 a.u. and 〈γ^env_cccc〉 = 1400 a.u.

In complement to the results published by Osted et al.,⁵⁹ we show that the second hyperpolarizability γ also fluctuates, but in a smaller extent than β: the standard deviations of the different components σ [γ^env] are about 20% of their average values 〈γ^env〉 – see Fig. 1 for γ_cccc, and ESI,† Table S4.

Fig. 1 moreover shows that the values of β and γ are not strongly correlated. Hence, in the following, we neglect the γ^env dispersion and attribute the same second hyperpolarizability 〈γ^env〉 to all the water molecules.

3.2 Using γ to calculate individual β in condensed phases

After describing the second hyperpolarizability γ^env in the liquid phase, we will question if it can be used to predict the first hyperpolarizability of water molecules in the liquid phase using different approximations inspired by eqn (2). We shall compare the β^env obtained by QM/MM to linear expansions relative to a perturbing electric field.

3.2.1 Limitations of the gas-phase reference model. Comparing eqn (2) and (3), one straightforward approach to include the effect of the environment on the first hyperpolarizability β is to use the second hyperpolarizability γ and the field created by the environment e^env:

β^GR = β^vac + γ⁰·e^env, (5)

where a linear perturbation is done relative to the Gas-phase Reference (GR), for which e^env = 0 and and the reference first and second hyperpolarizability would be the one of the vacuum for instance. For pedagogical reasons, we will first highlight here the limitation of this gas-phase reference model, that was already mentioned in the literature.^36,57

If the linear proportionality factor γ⁰ is known, and the same for all molecules, this approach would be very helpful: the liquid first hyperpolarizability depends only on e^env. An MD simulation can provide e^env for each molecule, and the liquid phase first hyperpolarizability is calculated using eqn (5), so that no expensive QM/MM approach is needed. Moreover, such an approximation is commonly used in the Surface-SHG community (for example ref. 18 and 19).

Liang et al.⁵⁷ have reported that eqn (5) using the vacuum value for the second hyperpolarizability γ⁰ = γ^vac is a poor estimator for β^env. We have thus applied eqn (5) using the average bulk second hyperpolarizability tensor γ⁰ = 〈γ^env〉 for all the molecules (values from Table 1 and from ESI,† Table S4), and the electrostatic field generated by the environment e^env. For liquid water, the distribution of the largest component of e^env, along the water molecular axis c is represented on Fig. 2: it has typical values around 3 × 10⁻² a.u., i.e. 1.6 V Å⁻¹. The other components average to zero, with approximately the same standard deviations (data not shown).

Fig. 2 Distribution of the electrostatic field generated by the environment e^env along the molecular c-axis. The electrostatic field is measured at 2 positions: at the oxygen atom position (O: orange squares), or at the negative charge one (T: blue dots) according to the TIP4P/2005 force field. The O and T positions are distant by 0.15 Å.

The predicted average values 〈β^GR〉 are compared to the QM/MM values 〈β^env〉 in Table 2. First, the addition of the linear correction in eqn (5) allows recovering signs for β^GR components in agreement with β^env. However, the β^GR model largely over-estimates the reference values: the correction γ⁰·e^env should be smaller. This overshooting is even worse when the vacuum value γ^vac is used for γ⁰, as Liang et al. had tested.⁹ Thus, even if it contains some relevant physics, eqn (5) describing a correction that is proportional to e^env is not workable: it predicts an average hyperpolarizability 〈β^GR〉 one order of magnitude too large.

Table 2 Comparison between relevant QM/MM components of the first hyperpolarizability, β^env, the vacuum ones, β^vac, and β^GR predicted by eqn (5) with γ⁰ = 〈γ^env〉 from Table 1. All values are in atomic units and are calculated for fundamental wavelength of 800 nm. The values of e^env obtained at the center-of charge position have been used

ijk	〈β^vac〉	〈β^env〉	〈β^GR〉
ccc	−15.3	4.1	28.5
caa	−12.5	−2.0	−0.6
aca	−12.4	−2.0	−0.8
cbb	−5.0	2.5	17.7
bcb	−7.4	2.2	16.9

---

Moreover, the value of e^env in eqn (5) is ill-defined because of its strong spatial heterogeneity. To highlight this fact, Fig. 2 presents the electrostatic field created by the environment along the molecular c-direction, e^env_c, calculated at 2 positions within the embedded molecule: either at the position of the negative charge of the TIP4P/2005 model, or at the position of the oxygen atom. While the 2 distributions look very similar, a difference on the average values of about 10⁻³ a.u. appears. Given the large value of γ_cccc (≈1400 a.u., see Table 1), the component β_ccc predicted by β^GR changes by about 40%, depending on where e^env is calculated. Therefore, we attribute the failures of eqn (5) to the strong intensity and large gradient of the embedding field e^env, which are incompatible with a simple linear expansion. Such limitations of perturbation field expansions to model solvation effects had already been discussed in the literature. When studying the nuclear magnetic shielding tensors of the atoms of N-methylacetamide solvated by liquid water, Kjaer et al.³⁶ have shown that the field experienced by the nuclei is very inhomogeneous, so that considering one averaged field from all configurations and/or across all nuclei is not a good approximation.

3.2.2 Separation into short and long-range neighbors. To go beyond previous works,^57,59 we propose to rewrite the pertubation field expansion so that the perturbing electric field and electric field gradients are weak, and the linear regime is verified. The evolution of the electrostatic field e(R_c) generated by the environment around the target molecule of the QM/MM calculations with the environment size R_c is presented on Fig. 3. The main contribution to this electrostatic field is created in the first and second solvation shells, at distances between 2 and 5 Å from the target molecule. This neighborhood region also creates the largest electrostatic field gradients, see ESI,† Fig. S6.

Fig. 3 Averaged norm of the electrostatic field generated by the electrostatic environment (e(R_c)) felt at the T-position, for increasing environment size R_c. In insert, the difference Δe(R_c) relative to the value at R_c = 40 Å, in logarithmic scale as a function of the size R_c.

Given that the closest shells are responsible for both the intensity and the heterogeneity of e^env, we divide the neighborhood into two parts, as illustrated on Fig. 4.

Fig. 4 Sketch of the electrostatic embedding procedure. The QM box is defined for only one molecule of water at the center. Direct area, until R_c: all the neighbors are included in a QM/MM calculation with discrete solvation procedure, such as the PE formalism. Long-range area, from R_c to R_f: the neighbors are included by the electrostatic field they generate at one specific point of the QM box.

The direct area contains the closest neighbors of the target molecule (at the center), up to a distance R_c: they create a strong and heterogeneous electric field, noted e(R_c). The long-range area, due to neighbors beyond R_c from the target molecules, creates a less intense and more homogeneous electric field on the target molecule.

The total embedding electrostatic field e^env is thus separated into a short-range contribution e(R_c), and a correction, Δe(R_c) due to the long-range electrostatic interactions:

e^env = e(R_c) + Δe(R_c). (6)

where all the fields are evaluated at the position of the negative charge in the TIP4P/2005 model (T-position, see Fig. 2 inset). The evolution of Δe (R_c) as a function of R_c is available in the insert of Fig. 3. At R_c = 10 Å, the electrostatic field generated by the long-range neighbors is very small compared to the one created by close neighbors: the average correction Δe (R_c = 10 Å) is close to 10⁻³ a.u. The finite field calculations reported in Section 3.1 have shown that for such amplitudes of the perturbation, the dependence of the QM/MM β^env relative to an external, homogeneous field is linear (see ESI,† Fig. S1). Moreover, we have investigated whether the correction Δe (R_c) is spatially homogeneous by measuring the correction to the electric field gradient (see ESI,† Fig. S1). Generally, the corrections to the gradients are small (typically 10⁻⁴ a.u) for neighbors further away than 10 Å. Since the correction Δe(R_c) is more homogeneous than e^env, it hardly depends on the point on which it is calculated. In the following, we use the value of Δe(R_c) evaluated at the T-position.

3.2.3 Long-range correction to β using γ. The further the neighbor is from the target molecule, the better its contribution to e^env can be modeled as a weak and homogeneous electrostatic field. A linear γ correction to β could be tested to describe the long-range electrostatic effects.

For the first hyperpolarizability, as for e^env, we describe the impact of the environment in two different ways depending on the neighbors distance. The direct neighbors within the sphere of radius R_c are included explicitly in the QM/MM calculation through the PE approach. The long-range neighbors, beyond R_c, are included implicitly via the homogeneous field they produce on the target molecule, Δe(R_c), with a correction proportional to 〈γ^env〉. This approximation can be viewed as a linear expansion relative to the value of first hyperpolarizability already perturbed by the direct neighbors, β^PE(R_c):

β^PE+L(R_c) = β^PE(R_c) + 〈γ^env〉·Δe(R_c). (7)

In the following, we compare this approximation of β, noted β^PE+L(R_c) and the usual QM/MM calculation where only the PE scheme up to a distance R_c is used, noted β^PE(R_c). To estimate the accuracy of the two approximations PE and PE + L as a function of R_c, we compare them to our reference values β^env.

As a typical example, we compare on Fig. 5 the β_ccc computed with R_c = 5 Å either using PE approximation (top) or the PE + L one (bottom) to the reference value β^env_ccc. The diagonal dashed lines represent the case where the approximation reproduces exactly the reference.

Fig. 5 Comparison of the reference individual molecular β^env_ccc component and the one obtained using PE (TOP) and PE + L (BOTTOM) approximation with R_c = 5 Å. The dashed line corresponds to the ideal case where a correction matches perfectly the expected value. Each dot represents a molecule.

For the usual calculation (β^PE), there is a systematic error leading to inaccuracy: at R_c = 5 Å, the average value of β^PE_ccc is smaller than the one of β^env_ccc. Moreover, the error fluctuates a lot, it can reach tens of a.u., depending on the molecule. On the contrary, the PE + L approach almost quantitatively reproduces the reference, even for R_c as small as 5 Å

To illustrate the effect of R_c, Fig. 6 depicts the evolution of 〈β_PE〉 and 〈β_PE+L〉 as a function of R_c for the ccc component (similar trends are obtained for other components). As expected, the inaccuracy of the PE and PE + L approaches decreases when R_c increases.

Fig. 6 Average of the β^X_ccc (R_c) component as a function of R_c for the PE (blue disks) and PE + L (orange squares, see eqn (5)) approximations. The dashed line is the reference value β^env obtained for R_f = 40 Å.

When the correction is added, the convergence is reached faster (R_c = 4 Å) compared to the usual β^PE calculation (R_c ≃ 12 Å). The effect of the neighbors further than 5 Å can thus be very efficiently included by a second hyperpolarizability corrective term. Beyond the description of the average first hyperpolarizability, Fig. 5 also illustrates that the PE + L approach improves the precision of the individual predictions. This can be quantified through the mean absolute error (MAE), averaged over all the N = 2400 molecular configurations and all the 27 components, noted Δβ_T:

where X stands for the approximation PE, or PE + L, and the exponent (n) indicates that the values of the n-th molecule are used. Fig. 7 reports the evolution of Δβ^X_T as a function of R_c.

Fig. 7 Mean absolute error of the individual first hyperpolarizability tensor defined by eqn (8) as a function of R_c for the PE (blue disks), PE + L (orange squares) and PE + LI (green triangles) approximations. The dashed line corresponds to an error of 0.1 a.u. at which we define convergence.

The long-range correction provides very good results and reaches an error below 0.1 a.u. at R_c = 12 Å. To obtain the same degree of accuracy, one would need to include explicitly neighbors up to R_c ≃ 30 Å in the usual PE scheme without long-range correction. For R_c larger than 10 Å, the error obtained with the long range correction (PE + L) is about one order of magnitude lower than the one with the traditional PE approach. This shows that – once the strong short-range effects are explicitly included – it is possible to use the knowledge of γ to predict accurately the long-range electrostatic effects on β. Note that the long-range correction is proportional to the average value 〈γ^env〉, i.e.Eqn (9) neglects γ^env fluctuations. As discussed further in ESI,† Section S4, including the fluctuations of γ^env in eqn (9) only slightly improves the results. However, it implies to compute γ^env for every molecule, which greatly increases the numerical cost. Therefore, using an average γ^env is a very good compromise between computational cost and accuracy, at least for bulk liquid water.

3.2.4 Implicit long-range correction. The linear correction is validated here for water molecules in bulk liquid water. We consider this result as a first step towards more complex systems, anisotropic or charged ones, where electrostatic effects become even more important. In these new environments, the second hyperpolarizability of water might differ from the one of the bulk phase, and its fluctuations may become larger. With these next steps in mind, we show here that our explicit long-range correction can also be added to the system implicitly, during the QM calculation, via the addition of the homogeneous field e(R_c) on the top of the heterogeneous MM embedding field. We note this approach PE + LI for Long-range Implicit correction, than can be noted

β^PE+LI (R_c) = β^PE (R_c; e(R_c)). (9)

Fig. 7 illustrates that the PE + LI values of the first hyperpolarizability are very close to the reference values. The PE + LI approach is even more precise that the PE + L because (1) the second hyperpolarizability fluctuations are neglected in PE + L – not in PE + LI, and (2) the second hyperpolarizability is calculated with finite field approach in PE + L – not in PE + LI. The PE + LI correction does not require the knowledge of the second hyperpolarizability. It can be seen as a simplification of more elaborated QM/MM developments taking into account long-range environment in non-periodic QM/MM methods.^60–63 Indeed, long-range effects can be approximated by boundary potential models, by cutoff method coupled with switching and shifting functions, by the so-called Long-Range Electrostatic Correction, or by additional charges on the edge of the MM embedding, or in the MM embedding. Our simpler approach with a homogeneous electrostatic field is intended to provide the impact of the remote MM environment on the non-linear optical property of the QM system. It is limited to sequential QM/MM approaches where the trajectory is obtained independently. The advantages of the PE + LI are (1) its insight in the physical origin of the correction and its link with usual equations used in the experimentalist community, and (2) its simplicity and low computational cost. Finally, the value of e(R_c) could include periodic boundary conditions, but it is beyond the present work.

4 Conclusions

In this work, we have computed the first and second hyperpolarizability of water in the liquid phase using state-of-the-art QM/MM calculations. With this set of values, several approaches were tested to take into account the electrostatic environment in the first hyperpolarizability.

The second hyperpolarizability is high in the bulk phase. The neighborhood induces a dispersion on the second hyperpolarizability γ(2ω, ω, ω,0) values, but less than for the first hyperpolarizability.

Concerning the first hyperpolarizability β, linear expansions are valid for a weak and spatially-homogeneous electrostatic field, while the embedding electrostatic field created by the liquid surroundings is not. To compute individual molecular hyperpolarizability in the liquid phase, the molecular structure of the nearest neighbors have to be described explicitly, and QM/MM models are very relevant.

But long-range effects, up to several nanometers, can become heavy to take into account using QM/MM approaches, since a very large number of neighbors have to be included. This problem is expected to become all the more pregnant when the environment is anisotropic, with charges or polarization. We have proposed to consider the effect of long-range electrostatic environment on the liquid first hyperpolarizability using the average second hyperpolarizability γ. This long-term correction increases the precision of the hyperpolarizability calculation, and speeds up the convergence of the average value relative to the QM/MM environment size R_c. For pure bulk water, our approximation makes sense for a long-range part defined beyond 5 to 10 Å, depending on the quantity required. From a numerical point of view, the correction can be powerful for charged environments because the size of the explicit embedding region can be reduced drastically.

This work has demonstrated the promising potentialities of using the second hyperpolarizability, explicitely or implicitely to predict first hyperpolarizabilities. The present study is applied to pure water bulk phase at 300 K. However, the correction developed here should be also relevant for many other system geometries or compositions, where the long-range electrostatic effects can modify optical responses. If necessary, an extension including electric field gradients would also be straightforward using the Dalton package. Precision tests remain to be performed for more complex systems such as salted aqueous solutions, or charged interfaces.

Author contributions

GLB: conceptualization, methodology, validation, writing. OB: supervision, validation, writing. EB: supervision, writing, review, editing. CL: supervision, validation, writing.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank P.-F. Brevet for insightful discussions. We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources.

Notes and references

1. P.-F. Brevet, Surface Second Harmonic Generation, 1997.
2. M. C. Goh, J. M. Hicks, K. Kemnitz, G. R. Pinto, K. Bhattacharyya, K. B. Eisenthal and T. F. Heinz, J. Phys. Chem., 1988, 92, 5074–5075.
3. Q. Wei, D. Zhou and H. Bian, Phys. Chem. Chem. Phys., 2018, 20, 11758–11767.
4. M. N. Nasir, E. Benichou, C. Loison, I. Russier-Antoine, F. Besson and P.-F. Brevet, Phys. Chem. Chem. Phys., 2013, 15, 19919–19924.
5. G. Licari, J. S. Beckwith, S. Soleimanpour, S. Matile and E. Vauthey, Phys. Chem. Chem. Phys., 2018, 20, 9328–9336.
6. E. Donohue, S. Khorsand, G. Mercado, K. M. Varney, P. T. Wilder, W. Yu, A. D. MacKerell, P. Alexander, Q. N. Van, B. Moree, A. G. Stephen, D. J. Weber, J. Salafsky and F. McCormick, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 17290–17297.
7. M. E. Didier, O. B. Tarun, P. Jourdain, P. Magistretti and S. Roke, Nat. Commun., 2018, 9, 1–7.
8. S. Bancelin, C. Aimé, I. Gusachenko, L. Kowalczuk, G. Latour, T. Coradin and M. C. Schanne-Klein, Nat. Commun., 2014, 5, 1–8.
9. C. Liang, G. Tocci, D. M. Wilkins, A. Grisafi, S. Roke and M. Ceriotti, Phys. Rev. B, 2017, 96, 1–6.
10. P. Kaatz, E. A. Donley and D. P. Shelton, J. Chem. Phys., 1998, 108, 849–856.
11. J. Yang, E. V. Barnat, S.-k Im and D. B. Go, J. Phys. D: Appl. Phys., 2022, 55(22), 225203.
12. B. Levine and C. Bethea, J. Chem. Phys., 1976, 65, 2429–2438.
13. T. N. Ramos, S. Canuto and B. Champagne, J. Chem. Inf. Model., 2020, 60, 4817–4826.
14. B. F. Levine and C. G. Bethea, J. Chem. Phys., 1975, 63, 2666–2682.
15. E. C. Yan, Y. Liu and K. B. Eisenthal, J. Phys. Chem. B, 1998, 102, 6331–6336.
16. T. Joutsuka, T. Hirano, M. Sprik and A. Morita, Phys. Chem. Chem. Phys., 2018, 20, 3040–3053.
17. S. Pezzotti, D. R. Galimberti, Y. R. Shen and M.-P. Gaigeot, Phys. Chem. Chem. Phys., 2018, 20, 5190–5199.
18. L. Dalstein, K.-Y. Chiang and Y.-C. Wen, J. Phys. Chem. Lett., 2019, 10, 5200–5205.
19. A. Marchioro, M. Bischoff, C. Lütgebaucks, D. Biriukov, M. Předota and S. Roke, J. Phys. Chem. C, 2019, 123, 20393–20404.
20. R. Hartkamp, A. L. Biance, L. Fu, J. F. Dufrêche, O. Bonhomme and L. Joly, Curr. Opin. Colloid Interface Sci., 2018, 37, 101–114.
21. L. B. Dreier, C. Bernhard, G. Gonella, E. H. Backus and M. Bonn, J. Phys. Chem. Lett., 2018, 9, 5685–5691.
22. E. Ma, P. E. Ohno, J. Kim, Y. Liu, E. H. Lozier, T. F. Miller, H.-F. Wang and F. M. Geiger, J. Phys. Chem. Lett., 2021, 12, 5649–5659.
23. H. Chang, P. E. Ohno, Y. Liu, E. H. Lozier, N. Dalchand and F. M. Geiger, J. Phys. Chem. B, 2020, 124, 641–649.
24. Y. Foucaud, B. Siboulet, M. Duvail, A. Jonchere, O. Diat, R. Vuilleumier and J.-F. Dufrêche, Chem. Sci., 2021, 12, 15134–15142.
25. G. Maroulis, Chem. Phys. Lett., 1998, 289, 403–411.
26. J. Kongsted, A. Osted, K. V. Mikkelsen and O. Christiansen, J. Chem. Phys., 2004, 120, 3787–3798.
27. M. Hidalgo Cardenuto and B. Champagne, Phys. Chem. Chem. Phys., 2015, 17, 23634–23642.
28. T. Giovannini, M. Ambrosetti and C. Cappelli, Theor. Chem. Acc., 2018, 137(6), 74.
29. P. Beaujean and B. Champagne, J. Chem. Phys., 2019, 151, 064303.
30. T. N. Ramos, F. Castet and B. Champagne, J. Phys. Chem. B, 2021, 125, 3386–3397.
31. L. Jensen and P. T. Van Duijnen, J. Chem. Phys., 2005, 123(7), 074307.
32. G. Le Breton, O. Bonhomme, P.-F. Brevet, E. Benichou and C. Loison, Phys. Chem. Chem. Phys., 2021, 23, 24932–24941.
33. A. Willetts, J. E. Rice, D. M. Burland and D. P. Shelton, J. Chem. Phys., 1992, 97, 7590–7599.
34. O. Bonhomme, L. Sanchez, E. Benichou and P. Brevet, J. Phys. Chem. B, 2021, 125, 10876–10881.
35. A. Pardon, O. Bonhomme, C. Gaillard, P.-F. Brevet and E. Benichou, J. Mol. Liq., 2021, 322, 114976.
36. H. Kjaer, S. P. A. Sauer and J. Kongsted, J. Chem. Phys., 2011, 134, 44514.
37. Y. Luo, P. Norman, I. Introduction and H. Agren, J. Chem. Phys., 1998, 109, 3589–3595.
38. K. O. Sylvester-Hvid, K. V. Mikkelsen, P. Norman, D. Jonsson and H. Ågren, J. Phys. Chem. A, 2004, 108, 8961–8965.
39. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
40. J. L. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
41. G. J. Martyna, D. J. Tobias and M. L. Klein, J. Chem. Phys., 1994, 101, 4177–4189.
42. R. E. Isele-Holder, W. Mitchell and A. E. Ismail, J. Chem. Phys., 2012, 137.
43. R. E. Isele-Holder, W. Mitchell, J. R. Hammond, A. Kohlmeyer and A. E. Ismail, J. Chem. Theory Comput., 2013, 9, 5412–5420.
44. J. M. H. Olsen and J. Kongsted, Advances in Quantum Chemistry, Academic Press, 2011, vol. 61, pp. 107–143.
45. E. R. Kjellgren, J. M. Haugaard Olsen and J. Kongsted, J. Chem. Theory Comput., 2018, 14, 4309–4319.
46. K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O. Christiansen, R. Cimiraglia, S. Coriani and P. Dahle, et al. , Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 269–284.
47. H. C. Georg and S. Canuto, J. Phys. Chem. B, 2012, 116, 11247–11254.
48. N. H. List, H. J. A. Jensen and J. Kongsted, Phys. Chem. Chem. Phys., 2016, 18, 10070–10080.
49. M. T. Beerepoot, A. H. Steindal, N. H. List, J. Kongsted and J. M. H. Olsen, J. Chem. Theory Comput., 2016, 12, 1684–1695.
50. C. Steinmann, P. Reinholdt, M. S. Nørby, J. Kongsted and J. M. H. Olsen, Int. J. Quantum Chem., 2019, 119, e25717.
51. A. Marefat Khah, P. Reinholdt, J. M. H. Olsen, J. Kongsted and C. Hattig, J. Chem. Theory Comput., 2020, 16, 1373–1381.
52. T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2008, 393, 51–57.
53. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
54. P. Sałek, O. Vahtras, T. Helgaker and H. Ågren, J. Chem. Phys., 2002, 117, 9630–9645.
55. P. Besalú-Sala, S. P. Sitkiewicz, P. Salvador, E. Matito and J. M. Luis, Phys. Chem. Chem. Phys., 2020, 22, 11871–11880.
56. P. Beaujean and B. Champagne, Theor. Chem. Acc., 2018, 137, 50.
57. C. Liang, G. Tocci, D. M. Wilkins, A. Grisafi, S. Roke and M. Ceriotti, Phys. Rev. B, 2017, 96, 1–6.
58. P. Beaujean and B. Champagne, J. Chem. Phys., 2019, 151, 064303.
59. A. Osted, J. Kongsted, K. V. Mikkelsen, P. O. Åstrand and O. Christiansen, J. Chem. Phys., 2006, 124, 124503.
60. G. A. Cisneros, M. Karttunen, P. Ren and C. Sagui, Chem. Rev., 2014, 114, 779–814.
61. T. Vasilevskaya and W. Thiel, J. Chem. Theory Comput., 2016, 12, 3561–3570.
62. X. Pan, K. Nam, E. Epifanovsky, A. C. Simmonett, E. Rosta and Y. Shao, J. Chem. Phys., 2021, 154, 024115.
63. J. Nochebuena, S. Naseem-Khan and G. A. Cisneros, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2021, 11, e1515.

---

Footnotes

† Electronic supplementary information (ESI) available: Methodological details; symmetry, averages and standard deviations of water second hyperpolarizability in the liquid phase; heterogeneity of the electrostatic field generated by the neighborhood; impact of second hyperpolarizability fluctuations on the γ long correction. See DOI: https://doi.org/10.1039/d2cp00803c

‡ This code is deposited on Zenodo platform https://doi.org/10.5281/zenodo.5998193, and available on demand to the authors.

---