Marc Alías Rodríguez

Simone Bonfrate

Woojin Park

Cheol Nicolas Ferré

Ho Choi
email: cchoi@knu.ac.kr

Miquel Huix-Rotllant
email: miquel.huix-rotllant@cnrs.fr

Marc Alías-Rodríguez
email: marc.alias-rodriguez@univ-amu.fr

Cheol Ho Choi

Solvent Effects and pH Dependence of the X-ray Absorption Spectra of Proline from Electrostatic Embedding Quantum Mechanics/Molecular Mechanics and Mixed-Reference Spin-Flip Time-Dependent Density-Functional Theory

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

X-ray spectroscopy is gaining renewed interest due to the recent developments of new laser sources with ultrashort pulses like X-ray free electron lasers (XFELs). Such techniques can provide sufficient resolution even at subfemtosecond timescales, opening the way for developing new experimental techniques based on attosecond pulses. Among the X-ray spectroscopic techniques, X-ray absorption (XAS) is perhaps one of the most widely used and wellestablished techniques, providing at the same time electronic information from the Near-Edge X-ray Absorption Fine Structure (NEXAFS) and geometric information from the extended X-ray absorption fine-structure (EXAFS). 1 Xray spectroscopies are an element-specific technique that is particularly useful for following up on the dynamics of chemical reactions. [START_REF] Kim | Tracking reaction dynamics in solution by pump-probe X-ray absorption spectroscopy and X-ray liquidography (solution scattering)[END_REF] Indeed, time-resolved (TR) XAS spectroscopy provides a wealth of information on (photo)chemical re-activity in organic and organometallic complexes. [START_REF] Schmidt | Time-Resolved Soft X-Ray Absorption Spectroscopy of Molecules in the Gas and Liquid Phases[END_REF] X-ray absorption directly probes core-excited states by promoting core electrons to unoccupied virtual orbitals, following the usual dipole selection rules for absorption. In the K-edge spectra of most common organic molecules, the 1s electron can be used to probe transitions to π * , σ * , or Rydberg orbitals. These states can be affected by the electrostatic interaction with the solvent as well as the structural fluctuations due to the temperature. Solvent effects in the XAS spectra are well-documented experimentally for many chemical systems, [START_REF] Schmidt | Time-Resolved Soft X-Ray Absorption Spectroscopy of Molecules in the Gas and Liquid Phases[END_REF][START_REF] Penfold | Simulations of X-ray absorption spectra: the effect of the solvent[END_REF][START_REF] Messer | Local Hydration Environments of Amino Acids and Dipeptides Studied by X-ray Spectroscopy of Liquid Microjets[END_REF][START_REF] Messer | pH Dependence of the Electronic Structure of Glycine[END_REF][START_REF] Zubavichus | Innershell Absorption Spectroscopy of Amino Acids at All Relevant Absorption Edges[END_REF] but only a few theoretical protocols exist in the literature capable of efficiently simulating solvent effects on XAS. [START_REF] Penfold | Simulations of X-ray absorption spectra: the effect of the solvent[END_REF][START_REF] Messer | pH Dependence of the Electronic Structure of Glycine[END_REF][START_REF] Cappa | Nature of the Aqueous Hydroxide Ion Probed by X-ray Absorption Spectroscopy[END_REF][START_REF] Weinhardt | Probing hydrogen bonding orbitals: resonant inelastic soft X-ray scattering of aqueous NH3[END_REF][START_REF] Folkestad | Equation-of-Motion MLCCSD and CCSD-in-HF Oscillator Strengths and Their Application to Core Excitations[END_REF][START_REF] Loh | Observation of the fastest chemical processes in the radiolysis of water[END_REF][START_REF] Reinholdt | Nitrogen K-Edge X-ray Absorption Spectra of Ammonium and Ammonia in Water Solution: Assessing the Performance of Polarizable Embedding Coupled Cluster Methods[END_REF][START_REF] Tsuru | Solvent Effects in the Ultraviolet and X-ray Absorption Spectra of Pyridazine in Aqueous Solution[END_REF][START_REF] De Santis | Environment Effects on X-Ray Absorption Spectra With Quantum Embedded Real-Time Time-Dependent Density Functional Theory Approaches[END_REF] Here we develop a new model for describing efficiently and accurately the solvent effects on the X-ray spectra. Our model is based on a combination of the electrostatic embedding QM/MM method, used to equilibrate a water droplet around the molecule of interest, and mixed-reference spin-flip time-dependent density functional theory (MRSF-TDDFT) coupled to QM/MM to extract the X-ray spectroscopy. MRSF-TDDFT has been recently shown to be able to efficiently and accurately extract XAS, with errors that are within 1 eV of the experimental spectra. [START_REF] Park | Mixed-Reference Spin-Flip Time-Dependent Density Functional Theory for Accurate X-ray Absorption Spectroscopy[END_REF] In addition, an MRSF-TDDFT/MM coupling was shown to lead to accurate excited-state non-adiabatic dynamics of thymine in water. [START_REF] Huix-Rotllant | Photochemistry of Thymine in Solution and DNA Revealed by an Electrostatic Embedding QM/MM Combined with Mixed-Reference Spin-Flip TDDFT[END_REF] As a test case, we compute the XAS spectrum of the aminoacid proline in water. Proline is the only standard amino acid with a secondary amino group. The interaction of zwitterionic proline with water is well documented, indicating that the carboxylate group accepts less than two hydrogen bonds and the amine hydrogen atoms donate less than one hydrogen bond to the surrounding water molecules. [START_REF] Mclain | Structure and Hydration of L-Proline in Aqueous Solutions[END_REF] Proline has a rich XAS structure for the oxygen, nitrogen, and carbon K-edges, [START_REF] Zubavichus | Innershell Absorption Spectroscopy of Amino Acids at All Relevant Absorption Edges[END_REF][START_REF] Plekan | An X-ray absorption study of glycine, methionine and proline[END_REF] and the corresponding spectra depend on the pH. [START_REF] Messer | pH Dependence of the Electronic Structure of Glycine[END_REF][START_REF] Meyer | Resonant Inelastic Soft X-ray Scattering and X-ray Emission Spectroscopy of Solid Proline and Proline Solutions[END_REF] Indeed, at neutral pH, proline exists in a zwitterionic form, in which the carboxylic group is deprotonated and the amino group is protonated. At very acidic pH (<2), the carboxylic group is protonated, the proline being in a cationic form, whereas at very basic pH (>11) the amino group is deprotonated and the proline is in anionic form.

This study is organized as follows: first, we describe the MRSF-TDDFT methodology for X-ray spectroscopy and its electrostatic embedding coupling to the solvent via electrostatic potential fitted (ESPF) charges; second, we describe a theoretical protocol to simulate the solvent effects in X-ray spectroscopy via the combination of QM/MM molecular dynamics and MRSF-TDDFT/MM method; finally, we describe the solvent effects on proline by comparing the gas phase spectra, the spectrum generated from molecular dynamics snapshots with and without electrostatics for the three protonation forms of proline.

Methods

Protocol for X-ray spectra simulations in solution

The simulation of solvent effects on X-ray spectra implies the description of the solutesolvent interactions as well as the inhomogeneous broadening due to the fluctuations of the solute in the solvent that acts as a thermal bath. In Fig. 1, we summarize the main steps of our protocol. This protocol makes use of QM/MM partitioning, describing accurately only transition densities that fit well within the spatial extension of the QM fragment and the employed basis set.

Step 0: Model preparation. The preliminary step is the construction of a model system containing the solute and the solvent. To this aim, the gas-phase ground state optimized geometry of the neutralized solute is placed in the center of a box of pure solvent, the temperature and the pressure of which have been equilibrated at MM level by two following dynamics propagated respectively in the NVT and NPT ensembles using the Berendsen thermostat and barostat. In both cases, Periodic Boundary Conditions (PBC) have been employed. A droplet of suitable radius, defined by testing the dependence of the XAS spectra on droplet size, is then extracted from the box.

Step 1: Equilibration and production. The droplet structure is firstly minimized at the QM/MM level, where the solute is described at the QM level while the solvent is described classically. Spherical boundary conditions are imposed on all the simulations on the droplet system to prevent diffusion of the outer shell of the solvent. The QM/MM electrostatic embedding method is based on electrostatic potential fitted (ESPF) operators, as briefly described in Sec. 2.2.

Once the system has been relaxed, a two-step equilibration is done. The droplet is firstly equilibrated in the NVT ensemble with QM/MM molecular dynamics (MD) in which the Berendsen thermostat is used until it reaches the desired temperature. Second, the equilibration is further refined using the Bussi thermostat. An NVT production run using the Bussi thermostat is then performed. If available, a comparison with experimental radial distribution functions is a good test of the convergence of the trajectories. The length of the production run will depend on the desired number of independent snapshots that would be extracted from the trajectory, which is imposed by the autocorrelation function.

Step 2: X-ray spectrum in solution. The lowest states involving core excitations are calculated at ∆MRSF/MM level on the different snapshots extracted from MD. The ∆MRSF/MM method is briefly described in Sec. 2.3. For such calculations, a mixed uncontracted atomic basis set is employed for converging the density at the nuclear region of the core of interest, and a contracted basis set that can describe correctly the valence states for the rest of atoms. The accuracy of ∆MRSF depends on the choice of exchange-correlation functional and basis set, which can be validated by comparing to experimental gas phase spectra. Finally, the XAS spectra are simulated assigning an established Gaussian broadening to each X-ray state.

Solute-solvent interactions via electrostatic embedding QM/MM

The MRSF-TDDFT has been recently coupled to an MM electrostatic environment via an electrostatic embedding QM/MM approach based on electrostatic potential fitted (ESPF) charges. [START_REF] Huix-Rotllant | Photochemistry of Thymine in Solution and DNA Revealed by an Electrostatic Embedding QM/MM Combined with Mixed-Reference Spin-Flip TDDFT[END_REF] The main advantage of ESPF electrostatic embedding is the reduced number of integrals, which only depends on the number of QM atoms and the numerical grid for computing the ESPF atomic charges. [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF] Here, we summarize the main features of the model. The total energy of the QM/MM system can be decomposed as a sum of three terms,

E = E QM + E M M + E int , (1)
where E QM is the energy of the QM subsystem in the gas phase, E M M is the energy of the MM subsystem plus the van der Waals interactions between QM and MM subsystems and E int consist in the energy of the electrostatic interactions between the two subsystems. Here, we restrict our model to non-polarizable MM force fields. In the ESPF approach, it is convenient to define the QM/MM interaction energy in terms of Coulombic potential as

E int = N QM A=1 q A Φ M M A . (2)
The Coulombic potential felt by a QM atom A at position r A is given by,

Φ M M A = N M M j=1 q j |r Aj | . (3)
in which r Aj = r A -r j . In the equation, r j and q j are respectively the coordinates and the charge of MM atom j, and N QM and N M M are respectively the number of QM and MM atoms.

The ESPF partial charges at the QM centers, q A = Z A -Q A , are obtained as the difference between the nuclear charge Z A , and the atomic charge population Q A = µν P µν Q A,µν which in turn is derived by tracing the one-electron density matrix P with the atomic charge operator matrix Q A . The elements of the charge operator matrix are obtained by solving a linear set of equations to fit the electrostatic integrals computed on the Lebedev atom-centered grid constructed around the QM subsystem, [START_REF] Ferré | Approximate electrostatic interaction operator for QM/MM calculations[END_REF][START_REF] Huix-Rotllant | Analytic Energy, Gradient, and Hessian of Electrostatic Embedding QM/MM Based on Electrostatic Potential-Fitted Atomic Charges Scaling Linearly with the MM Subsystem Size[END_REF]

N QM A Q A,µν |r kA | = V k,µν . (4)
In the equation, the electrostatic integrals on a grid are defined as

V k,µν = d 3 rχ * µ (r) 1 |r -r k | χ ν (r) , (5)
where r k are the coordinates of the grid points and χ are the atomic orbitals. A correction is applied to the charge operator to enforce the conservation of the total charge of the QM subsystem. [START_REF] Huix-Rotllant | Analytic Energy, Gradient, and Hessian of Electrostatic Embedding QM/MM Based on Electrostatic Potential-Fitted Atomic Charges Scaling Linearly with the MM Subsystem Size[END_REF] Deriving the interaction energy (Eq.

2) with respect to the density matrix P, we can obtain the QM/MM interaction in the form of an interaction operator, simply given by

h int µν = - N M M A Q A,µν Φ M M A , (6)
that is added as a one-electron Hamiltonian contribution to the gas phase Fock operator.

For the calculations of ∆MRSF, this accounts for adding h int µν in the restricted open-shell Kohn-Sham (ROKS) Fock operator, with no extra terms to be added in the response equations. [START_REF] Huix-Rotllant | Photochemistry of Thymine in Solution and DNA Revealed by an Electrostatic Embedding QM/MM Combined with Mixed-Reference Spin-Flip TDDFT[END_REF] It is further possible to compute the gradient of the total QM/MM energy (Eq. 1) as

E x = E x QM + E x M M + E int,x .
Hereafter, the superscript E x = ∂E/∂x refers to the derivative with respect to an x atom coordinate. The analytical derivative E x QM depends on the ab initio method employed, but it is generally available for the most common electronic structure methods. Likewise, analytical expressions for E x M M are available. Here we focus on the derivatives of the interaction energy with respect to QM and MM atom coordinates. [START_REF] Huix-Rotllant | Analytic Energy, Gradient, and Hessian of Electrostatic Embedding QM/MM Based on Electrostatic Potential-Fitted Atomic Charges Scaling Linearly with the MM Subsystem Size[END_REF] The derivative with respect to MM atom coordinates is simply given by

E int,x i = N QM A=1 q A Φ M M,x i A (7)
where the derivative of the Coulombic potential with respect to an MM atom is

Φ M M,x i A = -q i x Ai |r Ai | 3 . ( 8
)
The derivative of the interaction energy with respect to the QM atom is given by

E int,x B = N QM A=1 q A Φ M M,x B A -Q (x B ) A Φ M M A , (9)
where Φ M M A is the derivative of the Coulombic potential (Eq. 3) with respect to the same QM coordinate, which is given by

Φ M M,x B A = - N M M j=1 q j x Aj |r Aj | 3 δ AB . ( 10
)
On the other hand, Q

(x B ) A
is the derivative of the ESPF charge with respect to the QM coordinate x B . The parenthesis in the derivatives indicates that the term corresponding to the derivative of the density matrix is excluded (see Ref. 22 for further details). The construction of the ESPF charge operator derivative is straightforwardly obtained by solving a set of linear equations

N QM A Q x A A,µν |r kA | = V x A k,µν - N QM A Q A,µν 1 |r kA | x A . (11)

X-ray spectrum from mixedreference spin-flip timedependent density-functional theory

MRSF-TDDFT [START_REF] Lee | Eliminating spin-contamination of spinflip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix[END_REF][START_REF] Lee | Efficient implementations of analytic energy gradient for mixedreference spin-flip time-dependent density functional theory (MRSF-TDDFT)[END_REF] is a promising new variant of linear response (LR) theory, which overcomes the limitations of LR-as well as Spin-Flip (SF)-TDDFT. [START_REF] Shao | The spin-flip approach within timedependent density functional theory: Theory and applications to diradicals[END_REF] The essence of the MRSF-TDDFT is in the utilization of two references, which naturally expands the response space resolving issues such as spin-contamination of SF-TDDFT. [START_REF] Lee | Time-Dependent Density Functional Theory[END_REF] For its realization, a new mixed reference state is introduced that has an equiensemble of M S = +1 and M S = -1 com-ponents of a triplet state, e.g.,

ρ MR 0 (x) = 1 2 ρ M S =+1 0 (x) + ρ M S =-1 0 (x) , (12)
where the idempotency of ρ MR 0 (x) was recovered by the concept of hypothetical single reference. [START_REF] Lee | Eliminating spin-contamination of spinflip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix[END_REF] The responses from the two components allow exact external spin adaptations, which yields completely decoupled linear-response equations for singlet and triplet excited states as, [START_REF] Lee | Eliminating spin-contamination of spinflip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix[END_REF] rs

A (k)(0) pq,rs + A ′(k) pq,rs X (k) rs = Ω (k) X (k) pq , k = S, T (13)
where k = S, T labels singlet and triplet states,

A (k)(0)
pq,rs is an orbital Hessian matrix derived by the linear response, and

A ′(k)
pq,rs is a coupling matrix between configurations originating from different components, M S = +1 and M S = -1, of the mixed reference. [START_REF] Lee | Eliminating spin-contamination of spinflip time dependent density functional theory within linear response formalism by the use of zeroth-order mixed-reference (MR) reduced density matrix[END_REF][START_REF] Lee | Efficient implementations of analytic energy gradient for mixedreference spin-flip time-dependent density functional theory (MRSF-TDDFT)[END_REF]

X (k)
pq and Ω (k) are the amplitude vectors and the excitation energies with respect to the reference state, respectively. As a result, not only the spincontamination but also the ambiguity with the cumbersome identification of the resulting response states of SF-TDDFT are eliminated.

An important advantage of MRSF-TDDFT over the usual linear-response TDDFT [START_REF] Casida | ACS Symposium Series[END_REF][START_REF] Fiolhais | A primer in density functional theory[END_REF][START_REF] Casida | Progress in time-dependent density-functional theory[END_REF] is that MRSF-TDDFT enables the proper computation of the S 1 /S 0 conical intersections, as both states are treated at the same level of theory in MRSF-TDDFT. [START_REF] Lee | Conical Intersections in Organic Molecules: Benchmarking Mixed-Reference Spin-Flip Time-Dependent DFT (MRSF-TD-DFT) vs Spin-Flip TD-DFT[END_REF] Additionally, MRSF-TDDFT not only includes singly but also some important doubly excited configurations that are missing in the LR-TDDFT, [START_REF] Park | Internal Conversion between Bright (1 1 B + u ) and Dark (2 1 Ag ) States in s-trans-Butadiene and s-trans-Hexatriene[END_REF] significantly improving its prediction accuracy in the study of diradicals 32 and multi-configurational states. [START_REF] Horbatenko | How Beneficial Is the Explicit Account of Doubly-Excited Configurations in Linear Response Theory?[END_REF][START_REF] Park | Impact of the Dynamic Electron Correlation on the Unusually Long Excited-State Lifetime of Thymine[END_REF] The particular formulation of MRSF-TDDFT combined with maximum-overlap method (MOM) also provides a significant advantage in the computation of X-ray absorption. [START_REF] Park | Mixed-Reference Spin-Flip Time-Dependent Density Functional Theory for Accurate X-ray Absorption Spectroscopy[END_REF][START_REF] Gilbert | Self-Consistent Field Calculations of Excited States Using the Maximum Overlap Method (MOM)[END_REF] Although LR-TDDFT [START_REF] Casida | Progress in time-dependent density-functional theory[END_REF] has been the most popular method for valence excited state studies, the accuracy of TDDFT to predict experimental Xray states is very crude, with errors as large as ∼10 eV for second-period elements like C, N, O, and F, [START_REF] Besley | Timedependent density functional theory calculations of the spectroscopy of core electrons[END_REF] and even larger errors for heavier elements. These large shifts are attributed to the lack of description of core-hole orbital relaxation effects in the approximate exchangecorrelation functionals used in TDDFT, which cannot account for the drastic change in the Coulomb interaction of the remaining electrons when a core electron is excited. Properly tuned self-interaction corrected range-separated exchange functionals can remedy the large errors in core-excitation energies by DFT methods. [START_REF] Besley | Density Functional Theory Based Methods for the Calculation of Xray Spectroscopy[END_REF] However, a simple remedy to orbital relaxation is not available in the context of the TDDFT approach. These major challenges can be easily overcome by MRSF-TDDFT with the help of its high-spin triplet reference, which allows a core-hole particle (CHP) pair. The particular ∆CHP-MRSF flavor (referred to for simplicity as ∆MRSF in the following) was designed in our previous work, [START_REF] Park | Mixed-Reference Spin-Flip Time-Dependent Density Functional Theory for Accurate X-ray Absorption Spectroscopy[END_REF] that we briefly summarize here. In the ∆MRSF approach, the excitation energy is calculated as the energy difference between the X-ray state computed with MRSF-TDDFT using a core-hole triplet state as the reference and the singlet ground state computed at the restricted Kohn-Sham level. The orbital relaxation of the CHP pair configuration can be accomplished by replacing the two singly occupied open (O1 and O2) orbitals of restricted open-shell Hartree-Fock (ROHF), or ROKS applying DFT, reference with the particular core (the 1s orbital for XAS) and the virtual orbitals with MOM. [START_REF] Gilbert | Self-Consistent Field Calculations of Excited States Using the Maximum Overlap Method (MOM)[END_REF] In that way, the orbital relaxation due to the excitation of a core electron is embedded in the orbital energies, leading to a larger accuracy of the spectrum. Finally, to simplify the excitation space, a restricted excitation window (REW) is also applied to all calculations. This restricts MRSF-TDDFT excitation space from unwanted occupied orbitals. This simple strategy exhibited near-perfect predictions with RMSE ∼ 0.5 eV, featuring a median value of 0.3 and an interquartile range of 0.4 for Kedge excitation energies of 24 second-and 17 third-row molecules. [START_REF] Park | Mixed-Reference Spin-Flip Time-Dependent Density Functional Theory for Accurate X-ray Absorption Spectroscopy[END_REF]

Computational details

The model for proline in the gas phase has been constructed by performing minimum energy geometry optimization of the ground state. The geometries without counterions are optimized in GAMESS-US at BH&HLYP/6-31G* level. For the zwitterionic case, the N-H distance has been kept frozen during the process in order to avoid a proton transfer to the carboxylic group. The geometries including Na + and Cl -counterions are fully relaxed without constraints.

For including the solvent effects, the gas phase proline structures have been soaked in a pre-equilibrated TIP3P [START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF] water cubic box with an edge size of 50 Å. From such structures, we proceed to extract spherical droplets with a radius of 20 Å centered at C α of proline. These systems have been first relaxed with a QM/MM minimization at the HF/3-21G//Amber99 level, [START_REF] Wang | How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules?[END_REF] where the QM sub-region includes proline and counterions, whereas the MM region contains all the water molecules. A spherical solvent boundary has been applied to prevent the diffusion of the outer water molecules during the simulations. [START_REF] Kratky | Computational techniques for spherical boundary conditions[END_REF] After minimization, the systems were equilibrated at QM/MM MD level using the same conditions as the minimization (HF/3-21G//Amber99 in which the QM sub-region includes proline and counterions). A first equilibration was run for 150 ps in the canonical ensemble using the Berendsen thermostat with a target temperature of 298 K. The QM/MM MD equilibration was further refined at the same temperature for another 50 ps using the Bussi thermostat with spherical solvent boundary (see Fig. S1 in the supporting information.) Finally, we proceeded to the production run for 150 ps using the Bussi thermostat. A set of 50 independent snapshots has been extracted from each production trajectory. The statistical significance of the sampling has been tested by means of Z-test. The analysis of the p-values (reported in Tab. S1 in the supporting information) shows that in all cases but one (i.e. N -C 1 distribution in the anion trajectory) the p-values are greater than 0.05, indicating that the snapshot set constitutes a good sam-ple of the full trajectory in the most relevant proline distances.

For all the MD simulations we used the following parameters. The equations of motion have been integrated using the Beeman integrator with an integration time step of 1 fs. The coupling time (τ T ) of the thermostat was set to 0.1 ps for the Berendsen thermostat and 0.2 ps for the Bussi thermostat. An additional rescaling of the velocities has been applied at each integration step in order to enforce the translational and rotational invariance of the center of mass of the system.

Minimization and molecular dynamics were run with a development version of GAMESS-US, version 2021 (R2 Patch2), [START_REF] Barca | Recent developments in the general atomic and molecular electronic structure system[END_REF] coupled to a modified version of Tinker 8.10.1. [START_REF] Ponder | Software tools for molecular design[END_REF] The image of the droplet model was generated using the UCSF ChimeraX molecular visualization tool. [START_REF] Pettersen | UCSF ChimeraX: Structure visualization for researchers, educators, and developers[END_REF] The K-edge spectra have been simulated using ∆MRSF as implemented in GAMESS-US [START_REF] Barca | Recent developments in the general atomic and molecular electronic structure system[END_REF] with the scalar relativistic effects through second-order Douglas-Kroll-Hess Hamiltonian. [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF][START_REF] Hess | Relativistic electronicstructure calculations employing a two-component no-pair formalism with external-field projection operators[END_REF] The DFT exchange-correlation functional employed for all calculations is BH&HLYP. [START_REF] Becke | A new mixing of Hartree-Fock and local density-functional theories[END_REF] The basis set employed has been aug-pc-X2 for the core of interest and aug-pc-1 for the other QM centers. [START_REF] Jensen | Polarization consistent basis sets: Principles[END_REF][START_REF] Ambroise | Probing Basis Set Requirements for Calculating Core Ionization and Core Excitation Spectroscopy by the ∆ Self-Consistent-Field Approach[END_REF] For converging the core-hole triplet reference, the standard orbital rotation restriction based on the maximum overlap method has been used. From this reference, the lowest 14 X-ray states were extracted from linear-response MRSF-TDDFT. The overall spectra were obtained as a sum of Gaussian functions for each transition, whose position was determined by the energy at the ∆MRSF level, the height corresponded to the transition dipole moment and the width was fixed at 0.5 eV.

Results & Discussion

The X-ray absorption spectra for the carbon, nitrogen, and oxygen K-edge have been computed for a proline model for the ionic forms (see Fig. 2) in the gas phase and in the water solution using the methodology protocol de-scribed in the previous sections. The gas phase spectra of the neutral form of proline for all Kedges are found in the supporting information, showing an excellent agreement between the theoretical and experimental spectra. The calculations have been performed with and without counterions. Here, we only show the spectra of proline with counterions, in which the right amount of Na + and Cl -ions have been added so that the total net charge is zero (for the effect of the counterions in the XAS spectrum of proline, see Figures S12,S13 and S14 of the supporting information).

Proline can exist in three protonation forms, namely, cationic, anionic, and zwitterionic forms depending on the pH (see Fig. 2). The neutral form is only stable in the gas phase when there are no counterions present, and will not be considered in the current study. The focus of the following sections is on analyzing the proline-water interactions and describing the X-ray transitions. The energetic shifts due to the solvent for each of the protonation forms of proline are discussed independently for each K-edge, compared to the gas phase, and to the snapshot average deactivating the QM and MM electrostatic interactions (labeled "no electrostatics"), allowing us to determine the energetic shift due to the interaction with the solvent, and due to the fluctuations of proline. Finally, we reconstruct the pH-dependent spectra of proline.

Proline models in water

In order to understand better the proline-water interactions in the droplet models, we extracted the radial distribution function (RDF) from the production trajectories. In particular, we focused our attention on the interactions between water molecules and the amino and carboxylic groups of proline in the three different protonation states. The RDF plots obtained for the zwitterionic structure (in orange in Fig. 3) are in good agreement with the experimental RDF plots reported by McLain et al. as both the position of the main peaks and their relative intensity match well. [START_REF] Mclain | Structure and Hydration of L-Proline in Aqueous Solutions[END_REF] The water-amino group RDF plots (Fig. 3, The black lines correspond to the experimental spectra of the zwitterionic form, registered in a 1:10 proline/water solution, reported in Ref. 18. The dashed line corresponds to the carboxylic hydrogen and OW in the cationic form. a), instead, show that in the cationic and zwitterionic form, which present a peak of similar intensity around 2 Å, the hydrogens of the amine group act as hydrogen bond donors toward the oxygens of water. In contrast, the RDF of the anionic structure does not present such a peak suggesting that, in this protonation form, the only hydrogen of the amino group establishes preferentially an intramolecular hydrogen bond with one of the oxygens in the carboxylic group (see Fig. S7 of the supporting information). This hypothesis is further supported by the calculation of the average distance (r O 2 -HN = 3.28 ± 0.25 Å) between the two oxygens and the hydrogen atom.

In the case of water-carboxylic group RDF plots (Fig. 3,b), the presence of a peak around 1.75 Å for each protonation form, indicating clearly the formation of hydrogen bonds between water and the two oxygen of the carboxylic group (the latter acts as hydrogen bond acceptor.) The distances at which the peaks occur suggest also a correlation in all the protonation forms, although the lower intensity of the peaks for the cation model implies that the carboxylic group is accepting fewer hydrogen bonds than the other forms. The presence of a peak at a comparable distance in the RDF plot of the water oxygen and hydrogen of the carboxylic group in the cationic structure indicates that the latter acts also as a hydrogen bond donor.

As part of the geometrical analysis of the trajectories, we report also the behavior of the counterions along the trajectories that contain them, since it can affect the XAS spectrum, especially for the σ * peaks and their mixture with Rydberg-type orbitals (see Fig. S2 of the supporting information.) The sodium ion in both anionic and zwitterionic models exhibits interaction with the oxygens of the carboxylic group although, for the anion, the sodium ion interacts preferentially with one of the oxygens (r Na-O 2 = 2.17 ± 0.08 Å and r Na-O 1 = 4.14 ± 0.12 Å). In the zwitterion model instead, the interaction is equally shared between the two oxygens (r Na-O 2 = 2.40 ± 0.21 Å and r Na-O 1 = 2.33 ± 0.16 Å). The chloride ion does not form a strong interaction with the amine moiety of proline. In the cation model, in fact, the distribution of the distances between the counterion and the nitrogen atom, which is non-gaussian, has a median value of rCl-N = 11.51 ± 3.61 Å, indicating a weak electrostatic interaction with proline. Similarly, in the zwitterion case, the chloride-amine distance has a value of r Cl-N = 6.70 ± 0.33 Å exhibiting thus a larger interaction with the sodium ion. and(c) cationic in solution (solid line), gas phase (dashed), and no electrostatics (dotted). The gas phase XAS spectrum comes along with its analysis (shaded area) and corresponds to 1s → σ * N -H (blue) and 1s → σ * C-N (red). The spectra for the solvent model converged up to 406.96, 404.21 and 407.51 eV for the zwitterionic, anionic, and cationic structures, respectively.

N K-edge

The gas phase, solvent, and no electrostatics N K-edge spectrum for the three protonated structures are shown in Figure 4. Since there is only one nitrogen core, two types of electronic transitions are observed either to σ * N -H localized along one of the N-H bonds or to σ * N -C orbitals localized over both N-C 2 and N-C 5 bond (see Fig. S9 in the supporting information). In the gas phase, the spectrum of the zwitterionic structure shows a weakly intense peak at 404.1 eV corresponding to the σ * N -H transition and a broader more intense peak at 406.4 eV corresponding to the σ * N -C . A similar structure is observed for the cationic form, slightly blueshifted with respect to the zwitterionic form. Indeed, the σ * N -H peak is found at 404.4 eV, whereas the σ * N -C is located at 406.6 eV. Finally, the anionic spectrum is less intense and redshifted with respect to the zwitterionic or cationic forms. The first peak is located at 401.9 eV, 0.3 eV above the σ * N -H , whereas the σ * N -C is found at 403.2 eV. The N K-edge spectra of proline in water can be largely affected by the changes in the σ * orbitals, which depend directly on the N-C (see Fig. S4 in the supporting information.) and N-H distances (see Fig. S5 in the supporting information.) The zwitterion and cation distances exhibit similar average distances. The average N-C bond length (r N-C = 1.54 ± 0.04 Å) is 0.05 Å longer than in the gas phase. On the contrary, the average N-H bond distance in solution (r N-H = 1.03 ± 0.03 Å) is 0.07 Å shorter than in the gas phase. For the anionic form, the gas phase and solvent average distances are similar. The average N-C distance in solution (r N-C = 1.49 ± 0.03 Å) is 0.03 Å longer than in the gas phase, while the N-H distance in solution (r N-H = 1.01 ± 0.02 Å) is 0.01 Å longer.

The N K-edge spectra in solution show a strong impact both on the main peak corresponding to the 1s → σ * N -C transition and the σ * N -H . Indeed, for the zwitterionic and cationic forms, the σ * N -H intensity is largely decreased with respect to the gas phase, which only appears as a tail of the σ * N -C transition. In the anionic form instead, the first peak appears at the same position as in the gas phase. In the zwitterionic and cationic forms, a double peak structure is observed for the σ * N -C transition in both cases. The full width at half maximum of these peaks has a broader distribution of 2.2 eV for the zwitterion and 2.4 eV for the cation. This distribution is explained by the structural fluctuation of the N-C 2 and N-C 5 distances (see Fig. S4 in the supporting information), which largely affect the energy of the σ * N -C orbitals. For the zwitterion, the peak maximum is redshifted around 0.4 eV with respect to the gas phase, whereas in the cationic form, the peak maximum is blueshifted 0.4 eV. The anionic form, on the contrary, has a peak structure and width similar to the gas phase. The σ * N -C transition is blueshifted more than 1.0 eV.

Comparing the spectrum in solution with the spectrum without electrostatic interactions, we can unravel the effect of interactions of proline with water on the σ * transitions. For the anionic form, the spectra with and without electrostatics are identical, indicating that the effect of the solvent on neutral nitrogen is negligible. For the zwitterionic and cationic forms, the electrostatic interaction modulates slightly the intensity of the first σ * N -C peak, without affecting its excitation energy. For the second peak, on the contrary, the transition is redshifted by the solvent in both cases. and(c) cationic in solution (solid line), gas phase (dashed), and no electrostatics (dotted). The gas phase XAS spectrum comes along with its analysis (shaded area) and corresponds to oxygen 1 (blue) and 2 (red) according to the nomenclature of Fig. 2. The spectra for the solvent model converged up to 536.51, 536.61, and 537.58 eV for the zwitterionic, anionic and cationic structures, respectively.

O K-edge

The gas phase, solvent, and no electrostatics O K-edge spectrum for the three protonated structures are shown in Figure 5. The total gas phase spectra are decomposed in the K-edge of O 1 and O 2 , according to the nomenclature in Figure 2. For the equivalent decomposition of the solvent spectra, see Fig. S15 of the supporting information. The three forms show a similar spectrum, with a first main peak corresponding to the 1s(O)→ π * transition between 532 and 534 eVs. In the case, of the anion and zwitterion, the π * orbital is delocalized over the carboxylate group (see Fig. S10 supporting information). In the anionic forms, the O 1 and O 2 transitions are overlapping due to the symmetry between the two oxygens, leading to a maximum peak at around 532.9 eV. In the case of zwitterionic form, the 1s(O 2 )→ π * transition is slightly blueshifted with respect to O 1 due to the intermolecular hydrogen bond that is formed with the amine group. The absorption maximum is found at 533.2 eV. In the cationic form, the transition from the core-orbital of O 1 is to the σ * O-H , shifted by more than 3.4 eV due to the protonation of the oxygen. For the lower part of the spectrum, the absorption maximum is found at 532.7 eV. Above 534 eV, the spectrum corresponds to transitions to σ * or the lone pairs of oxygen, which are mixed with the Rydberg character. These transitions are less intense and structured than the π * transitions. In the case of the cationic form, they span a larger range of energies due to the shift caused by the protonation of O 1 .

The O K-edge spectra of proline in water can be largely affected by the changes both in the π * and the σ * orbitals, which depend on the fluctuation of C-O distances for all forms. The cationic form can also be affected by the O-H distance. The torsion of the carboxylic group could also affect the energetics of π * , but this motion is not observed in our dynamics (see Fig. S8 in the supporting information). From the analysis of the C-O distance along the trajectories, it is possible to observe that there is an almost perfect delocalization of the double bond across the two oxygens in the zwitterion model (r C-O = 1.26 ± 0.02 Å), 0.01 Å larger to the gas phase. The anion model presents a similar behavior although, in such a case, one bond is slightly longer than the other (r C-O 1 = 1.25 ± 0.02 Å and r C-O 2 = 1.27 ± 0.02 Å), close to the gas phase value of 1.26 Å. Instead, as expected in the cation model, one can distinguish clearly a C=O double bond (r C-O = 1.21±0.02 Å) and a single bond (r C-OH = 1.34 ± 0.03 Å). The double bond is 0.01 Å shorter in the gas phase, whereas the single bond is 0.02 Å shorter.

The O K-edge spectra for the three protonated structures in solution are shown in Fig. 5. In the zwitterionic structure, the main π * transition is slightly redshifted (around 0.2 eV) and its intensity is slightly reduced in solution compared to the gas phase one. A similar trend is observed for the lowest π * peak in the cationic form, in which energy is redshifted by around 0.4 eV and a similar decrease of intensity. The σ * O-H transition from O 1 core of the cationic form appearing at 536 eV is smeared out, overlapping with other σ * transitions. These transitions fall however at the water window region, [START_REF] Spielmann | Generation of Coherent Xrays in the Water Window Using 5-Femtosecond Laser Pulses[END_REF][START_REF] De Stasio | Feasibility tests of transmission X-ray photoelectron emission microscopy of wet samples[END_REF] making it impossible to record the proline spectra at such high energies. Finally, the anionic spectrum appears blueshifted with respect to the gas phase, with a maximum peak shift of 0.1 eV.

From the comparison between the gas phase, the solvent spectra, and the spectrum deactivating electrostatic interactions, we can see a cancellation effect between the structural fluctuations and the solvent electrostatic interactions. On the one hand, the structural fluctuations tend to redshift by as much as 1 eV the π * peaks, which are broadened and less intense than in the gas phase. This is probably due to the formation of hydrogen bonds with water. On the other hand, the electrostatic interaction blueshifts the π * due to an unfavorable electrostatic interaction, which brings back the π * states at almost the same excitation energy as in the gas phase.

C K-edge

The gas phase, solvent, and no electrostatics C K-edge spectrum for the three protonated structures are shown in Figure 6. The three forms show a similar spectrum, due to the fact that the protonation or deprotonation does not and(c) cationic in solution (solid line), gas phase (dashed), and no electrostatics (dotted). The gas phase XAS spectrum comes along with its analysis (shaded area) and corresponds to carbon 1 (blue) and carbon 2 to 5 (red) according to the nomenclature of Fig. 2. The spectra for the solvent model converged up to 290.16, 290.49 and 290.10 eV for the zwitterionic, anionic, and cationic structures, respectively. affect strongly the electronic structure of the carbons. The spectra have essentially two main contributions (see Fig. S11 for the supporting information for the molecular orbital pictures). A first main peak is found between 288 and 290 eV, and it corresponds to the 1s(C 1 )→ π * transition, in which the π * is centered on the carboxylic group. This peak is slightly blueshifted with respect to the other two forms. The second set of transitions forms a broader band found between 288 and 294 eV, corresponding to transitions of the other carbons to the σ * orbitals. The σ * orbitals are mixed with Rydberg-type orbitals, which explains the broad range of energies. For all forms, a second intense peak is found around 291 eV, followed by several shoulders at higher energies.

The C K-edge spectra of proline in water can be largely affected by the changes in the σ * orbitals, which depend directly on the C-C distances for all the carbons (see S6 in the supporting information), on C-O distances for C 1 (see Fig. S3 in the supporting information) and on N-C distances for C These values are between 0.01 and 0.02 Å larger than in the gas phase.

The solvent effect is small in the case of the C K-edge (for the decomposition of the solvent spectra, see Fig. S16 of the supporting information.) For the zwitterionic and anionic forms, the π * transition is blue-shifted by 0.1 and 0.3 eV and slightly more intense with respect to the gas phase, whereas the cationic form is redshifted by 0.1 eV and slightly less intense with respect to the gas phase. The effect of the solvent on the σ * peaks is small. Indeed, only the cationic form has a splitting into two peaks appearing at 290.6 and 292.2 eVs.

Comparing the gas phase, solvent, and noelectrostatics, we observe a similar trend as in the case of the oxygen K-edge in the cationic form. Indeed, the fluctuation of the structure blueshifts the π * transition, while the electrostatics redshifts it with respect to the gas phase. In the zwitterionic, this effect is more clearly observed for the σ * transitions. The anionic form seems to be largely unaffected by the presence of electrostatic interactions.

pH-dependent spectra

Proline exhibits two pKa values, one corresponding to the deprotonation of its carboxylic acid group (pKa 1 =1.95) and one corresponding to the deprotonation of its amine moiety (pKa 2 =10.6). [START_REF] Dawson | Data for Biochemical Research[END_REF] Here, we construct the pHdependent X-Ray spectra. For the sake of simplicity, we will focus the discussion on the N K-edge, for which experimental results are also available. [START_REF] Messer | pH Dependence of the Electronic Structure of Glycine[END_REF] The pH-dependent O and C K-edge spectra can be found in Fig. S17 of the supporting information.

The pH-dependent XAS spectra are constructed by using the computed spectra of the solvated anion, zwitterion, and cation weighted by each population extracted from the Henderson-Hasselbach equation for the zwitterionic (p Z ), anionic (p A ), and cationic (p C ) forms (see Fig. 2). Imposing that their populations add up to one, the populations are given by p Z = 1 1 + 10 pH-pKa 2 + 10 pKa 1 -pH (14)

p C =p Z • 10 pH-pKa 1 p A =p Z • 10 pH-pKa 2 . ( 15
)
The populations for the different pH values are plotted in Fig. 2. Fig. 7 (top) shows the computed N K-edge spectra at pHs 6.3 and 12.2 compared to the experimental spectra. [START_REF] Messer | Local Hydration Environments of Amino Acids and Dipeptides Studied by X-ray Spectroscopy of Liquid Microjets[END_REF] Since we simulated only the lower part of the spectrum, the higherenergy region (above 408 eV) is not reproduced in our theoretical model. Converging such a spectrum would require a convergence of the continuum and probably an extension of the QM region to better describe that quasiionization region and the electron transfer to water, which is beyond the scope of the present model.

At pH 6.3, the spectrum corresponds essen- tially to the zwitterionic form, whereas at pH 12.2 the spectrum corresponds to the anionic form. At both pH values, we observe an excellent agreement between the experimental and theoretical spectra, both in terms of intensity and excitation energies. We confirm thus the disappearance of the σ * N -H peak at pH=6.3, whereas this is still observed for pH=12.2. The broadening of the main σ * C-N transition at pH=7 is slightly more compact in the theoretical spectrum compared to the experimental spectrum. The double peak structure observed for that peak might is thus not observable experimentally.

In Fig. 7, we reconstruct the theoretical XAS spectrum of proline along a range of different pH values (0, 2, 7, 10.5, and 12). At pH=0, the population corresponds to the cationic form. At pH=2, the population is half-half between the cationic and zwitterionic forms. At pH=7, the population is essentially the zwitterionic state. At pH=10.5, a half-half population corresponds to the zwitterionic and anionic forms, and at pH=12 only the anionic form is present. As we have seen from the experimental spectra, it is possible to clearly distinguish the anionic and zwitterionic forms from the loss of intensity of the σ * C-N peak, the redshift of the spectrum and the appearance of the σ * H-N peak. The next question is whether the cationic form (appearing at acidic pH values) can be distinguished from the zwitterionic form through the N K-edge XAS spectrum. In principle, despite the double peak structure being experimentally not observed, the two forms could be distinguished by a change in the absorption maximum of the σ * C-N peak, which should be slightly blueshifted at acidic pH. For the oxygen Kedge, which could potentially be used for identifying the three protonation forms, most of the spectrum cannot be recorded since it falls outside the water window at which water oxygens absorb above 533 eV. For the carbon K-edge, the absorption spectrum is similar for the main transition to the π * , and not sufficiently structured for the σ * transitions to clearly identify the three protonation forms.

Conclusions

We have developed an efficient theoretical protocol for computing the solvent effects of X-ray spectra using an atomistic description of the solvent. The protocol is based on the combination of mixed-reference time-dependent densityfunctional theory (MRSF-TDDFT) with electrostatic embedding QM/MM based on electrostatic potential fitting charges.

We have applied this protocol to compute the nitrogen, oxygen, and carbon K-edges of three protonation forms of the amino acid proline. We describe the solvent effects in each of the edges by comparing the gas phase spectra with the spectra generated in water droplets with and without the electrostatic interaction with the solvent. We show that the 1s → π * peak, the most intense peak in the carbon and oxygen K-edges, is largely unaffected by the presence of the solvent, experiencing only a small shift due to the electrostatic interaction with water. At variance, the 1s → σ * transitions have an effect on both the electrostatic interaction with the solvent and the fluctuations of the structure. These fluctuations affect the mixed character of σ * orbitals with Rydberg-type orbitals. In the case of the N K-edge, in which only σ * transitions are observed, the largest effect of the solvent is observed. Indeed, in this core, the three forms of proline can be distinguished from the X-ray absorption spectra as shown from the reconstruction of pH-dependent spectra.

In summary, the ∆MRSF/MM protocol provides an efficient theoretical framework for describing solvent effects on X-ray spectroscopy. Currently, one limitation exists in this methodology, namely, that the SOMO→SOMO configuration is excluded from the excitation space. This implies that the converged core-hole reference must be in a state with a negligible transition dipole moment for constructing a reasonable X-ray absorption spectrum. This will be solved in the future by extending the excitation space included in the MRSF methodology.

Supporting Information Available

The supporting information contains the comparison of experimental and theoretical proline in the gas phase, statistical analysis of geometrical parameters for the QM/MM MD trajectories, the molecular orbitals representative of the different X-ray transitions, the effect of the counterions on the gas phase X-ray absorption spectra of proline, the decomposition of the solvent spectra in cores, and the pH-dependent oxygen and carbon K-edge absorption spectra.

Figure 1 :

 1 Figure 1: Schematic representation of the theoretical protocol for simulating solvent effects of X-ray spectra.

Figure 2 :

 2 Figure 2: (a) Cationic, zwitterionic, and anionic forms of proline depending on the pH. The atom numbering used throughout the article is indicated in the cationic form. The pH populations (lower panel) have been constructed using the Henderson-Hasselbach equation, with carboxylic pKa 1.99 and amino pKa 10.6. The pKa values in water have been extracted from Ref. 49. (b) Structure of the droplet model for the zwitterionic form of proline including sodium (in purple) and chloride (in green) counterions. Similar structures have been obtained for the anionic and cationic forms, with or without suitable counterions.

Figure 3 :

 3 Figure 3: Normalized radial distribution functions for the proline droplet models for the anionic, zwitterionic, and cationic forms during the 150 ps production QM/MM molecular dynamics. The distance shown in (a) corresponds to amine's hydrogen (H) with water oxygen (OW), and in (b) corresponds to the carboxylic oxygen (O) with water hydrogen (HW). The black lines correspond to the experimental spectra of the zwitterionic form, registered in a 1:10 proline/water solution, reported in Ref. 18. The dashed line corresponds to the carboxylic hydrogen and OW in the cationic form.

Figure 4 :

 4 Figure4: Nitrogen K-edge absorption of (a) zwitterionic, (b) anionic, and (c) cationic in solution (solid line), gas phase (dashed), and no electrostatics (dotted). The gas phase XAS spectrum comes along with its analysis (shaded area) and corresponds to 1s → σ * N -H (blue) and 1s → σ * C-N (red). The spectra for the solvent model converged up to 406.96, 404.21 and 407.51 eV for the zwitterionic, anionic, and cationic structures, respectively.

Figure 5 :

 5 Figure5: Oxygen K-edge absorption of (a) zwitterionic, (b) anionic, and (c) cationic in solution (solid line), gas phase (dashed), and no electrostatics (dotted). The gas phase XAS spectrum comes along with its analysis (shaded area) and corresponds to oxygen 1 (blue) and 2 (red) according to the nomenclature of Fig.2. The spectra for the solvent model converged up to 536.51, 536.61, and 537.58 eV for the zwitterionic, anionic and cationic structures, respectively.

Figure 6 :

 6 Figure6: Carbon K-edge absorption of (a) zwitterionic, (b) anionic, and (c) cationic in solution (solid line), gas phase (dashed), and no electrostatics (dotted). The gas phase XAS spectrum comes along with its analysis (shaded area) and corresponds to carbon 1 (blue) and carbon 2 to 5 (red) according to the nomenclature of Fig.2. The spectra for the solvent model converged up to 290.16, 290.49 and 290.10 eV for the zwitterionic, anionic, and cationic structures, respectively.

  2 and C 5 (see Fig. S4 in the supporting information.) Unlike the C-N and C-O bonds, the C-C bonds are less affected by the change in the protonation state of proline, with average values that are very close in the three different models (r C-C = 1.54 ± 0.03 Å for the zwitterion, r C-C = 1.55 ± 0.04 Å for the anion and r C-C = 1.54 ± 0.04 Å for the cation).

Figure 7 :

 7 Figure7: pH-dependent N K-edge spectra of proline. (top) comparison between experimental (dashed) and theoretical (solid) spectra at pH=6.3 (orange), and pH=12.2 (green); (bottom) theoretical spectra at pH=0 (pink), pH=2 (blue), pH=7 (orange), pH=10.5 (red), and pH=12 (green). The spectra converged up to 408.07, 405.19, and 408.55 eV for the zwitterionic, anionic, and cationic structures, respectively. Experimental data are reproduced from Ref. 5.

Acknowledgement M.H.R and N.F. acknowledge the support from "Agence Nationale de la Recherche" through the project MAP-PLE (ANR-22-CE29-0014-01). Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high-performance computing resources. This work was also supported by the NRF funded by the Ministry of Science and ICT ( 2020R1A2C2008246 and 2020R1A5A1019141 to C.H.C.).

TOC Graphic

Theory Experiment