Elise Duboué

Driven molecular dynamics DMSO Dimethyl sulfoxide DO Distributed origin ECD Electronic circular dichroism FF Force field FFMD

First I would like to thank my supervisors, Carine Clavaguéra and Florent Calvo

I would like to include Christian and Lou, for their motivation during their internships with me, Nitik for starting the muffin trend on that fateful day after a CESFO meal, the members of the Institut de Chimie Physique (ICP) laboratory and all the ThéoSim group for their kind interactions over the last 3 years.

Thank you to my parents for always being there to support me, and thank you dad for the last minute proofreading of my thesis. Rosalyn, thank you for putting up with your nerdy sister and for the champagne. Love you, dad, mom and Rosalyn. Thank you to friends and family, that I have not previously mentioned for being part of my life. Umber, you will not be able to read this, as your a horse, but thank you for all the amazing rides we've had together.

Last but definitely not least, Maxime, thank you for your profound support and love. The page of my thesis is turning, the next page I wish it to be with you. Je t'aime La forme des spectres VCD est très sensible aux changements mineurs de conformation et aux interactions moléculaires, ce qui en fait une sonde de l'isomérie conformationnelle et de la solvatation. L'attribution des spectres VCD expérimentaux nécessite en général une comparaison avec les calculs théoriques, ce qui implique également que le solvant soit précisément décrit.

Synthèse en français

Ce manuscrit décrit comment une modélisation basée sur des champs de force polarisables peut être exploitée pour interpréter les signaux VCD. Cette modélisation est appliquée à des molécules flexibles pour lesquelles les approches conventionnelles de chimie quantique reposant sur des approximations statiques sont limitées. En particulier, il est difficile de prendre en compte les anharmonicités et les effets de température finie, et le solvant est typiquement modélisé de manière simplifiée par un continuum ou par un nombre limité de molécules explicites.

Dans le travail présenté ici, le champ de force AMOEBA mis en oeuvre dans le logiciel Tinker a été utilisé pour modéliser les moments dipolaires électriques et magnétiques de diverses molécules pour des applications en phase gazeuse, ainsi que dans les phases condensées des systèmes solvatés ou cristallins. L'échantillonnage et les informations dynamiques ont été obtenus par simulations de dynamique moléculaire classique. Les spectres IR et VCD anharmoniques obtenus à partir des trajectoires ont été interprétés à l'aide d'une analyse des modes effectifs dans laquelle les coordonnées internes sont décomposées linéairement par projection sur des modes effectifs appropriés, les poids correspondants étant obtenus de telle sorte que le pic de chaque mode soit le plus localisé possible. Cette analyse permet de mettre en évidence les contributions des différents groupes fonctionnels de la molécule aux pics individuels du spectre.

Comme premier exemple, le cas du trans-1-amino-2-indanol solvaté dans le diméthylsulfoxyde (DMSO) a été examiné en combinant des explorations statiques au niveau de la théorie de la fonctionnelle de la densité, des trajectoires courtes de dynamique moléculaire ab initio (AIMD) et un large échantillonnage de la surface d'énergie potentielle à l'aide du champ de force. La comparaison avec les mesures existantes confirme l'importance de la flexibilité pour cette molécule, le rôle des molécules de solvant proches et leur contribution aux signaux VCD. La fluxionalité du système composé d'une molécule de trans-1-amino-2-indanol et d'une molécule de DMSO est illustrée dans la figure de courbes de contour suivante : La détermination des spectres VCD à partir du champ de force polarisable a ensuite été entreprise et ses performances ont été évaluées par comparaison avec des calculs de structure électronique. Cette comparaison a été effectuée sur la même trajectoire MD, en utilisant l'acide aminé alanine en phase gazeuse comme système de référence. Les différentes contributions aux moments dipolaires sont discutées. Après avoir validé la méthodologie, la convergence statistique des spectres IR et VCD a été quantifiée via des mesures d'erreurs dédiées. Ces erreurs montrent à leur tour que la convergence des spectres VCD est plus lente que celle des spectres IR et quelle n'est atteinte qu'après un échantillonnage bien plus important que celui réalisable par AIMD. L'approche par champ de force polarisable a ensuite été appliquée à l'alanine en phases condensées. Là encore, la comparaison avec les simulations AIMD confirme que les moments dipolaires électriques et magnétiques peuvent être évalués de manière fiable.

- Ce travail ouvre la voie à la détermination du spectre VCD de systèmes chimiquement complexes inaccessibles via des méthodes qui tiennent compte explicitement de la structure électronique.

Bibliography [1] Jasco, Spectra of alpha-pinene, 2022. https://jascoinc.com/applications/ measurement-of-vibrational-circular-dichroism-spectra-using-the-fvs-6000/ (accessed April 25th, 2022).

Chapter 1 Introduction

Chirality

A chiral molecule is a type of molecule that cannot be superimposed on its mirror image; like our hands. This resemblance is also shown in the etymology of the word chirality which is derived from the Ancient Greek χειρ (kheir) "hand". A chiral molecule has two enantiomers denoted R and S, represented in figure 1.1. These names come from rectus and sinister, the Latin for right (as in the sense of correct or virtuous) and left. Many occurences of chirality can be found in Nature. For example, mint and caraway seeds are composed of the same molecule: carvone. The mint contains the R enantiomer and the caraway seed the S enantiomer [1]. The shells of snails provide another case of chirality, with most shells turning anticlockwise when the head of the snail is to the left. A small minority have a shell of opposite chirality, with the shell turning clockwise [2]. Both of these examples are depicted in figure 1.2. Amino acids are organic molecules composed of an amino group (-NH 2 ) and a carboxyl group (-COOH). They are an essential building block of life as they form proteins. Most amino acids are chiral and have a dominant form just like snail shells. The nomenclature for amino acids is different, they are referred to as L (laevus) and D (dexter). To tell which enantiomer is shown, the CORN law can be used. CORN is an acronym for -COOH , the -R group (representation of side chain) and -NH 2 . According to the CORN law, looking along the hydrogen-carbon bond represented in figure 1.1, if the -COOH, R and NH 2 are clockwise then the molecule is in the L form. If they are anticlockwise then the molecule is in the D form. The L form corresponds to the R enantiomer and the D form to the S notation. This is independent of the (+) and (-) notation, where the (-) configuration rotates linear polarised light anticlockwise and the (+) configuration rotates this light clockwise. Most amino acids are of the L form, and this enantiomeric excess is also found in meteorites [3]. On the contrary, sugars are found mostly in the D form. One hypothesis for this dissimilarity is that the L amino-acids caused the preference for D sugars by acting as a catalyst [4].

Carvone

Deoxyribonucleic acid (DNA), which contains the genetic information of living organisms, forms a double helix, as shown in figure 1.3. This structure is chiral, as any kind of helix. There are 3 main configurations of DNA: A-DNA, B-DNA, and Z-DNA. The first two are right-handed and the last one left-handed. In living cells, most of the DNA is in a mixture of A-and B-DNA, with a few small regions capable of forming Z-DNA [5]. Like the helix, some mathematical objects also include chirality such as the Möbius strip. Finally, chiral molecules are present in medication. In this case, one enantiomer can cure but the other can be harmful. For example, the S enantiomer of ethambutol cures tuberculosis but the other enantiomer causes blindness [7]. It is the same for D-penicillamine, which is used for the treatment of rheumatoid arthritis, while Lpenicillamine is toxic as it inhibits the action of pyridoxine, an essential B vitamin [8]. It can thus be crucial to differentiate between the two enantiomers of the same chiral molecule.

Vibrational spectroscopy

Spectroscopy refers to the response of matter, in absorption or emission, to an excitation by an oscillating electromagnetic field as a function of the wavelength or frequency of this oscillating field. A light beam has its electric and magnetic component both perpendicular to the direction of propagation. A specific direction for these components can be selected with a polariser. A vertical polariser will give linear polarised light. The infrared (IR) light has a wavelength longer than visible light, as shown in figure 1. 4. The mid-IR-range corresponds to light with a wavelength between 16 667 nm and 2 500 nm, or between 600 and 4 000 cm -1 . Molecules absorb light in this range at wavelengths corresponding to the vibrations between its atoms. Thus, depending on how the atoms are assembled in the structure, different wavelengths will be absorbed.

Spectroscopy that studies molecular vibrations is called vibrational spectroscopy. It is a useful tool to determine the structural composition of complex molecules, including organic molecules. Current practice relies on the (successful) comparison between experimental measurements and theoretically calculated spectra, which is particularly useful to assign the peaks to specific molecular motions. Two methods of vibrational spectroscopy will be discussed here: IR and Raman spectroscopy.

Infrared spectroscopy

The most common vibrational spectroscopic technique relies only on the absorption of IR light. The frequency that is absorbed matches the frequency of the molecular vibration or vibrational mode of a sample. In order for a vibrational mode to be considered as IR active, it must be associated with changes in the electric dipole moment.

IR spectroscopy is exploited in organic and inorganic chemistry. It is used by chemists in the determination of molecular structures and functional groups and in the qualitative and quantitative identification of chemical species. This is possible as, usually, stronger bonds and lighter atoms will vibrate at a higher stretching frequency, thus the different groups will show various peaks at different places in the spectrum. IR spectroscopy is used in the investigation of matter in the solid, liquid, and gas phases. If a solvent, for example water which has a large electric dipole, is used, it must not absorb IR light in the same region as the molecule of interest as not to hide the signal.

Raman spectroscopy

Another vibrational spectroscopic technique uses light scattering, which was first predicted by Smekal in 1923 [12]. In this phenomenon, the incident light of energy ℏω i is scattered by a sample and subsequently some light of energy ℏω s is emitted. Most of the scattered light is at the same wavelength as the incident light source and does not provide useful information; this is called Rayleigh scattering. However a small amount of light (typically 0.0000001%) is scattered at different wavelengths, which depend on the chemical structure of the system; this is called Raman scattering. This effect is named after one of its discoverers, who observed it in organic liquids in 1928 together with Krishnan. Independently Landsberg and Mandelstam observed the effect in inorganic crystals [13,14]. The difference in wavelength in Raman scattering arises from a photon exciting the molecule into a virtual energy state for a short time before the photon is emitted. The emitted photon is of either lower or higher energy than the incident photon, as depicted in figure 1.5. While IR intensities depend on the changes in the electric dipole moment associated with the corresponding molecular motion, Raman scattering probes the higher-order component, namely the electronic polarisability.

Raman scattering is used to identify molecular structures, as IR spectroscopy does. Both spectroscopic methods are used for determining chemical bonds, with IR being sensitive to hetero-nuclear functional group vibrations and polar bonds, especially OH stretching in water, whereas Raman scattering is more sensitive to homo-nuclear molecular bonds. For example, it can distinguish between C-C, C=C and C≡C bonds. In solid-state physics, Raman spectroscopy is used to characterize materials, measure temperatures, and locate the crystallographic orientation of a sample.

The key advantage of Raman spectroscopy is that it requires little to no sample preparation while the IR method has constraints on sample thickness, uniformity and dilution to avoid saturation. Raman spectra can be collected from a very small volume (<1 µm in diameter) and also water does not interfere with the signal. However with Raman spectroscopy, fluorescence can obstruct the ability to record the spectra, the intense laser radiation can destroy the sample or swamp the Raman spectrum and Raman scattering is typically very weak. To enhance the weak scattering in the Raman process, nanostructured materials can be used as a support. This is called surface-enhanced Raman spectroscopy (SERS) and is depicted in figure 1.6. Figure 1.6: Raman scattering enhancement from a surface-enhanced Raman spectroscopy (SERS) substrate. This figure has been reprinted from reference [15].

Circular dichroism

While conventional spectroscopy in the IR or optical ranges provides a wealth of information about molecular structure, it is unable to distinguish between enantiomers. The incident electromagnetic field itself must carry some chirality in order to make the two enantiomers of a same molecule react differently, and this is achieved by using circularly polarised light.

With a quarter-wave plate, linear polarisation can be turned into circular polarisation. Such a plate dephases a given linear component of light relative to its orthogonal linear component. Circularly polarised light turns around the axis of propagation in both directions giving, left-and right-handed circular polarisations. Both of these polarisations form a helix which has a chiral form. Depending of the direction of polarisation, chiral molecules absorb this light differently. This is called circular dichroism (CD) (see figure 1.7).

Circular dichroism is mainly used in two ranges: ultraviolet/visible (UV-vis) and the IR range. The UV-vis range gives rise to electronic circular dichroism (ECD) and the IR range gives vibrational circular dichroism (VCD).

The rotation of polarised light by a molecule can also be induced by the application of a magnetic field. This perturbs the energy levels of the system, as was first observed by Faraday [16]. This is called the Faraday effect, and is the basis of magnetic circular dichroism (MCD) [17]. It can be used to probe both the electronic and geometric structures of metal-containing systems, such as metal centres in metalloproteins [18].

Circular dichroism can also be observed in photoelectron angular distributions following photoionization by left-and right-circularly polarised light [19]. This is called photoelectron circular dichroism (PECD). The photoelectron angular distributions show the distribution of directions where the electrons are emitted. A study of camphor, by Garcia and coworkers showed changes in the photoelectron angular distributions depending on the enantiomer and on the direction of the polarisation of light, as depicted in figure 1.8. figure has been reprinted from reference [19].

Electronic circular dichroism

Circular dichroism was first discovered by the French scientist Cotton in 1896, who measured CD in the near-ultraviolet (UV) spectral region [20,21]. UV-visible light is used to scrutinize electronic transitions of a system [22].

In the context of biopolymers, ECD is useful to give insight into the local planar amides which have conformational chirality due to their coupling in the polymer. This coupling gives the characteristic ECD bandshapes for helices and sheets. ECD requires lower concentrations than VCD and has little spectral interference with important solvents. However, important transitions of biosystems are present at wavelengths below 250 nm, which extend into a region (under 200 nm) where ECD is less sensitive and where there is more interference from both water and quartz. In addition, spectral simulation of ECD is challenging due to the fact that it involves transitions to excited electronic states, which are demanding for quantum chemistry methods [23].

Vibrational circular dichroism 1.3.1. Introduction to vibrational circular dichroism

A chiral molecule is considered with its two enantiomers. Each enantiomer absorbs the opposite amount of left or right circularly polarised light. The vibrational circular dichroism intensities arise from the difference between the absorbance of left-circularly polarised (LCP) and right-circularly polarised (RCP) light: ∆A = A L -A R , in the IR range. The intensities of the VCD spectrum for each enantiomer are opposite in intensity [24]. This is shown in figure 1.9 with the IR and VCD spectra of alpha-pinene taken as an example. The spectrum of the R enantiomer is represented in green and the one of the S enantiomer in orange. The IR spectra overlap but their VCD spectra are opposite.

VCD has become an invaluable tool in the determination of absolute configurations by assigning recorded spectra to structural information obtained from theoretical calculations [25]. It has been applied in the fields of nanoscience [26], catalysis [27], solid state organisation [28,29,30], and has also been employed to detect protein fibrils such as amyloids that are responsible for neurodegenerative disorders [31].

From a more fundamental point of view, VCD provides a very sensitive probe of conformational flexibility and molecular interactions [32,33,34,35,36,37,38,39,40]. At room temperature, flexible molecules can adopt many stable conformations corresponding to different local minima, each of them contributing to the VCD spectrum.

Non-covalent interactions have a strong effect on the VCD spectra, even in cases where IR absorption remains unchanged [37,38,39,41]. Consequently, effects like supramolecular chirality and chirality transfer can be addressed as well [33,42].

IR spectrum VCD spectrum

Figure 1.9: Spectra of alpha-pinene (neat, 50 µm pathlength BaF 2 liquid cell). This figure has been adapted from reference [24], with Abs. corresponding to the absorption.

Experimental VCD spectroscopy

While early experimental predictions of VCD were first conceptualised around 1970 [43,44,45,46,47], experimental work started on instruments capable of measuring VCD at the Stephens laboratory [48] and the Holzwarth laboratory [49]. VCD was first observed in the liquid state by Holzwarth et al. in 1974 [50], and by Nafie et al. in 1975 [51] on 2,2,2trifluoro-L-phenylethanol.

To make VCD a practical technique for determining the enantiomer of a chiral molecule, two essential tools needed to be developed. First, a spectrometer with a frequency limit lower than 1600 cm -1 and a higher sensitivity (i.e. the signal-to-noise ratio) since VCD signals are small (∼10 -4 ) compared to IR signals (∼10 0 ). This first goal was achieved around 1980 with Fourier transform infrared (FTIR) spectrometers.

The usefulness of the technique could be demonstrated for differentiating polypeptides. For example, poly-L-lysine can be converted from a disordered structure to α-helix and to β-sheet conformations by increasing the pH and the temperature [23]. This type of study in an aqueous solution is simplified by exchanging H 2 O to D 2 O, which involves an H/D exchange of the amides but enables measurements on more dilute solutions. In figure 1.10, the corresponding IR spectra are shown in the lower frames and the VCD spectra in the upper frames. These correspond to the spectra of the three configurations of poly-L-lysine and show that both helix conformations have nearly identical IR spectra but different VCD spectra. The peaks correspond to the amine I region of the spectra.

Theoretical VCD and the comparison to experimental VCD spectra

The second essential tool needed for the development of VCD was efficient theoretical models to interpret and assign the experimental results [52]. The VCD signal can change sign with a slight modification in conformation or non-covalent interactions, so VCD models need to be sufficiently accurate to be able to compare results reliably [53].

From the 1990s, the Stephens theory of VCD using density functional theory (DFT) [54] was implemented in the Gaussian suite of programs by Cheeseman and Frisch [21,Chapter 1]. More details about the DFT method will be given in section 2.3.

Example of the use of VCD spectroscopy to determine absolute molecular configurations

As an illustration of the state-of-the-art in the 2000s, the IR and VCD spectra determined using DFT for molecules of the pyrazole type, shown in Fig. 1.11 [55], are reproduced in Fig. 1.12. The pyrazole type molecules were separated into 2 diastereomers ( 2a and 2b), for this example, the absolute configuration of diastereomer 2a was examined. This figure is adapted from reference [55].

The DFT calculations were performed, for various configurations, with the B3LYP functional and 6-31G(d) basis-set in the Gaussian 03 software. For the experimental spectrum, a solvent of CDCl 3 was used.

The configurations 2a-A and 2a-B are the lowest in energy and their spectra are the most similar to the experimental spectrum. It was shown that a mixture of 60% of 2a-A and 40% of 2a-B gives the VCD spectrum closest to the experimental spectrum. The study accomplished in the article used for this example did not include solvation for the theoretical part. However in VCD experiments, the molecules of interest are usually solvated and one fundamental issue is to assess the role of the solvent on the molecular conformation and on the spectrum itself. This role can depend on the nature of the solvent and whether it has a strong interaction with the solute. Successful assignment of VCD signals based on computed spectra therefore requires that the solvent is correctly described.

Solvation and theoretical VCD

VCD is known to have a high sensitivity to bulk solvent effects. This was demonstrated by Dezhahang et al. [56], who used a system of the TOLBINAP (2,2'-Bis(di-p-tolylphosphino)-1,1'-binaphthalene) palladium complex Pd(TOLBINAP)Cl 2 with the structure depicted in figure 1.13. Figure 1.13: Structures of the TOLBINAP (2,2'-Bis(di-p-tolylphosphino)-1,1'-binaphthalene) (Ar=4-CH 3 -C 6 H 4 ) palladium complex, as taken from reference [56].

The goal of this study was to determine the conformational changes in the TOLBI-NAP ligands before and after their coordination with palladium [56]. The incorporation of an implicit polarisable continuum solvation model provided much better agreement between theory and experiment. This is shown in figure 1.14 where the spectra with the solvent described implicitly with a polarisable continuum (PCM) model show better agreement with the experimental spectrum. The regions showing noticeably different VCD patterns in the gas phase and in solution are highlighted with dashed-dotted lines.

In both the studies of Tur et al. [55] and Dezhahang et al [56], the calculations were performed assuming static configurations and without any explicit solvent. A way to include an explicit solvent is through the cluster-in-a-liquid model where the solute-solvent interactions are described explicitly through n solvent molecules that form a cluster with the solute. The surrounding bulk environment is accounted for with an implicit solvent [57]. To go beyond the static view, dynamic simulations such as molecular dynamics (MD) can be employed. MD is a powerful tool to sample the phase space of a molecular system at finite temperature in a realistic way [58] and more details about this technique will be given in chapter 2.

Thesis outline

After having introduced VCD spectroscopy in the previous sections, chapter 2 will illustrate various methods used to compute VCD spectra, pointing out the limitations of static approximations. This will bring us to the utility of the development of a method to obtain VCD spectroscopy through classical molecular dynamics based on force fields to include anharmonicites and finite temperature effects. It is this method that has been implemented during this PhD using the AMOEBA polarisable force field in the Tinker software package.

In chapter 3, a first illustrative case concerns trans-1-amino-2-indanol solvated in dimethyl sulfide which is examined using a combination of methods. This case demonstrates the importance of the sampling of the potential energy surface and of the ability to reach VCD spectra through classical molecular dynamics based on advanced force fields.

In chapter 4, this implementation using the AMOEBA polarisable force field is validated first with alanine in the gas phase, revealing the importance of extensive sampling of the system, which goes beyond the capacities of ab initio molecular dynamics. The various contributions of the electric and magnetic dipole moments are discussed, illustrating the importance of including electrostatics effects in the force field. Comparison with density-functional theory results along one MD trajectory confirms that electric and magnetic dipole moments and thus VCD spectra can be reliably evaluated using the polarisable force field. Information about the spectral assignment is not directly included in the molecular dynamics, thus a method to determine effective vibrational modes from the analysis of the MD trajectories is employed. In this method, the internal coordinates are linearly decomposed by projection onto suitable effective modes, the corresponding weights being obtained in such a way that the peak of each mode is as localized as much as possible. This effective mode analysis sheds light onto the various contributions of the different functional groups of the molecule to individual peaks in the vibrational spectrum.

In chapter 5, the force field combined to MD approach is extended and applied to alanine in crystalline and hydrated phases. The method is also validated with a comparison to DFT results and the effective mode analysis gives the vibrational modes corresponding to the peaks in the spectrum. Alanine solvated in dinitrogen is also examined, as this solvent binds more weakly to the alanine solute, the absence of hydrogen bonding between them limits the perturbation on the vibrational spectra. The final system studied is phenylcyclohexanediol in the crystalline phase and solvated in dimethyl sulfoxide. This relatively large molecule presents inter-molecular hydrogen bonding between molecules in the crystal and between the solvent and the solute in the solvated phase. The different environment of the hydroxyl groups are shown to have an impact on the spectra.

In chapter 6, the results are summed up in a conclusion and future work is detailed in the perspectives.

Chapter 2 Methodology

General theory of vibrational circular dichroism

Circular dichroism is the difference in absorption of left-and right-circularly polarised light. It originates from molecular interactions with the electric and magnetic fields of the light, which causes an oscillation of the electrons. These oscillations are in phase with the light's fields for one circular polarisation and out of phase for the other [1].

The interaction Hamiltonian, which expresses the molecule's energy in the lights electromagnetic field, is

H = -ì µ • ì E -ì m • ì B (2.1)
where ì µ is the molecule's electric dipole moment and ì m is the molecule's magnetic dipole moment. ì E is the externally applied electric field and ì B is the externally applied magnetic field. The electric field (E-field) and the magnetic field (B-field) of linear polarised light, propagating toward the positive x-axis vary as:

ì E(t) = E 0 êe i(κ x-ωt) ì B(t) = B 0 êe i(κ x-ωt) (2.2)
where E 0 and B 0 are the wave amplitudes, ê is a polarisation vector with unit magnitude, κ is the wavenumber, related to the wavelength λ through κ = 2π/λ, and ω is the temporal frequency [2]. These fields are represented in red for the E-field and in blue for the B-field in figure 2.1A.

For circularly polarised (CP) light, the E-field rotates around the propagation axis. For right-circularly polarised (RCP) light, the E-field rotates anti-clockwise as the viewer looks down the propagation axis toward the origin, as shown in figure 2.1B. For left-circularly polarised (LCP) light the rotation is clockwise. This convention is reversed if the viewer looks up the propagation axis toward the detector as in 2.1C and D. These rotations can be described by adapting the polarisation vector ê. Therefore the electric field of rightand left-circularly polarised light propagating toward the x-axis can be expressed as:

ì E(t) RCP = E 0 √ 2 (ŷ -iẑ)e i(κ x-ωt) , ì E(t) LCP = E 0 √ 2 (ŷ + iẑ)e i(κ x-ωt) (2.3)
where ŷ and ẑ are unit vectors on the y-axis and z-axis, respectively.

For CP light, the B-field rotates 90 • out of phase, behind for RCP and ahead for LCP as depicted in figure 2.1C and D. The B-field amplitude is equal to E 0 c , thus the magnetic field for left-and right-circularly polarised light can be written:

ì B(t) RCP = i c ì E(t) RCP , ì B(t) LCP = - i c ì E(t) LCP (2.4)
where c is the velocity of light. The molecular extinction coefficient (ε) conveys how strongly a molecule absorbs light at a given wavelength. This extinction coefficient is constructed from Fermi's golden rule which gives the transition probability between an initial state |i and a final state | f [3]:

ε = πN A ωδ(ω -ω fi ) ε 0 ℏc log10 1 E 2 0 | f | H |i | 2 (2.5)
where N A is Avogadro's number, ε 0 is the electric permittivity of free space, ℏ = h 2π with h being Planck's constant and δ(ωω fi ) a Dirac δ function.

The Hamiltonian can be adopted using the previous equations for the electric and magnetic dipole moments. Although these equations come from classical physics, they also apply in quantum mechanics via the replacement of the position and velocity vectors with their quantum operators. These operators can be written as:

ì µ = ì µ ii ì µ fi ì µ if ì µ ff ì m = ì m ii ì m fi ì m if ì m ff (2.6)
where each matrix element represents a bracket. In the electric dipole matrix, ì µ ii = i| ì µ |i is the static electric dipole moment of the molecule in its initial state, and ì After including the equations for electric and magnetic fields and the electric and magnetic dipole moment operators in equation (2.5) and following the rotational averaging procedure described by Andrews [4], the following expression for the extinction coefficient for circular dichroism is found:

µ ff = f | ì µ | f ,
ε CD = 4πN A ωδ(ω -ω fi ) 3ε 0 ℏc 2 log10 ( ì µ fi • ì m if ).
(2.7)

This expression corresponds to the Rosenfeld equation [5]. Different versions of the prefactor exist depending on the units chosen [6,7,8]. The magnetic dipole moment is taken into account for the VCD representation so we can consider the rotation around the axes of propagation, characteristic of circularly polarised light. The corresponding equation of the infrared (IR) signal [9,10] is:

ε IR = 4πN A ωδ(ω -ω f i ) 3ε 0 ℏc 2 log10 ( ì µ fi • ì µ if ) (2.8)
where ( ì µ fi • ì µ if ) corresponds to the dipole strength [11].

Harmonic approach

A transition from an initial to a final state occurs due to the absorption or emission of light. A transition has an associated frequency depending on the energy of the light absorbed or emitted ∆E, given by Planck's relation:

∆E = hω f i (2.9)
With light in the IR range, this transition corresponds to a vibrational transition, from one vibrational energy level to another, as shown in figure 2.2. The energy of a molecular system is the sum of its potential (V) and kinetic (T ) energies. For molecular conformations that reside near a local minimum of the potential energy surface, a Taylor expansion can be used to approximate the potential energy as:

V = V 0 + 3N j=1 ∂V ∂p j 0 p j + 1 2 3N j,k ∂ 2 V ∂p j ∂p k 0 p j p k + • • • (2.10)
where p j are the mass-weighted coordinates (MWC) of atom j. MWCs are given relative to the equilibrium position (x j,0 , y j,0 , z j,0 ) and are multiplied by the square root of the mass of the corresponding atom:

p 1 = √ m j (x j -x j,0 ), p 2 = √ m j (y α -y j,0 ), p 3 = √ m j (z α -z j,0 ), ...
As we consider an equilibrium geometry, the geometry corresponds to a minimum of the potential energy surface, where the gradient term vanishes. The consideration of only equilibrium geometries is indicated via () 0 . The harmonic approximation ignores the higher terms of the expansion. Thus the potential is reduced to:

V harmonic = V 0 + 1 2 3N j,k ∂ 2 V ∂p j ∂p k 0 p j p k = V 0 + 1 2 3N j,k f MWC i j p j p k (2.11)
where

f MWC i j = f MWC ji
is the i j component of the mass-weighted Hessian denoted H MWC .

To obtain IR and VCD spectra, the vibrational frequencies are computed by determining the second partial derivatives of the potential energy V with respect to displacement of the atoms in MWC, at the equilibrium geometry. This corresponds to the ij component of the mass-weighted Hessian, f MWC i j

, mentioned in equation (2.11):

f MWC i j = ∂ 2 V ∂p i ∂p j 0 (2.12)
The mass-weighted Hessian matrix, also called the dynamical matrix, is diagonalised to obtain the 3N eigenvalues corresponding to the harmonic frequencies of the system and the 3N eigenvectors corresponding to the atomic displacements associated with each frequency. These eigenvectors are named normal modes and correspond to linear combinations of atomic displacements for which the dynamical matrix is diagonal. The normal modes corresponding to global rotations and translations of the molecule are associated with vanishing eigenvalues. Removing these modes leaves N vib = 3N -6 independent normal modes, or 3N -5 if the molecule is linear.

The IR intensities arise from the dipole strength D i and VCD intensities come from the rotational strength R i [12]. As shown in equations (2.8) and (2.7), the vibrational in- tensities arise from the determination of the electric ( ì µ fi ) and magnetic transition dipole moments ( ì m if ) [13].

A second harmonic approximation is also taken where fundamental properties describing the spectroscopic intensities, in our case the dipole moment surfaces, are treated as harmonic [14].

Density-functional theory

Basic principles of DFT

To describe an atomic or molecular system and to determine its properties, various computational modelling methods can be used. The first one explored here is densityfunctional theory (DFT).

Density-functional theory is based on the first and second Hohenberg and Kohn theorems [15]. The first theorem declares that all physical properties of a molecular system can be obtained from its electronic density ρ(ì r). Thus the energy is a density functional that depends only on the three spatial coordinates :

E[ρ] = T[ρ] + V ee [ρ] + V N e [ρ]
(2. 13) where T[ρ] is the kinetic energy of the electrons, V ee [ρ] is the electron-electron repulsion energy and V N e [ρ] the electron-nuclei attraction energy. This dependence of the energy on the electronic density differs from the Hartree-Fock (HF) method, where the energy is a function of the wavefunction that depends on 3N space coordinates and N spin coordinates (for N electrons).

The second Hohenberg and Kohn theorem states that the variational principle can be applied. Thus the ground-state energy, E 0 , is always less than or equal to the expected value of the energy E:

E[ρ] ⩾ E[ρ 0 ] = E 0 .
(2.14) However the exact analytical expression of T[ρ] and V ee [ρ] is not known. To attempt to access the energy, the real system of interacting electrons is substituted by a fictitious system of independent electrons having the same density, following the Kohn-Sham approach [16]. One-electron wavefunctions χ n are introduced and obtained by solving the Kohn-Sham equation, where v eff is an effective potential:

-1 2 ∇ 2 + v eff (r) χ n = ε n χ n .
(2.15)

Hence, the corresponding non-interacting Hamiltonian does not contain the explicit electron-electron repulsion term but rather an effective potential. This equation is solved iteratively until self-consistency is achieved.

The effective potential contains an exchange-correlation potential, V XC . The exchangecorrelation can be decomposed into the exchange term V X , which represents the interaction of electrons with the same spin, and into the correlation term V C , which represents the interaction of electrons with opposite spins. The exchange-correlation energy functional corresponds to the integration of the exchange-correlation potential over the electronic density:

E XC [ρ] = ∫ ρ(ì r)V XC dì r.
(2.16)

Unfortunately, the exact exchange-correlation functional is not known either. Thus, approximate exchange-correlation functionals are employed, leading to approximate densities and energies. Over the last decades, many exchange-correlation functionals have been proposed and tested [17,18,19,20,21,22,23,24]. They are usually distinguished as the rungs of a ladder, going from basic/simple expressions to improved ones [25], as shown in figure 2.3. The family of least accurate exchange-correlation functionals are functionals that use the local density approximation (LDA), they depend only on the electronic density. Then there are the functionals that employ the generalized gradient approximation (GGA) that depend on the density and its gradient. If the functional also depends on the Laplacian of the density, it is a meta-GGA. The functionals with the highest accuracy exploit exact exchange where the functional can be expressed in terms of occupied orbitals. During this thesis, if nothing is mentioned, the functionals B3LYP and the basis-set 6-311++G(d,p) were used. The comparison of the quality of mid-IR spectra of 4-methyl-2-oxetanone predicted using various density functionals demonstrates the quality of the B3LYP functional [26], though this reference is dated. More recent work on on vibrational spectroscopy also used this functional [27,28].

For the solvated systems, the D3BJ correction was added to the B3LYP functional to include dispersion effects. The term D3BJ corresponds to the D3 version of Grimme's dispersion with Becke-Johnson (BJ) damping [29]. Grimme's dispersion is an atom-pair wise potential for the treatment of the dispersion energy which helps to reach the so-called chemical accuracy, currently assumed to be 1 kcal/mol for the energies [30,31]. Including the BJ damping shows better results for nonbonded distances, which is important when considering the interactions between the solute and the solvent.

Application to infrared and vibrational circular dichroism spectroscopies

An early application of density-functional theory to obtain IR and VCD intensities has been reported by Cheeseman et al. [32]. The vibrational dipole and rotational strengths of 2,3-trans-d 2 -oxirane were determined within the harmonic approximation and compared to experimental values. The structure of the molecule is shown in figure 2 The two modes shown with a red arrow are the main contributors to the reduction of the mean absolute deviation from experiment, as their correlation with experiment is excellent for DFT. The arrows point to the HF values and the corresponding DFT point are on the same vertical line. These modes are assigned experimentally to hydrogenic bending modes [33]. In contrast, the DFT agreement with experimental results is worse for 3 of the 4 hydrogenic stretching modes, shown with black arrows. These modes hinder the reduction of the mean absolute deviation from experiment. Ignoring these modes, the mean absolute deviation drops from 9.7 for HF to 4.0 for the DFT calculations.

By comparing frequencies and dipole strengths computed at different levels of theory and comparing theory and experiment for other achiral isotopomers of oxirane, the authors also deduced that an important source of error of the calculated results compared to experimental results is the neglect of anharmonicity. Thus, in the next section, anharmonicity will be discussed.

Anharmonicity effects

Within the harmonic approximation, the potential energy between 2 atoms can be described as the following potential:

V N=2 harmonic = 1 2 kr 2 , with k = ∂ 2 V ∂r 2
(2.17) and r 2 = p 1 p 2 , with p i the mass-weighted coordinates of atom i [34].

A more accurate model involves the use of the Morse potential, which is described with the following equation [35]:

V N=2 Morse = D e 1 -e -a(r-r e ) 2
(2.18)

Here r is the distance between the atoms, r e is the equilibrium bond distance, D e is dissociation energy, and a controls the width of the potential, the smaller a is, the larger the well.

The harmonic and Morse potentials are compared to each other, together with their vibrational levels, in figure 2.6. The energy spacing between vibrational levels and thus the frequency, decreases with increasing energy of the vibrations. The harmonic and Morse potentials are close to each other only for the ground vibrational level (v = 0). Also, the further the value of the internuclear separation r is from the value of equilibration r c , the more the two potentials differ. Hence the harmonic approximation is only valid for small amplitudes of vibration. has been reprinted from reference [36].

For the value of the vibration levels higher than 1, there is a shift between the corresponding energies for the Morse and harmonic values. This leads to an overestimation of the harmonic frequencies which is also found in larger systems where the harmonic frequencies are higher compared with experimental results. A first solution to go beyond the harmonic approximation is to apply a scaling factor that shifts the values of the frequencies to correct for anharmonicities. Scaling factors also empirically correct for the various deficiencies of the underlying quantum chemical methods, including basis set incompleteness. Their ideal values have been tabulated for many methods. [37,38,39,40,41]. Depending on the functional and basis set used, these scaling factors range between 0.8 and 0.99, for common organic molecules. For the functional B3LYP and the basis-set 6-311++G(d,p), a scaling factor of 0.98 was used in accordance with similar theoretical methods [28].

The potential surfaces of fluxional molecules, namely molecules that have multiple conformations, are composed of multiple wells that can hardly be described with a single harmonic potential. Also, at sufficiently high temperature, the system has enough energy to access higher vibrational modes which are misrepresented with the harmonic potential, as shown in figure 2.6.

In a second step, approaches that improve over the harmonic approximation are discussed. Two complementary methods will be considered, namely second-order perturbation theory [42] and the so-called cluster-in-a-liquid approach [43].

Perturbation approaches

Within second-order vibrational perturbation theory (VPT2), the potential energy can be expanded up to the next non-vanishing term (order 4) as:

V anharm = 1 2 3N j,k ∂ 2 V ∂p j ∂p k 0 p j p k + 1 6 3N j,k,l ∂ 3 V ∂p j ∂p k ∂p l 0 p j p k p l (2.19) + 1 24 3N j,k,l,m ∂ 4 V ∂p j ∂p k ∂p l ∂p m 0 p j p k p l p m + ...
This equation can be reduced to :

V anharm = 1 2 3N j,k f MWC i j p j p k + 1 6 3N j,k,l k jkl p j p k p l + 1 24 3N j,k,l,m k jklm p j p k p l p m + ... (2.20)
where k jkl and k jklm are the anharmonic force constants.

The idea of VPT2 is to solve the vibrational Hamiltonian within a 2nd order perturbation approach. The vibrational energy [44] of mode n out of a total of M modes, about an equilibrium geometry, is found as a so-called Dunham expansion expressed as:

E n = E (0) n + ℏ M i ω i (n i + 1 2 ) + ℏ M i ≤ j χ i j (n i + 1 2 )(n j + 1 2 ) (2.21)
where the ω i 's are the harmonic wavenumbers and the χ i j 's are the anharmonic contributions that can be given as a function of the anharmonic force constants [45].

The fundamental bands ν i , overtones [2ν i ], and combination bands [ν i ν j ] are then expressed by:

ν i = ω i + 2 χ ii + 1 2 N j=1 j i χ i j (2.22) [2ν i ] = 2ω i + 6 χ ii + N j=1 j i χ i j = 2ν i + 2 χ ii (2.23) [ν i ν j ] = ω i + ω j + 2 χ ii + 2 χ j j + N k=1 j j,i ( χ ik + χ jk ) (2.24)
Foldes et al. used VTP2 to show that hydrogen bonding has a profound effect on anharmonicity [46]. For example, the formation of the hydrogen bond induces a red-shift in the OH and NH vibrations.

Boltzmann weights and the cluster-in-a-liquid model

The VPT2 method generally improves the positions of vibrational bands, however, it still remains unable to account for fluxionality. In this section, the anharmonicity due to the fluxional character is accounted for through a large number of conformations of the system. Each conformation represents a well of the potential energy surface. For each one, a spectrum and a Boltzmann weight, depending on the energy and temperature of the configuration, can be determined. Thus an average spectrum can be obtained by multiplying each individual harmonic spectrum by its Boltzmann weight and summing all of the contributions. In the examples given in this work, the properties of each conformation are determined using DFT.

Moreover, hydrogen bonding is poorly described with second-order perturbation theory. To consider hydrogen bonds and solvation, Perera et al. [43] introduced the use of a cluster-in-a-liquid model where the hydrogen bonds are described explicitly through solute-(solvent) n clusters. Thus the explicit system includes a solute molecule surrounded by n solvent molecules. The surrounding bulk environment is accounted for with an implicit solvent. Including the whole solvent as a explicit bulk would be too expensive for DFT calculations.

Continuum models of solvation represent the solvent as a continuous medium, with the same dielectric constant as the pure solvent, surrounding the solute, instead of individual explicit solvent molecules [47]. In the polarisable continuum model (PCM), the effective solute-solvent interaction potential is determined from the dielectric medium polarised by the charge distribution of the molecular solute [48]. An integral equation formalism variant of the implicit polarisable continuum model (IEFPCM) also exists [49] and is the current version of PCM applied in common quantum chemistry packages.

Merten et al. [50] demonstrated the importance of hydrogen bonding for VCD in the case of methylbenzyl-amine-borane in CDCl 3 . In this molecule there is a special kind of hydrogen bond, named dihydrogen bond (DHB): N-H δ+•••δ-H-B. This bond was identified to be responsible for the unusual stability of the system [51] and is shown by the dotted lines in the left of figure 2.7.

For the VCD spectra, a better agreement is achieved between the experimental spectrum and the spectrum of the dimer than the one of the monomer. The features at ∼1585 cm -1 and 1400-1350 cm -1 are only reproduced in the dimer spectrum, when the DHB is present. Both the monomer and dimer are likely to coexist in solution, though the stronger agreement of the dimer's spectrum to experimental results supports the proposal that a majority of the dimer configuration exists in solution. This confirms that VCD is sensitive to the effect of hydrogen bonding and thus this effect needs to be accounted for when calculating VCD spectra.

Weirich et al. [52] employed the cluster-in-a-liquid model to study how hydrogen bonding affects the VCD spectral signatures of chiral alcohols. They found that hydrogen bonding to the OH-group alters mainly the vibrational modes that have contributions of the H-O-C bending and C-O related deformation motions. They also showed that the farther away the OH-group from the stereocenter, the smaller the effect on the VCD spectrum. However VCD spectra obtained by Boltzmann-averaging of the different minima contributions do not show optimal agreement with experiment for floppy molecules [52,53]. Moreover Ghidinelli et al. [54] claimed that solvents with donor capability like methanol seem particularly difficult to treat satisfactorily with the cluster-in-a-liquid model. According to these authors, the inclusion of anharmonicity and temperature effects through the Boltzmann weights provides a flimsy result. It relies on determining enough structures to describe the system fully and the temperature is considered after the determination of the individual spectrum. Finally the cluster-in-a-liquid model does not include long-range interactions with an explicit solvent.

Classical molecular dynamics

In the previous section, various methods were discussed that offer an improvement on the harmonic approximation, though they can be inaccurate for very floppy molecules which have a very large number of local minima in their potential energy surface. In the cluster-in-the-liquid approach, the temperature is considered through the Boltzmann weights, subsequently to the determination of the frequencies and spectral intensities.

Atomistic simulations such as molecular dynamics (MD) or Monte Carlo are powerful tools to sample the phase space of a molecular system at finite temperature, realistically.

Principle of molecular dynamics

Molecular dynamics simulation methods naturally account for anharmonic and temperature effects and thus the fluxionality of the system. In MD, successive configurations of the system are generated to give a classical trajectory that specifies how the positions and velocities of the particles vary with time. This is completed by integrating Newton's second law of motion, shown here for only one spatial coordinate x:

d 2 x i dt 2 = F x i m i (2.25)
where m i is the mass of particle i, d 2 x i dt 2 is its acceleration along the x-axis and F x i is the force on the particle along the x direction. To integrate Newton's equation of motion, the positions and the dynamical properties (velocities, acceleration, etc.) are approximated as Taylor expansions:

ì r(t + δt) = ì r(t) + δtì v(t) + 1 2 δt 2 ì a(t) + 1 6 δt 3 ì b(t) + 1 24 δt 4 ì c(t) + ... (2.26) ì r(t -δt) = ì r(t) -δtì v(t) + 1 2 δt 2 ì a(t) - 1 6 δt 3 ì b(t) + 1 24 δt 4 ì c(t) + ... (2
ì b(t + δt) = ì b(t) + δt ì c(t) + ... (2.30)
where ì r = (x, y, z) is the position, ì v the velocity, ì a the acceleration and ì b and ì c the third and fourth derivatives of the position with respect to time [55]. To propagate the particles and determine their positions and velocities at the subsequent time step t + δt from the knowledge at time t, the Verlet algorithm is adopted. In this algorithm, the sum of equations (2.26) and (2.27) is used to determine the positions at the next step :

ì r(t + δt) = 2ì r(t) -ì r(t -δt) + δt 2 ì a(t) (2.31)
Then to determine the velocities a simple approach is to divide the difference of the positions at time t + δt and time t -δt by 2δt:

ì v(t) = ì r(t + δt) -ì r(t -δt) 2δt (2.32)
However with this algorithm, which is straightforward and for which the storage requirements are modest, the velocities are only available when the positions of the next step have been computed. In this work, it is the Beeman algorithm [56] that was chosen, as it is naturally implemented in the Tinker software package on which our work relies. The key equations for the Beeman method are given in equations (2.33) and (2.34) for the positions and velocities, respectively:

ì r(t + δt) = ì r(t) + δtì v(t) + 2 3 δt 2 ì a(t) - 1 6 δt 2 ì a(t -δt) (2.33) ì v(t + δt) = ì v(t) + 1 3 δt ì a(t) + 5 6 δt ì a(t) - 1 6 δt ì a(t -δt) (2.34)
This expression of the velocities is also more precise than the Verlet algorithm, which is important when determining the temperature of the system, as explained in the following.

Thermostatted molecular dynamics

To control the temperature, an NVT or canonical ensemble is employed. In this ensemble the number of particles N and the volume V are constant and the system is maintained at a temperature T. This is useful to equilibrate our systems at the desired temperature, generally 300 K. The temperature of the system is related to the time average of the kinetic energy [55] by:

E kin = 3 2 N k B T (2.35)
The average kinetic energy is fixed at a certain temperature, the instantaneous kinetic energy of an N-body system fluctuates with the velocity:

E kin (t) = 1 2 N i=1 m i v 2 i (t) (2.36)
Thus the instantaneous temperature [55,57] can be given by:

T(t) = N i=1 m i v 2 i (t) 3N k B (2.37)
Hence by modifying the velocities is possible to control the temperature. This influence on the MD trajectory to keep the temperature constant is called a thermostat. Three main thermostats exist: Berendsen, Andersen and Nosé-Hoover.

The Berendsen thermostat [58] couples the system to an external heat bath fixed at the desired temperature. At each step, the velocities are scaled such that the rate of change in temperature is proportional to the difference in instantaneous temperature between the bath and the system. However the kinetic energy converges towards a constant with this thermostat, thus the fluctuations of the instantaneous kinetic energy become zero. This absence of fluctuation does not respect the canonical ensemble.

The canonical sampling through velocity rescaling (CSVR) thermostat [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] achieves sampling of the canonical ensemble by propagating the positions and velocities using an areapreserving NVE integrator, then updates the kinetic energy for one time step using some stochastic dynamics in the canonical ensemble, and finally rescales all velocities using the updated kinetic energy.

The next two thermostats respect the canonical ensemble as the instantaneous kinetic energy fluctuates around its average value. The Andersen thermostat [START_REF] Andersen | Molecular dynamics simulations at constant pressure and/or temperature[END_REF] chooses a particle at random and changes its velocity to a value randomly selected from a Maxwell-Boltzmann distribution. This is equivalent to random particles of the system undergoing stochastic collisions with a fictitious particle. The Nosé-Hoover thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] corrects the velocities by exchanging energy with an extra degree of freedom. This method produces correct fluctuations in the velocities and reproduces well the desired canonical ensemble. It is this thermostat that will be adopted in this work.

NVT ensembles influence the MD trajectory via their thermostat. This influence on the trajectory could interfere with the time series of dipole moments computed. Thus to explore the potential energy surface in the goal of retrieving the electric and magnetic dipole moments, the NVE, or microcanonical ensemble is used once the system is equilibrated at the wanted temperature.

Another way to explore the potential energy surface is through the Monte Carlo method where one or more coordinates of the structure are randomly modified at each step. Generally the Metropolis algorithm [START_REF] Metropolis | The beginning of the Monte Carlo method[END_REF] is used where if the new structure has a lower energy that the previous one, it is automatically accepted. If the energy is higher, the new structure is only accepted with a Boltzmann probability depending on the temperature of the system. This method does not permit to access dynamical properties of the system and thus will not be exploited here. The next section will develop the acquiring of IR and VCD spectra with time-dependent properties.

Infrared and VCD absorption spectroscopies from time-dependent properties

While IR spectroscopy measures the response of the electric dipole moment to an external excitation, the VCD signal comes from simultaneous changes in electric and magnetic dipole moments of the molecule due to nuclear motion [START_REF] Berova | Circular dichroism: principles and applications. 605 Third[END_REF]Chapter 4].

To represent the dynamical changes of these two moments, a time-correlation function (TCF) C(τ) can be used. This function provides a statistical description of the time evolution of one or two random variables at different points in time (t i and t j , with τ = t jt i ).

If the random variables represent the same quantity: Abbate and co-workers [START_REF] Abbate | The use of cross-correlation functions in the analysis of circular dichroism spectra[END_REF] introduced the concept of cross-correlation functions of time-dependent electric and magnetic dipole moments for the computation and analysis of circular dichroism spectroscopy, following the work of Gordon [START_REF] Gordon | Correlation functions for molecular motion[END_REF] on autocorrelation functions for IR spectroscopy.

C(τ) = A(t i ) • A(t j ) ),
According to the Wiener-Khintchine theorem, [START_REF] Mcquarrie | Statistical mechanics. Harper's chemistry series[END_REF]Chapter 22] the correlation function and the spectral density S(ω) are the Fourier transform of each other (see Equation 2.38). This theorem is taken within the ergodic hypothesis, which states that for a stationary random process, a large number of observations made on a single system at n arbitrary instants of time have the same statistical properties as observing n arbitrarily chosen systems at the same time from an ensemble of similar systems:

S(ω) = ∫ +∞ -∞ C(τ)e -iωτ dτ (2.38)
Thus for describing the absorption intensity line shape of the IR and VCD signal [8,[START_REF] Abbate | The use of cross-correlation functions in the analysis of circular dichroism spectra[END_REF][START_REF] Calvo | Atomistic modeling of IR action spectra under circularly polarized electromagnetic fields: Toward action VCD spectra[END_REF], we can write:

I IR (ω) = 2π βω 2 3V cn(ω) ∫ +∞ -∞ ì µ(0) • ì µ(τ) e -iωτ dτ
(2.39)

I VCD (ω) = 4π βω 3V cn(ω) ∫ +∞ -∞ ì µ(0) • ì m(τ) e -iωτ dτ (2.40)
However the IR spectral line shape can also be obtained from the autocorrelation function of the time derivative of ì µ, and the VCD spectral line shape from the cross correlation function of the time derivative of ì µ with ì m, which offers a number of practical advantages, including origin independence [8,12]. Various methods were used to apply Gordon and Abbate's concepts to the determination of IR and VCD spectroscopy from a classical trajectory, starting with MD simulations that employ ab initio harmonic force fields to determine the value of the vibrational frequencies of α-pinene, D-gluconic acid, formaldehyde dimer and acetylprolineamide [9].

In
The quantum mechanics/molecular mechanics (QM/MM) method was also used to study alanine dipeptide in water [6,[START_REF] Kwac | Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: Comparisons with IR and vibrational circular dichroism spectra[END_REF], (R)-methyloxirane and (L)-alanine in aqueous solutions [START_REF] Giovannini | Effective fully polarizable QM/MM approach to model vibrational circular dichroism spectra of systems in aqueous solution[END_REF] and (1S)-(-)-β-pinene in CCl 4 [START_REF] Choi | Direct calculations of mid-and near-IR absorption and circular dichroism spectra of chiral molecules using QM/MM molecular dynamics simulation method[END_REF]. Subsequently ab initio molecular dynamics was adopted to study gas [START_REF] Scherrer | Nuclear velocity perturbation theory of vibrational circular dichroism[END_REF][START_REF] Thomas | Classical magnetic dipole moments for the simulation of vibrational circular dichroism by ab initio molecular dynamics[END_REF][START_REF] Ditler | Implementation of nuclear velocity perturbation and magnetic field perturbation theory in CP2K and their application to vibrational circular dichroism[END_REF] and liquid phases [8]. The last method employed is force field-based molecular dynamics to study camphor [START_REF] Calvo | Atomistic modeling of IR action spectra under circularly polarized electromagnetic fields: Toward action VCD spectra[END_REF].

Ab initio molecular dynamics (AIMD) and force field-based molecular dynamics (FFMD) will be the subject of the next two sections.

Ab initio molecular dynamics

In an ab initio molecular dynamics (AIMD) calculation, to include the electronic structure in MD simulations, electrons are treated quantum mechanically and the nuclei which are heavier are modelled as classical particles. Thus the nuclei are propagated via Newton's second law, equation (2.25) and the energies and the forces are computed at the quantum level.

Born-Oppenheimer molecular dynamics

In the Born-Oppenheimer molecular dynamics (BOMD) method, the static electronic structure problem is solved at each molecular dynamics step given the set of fixed nuclear positions at that instant of time. Thus the electronic structure part is treated within the time-independent Schrödinger equation. This corresponds to an adiabatic treatment where the electronic wavefunction is considered to adapt quasi instantaneously to the nuclei movement.

Car-Parrinello molecular dynamics

Though in BOMD the electrons are treated statically, Car and Parrinello [START_REF] Car | Unified approach for molecular dynamics and densityfunctional theory[END_REF] proposed the use of an electronic wavefunction as a variable of the dynamical trajectory, which enables a dynamically evolving electronic subsystem.

While BOMD separates electron and nuclei via time, Car-Parrinello molecular dynamics (CPMD) divides nuclei and electrons through energy to maintain adiabaticity and maps both the electronic and nuclear component classically. Car and Parrinello introduced a new class of Lagrangians that takes into account the energy of the electronic subsystem Ψ 0 | H e |Ψ 0 and the possible constraints that can be imposed on the set of orbitals such as orthonormality:

L C P = I 1 2 M I . R 2 I + i 1 2 µ i . ψ i . ψ i kinetic energy -Ψ 0 | H e |Ψ 0 potential energy + constraints orthonormality (2.41)
with M I and R I the nuclear masses and positions, respectively. The equations of motion of Car-Parrinello MD are:

           M I .. R I = - ∂ ∂R I Ψ 0 | H e |Ψ 0 + ∂ ∂R I constraints σ .. ψ i = - δ δψ * i Ψ 0 | H e |Ψ 0 + δ δψ * i constraints (2.42)
where σ is a fictitious mass associated with the electronic wavefunction. It ensures an integration of the equations of motion within a reasonable time step and adiabatic separation of the electrons and nuclei. The unit of this mass parameter is energy times a squared time for dimensional balance [START_REF] Iftimie | Ab initio molecular dynamics: Concepts, recent developments, and future trends[END_REF][START_REF] Marx | Ab initio molecular dynamics: basic theory and advanced methods[END_REF].

IR and VCD spectra from nuclear velocity perturbation theory

The adiabatic treatment is inadequate for predicting VCD [START_REF] Nafie | Vibronic coupling theory of infrared vibrational transitions[END_REF][START_REF] Nafie | Adiabatic molecular properties beyond the Born-Oppenheimer approximation. complete adiabatic wave functions and vibrationally induced electronic current density[END_REF][START_REF] Scherrer | Nuclear velocity perturbation theory for vibrational circular dichroism: An approach based on the exact factorization of the electron-nuclear wave function[END_REF]. To circumvent this issue, nuclear velocity perturbation theory (NVPT) can be used. In this approach, the total electron-nuclear wavefunction Ψ(ì r, ì R, t) is factorized in the nuclear part χ( ì R, t) and the electronic part Φ ì R (ì r, t), which depend on the electronic coordinates ì r and on the nuclear coordinates ì R:

Ψ(ì r, ì R, t) = Φ ì R (ì r, t) χ( ì R, t) (2.43)
In NVPT, the electronic part is taken as the Born-Oppenheimer (BO) ground state Φ (0) ì R (ì r), corrected by a time-dependent, purely imaginary perturbative contribution iΦ (1) ì R (ì r, t):

Φ ì R (ì r, t) = Φ (0) ì R (ì r) + iΦ (1) ì R (ì r, t) (2.44)
The electronic Hamiltonian can be partitioned into the standard BO Hamiltonian ĤBO with eigenfunctions Φ (0) ì R (ì r) and a small NVPT perturbation:

Ĥel (ì r, ì R) = ĤBO + N n ν=1 λ ν ( ì R, t) • (-iℏ∇ ν ) (2.45)
with the nuclear gradient operator acting on the parametric dependence of the electronic wavefunction. Here, the sum runs over the N n nuclei. In the semiclassical limit for the nuclei, the perturbation parameter λ ν ( ì R, t) corresponds to the nuclear velocity ì V ν ,

λ ν ( ì R, t) = 1 M ν -iℏ∇ ν χ( ì R, t) χ( ì R, t) = ì V ν (2.46)
with the ionic mass M ν . The total perturbative correction to the electronic wavefunction iΦ (1) ì R (ì r, t) is obtained as a weighted sum over the corrections for each component:

Φ (1) ì R (ì r, t) = N n ν=1 3 α=1 λ ν α ( ì R, t)Φ (1) ì R,να (ì r) (2.47)
Using the NVPT treatment, it is thus possible to access the VCD properties through the following expressions:

ì µ = ì µ e + ì µ n = - N e i=1 q ì r i + N n ν=1 Z ν q ì R ν (2.48)
and

ì m = ì m e + ì m n = - N e i=1 q 2c ì r i × ì v i + N n ν=1 Z ν q 2c ì R ν × ì V ν (2.49)
Here, q is the electronic charge, Z ν q is the nuclear charge, and c is the velocity of light in vacuum. The position and velocity operators for the electronic subsystem are indicated as ì r i and ì v i , respectively, and similar symbols are used for the nuclear operators ì R ν and ì V ν .

An implementation of nuclear velocity perturbation theory (NVPT) was carried out in the plane-wave codes CPMD [START_REF] Cpmd | Copyright IBM Corp 1990-2019[END_REF] and in CP2K [START_REF] Ditler | Implementation of nuclear velocity perturbation and magnetic field perturbation theory in CP2K and their application to vibrational circular dichroism[END_REF] and tested in the gas [START_REF] Scherrer | Nuclear velocity perturbation theory of vibrational circular dichroism[END_REF] and liquid [8] phases.

Example of propylene oxide: evidence for chirality transfer

The NVPT approach was applied to (R)-propylene-oxide (RPO) in water, for which chirality transfer was found from a chiral solute to an achiral solvent. With chirality transfer, an achiral species gains some of the chiral properties of the chiral system. This transfer has been theoretically demonstrated for RPO in water but also methyl lactate in water [START_REF] Losada | Solvation of propylene oxide in water: Vibrational circular dichroism, optical rotation, and computer simulation studies[END_REF][START_REF] Losada | Chirality transfer through hydrogen-bonding: experimental and ab initio analyses of vibrational circular dichroism spectra of methyl lactate in water[END_REF]. The system of study is composed of 1 RPO molecule, represented in figure 2.8 and 53 water molecules. According to Scherrer et al. [8] the most relevant contributions to the VCD spectra arise from the molecules in the first solvation shell, for RPO in water. Scaled molecular moments were introduced which include only the molecular moments of a finite region around a chosen center. For a chiral molecule K, the surrounding molecules J are included as:

ì µ scale K (t) = ì µ K (t) + J K P K J (t) ì µ J (t) (2.50) ì m scale K (t) = ì m K (t) + J K P K J (t) ì m J (t) (2.51)
with the damping function being

P K J (t) = 1 + e ( | ì R K J (t) |-R 0 )/D -1
, for intermolecular contributions. R 0 corresponds to the distance included from the center of mass of the solute, it will be referred to as the cutoff. D is the sharpness parameter and it is chosen equal to 0.25 Å as in reference [START_REF] Heyden | Dissecting the THz spectrum of liquid water from first principles via correlations in time and space[END_REF]. It is this cutoff that will be used in chapter 5 for hydrated alanine.

The IR signal of water and VCD signals of RPO in water are shown in figure 2.9. The blue lines correspond to the experimental data [START_REF] Losada | Solvation of propylene oxide in water: Vibrational circular dichroism, optical rotation, and computer simulation studies[END_REF], the red lines to the VCD signal considering only the moments of the RPO molecule. The dashed line is the VCD signal of the RPO molecule and the molecular contributions of the water molecules near to the solute. indicate the standard deviation of the statistically independent microcanonical simulations. The experimental data are taken from reference [START_REF] Losada | Solvation of propylene oxide in water: Vibrational circular dichroism, optical rotation, and computer simulation studies[END_REF] and are in arbitrary units. The computed spectra are divided by the refractive index of bulk water [START_REF] Bertie | Infrared intensities of liquids XVI. accurate determination of molecular band intensities from infrared refractive index and dielectric constant spectra[END_REF] to compare with the experiment. This figure has been reprinted from Scherrer et al. [8].

At 1650 cm -1 the signal corresponding to the bending mode of water is expected. The presence of this VCD signal of achiral water shows the existence of chirality transfer as the VCD spectra of an achiral molecule is expected to be vanishing.

If we consider a cutoff of R 0 = 3.75 Å, only a slight contribution of the water is seen, though if the cutoff is increased to 4.0 Å, the signal becomes more intense and noisy. The disorder can be measured through the mean standard deviation of the VCD spectra:

σ(R 0 ) = N I VCD i (R 0 ) -I VCD av (R 0 ) 2 N (2.52)
The spectra from N=8 microcanonical simulations (I VCD As the value of the cutoff increases, the VCD spectra get more and more disordered, thus considering a cutoff of R 0 = 3.75 Å is a good compromise as it also includes the first solvation shell.

In BOMD or CPMD studies of vibrational spectroscopy, the solvent is typically modeled by a continuum or very few solvent molecules [START_REF] Scherrer | Nuclear velocity perturbation theory for vibrational circular dichroism: An approach based on the exact factorization of the electron-nuclear wave function[END_REF]. This method is also expensive which leads to a short exploration time, a few dozen picoseconds for a system with around 100 atoms.

Classical force fields for molecular simulations

In the previous section, the electronic part was treated with quantum chemistry methods. Here the whole system is treated with classical dynamics. Thus the potential energy and dipole moment surfaces are calculated using a collection of equations and associated constants describing the dependence of the energy of a system on the coordinates of its particles. This is called a force field (FF) thus the name force field-based molecular dynamics (FFMD) when the collection of equations is used within MD. The use of a force field leads to the ability to explore for longer periods of time as the calculation for each step is less expensive than with AIMD. Larger systems can also be considered to include the solvent explicitly.

Nonpolarisable contributions

Force fields for organic systems can be divided into four key components: bond stretching, angle bending, torsional terms and non-bonded interactions. The potential energy V of a system of N atoms can be written as a sum of these four contributions:

V(ì r 1 , ì r 2 , ..., ì r N ) = V bond + V angle + V torsion + V non-bonded (2.53) V(ì r N ) = bond k r 2 (r -r 0 ) 2 + angle k θ 2 (θ -θ 0 ) 2 + torsions V n 2 (1 + cos(nω -γ)) + N i=1 N j=i+1 4ε i j σ i j r i j 12 - σ i j r i j 6 + q i q j 4πε 0 r i j (2.54)
where (rr 0 ) and (θθ 0 ) are the deviation of the bond length and of the angle from their reference values, k r and k θ are corresponding force constants. Thus each bond and angle is described by a force constant and a reference value. The torsions are described by a constant V n , the torsion angle ω and the phase factor γ. The non-bonding part includes a van der Waals term and the electrostatic interactions [55].

Bond length term

This quadratic description of the potential due to the variations of the bond length is appropriate for large systems (proteins, polymers...) and is implemented in the CHARMM [START_REF] Foloppe | All-atom empirical force field for nucleic acids: I. parameter optimization based on small molecule and condensed phase macromolecular target data[END_REF], AMBER [START_REF] Cornell | A second generation force field for the simulation of proteins, nucleic acids, and organic molecules[END_REF], OPLS [START_REF] Jorgensen | Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids[END_REF] and GROMOS [START_REF] Scott | The GROMOS biomolecular simulation program package[END_REF] force fields. For smaller molecules, and to more accurately model the Morse curve, terms with higher orders can be added to V bond , according to a Taylor expansion:

V bond (r) = k r 2 (r -r 0 ) 2 1 -k (r -r 0 ) -k (r -r 0 ) 2 -k (r -r 0 ) 3 ... (2.55)
This expression of V bond is used in the MM1 [START_REF] Allinger | Calculation of molecular structure and energy by force-field methods[END_REF], MM2 [92], MM3 [93] and MM4 [94] models developed by Allinger and dedicated to organic molecules. As the cubic term passes through a maximum, this potential must only be used close to the reference value, where it describes the Morse potential well. To prevent the overly large lengthening of bonds, a quartic term can also be added, which is the case in MM3. The AMOEBA force field used in this thesis, and described below, is based on the MM3 formalism [START_REF] Ponder | Current status of the AMOEBA polarizable force field[END_REF].

Figure 2.12: A cubic bond-stretching potential passes through a maximum but gives a better approximation to the Morse curve close to the equilibrium structure than the quadratic form. This figure has been reprinted from reference [55].

Angular term

A Taylor expansion can also be applied to the angle bending, through the following expression:

V angle (θ) = k θ 2 (θ -θ 0 ) 2 1 -k (θ -θ 0 ) -k (θ -θ 0 ) 2 -k (θ -θ 0 ) 3 ... (2.56)
As the Taylor expansion of the bond length term, this expression of the angle terms is used in the MM1, MM2, MM3 and MM4 series of force fields and in the AMBER, CHARMM and AMOEBA force fields.

Torsional term

Most force fields for organic molecules include the torsional term explicitly. This is mostly achieved through a cosine series expansion:

V torsions (ω, γ) = torsions V 1 2 (1 + cos(ω -γ)) + V 2 2 (1 -cos(2ω -γ)) + V 3 2 (1 + cos(3ω -γ)) (2.57)

Non-bonding terms

The non-bonding terms describe inter-and intra-molecular interactions which are described with a van der Waals (vdW) term and an electrostatic term. The van der Waals term includes the repulsion (∝ r -12 ) and attraction (∝ r -6 ) terms are often described by the 12-6 Lennard-Jones potential:

V vdW = 4ϵ σ r 12 - σ r 6 (2.58)
where σ is the distance at which E vdW = 0, ϵ is the depth of the well represented in figure 2.13 and r is the distance between two atoms. The repulsion term accounts for Pauli's repulsion between electrons which prevents the overlapping of the electron clouds. This term dominates when r is lower than the equilibrium value r 0 . At large values of r, it is the attractive term that dominates. This term accounts for the weak dispersion attraction between two atoms at long range. The idea of the electrostatic term is to represent the electronic response of the system to its environment, it is represented by Coulomb's law:

Intermolecular potential V(r)
V Coulomb = q A q B 4πε 0 r AB (2.59)
This law describes the interaction between two point charges on two atoms (A and B). The corresponding force F Coulomb (r) = -dV Coulomb /dr can be represented by the following, as in figure 2.14. This description of the vdW term and of the electrostatic interaction terms is used in conventional force fields. However these force fields are not accurate enough for spectroscopy. A reason for this is the neglect of molecular polarisation. This has led to recent FF development to improve the treatment of electrostatics [START_REF] Jing | Polarizable force fields for biomolecular simulations: Recent advances and applications[END_REF][START_REF] Huang | Mapping the drude polarizable force field onto a multipole and induced dipole model[END_REF][START_REF] Lopes | Current status of protein force fields for molecular dynamics simulations[END_REF]. These improvements include the use of multipoles centered on each atom and permit the introduction of polarisation and charge transfer effects [START_REF] Semrouni | Finite temperature infrared spectra from polarizable molecular dynamics simulations[END_REF][START_REF] Thaunay | Strategy for modeling the infrared spectra of ion-containing water drops[END_REF].

r AB

F A +e A F B +e B F A +e A F B -e B

Accounting for polarisation

Polarisation of electronic clouds is the deformation experienced by the charge on atoms in response to the induced electric field from the charges of the other atoms. To describe this effect, three models are commonly employed: the fluctuating charges model [START_REF] Rick | Dynamical fluctuating charge force fields: Application to liquid water[END_REF], the Drude model [START_REF] Rick | Potentials and algorithms for incorporating polarizability in computer simulations[END_REF] and the induced dipoles model [START_REF] Jing | Polarizable force fields for biomolecular simulations: Recent advances and applications[END_REF][START_REF] Cieplak | Polarization effects in molecular mechanical force fields[END_REF]. The manner in which polarisation is accounted for in these models is schematised in figure 2.16.

(a) (b) (c) q q + Q q +δe -Q μ ind

Fluctuating charges

In the fluctuating charges model, the fixed charges are replaced with charges that can vary by incorporating an additional electrostatic energy term. This term represents the energy required to modify the charge of an individual atom. The idea is to redistribute the atomic charges to equalize the electronegativity at each site, according to the electronegativity equalization principle [START_REF] Bret | A chemical potential equalization model for treating polarization in molecular mechanical force fields[END_REF]. The total electrostatic energy is minimized, for each time step, with respect to the atomic charges. For a system of N atoms with charges q 1 , q 2 , ..., q N , the electrostatic energy is written as:

V Fluct = N i=1 E i0 + χ 0 i q i + 1 2 N i=1 N j i q i q j J i j (2.60)
where E i0 is the energy of atom i at zero charge, χ 0 i is the electronegativity and J i j is a term that mimics the Coulomb interaction at long distance but, typically, smoothly reaches finite values at vanishing distances.

To obtain the charges of the system, equation (2.60) is minimized with respect to the sum of the charges q 1 , q 2 , ..., q N which are fixed for the whole system. This yields a (N + 1) × (N + 1) linear system to solve:

∂V Fluct ∂q 1 = ∂V Fluct ∂q 2 = ... = ∂V Fluct ∂q N and i q i = Q tot = const.
(2.61)

This model permits the calculation of polarisation without the addition of any new interactions. However one problem with the fluctuating charges method is that it does not reproduce out-of-plane polarisation for planar or linear chemical moieties. This is understandable because electronegativity equalization proceeds along the bonds between atoms.

Drude model

In the Drude model (or shell model), the fluctuation of the charge is represented by a charge (Q) on a harmonic oscillator centred on the atom. The displacement of this charge, by a distance ì d, in response to an electric field creates a dipole moment ì µ Drude on atom i:

ì µ Drude i = -Q i ì d i (2.62)
The total energy of a system described with the Drude model can be written:

V Drude (ì r i , ì d i ) = N i=1 1 2 k i d 2 i + q i ì r i • E 0 i -(r i + d i ) • E 0 i + N i=1 1 2 N i=1 N j i 1 r i j - 1 |r i j -d j | - 1 |r i j + d i | + 1 |r i j -d i + d j | (2.63)
where k i is the spring constant, E 0 i is the static field at the location of the core charge, r i , and E 0 i is the static field at the location of the shell charge,

r i + d i . Note that E 0 i E 0 i
in general. This model uses particle-particle interactions, thus keeping the computational cost low.

Induced dipoles

Most polarisable models developed nowadays employ the induced dipole model approach [START_REF] Caldwell | Structure and properties of neat liquids using nonadditive molecular dynamics: water, methanol, and n-methylacetamide[END_REF][START_REF] Cieplak | A new water potential including polarization: Application to gas-phase, liquid, and crystal properties of water[END_REF][START_REF] Dang | Importance of polarization effects in modeling the hydrogen bond in water using classical molecular dynamics techniques[END_REF][START_REF] Ponder | Force fields for protein simulations[END_REF][START_REF] Ren | Polarizable atomic multipole water model for molecular mechanics simulation[END_REF][START_REF] Walsh | A multipole-based water potential with implicit polarization for biomolecular simulations[END_REF]. The induced dipole model adds induced dipoles, centered on each atom i, to the permanent charge:

ì M i = ì M 0 i + ì M ind i (2.

64)

M i includes fixed point charges, as well as possibly permanent dipoles and quadrupoles:

ì M i = [q i , µ ix , µ iy , µ iz , Q ixx , Q ixy , Q ixz , . . . , Q izz ] (2.65)
with q i the point charge on atom i, µ the dipole and Q the quadrupoles.

The induced dipole is proportional to the electric field, ì E i on site i. The electric field arises both from the permanent charges of the system and from the other induced dipoles:

ì µ i ind = α i ì E i = α i j T 1 i j ì M 0 j + j T 11 i j ì µ ind j (2.66)
with α i the isotropic atomic polarisability, T 1 i j a 3 x 13 matrix representing the second to the fourth rows of the multipole-multipole interaction matrix T i j and T 11 i j is a 3 x 3 submatrix consisting of elements in T i j corresponding to the dipole moments. The sum over j represents all atomic sites outside the polarisation group of atom i and the sum over j corresponds to all atomic sites other than i itself. T i j is the interaction matrix between sites i and j:

T i j =             1 ∂ ∂x j ∂ ∂y j ∂ ∂z j . . . ∂ ∂x i ∂ ∂x i x j ∂ ∂x i y j ∂ ∂x i z j . . . ∂ ∂y i ∂ ∂y i x j ∂ ∂y i y j ∂ ∂y i z j . . . ∂ ∂z i ∂ ∂z i x j ∂ ∂z i y j ∂ ∂z i z j . . . . . . . . . . . . . . . . . .             1 r i j (2.67)
Since the induced dipoles also contribute to the electric field, they are solved iteratively through a self-consistent field.

The polarisation energy due to the induced dipoles can be expressed as:

V ind = - 1 2 N i=1 ì µ i ind t • ì E 0 i (2.68)
where ì E 0 i is the electric field due to the permanent charges and () t the transpose.

Most developments in the polarisable force fields for biomolecules relate to explicitly induced dipoles and somewhat less to fluctuating charges and Drude oscillators. This is despite induced dipoles requiring iterative computation, which is a much slower process compared with the other additive models. The main advantage of this method is it only requires the addition of a few parameters to the force field, namely the atomic polarisabilities. It is the induced dipoles model that is used by the AMOEBA force field.

Periodic boundary conditions

To model liquid and solid phases, periodic boundary conditions (PBC) [55] are needed to mimic infinitely extended systems. For a simulation of a system using PBC, the molecules are represented in a unit cell, any molecule that crosses one boundary of the cell reappears on the opposite side, as represented in figure 2.17.

Figure 2.17: Visual representation of the idea of periodic boundary conditions. This figure has been reprinted from reference [START_REF]Model Box Periodic Boundary Conditions[END_REF].

In PBC and with a globally neutral system, long distance interactions can been described with the Ewald summation. In this method, the long-range energy contributions (usually selected by a cutoff function) are calculated in Fourier space, where they converge faster than in real space [START_REF] Toukmaji | Ewald summation techniques in perspective: a survey[END_REF]. In the program used for the MD simulations, Ewald summation is used for the polarisable atomic multipoles.

The non-bonded van der Waals interaction is typically of short range relative to, for example, electrostatic forces. This interaction is truncated at a cutoff distance. The value of the corresponding van der Waals potential is brought to zero for distances larger than the cutoff [START_REF] Roth | Van der waals interactions involving proteins[END_REF].

The AMOEBA force field

The AMOEBA (atomic multipole optimized energetics for biomolecular simulation) polarisable force field [START_REF] Ponder | Current status of the AMOEBA polarizable force field[END_REF] has been used in this thesis to describe the molecular systems of interest through molecular dynamics simulations. AMOEBA was originally developed for biomolecules such as proteins and DNA, and some organic molecules. The main advantage of the force field is its ability to reproduce relative conformational energies [START_REF] Rasmussen | Force field modeling of conformational energies: importance of multipole moments and intramolecular polarization[END_REF].

Conventional force fields, that do not include polarisation effects, only reproduce ∼50 % of the conformations found at the MP2/aug-cc-pVTZ level of theory whereas AMOEBA reproduces ∼80 % of the conformations [START_REF] Kaminskỳ | Force field modeling of amino acid conformational energies[END_REF]. This leads to a good precision for thermodynamic properties such as free energies of solvation and complexation. The AMOEBA force field has also been extended to IR spectroscopy very satisfactorily [START_REF] Semrouni | Finite temperature infrared spectra from polarizable molecular dynamics simulations[END_REF][START_REF] Thaunay | Vibrational mode assignment of finite temperature infrared spectra using the amoeba polarizable force field[END_REF].

The good reproduction of the relative energy is due to the use of permanent electrostatic multipolar moments on each atom by this force field. These moments include point charges as well as dipole and quadrupoles moments and are derived from high-level gas phase quantum mechanical calculations [START_REF] Shi | Polarizable atomic multipole-based amoeba force field for proteins[END_REF]. This derivation is performed via Stone's distributed multipole analysis (DMA) where the main idea is to distribute the charge density from the quantum wavefunction into local multipoles moments [START_REF] Stone | Distributed multipole analysis, or how to describe a molecular charge distribution[END_REF].

In practice, AMOEBA combines the traditional ingredients of nonpolarisable models, including bonded and non-bonded terms, with a multipolar description of electrostatics on each atom up to the quadrupole, and induced dipoles reacting self-consistently to the presence of local electric fields on each polarisable site.

Previously a general overview of the contribution to the energy of the system was discussed, then, the specificities of AMOEBA are now presented [START_REF] Ponder | Current status of the AMOEBA polarizable force field[END_REF].

V = V bond + V angle + V bθ + V oop + V torsion + V vdW + V ElectroStat (2.69)
where the first five terms correspond to bond stretching, angle bending, bond-angle cross term, out-of-plane bending, and torsional rotation. These are the short-range valence interactions. The last two terms correspond to the non-bonded van der Waals and electrostatic contributions, including Coulomb and polarisation interactions.

Covalent interactions

The bond stretching V bond , angle bending V angle and the coupling between the stretching and bending V bθ are represented by equations (2.70) to (2.72). These correspond to those of the MM3 force field [START_REF] Allinger | Molecular mechanics. the MM3 force field for hydrocarbons[END_REF] in which anharmonicity is accounted for by the terms of order 3 and 4.

V bond = K b (r -r 0 ) 2 1 -2.55(r -r 0 ) + 7 12 2.55 2 (r -r 0 ) 2
(2.70)

V angle = K θ (θ -θ 0 ) 2 1 -0.014(θ -θ 0 ) + 5.6 × 10 -5 (θ -θ 0 ) 2 -7 × 10 -7 (θ -θ 0 ) 3 + 2.2 × 10 -8 (θ -θ 0 ) 4
(2.71)

V bθ = K bθ (r -r 0 ) + (r -r 0 ) (θ -θ 0 ) (2.72)
where K b , K θ and K bθ are force constants.

The out-of-plane bending term V oop is controlled by a Wilson-Decius-Cross function [START_REF] Wilson | Molecular vibrations: the theory of infrared and Raman vibrational spectra[END_REF] for the sp 2 -hybridized trigonal centers only:

V oop = K χ χ 2 (2.73)
where K χ is the force constant and χ the dihedral angle between the 4 atoms that are considered.

The torsion terms here originate from a Fourier expansion of a 1-fold through 6-fold trigonometric form:

V torsion = 6 n=1 K nω [1 + cos(nω ± γ)] (2.74)

57

The other coupling terms (torsion-bond, torsion-angle and angle-angle) incorporated in the MM3 force field have not been conserved in AMOEBA.

Non-bonded interactions

The first non-bonded interaction is the van der Waals interaction term V vdW which takes the buffered 14-7 [START_REF] Halgren | The representation of van der waals (vdW) interactions in molecular mechanics force fields: potential form, combination rules, and vdW parameters[END_REF] functional form instead of the standard 12-6 Lennard-Jones form:

V vdW (i j) = ϵ i j 1.07 ρ i j + 0.07 7 1.12 ρ 7 i j + 0.12 - 2
(2.75)

with ρ i j = r i j /r 0 i j (2.76)
r 0 i j = (r 0 ii ) 3 + (r j j ) 3 (r 0 ii ) 2 + (r 0 j j ) 2
(2.77)

ϵ i j = 4ϵ ii ϵ j j ( √ ϵ ii + √ ϵ j j ) 2 (2.78)
where ϵ i j is the potential well depth, r i j is the distance between atoms i and j and r 0 i j is the minimum energy distance. The equation applies to heterogeneous atom pairs. The parameters involved in these expressions are fitted to both gas phase and bulk phase experimental properties. Halgren [START_REF] Halgren | The representation of van der waals (vdW) interactions in molecular mechanics force fields: potential form, combination rules, and vdW parameters[END_REF] showed that the conventional 12-6 Lennard-Jones potentials are too repulsive at short contact distances and the buffered 14-7 potentials reproduced well ab initio results of noble gases.

The electrostatic term has two contributions. The first arises from the permanent mutipolar moments on each atom, and replaces the Coulomb term described previously. The second contribution emerges from the induced dipole moments.

In the induced dipoles model, the set of induced dipoles is solved iteratively and is accelerated in the AMOEBA force field via successive overrelaxation (SOR) [START_REF] Young | Iterative solution of large linear systems[END_REF] using a default value of ω = 0.7 in:

ì µ i ind (n + 1) = (1 -ω) ì µ i ind (n) + ωα i j T 1 i j ì M 0 j + j T 11 i j ì µ ind j (n) (2.79)
In this expression of the induced dipoles, the atomic polarisabilities α i are present, for which the values have been previously determined through Applequist, Carl, and Fung's model [START_REF] Applequist | Atom dipole interaction model for molecular polarizability. Application to polyatomic molecules and determination of atom polarizabilities[END_REF].

However Thole [START_REF] Thole | Molecular polarizabilities calculated with a modified dipole interaction[END_REF] pointed out that at small distances, r 12 ≤ (4α 1 α 2 ) 1/6 , the polarisation energy becomes infinite. To circumvent this "polarisation catastrophe", he proposed a damping function to this energy:

ρ = 3a 4π exp(-au 3 ) (2.80)
with u = r i j /(α i α j ) 1/6 and a a dimensionless parameter that controls the damping strength.

In the AMOEBA force field this distance-dependent damping function is combined with the exclusion of some of the 1-2 bonded polarisation interactions [START_REF] Ren | Consistent treatment of inter-and intramolecular polarization in molecular mechanics calculations[END_REF].

Summing all the energy contributions gives the energy of the system as described by the AMOEBA force field.

The AMOEBA force field is distributed with the Tinker software package [START_REF] Ponder | TINKER -Software Tools for Molecular Design (version 8)[END_REF], developed by Ponder and coworkers. It is this software that was modified to incorporate the ingredients needed to calculate VCD spectra through force field molecular dynamics.

The Tinker software [START_REF] Ponder | TINKER -Software Tools for Molecular Design (version 8)[END_REF] provides a good framework to perform molecular dynamics as well as other energy-based and structural manipulation calculations such as local energy minimisations. It is a code composed of relatively small programs that inter-operate to perform complex computations. Both these elements ease the addition of new code.

Implementation of VCD spectroscopy within AMOEBA

In section 2.5.3, the absorption intensity line shape of an IR spectrum is shown to be obtained via the Fourier transform of the autocorrelation of the electric dipole moment. Likewise, the line shape of a VCD spectrum can be obtained via the Fourier transform of the cross-corelation of the electric and magnetic dipole moments. Thus to reach these spectra, the electric and magnetic dipole moments must be acquired throughout time.

The dynamics of our molecular system is described by an atomistic force field, in which the potential energy and gradient assume a particular form that explicitly accounts for multipolar electrostatics, including polarisation forces and induced dipoles. For any configuration R of the system, the total electric dipole moment ì µ contains a contribution at lowest order from the partial charges {q i } distributed at the respective positions ì r i , and a contribution of the induced dipoles ì µ ind i on each polarisable site i:

ì µ(t) = N i q i ì r i (t) + ì µ ind i (t) (2.81)
where N denotes the number of atoms. Once the electric dipole moment is determined for a molecule, it's time derivative is calculated.

However the calculation of the magnetic dipole moment ì m(t) is not included in the force field and software adopted here. Thus we define our magnetic moment by basing the expression on the magnetic moment used for conventional force fields [START_REF] Calvo | Atomistic modeling of IR action spectra under circularly polarized electromagnetic fields: Toward action VCD spectra[END_REF]:

ì m ConvFF (t) = 1 2c N i q i ì r i (t) × ì v i (t) (2.82)
where ì v i (t) denotes the velocity vector of atom i at time t, q i its charge and c the velocity of light.

The transformation from conventional force field to polarisable FF involves the addition of the induced dipole moment that needs to participate in the magnetic moment.

Thus a second term is added, corresponding to the moving induced dipole as given by Hnizdo [START_REF] Hnizdo | Magnetic dipole moment of a moving electric dipole[END_REF]:

ì m MoveInd (t) = 1 c ì µ ind i (t) × ì v i (t) (2.83)
The effects of the time varying induced dipoles moments due to the motion of neighbouring charges also need to be considered. Therefore a last term is added to the magnetic moment:

ì m Fluc (t) = 1 2c ì r i (t)× ì µ ind i (t) (2.84)
This contribution is also present in the description of the total magnetic dipole moment in equation ( 41) of reference [8].

Thus we define the magnetic dipole moment as the sum of the moving point charges, the moving induced dipoles and the time varying induced dipoles moving around fixed positions:

ì m(t) = N i 1 2c q i ì r i (t) × ì v i (t) + 1 c ì µ ind i (t) × ì v i (t) + 1 2c ì r i (t)× ì µ ind i (t) (2.85)
It is this expression of the magnetic moment that was added to the Tinker software package, for this work.

Returning to the expression of the electric dipole moment, the velocity form of this

moment ì µ(t) = d ì µ(t)
dt is exploited to avoid the origin dependence that arises through the positions [12,Chapter 4]. This origin dependence also arises within a periodic cell, which will be developed in the next section.

Specific difficulties under periodic boundary conditions

In the case of PBC, the minimum image convention is used. This ensures that the distance between the origin O and any atom i at a position ì r i is the lowest, taking the periodic image if necessary. This convention is expressed by translating ì r i by ∆ PBC Oi , a conditional lattice translation, which is a multiple of the size of the PBC cell.

Thus, depending on where the origin is placed, the original value of ì r i or a periodic image value will be used. This changes the values of ì r i and therefore ì m according to the origin, as shown by the black lines in figure 2.18. To address this issue, Jähnigen et al [7] have recently proposed an alternative formalism that establishes an origin invariant version of the TCF function. This is achieved by considering the nearest periodic image directly between atoms instead of through the origin, as represented by the red lines in figure 2.18 and by the following equation:

C indep (t) = ì µ(0) • ì m(t) + 1 2c i ì µ i (0) j ì r j (t) + 1 2 ∆ PBC i j (0) × ì µ j (t) (2.86)
where

∆ PBC i j
is the conditional lattice translation used for the nearest-image convention between i and j. The first term is the sum of the correlation of moments within and between locally finite cells. These cells, generally corresponding to one molecule, each have a local origin. Partitioning the space as such corresponds to the distributed origin (DO) gauge [12]. The second term is the transformation from DO gauge to the common origin gauge, where all the particles have the same origin.

This equation was implemented in the ChirPy program [START_REF] Jähnigen | Chirpy -A python package for analysing supramolecular and electronic structure, chirality and dynamics[END_REF], which was used to proceed from the time correlation functions to the corresponding spectra, as discussed in the next section.

ChirPy is a python package for analysing supramolecular and electronic structure, chirality and dynamics [START_REF] Jähnigen | Chirpy -A python package for analysing supramolecular and electronic structure, chirality and dynamics[END_REF]. This package can, among other things, compute timecorrelation functions (TCF) of moments obtained from MD trajectories. This package is used to compute the IR and VCD spectra from FFMD trajectories as the TCF function and the origin invariant version of the TCF function (C indep ) for the PBC case is already implemented.

From the time correlation functions to the spectra

The moments accumulated over the course of the trajectories were processed using the ChirPy program [START_REF] Jähnigen | Chirpy -A python package for analysing supramolecular and electronic structure, chirality and dynamics[END_REF], from which the time correlation functions and resulting IR and VCD spectra were determined. So far the equations given for the determination of the IR and VCD spectra through MD have been given in arbitrary units, as proportional to the Fourier transform of the appropriate TCF [equations (2.39) and (2.40)]. These correlation functions are classical, thus they are real and even, hence they obey C cl (-t) = C cl (t). The Fourier transform conserves parity, giving the following relation for the classical spectral line shape:

I cl (-ω) = I cl (ω).
However, the spectral line shape should satisfy the detailed balance condition:

I(ω) = e βℏω I(-ω) (2.87) with β = 1 k B T
, T being the temperature. To restore this condition, a frequency and temperature dependent prefactor is added: -ω). This prefactor is often referred to as a quantum correction factor. Various corrections exist in the literature, [START_REF] Egorov | Semiclassical approximations to quantum time correlation functions[END_REF][START_REF] Ramırez | Quantum corrections to classical time-correlation functions: Hydrogen bonding and anharmonic floppy modes[END_REF] the correction implemented in ChirPy is called the harmonic approximation :

I(ω) = Q QC (ω) e βℏω I(
Q HA QC (ω) = βℏω 1 -e -βℏω
(2.88) Hence the line shape determined by ChirPy has the following expression [7], using Gaussian-centimetre-gram-second (Gaussian-cgs) units:

I VCD (ω) = 8π(1 -e -βℏω ) 3Nℏcn(ω) βℏω (1 -e -βℏω ) ∫ +∞ -∞ ì µ(0) • ì m(τ) e -iωτ dτ (2.89) = 8π βω 3Ncn(ω) ∫ +∞ -∞ ì µ(0) • ì m(τ) e -iωτ dτ
with N being the number of units, n(ω) the refraction index of the medium and c the velocity of light. The prefactor is valid for isotropic systems at thermal equilibrium [8] and is constructed through Fermi's golden rule for a transition probability between two states. The time correlation function is taken at the classical limit and Kubo-transformed [7,8] to give equation (2.90).

It is now possible to obtain the line shape of a VCD spectrum from the electric and magnetic dipole moments determined during a force field molecular dynamics trajectory. The next step is to assess the conditions of the MD simulations needed to obtain correct spectra.

Assignment of simulation parameters: the case of trans-1-amino-2 indanol

AMOEBA force field parameters for trans-1-amino-2 indanol (trans-AI) had already been determined and tested previously [START_REF] Barbu-Debus | Assessing cluster models of solvation for the description of vibrational circular dichroism spectra: synergy between static and dynamic approaches[END_REF]. It is thus a good molecule to try and implement the present methodology. The first test was to compare an IR spectrum obtained from a trajectory with a long sampling time and one obtained with an average over trajectories with a shorter sampling time. The spectra from a simulation of 1 ns were compared to one from 10 simulations of 100 ps. The MD trajectories were performed with a time step of 0.2 fs with a temperature of 150 K. The IR spectra obtained from these trajectories are shown in figure 2.20.

We find that the the bands are generally defined with a higher number of peaks when averaging over multiple simulations especially just below 1600 cm -1 . This is because with multiple simulations multiple starting points can be chosen and also velocities are chosen at random in the beginning which improves sampling. So the simulation can explore a wider variety of configurations. The peak of the OH stretch, which is the vibration with the highest frequency for this molecule, is found experimentally at 3616 cm -1 . The blue spectrum using a time step of 0.2 fs reproduces this value. However as the value of the time step increases, a larger shift from this value appears. We chose the value of 0.5 fs as it is a good compromise between accuracy and computational time. Horníček et al [9] reached the same conclusion, explaining that the error between the expected and actual values of the peak corresponding to the OH stretch for long time steps, originates from the Beeman integration during the molecular dynamics trajectory.

The molecular dynamics simulations and thus the dipoles moments are calculated with a time step of 0.5 fs however the time correlation functions sample the dipole moments with a time step of 4 fs. Using a time step of 4 fs for the time correlation functions gives the same line shape as with a time step of 0.5 fs but with a reduced computational cost.

Such a short time step is only needed for determining spectra. For exploring the potential energy surface a longer time step (1 fs) is used as it permits a longer simulation time and thus a better exploration of the surface.

The next goal is to determine the vibrational modes that correspond to the bands in the spectrum.

Assignment of spectral modes

In section 2.3.2, the assignment of the frequencies in the harmonic approximation through the diagonalisation of the mass-weighted Hessian matrix was discussed. However, the interest here is for time-dependent, anharmonic systems. The correlation functions do not provide any information about the assignment of the frequencies, namely spectral modes. This information is essential to compare theoretical and experimental spectra in order for reliable conclusions to be drawn.

Methods to obtain vibrational modes

Within the harmonic approximation, it is straightforward to assign all modes through the diagonalisation of the mass-weighted Hessian matrix of static systems. For dynamical systems, that do not obey the harmonic approximation, various methods exist to determine normal modes that extend the concept of eigenmodes in harmonic systems.

The first method developed by Wheeler et al. [START_REF] Wheeler | Quasiharmonic vibrations of water, water dimer, and liquid water from principal component analysis of quantum or QM/MM trajectories[END_REF] is called principal modes analysis (PMA). The displacement of the mass-weighted molecular coordinates from their equilibrium position through a MD simulation is described in a covariance matrix. The eigenvalues and eigenvector of the inverse of this matrix are used to provide the frequencies and the normal modes.

A second method is the instantaneous normal mode analysis (INMA) where the massweighted Hessian matrix is diagonalised at each time step. However this can get expensive especially if the system is large as the diagonalisation, which is repeated, scales approximately as O(N 3 ) with N the number of atoms. Futhermore, INMA must deal with multiple negative eigenvalues. [START_REF] Nonella | Structures and vibrational frequencies of the quinones in Rb. sphaeroides derived by a combined density functional/molecular mechanics approach[END_REF][START_REF] Pejov | Using MD snapshots in ab initio and DFT calculations: OH vibrations in the first hydration shell around Li + (aq)[END_REF] Bowman et al. [START_REF] Bowman | Normal-mode analysis without the hessian: A driven molecular-dynamics approach[END_REF] proposed an approach called driven molecular dynamics (DMD). This method becomes competitive compared to INMA for systems with more than 1 000 atoms. In this method a sinusoidal driving force V DMD is added to the Hamiltonian. This term is determined for each step of the dynamics as a function of the distance between the nuclei:

V DMD (t) = N i, j λ i, j r i, j sin(ω n t) (2.90)
where r i, j are internuclear distances, ω n , called the driving frequency, is a fixed exciting frequency and λ i, j are the coupling constants. The driving frequency can be varied through a frequency range and energy absorption should maximize locally at the normal mode frequencies. An examination of the molecular motion at a given resonance provides the corresponding normal mode. DMD was coupled with the determination of the IR spectra through the Fourier transform of the autocorrelation function by Thaunay et al. [START_REF] Thaunay | Vibrational mode assignment of finite temperature infrared spectra using the amoeba polarizable force field[END_REF]. Hence it was possible to give the corresponding normal modes of a spectrum that includes anharmonicity and temperature effects from the trajectory. This is demonstrated in figure 2.22 where the displacements of C=O and N-H of N-methyl-acetamide are shown for two resonant driving frequencies. At ω n = 1723 cm -1 the normal mode corresponding to the C=O stretch is shown and at ω n = 3518 cm -1 it is the N-H stretching mode. The Kohanoff procedure [START_REF] Kohanoff | Phonon spectra from short non-thermally equilibrated molecular dynamics simulations[END_REF] extracts vibrational frequencies and eigenvectors from short non-thermally equilibrated molecular dynamics trajectories. The frequencies are estimated, then the eigenvectors are determined through a least squares fit of the trajectory including orthogonality constraints. The trajectory is projected onto each normal mode corresponding to one frequency.

Effective mode analysis

For DMD and the Kohanoff procedure, short MD trajectories need to be performed for each frequency, which is not the case for the last method presented here. It corresponds to an effective mode analysis [START_REF] Martinez | Extracting effective normal modes from equilibrium dynamics at finite temperature[END_REF][START_REF] Agostini | Infrared spectroscopy of small protonated water clusters at room temperature: an effective modes analysis[END_REF]. In this method, for each mode, an IR spectrum and the corresponding vibrations are determined from a trajectory of internal coordinates and velocities of the system. This method is useful if there is some prior knowledge of the physical chemistry of the system, which can be used to define the internal coordinates.

The total infrared spectrum I n (ω) can be written as:

I n (ω) = βN A 6cϵ 0 3N -6 i=1 ∂ ì µ ∂q i (0) 2 P (q)
ii (ω)

(2.91)
where, aside from the temperature-dependent prefactor, the index i runs over 3N -6 internal degrees of freedom of the molecule, ì µ is the dipole moment of the whole molecule, q i is the i-th effective mode (to be determined), and

P (q)
ii (ω) is the power spectrum, i.e., the velocity-velocity correlation function of the i-th effective mode

P (q) ii (ω) = ∫ +∞ -∞ dt e -iωt q i (t) q i (0) (2.92)
In equation (2.91) the IR spectrum is decomposed as the sum of independent contributions, each arising from a single effective mode i; in addition, the effective modes are constructed such that their power spectra are as localized as possible in frequency. This feature is achieved by minimising the functional

Ω (n) = 3N -6 i=1 β 2π ∫ +∞ -∞ dω ω 2n P (q) ii (ω) - β 2π ∫ +∞ -∞ dω ω n P (q) ii (ω) 2 (2.93)
with respect to the elements of the transformation matrix Z -1 from internal coordinates

ζ to effective modes q, namely q = Z -1 ζ or q i = 3N -6 k=1 Z -1 ik ζ k (2.94)
The functional Ω (n) represents the variance of the n-th moment of P (q)

ii (ω) when interpreting P (q) ii (ω) as a probability distribution.

In order to derive equation (2.91) for the IR spectrum starting from the Fourier transform of the dipole-dipole time correlation function, the following hypotheses are made:

1. Positions and velocities are uncorrelated all along the trajectory; 2. Rotations and vibrations are uncorrelated, and the rotational correlation function decays much slower than the vibrational correlation function;

3. The correlation function for the derivatives of the dipole moment decays much slower than the vibrational correlation function, thus we can only consider its value at t = 0;

4. The derivatives of the dipole moment with respect to q i can be determined from the atomic polar tensor in Cartesian coordinates, then transformed to internal coordinates, then transformed to effective modes;

5. The off-diagonal elements of the q-power spectrum are negligible if compared to the diagonal elements.

The VCD spectrum has peaks present at the same frequencies as the IR spectrum as there can be no VCD activity if there is no IR absorption. Thus by determining the modes for IR spectroscopy, the information can be transferred to VCD.

Concluding remarks

Throughout this chapter, various methods to compute vibrational circular dichroism spectra have been presented as well as procedures to assign molecular vibrations to spectral peaks. After giving the general theory of VCD, first a static approach was discussed using DFT and ways to include anharmonicity were presented. Then molecular dynamics, which naturally accounts for anharmonic and temperature effects, was presented and how to achieve VCD spectroscopy from this method was shown. Subsequently two approaches using MD were examined. The first approach using a quantum chemical description of the electrons, and the other describing the system with a force field. The implementation of VCD spectroscopy into force field MD was then explained, including the difficulties in periodic boundary conditions, how the spectra were retrieved from electric and magnetic dipole moments, and the importance of the time step of the MD simulation. It is this magnetic dipole moment that was implemented in a force field MD program during this work and that will be assessed in chapters 4 and 5, dealing with systems in the gas and condensed phases, respectively. Finally methods to determine spectral modes from MD simulations were presented.

Chapter 3

Trans-1-amino-2-indanol: a case study for VCD and fluxionality

Introduction

In experimental VCD, the molecules of interest are usually solvated and one fundamental issue is to assess the role of the solvent on the molecular conformation and on the spectrum. Successful assignment of VCD signals based on computed spectra therefore requires that the solvent is correctly described, and that the spectra are calculated with sufficient accuracy.

Sufficient sampling of the conformational space is also important to capture the effects on the spectra from the various conformations. These conformations can be found by chemical intuition, a random-search method for finding low-energy conformations like the Monte Carlo Multiple Minimum method [1] or through molecular dynamics simulations.

Systems of interest

Only little attention has been paid so far to the influence of solvation on the VCD spectrum of 1,2-amino-alcohols. Tarczay and coworkers demonstrated the suitability of matrix isolation VCD spectroscopy for determining the absolute configuration of fluxional molecules that have strong intermolecular interactions, especially hydrogen bonds [2]. Matrix isolation involves condensing the substance to be studied with a large excess of inert gas (usually argon or nitrogen) at low temperature to form a rigid solid (the matrix) [3].

The aim of this chapter is to unravel the interplay between the intramolecular structure and the solvation network, and explore how the solvation dynamics affects the VCD spectrum of 1-amino-2-indanol, depicted in figure 3.1.

The cis isomer (cis-AI) was explored employing a Monte Carlo Multiple Minimum method and the cluster-in-a-liquid model. This isomer was studied in reference [4]. For the trans isomer (trans-AI), polarisable force field molecular dynamics (FFMD) was used to extensively explore the potential energy surface of the monomer as well as that of its complexes with one and two DMSO molecules. Being designed for bulk DMSO, the AMOEBA FF used is expected to reproduce satisfactorily the solvation trends observed in the system. This isomer was studied during the present work. Representative structures were extracted and further optimised at the density-functional theory level to compute IR absorption and VCD harmonic spectra, in the frame of the cluster-in-a-liquid model, with the contribution from our collaborators Anne Zehnacker and Katia Le Barbu-Debus. Special attention was paid to the factors that influence the reconstructed spectra, in particular through Boltzmann weights associated with each cluster.

Furthermore, FFMD at fixed temperature was used to explore the potential energy surfaces of the trans-AI molecule solvated by one, two and five DMSO molecules, and compare them to the bulk limit, as part of this thesis. A picture of the solvent fluxionality is provided, which aims at understanding full solvation from the interaction of trans-AI with a limited number of DMSO molecules. The most stable trans-AI monomer and its 1:1 solvent-solute clusters with a single DMSO molecule were used as starting points for ab initio molecular dynamics simulations and calculations of the VCD spectrum by means of the nuclear velocity perturbation theory approach.

An example of the study of VCD spectra calculated using a cluster-ina-liquid model: cis-1-amino-2-indanol

In the following example, the influence of the solvent on the spectra is studied for (1 S ,2 R)-(-)-cis-1-amino-2-indanol (cis-AI) in DMSO [4]. The solvent is described with the cluster-in-a-liquid model, as described in the previous chapter.

The structures of cis-AI with two molecules of explicit solvent, namely the 1:2 complex, were optimised at the B3LYP-D3BJ/6-311++G(d,p) level of theory and VCD spectra were computed at the same level of theory. Both the structure and the corresponding VCD spectra are represented in figure 3.2. Two molecules of explicit solvent were added to explore their interaction with the two functional groups of cis-AI. All structures were calculated in a DMSO continuum solvent with a polarisable continuum model (PCM). Special attention was paid to the modes that are possibly involved in hydrogen bonding interactions with the solvent. These modes are represented with β and, in parenthesis, the group in the solute that is hydrogen bonded to the solvent. For example, the mode β(OH) is found at ∼1400 cm -1 and the mode β(NH 2 ) is represented at ∼ 1600 cm -1 .

The individual spectra, which are the same color as the frames around the corresponding structure, show varying spectral shapes at 1400 cm -1 and 1600 cm -1 depending on the solvent-solute interaction. The blue individual spectrum bears a strong resemblance to the experimental spectrum shown above. Its structure is also the lowest in energy. Thus it can be confirmed that the most likely structure of cis-AI corresponds to the structure shown surrounded by a blue frame in figure 3.2. This structure has an hydrogen bond between the OH group of cis-AI and the oxygen of a DMSO molecule. However, the cluster-in-a-liquid model involves a limited number of explicit solvent molecules and the anharmonic and temperature effects are only included through Boltzmann averaging. The next system, trans-1-amino-2-indanol, will be explored with the cluster-in-a-liquid method, as previously, though it will also be analysed with molecular dynamics methods.

Cluster-in-a-liquid approach

In a first step, the system is studied within the cluster-in-a-liquid model, using the same level of theory as the previous system. This level of theory has been shown to reproduce the VCD spectrum of similar molecules in solution, as well as their structure and vibrational spectrum in the gas phase [5,6,7]. It has also been used successfully for calculating the VCD spectrum of amino-acids in aqueous solution [8]. Solvent effects were taken into account by the IEFPM implicit polarisable continuum model [9]. Vibrational frequencies were computed at the same level of theory and scaled by 0.98 to correct for anharmonicity and basis set incompleteness. This value is close to that used for similar theoretical methods [10,11]. The final vibrational spectra were obtained by convoluting the harmonic intensities with a Lorentzian line shape (FWHM 4 cm -1 ).

The individual spectra of the most relevant 1:1 complexes are shown in figure 3.4, where one configuration from each family is also depicted. The structures resulting from the exploration of the potential energy surface can be classified into three fami-lies, namely, bidentate, monodentate and non-hydrogen-bonded complexes. Each family encompasses several conformations of trans-AI that will be included in the Boltzmannaveraged contributions when simulating the IR absorption and VCD spectra. In the bidentate family, the DMSO oxygen atom interacts with both OH and NH 2 groups of trans-AI. In the monodentate families, the DMSO oxygen interacts with either the OH (mono OH ) or the NH 2 (mono NH ) group of trans-AI. The last family is called non-hydrogen bonded (noHB) and accordingly contains structures lacking any intermolecular hydrogen bond. However, this type of structure with a loose DMSO molecule is not likely to be realistic. find their counterpart in the experimental spectrum. The weak positive doublet calculated at 1447/1431 cm -1 corresponds to that observed at 1462/1443 cm -1 . An intense band is predicted at 1384 cm -1 , which meets the position of the intense experimental feature. The negative band calculated at 1331 cm -1 may correspond to that observed at 1318 cm -1 . The last two calculated negative bands at 1216 and 1126 cm -1 are in particularly good agreement with those observed at 1217 and 1123 cm -1 . Other complexes show partial overlap with the experimental spectrum and can contribute to it as well. The VCD spectrum of the most stable bi-eqg + G -complex displays, like the eqg + G -monomer, the bisignate transition centred at 1616 cm -1 . Two positive features calculated at 1387 and 1345 cm -1 are close to the intense and broad experimental band at 1380 cm -1 and its shoulder at 1355 cm -1 . The negative signal at 1320 cm -1 is in good agreement with that observed at 1318 cm -1 . However, the negative features observed at 1228 and 1125 cm -1 are not predicted for the bidentate complex. They are due to coupled motions of the NH 2 and OH groups, which are strongly modified by the intermolecular hydrogen bond network and lose intensity in the complex. The mono OH -ax complexes also display satisfactory overlap with the experimental spectrum; in particular, mono OH -axtG -shows the bisignate signature of the G -forms and it becomes one of the most stable complexes when dispersion corrections are considered. The axial complexes all have negative VCD signals in the region of ∼1200 cm -1 , as experimentally observed. While the experimental features find at least partial counterpart in the spectra of the mono OH or the bidentate families, this is not the case for mono NH , which lacks in particular the strong positive band between 1340 and 1390 cm -1 . This poor match, together with very high relative energy, allows discarding the mono NH complex. bi-eqg + T cannot contribute to the spectrum due to the wrong signs of the bisignate at 1620 cm -1 and of all the bands between 1100 and 1250 cm -1 . Lastly, mono OH eqtG -exhibits a wrong sign for the doublet at ∼1450 cm -1 .

The Boltzmann-averaged VCD spectra of the monomer and the 1:1 complex are compared to the experimental spectra in figure 3.5. The calculated VCD spectrum of the 1:1 complex is in good agreement with the experimental spectrum, in particular in the 1350 cm -1 range. The calculated positive band at 1624 cm -1 corresponds to that at 1614 cm -1 in the experimental spectrum. The 1462 and 1443 cm -1 bands and the intense positive band at 1380 cm -1 with the shoulder at 1355 cm -1 are well reproduced. The negative features at 1217 and 1123 cm -1 also have their counterpart in the calculated spectrum. However, the band calculated at 1281 cm -1 has no equivalent in the experimental spectrum. The influence of the solvent can be seen around 1380 cm -1 with a negative then positive signal present in the experimental spectrum and the 1:1 complex spectrum but not in the monomer spectrum. This VCD band contains contributions of β(OH) and β(NH) bends.

Although the hydrogen bond pattern does not change upon inclusion of a second DMSO molecule, the energetic ordering is strongly modified compared to the 1:1 complexes. This observation points to the difficulty of defining the relative contribution of the different structures by their Boltzmann weights alone and the intrinsic limitation of describing solvation by a finite number of molecules.

In summary for the 1:1 complexes, an important role of the solvent is to stabilise the structures selectively, therefore changing the contributions of the different conformers in the spectra. In this respect, it is worth noting that axial forms are stabilised in DMSO complexes, in particular mono OH -axtG -can contribute to the final spectrum. Although almost no intra-molecular structural modification happens upon solvation, hydrogen bonding induces some band shifting and changes in the intensity or sign of the VCD spectrum in the corresponding regions. This is especially significant in the region of the β(NH) and β(OH) modes, just below 1400 cm -1 . Not all interaction types perturb the spectrum to the same extent. Non-hydrogen bonded complexes show the same spectrum as the corresponding monomer. bi-eqg + G -shows strong modification of the VCD Figure 3.5: Structures of the 1:1 complex, experimental VCD spectrum and calculated spectra of the monomer and the 1:1 complex. Adapted from Le Barbu-Debus et al. [12]. spectrum compared to eqg + G -in the 1200 cm -1 region, as mentioned above. The interaction with NH 2 perturbs the spectrum much less than the interaction with OH does: the VCD spectrum of mono NH -eqg + G + is almost identical to that of the corresponding monomer, while that of mono OH complexes is more affected by solvation.

Beyond cluster-in-a-liquid: exploration of the potential energy surface with molecular dynamics

MD simulations were performed with the AMOEBA force field and with AIMD to explore the potential energy surface for trans-AI with a varying number of DMSO molecules.

To characterise the position of the DMSO molecules with respect to trans-AI, the distances between the nitrogen or oxygen of trans-AI and the oxygen of DMSO were evaluated throughout the trajectory and drawn as 2D contour plots using the Plotly library in the R software [13,14]. The distance between the nitrogen atom of trans-AI and the oxygen atom of the DMSO is represented on the y-axis and the distance between the oxygen atom of trans-AI and the oxygen atom of the DMSO is on the x-axis. Then, the z-axis represents the occurrence of the distance observed over the trajectory in a false colour scale. Each graph is divided in four regions: (i) a spot at (x 3 Å, y 3 Å ) corresponds to the bidentate family; (ii) a spot around (x 3 Å, y 5 Å ) represents a monodentate configuration with the DMSO hydrogen-bonded to the OH group of trans-AI; (iii) a spot at (x 4 Å, y 3 Å ) corresponds to a monodentate configuration with the DMSO hydrogen bonded to the NH 2 group of trans-AI; (iv) finally, for a spot at (x ≥ 5 Å, y 8 Å ), the DMSO is in the vicinity of trans-AI without being hydrogen bonded.

AIMD approach

For isolated trans-AI, the 1:1 as well as the 1:2 complexes with DMSO, one, three, and one starting structures were created, respectively, based on the optimised geometries found from the static calculations. AIMD calculations based on DFT were carried out, by our collaborator S. Jähnigen, with the Quickstep module of the CP2K software package [15,16].

The simulations were of Born-Oppenheimer MD type with a time step of 0.5 fs, using the BLYP or B3LYP exchange-correlation functional [17,18], Grimme's dispersion correction (D3) [19], GTH pseudopotentials [20,21,22], and the Gaussian and plane wave basis TZVP-MOLOPT-GTH with an energy cutoff of 400 Ry [23]. All simulations were performed in the canonical ensemble using the CSVR thermostat [24].

The simulations were performed at a slightly elevated temperature of 340 K to counterbalance the underestimation of temperature by the chosen functional [25]. However, the effect of this temperature on the intermolecular sampling is marginal, as shown in figure 3.6.

The AIMD simulations were carried out in a vacuum supercell, that is each setup, either the isolated molecule or the complexes (1:1 or 1:2) with DMSO, was placed in a box of predefined size: 16 3 Å 3 and 20 3 Å 3 for the isolated molecule and the complexes, respectively. Each sample underwent geometry optimisation followed by a 5 ps equilibration performed under massive thermostatting with a coupling constant in the 10-500 fs range. The production trajectory of 30 ps was carried out under global thermostatting with a coupling constant of 500 fs.

AIMD trajectories were out of reach for 1:5 complexes, but could be carried out at 340 K for three 1:1 (trajectories b1 to b3) and one 1:2 complexes using the BLYP functional and at 320 K for two 1:1 complexes (trajectories d1 to d2) using the B3LYP functional. The starting structures contained DMSO interacting either in a bidentate manner or via OH•••O interaction. The BLYP-AIMD trajectories are short and produce results that depend more strongly on the starting structure, their analysis indicates that trans-AI spends 28% and 30% in the axial position for the 1:1 and 1:2 systems, respectively.

The influence of the functional on the trajectories is illustrated by comparing the results described above with similar MD trajectories obtained with the B3LYP functional (figure 3.7). For the isolated molecule, the B3LYP trajectories also evidence contributions from both axial and equatorial configurations. The main difference compared to BLYP results is an increased contribution of G -structures. For the 1:1 complex also, a balanced contribution of axial (41%) and equatorial configurations (59%) is obtained, as well as a larger contribution of G -structures relative to BLYP trajectories. In the 1:1 complex, the contour plot on figure 3.7 shows that both the bidentate and monodentate (OH and NH 2 ) conformations are visited. The DMSO molecule remains close to the OH group of trans-AI. With two DMSO molecules, conformational exploration is very limited, as manifested by the contour plots of figure 3.8 where the simulation mainly explores one type of interaction with DMSO. However, even under this short time scale the dynamics shows a partial mobility away from the initial conformation of trans-AI.

Force field based MD

Polarisable FFMD simulations resting on the AMOEBA force field were carried out to model trans-AI in interaction with one, two and five DMSO molecules on time scales much longer than affordable with AIMD. A set of multipoles was generated for the trans-AI isolated monomer using the distributed multipole analysis of the MP2/ cc-pVTZ electron density [26,27].

The temperature was set to 300 K or 150 K and was controlled using the Nosé-Hoover thermostat. The additional trajectories at 150 K were performed to constrain the dynamics of the system to remain around a selected local potential well of trans-AI. At 300 K, the duration of the trajectories was 1 ns and a time step of 0.5 fs was used. For the clusters with 2 and 5 DMSO molecules, a spherical van der Waals potential of 17 Å of diameter was used to avoid evaporation. At 150 K, the duration of the simulations was increased up to 3 ns to improve potential energy surface exploration but a time step of 1 fs was used to reduce the computational time. MD simulations of bulk systems were carried out using a cubic DMSO box with an edge length of 19.20 Å that contains 56 DMSO molecules. This box was pre-equilibrated before soaking the trans-AI solute. Simulations were performed at constant volume and 300 K using periodic boundary conditions, the Nosé-Hoover thermostat, and a 0.5 fs time step. Ewald summation was used for the longrange electrostatic interactions. The Ewald real-space cutoff was 7 Å, the van der Waals cutoff was 7 Å. Several simulations were carried out for a total simulation time of 4 ns.

Trajectories at 150 K

A local and detailed description of the conformational landscape around a given structure in a potential depth was obtained using 150 K FFMD trajectories.

Under such conditions, the temperature is low enough to catch the trans-AI molecule in individual potential wells, allowing the intrinsic stability of its different conformations to be explored. Four trajectories were carried out for trans-AI in the equatorial position: one for the 1:1 bidendate conformation that corresponds to the most stable structure in the static sampling, two for the 1:2 complex (one with one bidentate DMSO and the second DMSO in the vicinity, one as mono-OH and mono-NH conformation) and one with five DMSO molecules spread around trans-AI with two molecules hydrogen-bonded to the OH and NH 2 groups. With trans-AI in the axial position, six trajectories were simulated, including two with one DMSO (mono-OH and mono-NH, respectively), three with two DMSO (mono-OH and mono-NH, mono-OH and in the vicinity, mono-NH and in the vicinity, respectively) and one with five DMSO (with similar configuration of the solvent molecules to that of the equatorial system).

The contour plots (figure 3.9 and 3.10) show that the complexes remain near their starting points, keeping the hydrogen-bond interactions with one or two DMSO molecules. One exception is the simulation starting from trans-AI in the equatorial position with one DMSO molecule. In this trajectory, trans-AI switches to the axial position after 1.5 ns and remains axial during the rest of the simulation. This result contrasts with the static DFT results, which give only a few stable structures in axial configuration for the 1:1 complex. It could be explained by an artifact of the force field, which may excessively favour the axial configuration. In figure 3.11, the contour plots of the 1:5 equatorial complex show the five DMSO molecules staying mainly around their initial position. One molecule stays hydrogen bonded to the OH group, two are in interaction with the NH 2 group with one closer to the OH, and the remaining two molecules are in the vicinity of the trans-AI solute. These results are consistent with the main intermolecular interactions expected in this complex, i.e. OH• • •O and NH• • •O hydrogen bonds, and solvent-solvent interactions. The dynamics of the third DMSO molecule is more diverse. Lastly, the bidentate conformation is present but in a smaller amount and does not correspond to a well-defined spot in figure 3.11. Still, the superimposed contour plots for the five DMSO molecules show that all the expected interaction sites of trans-AI (protic groups and aromatic ring) interact with DMSO, which demonstrates that 1:5 complexes give a realistic description of solvation for this system.

Trajectories at 300 K

FFMD simulations at 300 K were performed to assess the sampling of a few solvent molecules around trans-AI. At 300 K, the FFMD trajectories are less impacted by the starting configuration. Four trajectories were produced, namely two with one DMSO, and two with either two or five DMSO molecules, respectively. The dihedral angle of the 5membered ring was monitored to characterise the equatorial and axial conformations (see figure 3.12). In the FFMD trajectory for the 1:1 complex, trans-AI spent most the time in axial position (93%). However the results are more balanced for the system with two (80%) and even more so for five DMSO (45%). For the amino group, the structure with the G + orientation is clearly favoured with FFMD, however at the DFT level it is the one with the G -orientation that has the lowest energy with explicit solvation. distances in the 1:1 complex at 300 K in the FFMD simulations.

In the 1:1 complex (figure 3.14), the DMSO molecules mainly takes three positions, near the OH and NH 2 groups of trans-AI (red spot) and non-hydrogen bonded to trans-AI. Contrary to the static DFT calculations, the bidentate interaction is not favoured because trans-AI is mainly in the axial position during the simulation. Two starting points were exploited, one with the solvent hydrogen-bonded to the OH group and the other with the solvent on the NH 2 group. Yet, the two trajectories give similar contour plots.

With two DMSO molecules as explicit solvent, one DMSO molecule stays mainly hydrogen-bonded to the OH or the NH 2 group of trans-AI (figure 3.15, right). The second DMSO molecule takes mainly three positions, hydrogen-bonded to the OH or the NH 2 group of the trans-AI, or non-hydrogen bonded in the OH region (figure 3.15, left). At 300 K, the DMSO molecules are much more mobile than at 150 K, moving around the oxygen and the nitrogen atoms of trans-AI. Figure 3.16 provides the distances between the OH and NH 2 groups of trans-AI and DMSO as a function of time, showing the fluctuations of the two DMSO molecules moving close, then far away from trans-AI. Exchanges between the two molecules happen several times on the NH 2 group with an approximate residence time, i.e. lifetime of the hydrogen bond, between 100 and 200 ps. Nevertheless, the contour plots indicate that there is a tendency of each DMSO molecule to remain in certain spots throughout the simulation. FFMD simulations were also performed with five DMSO molecules as explicit solvent. The solvent moves around the trans-AI with an extensive sampling for all the molecules. The contour plots show similar preferred spots as for the 1:1 and 1:2 complexes: near the OH and NH 2 groups of trans-AI and non-hydrogen bonded in the OH region (figure 3.17). The bidentate interaction for one DMSO molecule is also found and can be correlated to a larger proportion of the equatorial conformation. Two main conclusions can be drawn from these graphs: (i) at short distances, one DMSO molecule is always hydrogen-bonded to the OH and NH 2 groups, although this interaction is dynamical with frequent exchanges, leading to a residence time of ca. 70 ps and 100 ps, respectively; (ii) the other molecules (three or four, depending on the observed conformation) solvate the trans-AI solute and are also in its vicinity. These observations point to the mobility of the solvent at room temperature. The RDF between trans-AI and the oxygens of 56 DMSO molecules are shown in figure 3.20. The RDF between the oxygen of trans-AI and the DMSO molecules (red) show an intense peak at 2 Å. The RDF between the nitrogen of trans-AI and the DMSO molecules (blue) show an intense peak at 3 Å. The first plateaux of the integration of the RDF are at 0.9 and 1.7 for the oxygen and the nitrogen of trans-AI, respectively. This confirms the hydroxyl group having nearly always an hydrogen bond with a DMSO molecule at 2 Å and the amino group mostly having two hydrogen bonds with a DMSO molecule at 3 Å. It also suggests the bidentate family being in a small minority. The second solvation shells are around 6 Å from the oxygen of trans-AI and 5.75 Å from the nitrogen of trans-AI. In terms of solvation, the results of the AIMD simulations for the 1:1 and 1:2 complexes parallel those of FFMD. They also highlight the mobility of the DMSO molecules (figure 3.7) and the preference of the DMSO for the OH group of trans-AI, which parallels the conclusion of static calculations where the mono OH or bidentate complexes are the most stable structures.

Distance between Nitrogen of the trans-AI and the

Locality of vibrational absorption and VCD and the role of the solvent

To further explore the effect of the solvent, a radial cutoff function for the moments was applied as it has been used and introduced in the previous chapter and used in earlier studies [5,28,29]. The electric and magnetic dipole moments entering the TCF were thus evaluated according to their originating position in space. The spatially resolved IR and VCD spectra derived from AIMD-NVPT calculations based on the B3LYP trajectories are shown in figure 3.21. This study of the dipole moment was carried out by our collaborator S. Jähnigen. In figure 3.21a, the radius r corresponds to the distance from the centre of mass of trans-AI as it is scanned starting from 0, where only local contributions to the IR/VCD spectrum appear, i.e. those contributions that stem from trans-AI itself. By moving farther from the chiral solute, non-local signatures can be captured; for trans-AI this corresponds to contributions stemming from DMSO. Hence, if a non-local signal can be found, the polarisation of the centre (trans-AI) is coupled with the environment (DMSO), which can be due to either coupled oscillations, or induced polarisation [30]. It becomes clear that the spatial region above 2 Å, denoted as "solvent effect", hardly returns any signal that may account for intermolecular IR/VCD; the main contributions are indicated as being of local origin, that is, stemming from trans-AI itself. Consequently, DMSO as a solvent does not significantly contribute to neither the IR nor the VCD spectrum of trans-AI in solution. There are negligible traces at about 1360 cm -1 and 1400 cm -1 in the IR spectrum and, additionally, at 1425 cm -1 in the VCD spectrum, where the large dipole moment of the solvent slightly polarises the solute. This does not come as a surprise since these frequencies correspond to modes localised at the OH and NH 2 groups, respectively, which bind the molecule to the solvent. Although in these regions the non-local contributions sum up to 20 %, they do not change qualitatively or quantitatively the overall shape of the spectra (figure 3.21b). Consequently, what matters is the conformational sampling of trans-AI and this is where solvent effects are discernible.

Such a conclusion is important as it means that in the case of DMSO electrostatics do barely influence the VCD response of trans-AI directly. This is at least true for the case of a 1:1 cluster of trans-AI and DMSO, on which the AIMD results rely. Yet, as the solvation study suggests, it can be assumed that this small-scale image generalises to full solvation. The number of solvent molecules surrounding the solute may ultimately lead to a rise of the small contributions shown in figure 3.21, but since the nature of the interaction does not change significantly with the number of DMSO molecules, it is unlikely to alter the locality of the VCD spectrum. Consequently, only the electronic wavefunction of the chiral solute is needed, which encourages the use of theoretical models like QM/MM approaches for IR and VCD determination.

Concluding remarks

The chosen system, (1S,2S)-trans-1-amino-2-indanol in DMSO, is especially challenging for VCD spectroscopy. Because of its flexibility, many structures coexist within tiny energy differences that strongly depend on the method used. Static calculations based on the cluster-in-a-liquid model capture the main features of the experimental VCD spectra for the 1:1 complex. Conformations higher in energy that are explored in FFMD or AIMD, like axial conformations, may contribute to the spectrum more than their static energy suggests.

The solvent influences only indirectly the calculated spectra through structural effects because of the locality of the VCD signal. This locality is an important finding that further supports the interpretation of vibrational spectra from the sole perspective of the solute conformational space and justifies the use of QM/MM calculations, with only a cluster of limited size being treated quantum mechanically. As with any MD method, ergodicity is an important issue. The limited duration of the AIMD trajectories, due to their computational cost, makes it important to define properly the starting structure. In this respect, combining FFMD for an exhaustive exploration of the PES with AIMD for the calculation of the VCD spectra proved to be efficient. The results obtained here stress the importance of solvent mobility onto the results, which is hardly taken into account by small finite-size clusters. Therefore, a thorough FFMD study seems essential to fully understand VCD for solvated compounds, as a complementary tool to the static cluster picture. This is why development of polarisable force fields amenable to very large sizes, close to full solvation, is necessary. No spectrum was obtained here with such a method as it was still in development when this study was in progress. It is the object of the next chapters where the determination of such spectra is assessed and demonstrated.

Chapter 4

Validation of the polarisable force field approach for VCD spectroscopy: alanine in the gas phase Classical molecular dynamics performed with a force field (FFMD) usually achieves a much better sampling than ab initio methods (AIMD), because trajectories can be propagated over significantly longer times and a larger number of initial conditions. This ability stems from the full classical description in a force field. This differs from AIMD, in which the electronic structure problem is solved quantum mechanically through DFT.

The present chapter aims to discuss and validate our implementation of the strategy described in Chapter 2 for the computation of VCD spectra with FFMD. The AMOEBA force field is designed to simulate organic molecules and biomolecules. Thus it can describe amino acids which are used as building blocks to characterise proteins. The simplest amino acid which is chiral is alanine. Thus isolated alanine was investigated as our benchmark system in the zwitterionic and neutral forms (see fig. 4.1 a and b). The atom numbering will be used in section 4.4. In the gas phase, alanine has been reported to be more stable in the neutral form though once it is hydrated or in a crystal, the zwitterionic form becomes more stable [1,2]. Initial tests were performed on the zwitterionic form as the corresponding force field parameters were already available. The calculations in the neutral form were performed later, once the parameters for this form had been determined. Various results for the vibrational spectroscopy of isolated alanine have been reported by other groups, from the theoretical [1] and experimental [3] perspectives.

In this chapter, the IR and VCD spectra of isolated alanine were determined with FFMD and the various contributions to the magnetic dipole moment were examined, as well as the time derivative of the electric dipole moment. The quality of sampling was quantified through the determination of statical errors and convergence rates of the spectra. Various configurations were determined from an FFMD simulation and reoptimised in DFT to assess the correctness of the FFMD results obtained. The DFT results were also used to evaluate the importance of anharmonicity. Moreover, to validate the dipole moments acquired from the FFMD, they were compared to dipole moments determined with DFT at each step of the same trajectory. Finally effective modes for neutral and zwitterionic alanines were extracted from respective FFMD trajectories. The various contributions from each functional group were also examined.

Vibrational spectra obtained from force field molecular dynamics (FFMD)

The idea of this section is to obtain IR and VCD spectra from FFMD trajectories. To ensure that the configurations used to determine the electric and magnetic dipole moments are as diverse as possible, different configurations were used to initiate a large number of trajectories. As shown in chapter 2, the resolution of the spectrum is improved by sampling a lot of smaller trajectories rather than one long trajectory. To collect various configurations to be employed as starting points, a long trajectory of 1 ns was achieved in the NVT ensemble. From this trajectory, 40 starting configurations were taken every 10 ps for the last 400 ps. For each starting point, a 5 ps NVT simulation was performed to equilibrate the system at the correct temperature then a 200 ps NVE simulation was carried out, from which the electric and magnetic dipole moments were collected. For each NVE trajectory, an individual spectrum was calculated through the Fourier transform of the time correlation function of the dipole moments, as described in chapter 2. Then the individual spectra were averaged to obtain a spectrum of an hopefully well sampled system.

For neutral alanine, the multipoles were not available in the AMOEBA force field. Thus the atomic multipoles were extracted from ab initio calculations thanks to the distributed multipoles analysis developed by Stone [4]. These multipoles were added to the bonded, van der Waals and polarisation parameters all ready included in AMOEBA. However the lowest configuration shown in literature [5] was not represented in an initial NVT ensemble FFMD trajectory. The force field parameters were not describing the NH 2 group with sufficient accuracy. Thus the bonded parameters of the NH 2 group were changed to those of similar functional groups already present in the force field. These trials yielded only a slight improvement. Furthermore, the correlation of the relative energy between DFT and FF results was improved by varying the charge on the oxygen of the OH group. The other charges of the system were reequilibrated to maintain the molecular charge at zero.

Calculating the electric and magnetic dipole moments for the IR and VCD spectra

The equations of the electric ì µ and magnetic ì m dipole moments are recalled below. The first term of electric dipole moment ì µ is expressed here as a function of the position, thus this electric dipole moment corresponds to the position form. If the electric dipole moment is differentiated with respect to time, ì µ = ∂ ì µ/∂t, the first term of the result becomes a function of the velocity, leading to the so-called velocity form of the electric dipole moment. The magnetic dipole moment is composed of 3 contributions. By ignoring any term with the induced dipole, ì µ ind i , the polarisation of the atoms is neglected. The term corresponding to the time varying induced dipoles moving around fixed positions can also be removed to explore the influence of this term. Second, the VCD spectra including a varying number of contributions from the magnetic dipole moment are shown in figure 4.2a. The blue spectrum depicts the simplest expression of the magnetic dipole moment combined with the electric dipole moment which is only composed of its first term. If the induced dipole is added to the electric dipole moment and the second term is added to the magnetic dipole moment, the orange spectrum is obtained. This curve differs from the blue one by peak inversions at 1150 cm -1 and 1700 cm -1 . The black spectrum arises from the cross-correlation from the full expressions of the electric and magnetic dipole moments. Compared to the orange spectrum, there are peak inversions at 1150 cm -1 , 1500 cm -1 , 1590 cm -1 , 1700 cm -1 , 1850 cm -1 and 1925 cm -1 .

ì µ(t) = N i q i ì r i (t) + ì µ ind i (t) (4.1) ì m(t) = N i 1 2c q i ì r i (t) × ì v i (t) + 1 c ì µ ind i (t) × ì v i (t) + 1 2c ì r i (t)× ì µ ind i (t)
Moreover, the position and velocity forms of the electric dipole moment are compared for the VCD spectra in figure 4.2b. They confirm that there is no influence on the spectrum between the two forms. The velocity form is chosen as the time correlation function (TCF) involving µ decays faster than that involving µ, thus facilitating convergence [6]. In the rest of the thesis, the velocity form of the electric dipole moment and the term µ × r i will be used.

These results show the importance of the inclusion of the polarisation, which is performed here through the induced dipole moments within the AMOEBA force field.

By comparing our spectra to the results obtained in the literature for the zwitterionic and neutral forms [1,7] the peaks above 1600 cm -1 could correspond to the C O stretch of the carboxylate or carboxylic acid and the bend of the amino or amine groups. The bending of the CH 3 is expected around 1450 cm -1 and the bending of the CH is predicted at 1300-1200 cm -1 . Below 1200 cm -1 the presence of the C-N-H bend can be reasonably assumed.

Statistical errors and convergence rates

The spectra of alanine shown above are the result of the average of 40 trajectories. To show the convergence of both forms with the sample size, the averages of an increasing number of trajectories (N T ) are plotted. Convergence is achieved when the intensities of the peaks stabilise at a certain value after a certain sampling time. Starting from the top of figure 4.3, the spectrum arises from the average of 4 trajectories (N T = 4), then the number of trajectories used increases by 4. Each spectrum is compared to the average of 40 trajectories shown in grey.

The spectra for the zwitterionic form converge faster as the charges create stronger electrostatic interactions. Also the IR spectrum converges faster than the VCD spectrum as VCD is more sensitive to a change of configuration of the molecule. Hornivcek et al. [8], found that 100 MD simulations of 20 ps are needed at 300 K to obtain convergence for α-pinene. Similar results are found here with the convergence starting at 10 trajectories of 200 ps. The average of the spectra from 40 FFMD trajectories is considered converged, thus the spectra of N T < 40 are compared to this averaged spectrum so as to determine their normalised error. The idea is to represent quantitatively the convergence with an increasing number of trajectories. For the IR spectra the standard error method is employed to determine the error between the converged spectrum and the others. For the VCD spectra, the SimVCD method [9] is used, which computes a quantitative similarity between two VCD spectra. With these methods the normalised error can be determined with the following equations:

E IR = A 12 (A 11 + A 22 ) (4.3) E VCD = 1 - Z 12 (Z 11 + Z 22 -|Z 12 |) = 1 -SimVCD (4.4)
where A i j is the quadratic difference between the IR spectra i and j, A ii their quadratic norm, Z i j the overlap integral between the VCD spectra i and j (also for i = j):

A i j = ∫ |I i (ω) -I j (ω)| 2 dω, ( 4.5
)

A ii = ∫ I i (ω) 2 dω, ( 4.6
)

Z i j = ∫ VCD i (ω) × VCD j (ω)dω, (4.7
)

Z ii = ∫ VCD i (ω) 2 dω (4.8)
Where I i (ω) is the IR spectrum intensity at frequency ω and VCD i (ω) the VCD spectrum intensity. Here the averaged spectrum that is considered converged, the so-called complete spectrum and the average of N T randomly selected trajectories, the so-called partial spectrum are compared. Y C (ω) denotes the complete spectrum at frequency ω and Y P (ω) denotes the partial spectrum at frequency ω. As a discrete number of frequencies is taken, the integrals are converted into the sums shown in equations (4.9) and (4.10).

As mentioned previously, N T trajectories are selected randomly for the partial spectrum. This allows us to use a Monte Carlo process to determine the error of the normalised errors. This Monte Carlo process relies on repeated random selections to obtain statistical results. Thus these normalised errors also have their own statistical error calculated with equation (4.11), giving the fluctuations over m = 10 000 trials.

E IR = ω i (Y C (i) -Y P (i)) 2 ω i Y 2 C (i) + ω i Y 2 P (i) (4.9) E VCD = 1 - ω i (Y C (i) × Y P (i)) ω i Y 2 P (i) + ω i Y 2 C (i) - ω i (Y C (i) × Y P (i)) (4.10) Fluc(E = E IR or E VCD ) = m i E 2 m - m i E m 2 (

Comparison with electronic structure calculations 4.2.1. Structures and relative energies

Each configuration of an FFMD simulation of 1 ns of alanine in the neutral form was optimised by the force field and 6 main configurations could eventually be identified from this analysis. A 7 th configuration was added, which corresponds to the lowest energy configuration from the article of Blanco et al. [5]. This configuration turns out not to be locally stable with AMOEBA, despite efforts to improve the parameters of the NH 2 and COOH groups. This is probably due to the NH 2 parameters not being fully optimised. Only the neutral form was studied as the zwiterrionic form is unstable: during the optimisation, the hydrogen moves back to the oxygen thus the neutral form is recovered. Rahmani et al. [1] argued that the change of structure is due to intramolecular hydrogen bonding between the hydrogen of the amino group and the oxygen of the carboxylic group followed by a rotation of the COO -group. These configurations were further classified into families based on their structural similarities, see Table 4.1. The configurations within each family differ mainly through the rotation of the CH 3 group. The families differ in the position of the OH group compared to the NH 2 group and the positions of the nitrogen hydrogens compared to the main carbon chain.

DFT calculations at the B3LYP/6-311++G(d,p) level of theory were performed for each of the 7 configurations using the Gaussian 16 software [10]. The energy of each configuration was then recalculated with the force field to compare the values of the relative energy and the electric dipole moments.

For each family, the root-mean-square deviation (RMSD) was calculated taking the configuration of family 1 as a reference. For the other families, the structure used to determine the values of the RMSD, the relative energy and the permanent dipole moment corresponds to the structure with the lowest energy of the family. Families 1 and 3 have a similar orientation of the COOH group, thus their RMSD are closer than for the families 1 and 2. The energy of each family is close, this is mainly due to the hydrogen bond always being present between NH 2 and an oxygen. The charges of the force field are described with a set of point charges, dipoles and quadrupoles centered on each atom, which give the electric dipole. In DFT, the dipole moment is defined through the electronic density. Thus the dipole moments differ as the distribution of charges are not considered in the same way. However the values of the dipole moments follow the same trend for DFT and the force field, which gives a first validation of the dipole moments. Table 4.1: The 3 families of configurations found for neutral alanine at the DFT/B3LYP/6-311++G(d,p) level of theory. Also included is the RMSD value compared to the 1 st configuration, the relative energy (∆E) calculated with DFT and with the force field (in kcal/mol) and the permanent electric dipole moment (µ tot ) calculated by the force field and DFT (in Debye).

IR and VCD spectra from static DFT calculations

For each of the families shown previously, IR and VCD spectra were determined. The harmonic frequencies were scaled by 0.98 to correct for anharmonicity and basis set incompleteness [11,12]. The anharmonic spectra were obtained using the second order perturbation theory (VPT2) described in chapter 2. No scaling factors were applied to anharmonic frequencies. The final vibrational spectra were obtained by convoluting the intensities with a Lorentzian line shape (FWHM 10 cm -1 ).

The comparison of spectra from the harmonic and anharmonic frequencies, in figure 4.6, shows that anharmonicity mainly affects the spectra at 1300 cm -1 with a band inversion in the VCD spectrum for all 3 configurations. The three VCD spectra show different signals at 1800 cm -1 , this is due to the different positions of the C=O group compared to the CH 3 group. The bands of the NH 2 scissor vibration have the same sign but the relative intensity is higher when the NH 2 group is aligned with the C=O group. Families 1 and 3 have the same sign for the bands around 1100 cm -1 .

The values of the frequencies of the most intense bands in the IR spectra and their assignments are reported in Table 4.2. The assignments were determined by Rahmani et al. [1]. For each value of frequency reported here, the vibrational modes were examined and similar values of frequencies are found for each type of vibration. This is not surprising as very similar methods are used to determine the frequencies and their assignments. Table 4.2: A comparison of harmonic scaled frequencies (in cm -1 ) with different basis sets and experimental frequencies (in cm -1 ) and their assignments. Rahmani et al. [1] used a 0.97 scaling factor whereas a scaling factor of 0.98 is used here as the level of theory differs. As the energy difference between each family is low, the assumption is made that they coexist. Thus each configuration needs to be included to form a realistic spectrum. For this purpose an average over the spectra from the 7 conformations was performed. To make the spectra more realistic, weights were used to give more importance to the spectra from configuration that could exist with a higher probability. These weights can be determined in two ways. The first uses the MD data, where the number of times each configuration was present in the 1 ns trajectory is considered. These will be called the MD weights. The second way relies on Boltzmann weights and requires, for each configuration α, the values of the partition function Z representing the thermal populations at canonical equilibrium. These will be called the DFT weights. In the harmonic approximation, the partition function Z is defined as the following:

Z α = e -βE α κ i=1 e -βhν i α /2 1 -e -βhν i α (4.12)
To account for intrinsic anharmonic effects, in equation (4.12) the (unscaled) anharmonic frequencies, ν i α , i = 1,..., κ = 3N -6, determined from our VPT2 calculations is used. β corresponds to 1/k B T and h is Planck's constant.

The average of the 7 spectra using the DFT and MD weights are presented in figure 4.7. The FFMD spectrum shown previously is also represented.

The average from the MD data corresponds mainly to the spectrum of family 3 (94.3%) and family 2 (5.7%) is present as a minority in the FFMD simulation. The weight of family 1 is zero as its structure is not locally stable with AMOEBA. The average from the DFT data has family 1 being the main contributor (49%) followed by family 2 (30.5%) and finally family 3 (20.5%). This implies that the force field and DFT favour the hydrogen from the carboxyl group positioned far from the nitrogen. The amino group is favoured by DFT to be perpendicular to the carboxyl group, whereas the force field favours an aligned configuration. The parameters of the NH 2 group are probably not optimal. Trials were performed to improve these parameters by varying the values corresponding to the angle, torsion and multipole terms. However only a slight improvement was noticed.

In the region 1600 cm -1 -1800 cm -1 , the bands from the DFT spectra are assumed to correspond to the peaks at 1700 cm -1 from the FFMD results. The sampling from the FFMD leads to broader peaks due to anharmonicity and temperature effects. The bands from the DFT results have been artificially widened through a convolution of the intensities with a Lorentzian line shape. However the FFMD and DFT are still quite different. This originates from the exploration of different configurations. As shown in figure 4.6, the spectra from each of family are very different. Neutral alanine is very flexible and is composed of a large number of configurations that each method explores differently.

Assessing moments and spectra from one trajectory

The electric and magnetic dipole moments determined with force field molecular dynamics are compared to the predictions of ab initio MD. For a robust comparison, a single MD trajectory (produced by the force field) is exploited. To be able to compare the electric and magnetic dipole moments of the two methods, they need to be determined from the same trajectory. The trajectory includes the positions and velocities throughout the MD simulation. Thus for both the zwitterionic and neutral forms, a trajectory of 200 ps at 300 K was determined with FFMD giving the dipole moments. The dipole moments were then computed again using now DFT within the CPMD software, by our collaborator S. Jähnigen. As the trajectory has already been determined, the DFT method using CPMD is equivalent to performing a DFT calculation for each time step, where the Troullier-Martins pseudopotentials [13] are used and the plane wave basis has a cutoff of 70 Ry. The trajectory was sampled with a time resolution of 4 fs.

From these dipole moments, the IR and VCD spectra can be determined using the Fourier transform of the corresponding time correlation functions. The spectra from both of these methods are shown in figure 4.8. For the zwitterionic form (Fig. 4.8a), the force field gives rise to more activity in the 1410-1450 cm -1 range and a lower intensity at 1925 cm -1 . The change in intensity can come from the difference in the description of the charges on the carboxylic acid. In the VCD spectra there is an inversion of sign at 1050 cm -1 . For the neutral form (Fig. 4.8b), the VCD signs are all the same, only the intensities vary in the IR and VCD spectra. For the electric dipole moment of the zwitterionic form of alanine, the signal is more disperse along the y axis. This is because there is a slight rotation of the molecule in the trajectory. This dispersity does not appear for the time derivative of the electric dipole moment and the magnetic dipole moment as the distribution of these dipole moments are centrosymetric. No rotation is shown in the dipole moments of the neutral form of alanine.

The value of r is high for the electric dipole moment and satisfactory for its time derivative and the magnetic dipole moment of the zwitterionic form. These values are improved for the neutral form.

-1 The results are quite similar between the DFT and FF approaches, which is promising for our new method. The neutral form shows a better correlation with DFT than the zwitterionic form, this is due to the description via point charges of the force field. Despite the charges being expressed by multipoles and by induced dipoles, this description is disadvantaged for the -COO -and -NH 3 + groups. Keeping in mind that the force field was parametrized for solvated molecules, a better description can be expected for the hydrated system compared to the isolated molecule.

0 1 2 3 µ DFT α (a.u.) -1 0 1 2 3 µ FF α (a.u

Assignment of the IR spectra into effective modes

The dipole moments and the vibrational spectra were obtained with good accuracy. In this section, the IR spectra are decomposed with the effective modes analysis method explained in chapter 2. The information of the vibrational modes is not included in the time-correlation functions, thus this analysis is needed to attain information about the origin of each peak. The effective modes were determined by analysing one trajectory (positions and velocities) of 200 ps at 300 K from FFMD. A quick attribution of peaks was fulfilled in sections 4.1.1 and 4.2.2 by comparing our results with the DFT calculations of alanine and with literature [1,7]. Here a more detailed description is provided.

Analysis of the zwitterionic and neutral forms of alanine

Using the configuration and atom numbering of figure 4.1, the internal coordinates reported in Table 4.3 were used to initialise the effective modes analysis.

bond lengths valence angles dihedral angles

N1-C2 C3-C2-N1 O4-C3-C2-N1 C2-C3 O4-C3-C2 H5-N1-C2-C3 C3-O4 C2-N1-H5 H8-C2-N1-H7 N1-H5 C2-N1-H6 C9-C2-N1-H7 (neutral: O4-H5) N1-H6 C2-N1-H7 H10-C9-C2-N1 N1-H7 N1-C2-H8 (neutral: H5-O4-C3-O13) C2-H8 C9-C2-N1 C2-C9 H10-C9-C2 C9-H10 H11-C9-C2 C9-H11 H12-C9-C2 C9-H12 O13-C3-C2 C3-O13
H5-N1-H6 (neutral: none) H7-N1-H6 H12-C9-H11 H10-C9-H11 O13-C3-O4 Table 4.3: List of the 33 internal coordinates used to initialise the effective modes analysis of alanine in the zwitterionic form (12 bond lengths, 16 valence angles, 5 dihedral angles). The labels of the atoms refer to Fig. 4.1(b). For the neutral form, the bond N1-H5 is replaced by the bond O4-H5, the valence angle H5-N1-H6 is suppressed, and the dihedral angle H5-O4-C3-O13 is added.

In figure 4.11, the IR spectrum from the previous FFMD simulations is shown and the peak for each effective mode is represented underneath, these are the power spectra of the effective modes. The effective modes analysis gives, in addition to a local spectrum for each mode, a representation of the corresponding vibration using vectors that is shown in figures 4.12 to 4.14. Note that the power spectra are shown without considering the | ∂ ì µ ∂q i (0) | 2 prefactor, from equation (2.91) in chapter 2, which would alter the relative intensity of the peaks.

The functional Ω (n) was tested with n = 2 and n = 4 and no major differences were found, thus the results are reported only for n = 2 in both zwitterionic and neutral forms of alanine. Note that the minimisation procedure may depend on the choice of initial internal coordinates ζ. Also in this case various options were tested and only one set of results is reported below (no major differences were found among the various tests performed).

Analysis of the zwitterionic form

In the low frequency region of Fig. 4.11, four bands can be clearly identified: two overlapping bands in the region 1000-1100 cm -1 , one band in the region 1100-1200 cm -1 , and one band localised at 1300 cm -1 . The four bands see the participation of the CH 3 , NH 3 and CO 2 groups, together with contributions from the vibration of the CH bond. This finding is confirmed in figure 4.12, where the COO motion is seen to participate to modes 1 and 4. The region 1350-1500 cm -1 of the spectrum in figure 4.11 is populated by four main bands, associated to the 4 modes shown in figure 4.13. This region of the spectrum sees the main participation from the bending and umbrella modes in the two fragments CH 3 and NH 3 . Some distorted COO bending motions contributes as well to the higher frequency modes (modes 7 and 8 in figure 4.13), giving rise to the overlapping band slightly above 1450 cm -1 in figure 4.11. In the remaining part of the spectrum of zwitterionic alanine in the gas phase, five modes are identified between 1500-1900 cm -1 . The high energy band at 1900 cm -1 is easily assigned to the asymmetric stretching motion in the COO fragment. This is mode 13, represented in figure 4.14.

ω5=1425 cm -1 ω6=1441 cm -1 ω8=1465 cm -1 ω7=1467 cm -1
At 1700 cm -1 , the band corresponding to mode 12 in figure 4.14 is identified. This band is produced by the excitation of the bending and the umbrella motion in NH 3 (see figure 4.18 for comparison), which is now coupled to the bending motion in the CH 3 fragment as well as to the vibration of the CH bond which is not involved in CH 3 . The three remaining bands are located in the region 1500-1600 cm -1 , and are mainly produced by the bending and umbrella vibrational modes of the NH 3 fragment. These bands and these effective modes are labelled 9, 10 and 11 in figure 4.14.

ω9=1538 cm -1 ω10=ω11=1580 cm -1 ω12=1702 cm -1 ω13=1919 cm -1
The IR data reported in reference [3] on the zwitterionic and neutral forms of alanine is compared to the results of our effective modes analysis in terms of frequencies and bands assignments. Note, however, that in reference [3] the bands assignments were performed based on harmonic normal modes calculations with electronic energies calculated at the Hartree-Fock level. The modes assignment is presented in the central column, where the notation for the modes used in reference [3] is adopted. In Table 4.4, the results for the zwitterionic form of alanine are reported.

Experimental frequencies (cm -1 )

Modes

Theoretical frequencies (cm -1 ) 1014 4.4: Left column: experimental frequencies collected from Ref. [3] in the spectral region 1000-2000 cm -1 for the zwitterionic form of alanine. Central column: bands assignment using the abbreviation from reference [3] (ρ for rocking, δ for bending, ν for stretching, ω for wagging; s for symmetric and as for asymmetric). Right column: theoretical frequencies from this work.

Most of the modes are represented in both experimental and theoretical results with similar frequencies. However for a few modes, the experimental and theoretical frequencies differ from one another. These modes are the NH 3 symmetric bend, the COO symmetric stretch, the CH 3 antisymmetric bend and the COO antisymmetric stretch. The CH 3 and NH 3 antisymmetric bend and wag are not described by theory.

Analysis of the neutral form

As expected from the different structures, the IR absorption spectrum of neutral alanine in the gas phase is very different from that calculated for the zwitterionic form. By comparing figure 4.11 and figure 4.15 some major differences can be pointed out: (i) in the region 1000-1300 cm -1 , the low intensity band at 1150 cm -1 in zwitterionic alanine becomes more intense and splits into two contributions for the neutral form; (ii) in the region 1350-1500 cm -1 , the large absorption band of zwitterionic alanine becomes less intense in the neutral form, and four main peaks appear better resolved; (iii) the region 1500-1900 cm -1 becomes more IR active for the neutral form of alanine than for the zwitterionic form, with the complete disappearance of absorption above 1800 cm -1 .

The effective modes of neutral alanine in the gas phase have been determined in order to identify the modes attesting to these differences between the IR spectra of the two forms of the molecule. For the sake of clarity, the 14 modes determined in the spectral region 1000-2000 cm -1 are not described, but the effective modes mainly responsible for the differences (i), (ii) and (iii) discussed above, are commented.

In figure 4.16, the selected effective modes responsible for the different features of the spectrum of neutral alanine with respect to zwitterionic alanine involve some vibrational motions of the fragments that mainly distinguish the two species: the OH bond in the COOH fragment and the NH 2 fragment itself. Mode 2 of figure 4.16 centered at around 1150 cm -1 is a seemingly wagging motion of the NH 2 moiety, which is absent in the zwitterionic form. Similarly, the wagging motion of NH 2 at 1200 cm -1 , which is coupled to the in-plane oscillations of the OH bond (in plane with COO), is totally absent in the absorption spectrum of zwitterionic alanine of figure 4.11. Mode 8, at around 1480 cm -1 , is similar to mode 3, just described, but it is coupled to a vibration in the CH 3 fragment that looks like an umbrella mode. Finally, that mode 12, slightly above 1700 cm -1 , clearly involves solely the HNH bending motion in NH 2 , and was not so well defined in the analysis reported in figure 4.14.

Similarly to the results reported in Table 4.4, Table 4.5 shows the comparison between the calculated frequencies and bands assignment of the effective modes analysis and the experimental results of reference [3]. 4.5: Left column: experimental frequencies collected from reference [3] in the spectral region 1000-2000 =cm -1 for the neutral form of alanine. Central column: bands assignment using the abbreviation from reference [3] (ρ for rocking, δ for bending, ν for stretching, ω for wagging; s for symmetric and as for asymmetric). Right column: theoretical frequencies from this work.

Experimental frequencies (cm

The C=O stretch is a lot closer to the experimental value here than the COO antisymmetric stretch for the zwitterionic form. The force field must better represent the C=O bond rather than the CO bond in the COO -group.

Thus the IR spectra of alanine in the zwitterionic and neutral forms have been analysed, and the corresponding molecular vibrations were given. For a finer analysis of these results for the zwitterionic form of alanine, the effective modes analysis will also be performed on the 3 fragments that constitute the molecule.

Fragment analysis

In this section, some of the molecular fragments of alanine in its zwitterionic form are analysed. Since the spectral region of interest, 1000-2000 cm -1 , is densely populated by overlapping and not well-defined absorption bands, this secondary study aims to clarify the results of the previous sections on the full molecule. From the entire FFMD trajectory, positions and velocities of three fragments, i.e., CH 3 , NH 3 , and COO are extracted in order to identify their effective modes and the bands of the full IR spectrum arising from their vibrations. Note that, while the analysis is performed on the isolated fragments, their dynamics accounts for the presence of the remaining degrees of freedom of the molecule.

IR analysis of the methyl group

To determine the effective-mode power spectrum of the methyl group, the following internal coordinates to determine the effective modes are introduced. The IR absorption bands of the CH 3 fragment are mainly located in the region of 1400-1500 cm -1 , as shown by the coloured bands in the lower panel of figure 4.17, and thus contribute to the large absorption region of zwitterionic alanine in gas phase in the region 1300-1600 cm -1 (black line in the upper panel of figure 4.17). The effective modes identified by this analysis are, unsurprisingly, combinations of bending modes (mode 1 and mode 3, depicted as arrows in the insets of figure 4.17

IR analysis of the carboxyl group

The following internal coordinates are used to determine the effective modes of the carboxyl group.

modes atoms stretch C1-O2 stretch C1-O3 bend O2-C1-O3
Table 4.8: List of COO (second column) internal coordinates in zwitterionic alanine used for the specific analysis of the carboxylic acid groups power spectrum.

The IR absorption of the COO fragment is composed of two main bands: one band appears essentially in the region 1500-1600 cm -1 , but it is coupled to smaller bands at (a bit less than) 1100 cm -1 and at 1300 cm -1 ; one band is very well defined and localised around 1900 cm -1 . These bands due to the power spectra of the effective modes are shown in the lower panel of figure 4.19. Therefore, also for this fragment, the bending motion in COO (mode 1, depicted as arrows in the inset of figure 4.19) is found to participate to the large absorption region of zwitterionic alanine in the gas phase in the region 1300-1600 cm -1 (black line in the upper panel of figure 4.19), while the asymmet-ric stretching in COO (mode 2, in the inset of figure 4.19) produces an absorption peak in the high frequency region of the analysed spectrum.

Alanine_gaz_plus_CO2

Frequency (cm -1 ) IR 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 IR spectrum mode1 mode2 ω1=1538 cm -1 ω3=1918 cm -1 CO2 in zwitterionic alanine The reported decomposition of the spectra of alanine confirms the capability of the effective modes analysis to unravel absorption spectra that can appear at first glance very broad and not well resolved. Since alanine is quite flexible and its dynamics was simulated at room temperature, the normal modes analysis would not be sufficient to satisfactorily characterise the complex absorption spectra discussed in the present work. The comparison between the two gas phase forms of alanine allowed us to identify the modes responsible for the changes in the IR spectra: even though it is expected that the changes in the neutral form are produced by the fragments of the molecule that are not present in the zwitterionic form, it is interesting to find that those changes are confirmed based on a rigorously-derived numerical procedure. The results from this effective modes analysis can be used to check the value of the frequency corresponding to each vibration. If the values of the frequencies are shifted from their experimental value, the parameters from the force field can be adjusted in consequence.

Concluding remarks

The results using our implementation of polarisable force field molecular dynamics for vibrational circular dichroism have been scrutinised on the example of the alanine molecule in the gas phase, in its neutral and zwitterionic forms. The importance of accounting for polarisation effects through the induced dipoles moments was demonstrated. The fluctuations of the charges need to be accounted for to represent com-pletely the VCD spectrum. Through scrutinising the convergence and the error to this convergence, the need to sample over a certain period of time was demonstrated. This sampling time (20x200 = 4 ns) is much longer than what is usually achievable with AIMD simulations. Harmonic and anharmonic DFT results were compared, giving a small difference between both spectra, showing the influence of anharmonicity on the IR and VCD spectra.

To confirm the FFMD results, a comparison with ab initio molecular dynamics for one trajectory was accomplished to determine the moments and spectra. This comparison permits the validation of the implementation in the gas phase. This chapter also discussed the molecular origin of the spectroscopic signals of alanine in the gas phase based on an effective modes analysis aimed to decompose the IR absorption spectra of alanine in zwitterionic and (stable) neutral forms. The decomposition relies on the identification of well-defined, maximally localised in frequency, absorption bands calculated as the power spectra of the effective modes, whose sum yields a good approximation of the full IR spectrum for the molecule at hand. The comparison to experimental frequencies displays a small number of discrepancies which hopefully should be solved in the future by adjusting the force field parameters.

Experimentally, VCD spectra are recorded in dense media where the signals are far more intense. The computational method presently assessed in the gas phase must therefore be extended to the relevant liquid and crystalline cases. This is the subject of the next chapter.

Chapter 5

Applications in condensed phases

The previous chapter was aimed at assessing the validity of force field molecular dynamics for VCD spectroscopy, using alanine in the gas phase as our testing ground. Here the scope of the study is extended by focusing on the same molecule but in solid (crystalline) and liquid (solvated) phases. Owing to the usually low intensities of VCD signals, a fair amount of material is needed to actually record them, making condensed phases the natural setting in experiments. As mentioned previously, the correct description of the environment is important in VCD [1].

The VCD spectroscopy of alanine in various condensed phases has been the subject of several earlier works. Jähnigen et al. showed that for crystalline alanine, supramolecular interactions need to be considered to represent properly the VCD spectrum [2]. This work explored the structure of the alanine crystal, depicted in figure 5.1a, where the amino groups are arranged in a helical fashion leading to a chiral supramolecular packing. The supramolecular interactions within a range of 4 to 5 Å give characteristic features of the vibrational spectra. In the hydrated phase, various studies examined the number of explicit water molecules needed to embed the alanine solute. With DFT methods, Tajkhorshid et al. demonstrated the need for 4 water molecules [4] whereas Frimand et al. argued the need for 9 explicit water molecules to represent the vibrational spectra [5]. From AIMD simulations Jalkanen et al. found that 20 water molecules are necessary to fully hydrate the zwitteri-onic species, including 11 solvent molecules on the methyl group, depicted in figure 5.1b [3]. However, they agreed on the importance of an explicit water solvent to capture the effect of the solvent and the chirality transfer effect. In this chirality transfer effect, the chiral properties of the solute are transferred to the solvent [6].

In the present chapter, the assets of FFMD for sampling time and system size are demonstrated for the alanine crystal and for hydrated alanine, both under periodic boundary conditions. An effective modes analysis was also performed to identify the vibrational contributions to the spectrum. To disentangle the chemical effects of hydrogen bonding from those of a physical nature, due to the presence of the embedding medium, alanine in a liquid dinitrogen solvent was also considered. This system was studied under the conditions of temperature and density that are more appropriate to this cryogenic medium. Finally, phenylcyclohexanediol was also studied as a crystal and solvated in dimethyl sulfoxide. This system is of chemical interest owing to the diol group that displays strong hydrogen bonding.

Crystalline alanine

Vibrational spectra obtained from force field molecular dynamics

In accordance with experimental crystallographic data, an alanine crystal was prepared for simulations. Given the chiral space group of alanine (P2 1 2 1 2 1 ), a periodic box of 12.073 × 12.342 × 11.57 Å was constructed containing 16 alanine molecules, as shown in figure 5.2. It is this system that was used to determine the electric and magnetic dipole moments throughout MD simulations. In the crystal, alanine is in its zwitterionic form so the molecules are closely held together by electrostatic interactions between the amino and carboxylate groups, represented by the black dotted lines for one amino group.

Comparing crystalline and gas phase vibrational spectra

Like in the gas phase, a certain number of trajectories was performed to sample the system. Each of these trajectories was initiated from a different starting point, chosen along a trajectory previously performed in the canonical ensemble. The electric and magnetic dipole moments were retrieved from each trajectory, then the Fourier transform of the corresponding time correlation function was determined. This gave an individual spectrum for each trajectory, all spectra being averaged to obtain a spectrum considered as converged. Under PBC, as explained in chapter 2, the VCD spectrum becomes origin-dependent through the minimum image convention, which changes the value of the position depending on the origin. Thus the nearest periodic image is directly chosen between particles to make the VCD spectrum origin-independent [7]. For the crystalline phase, 30 trajectories of 40 ps for a total simulation time of 1 200 ps at 300 K were simulated. In figure 5.3 the average spectra in both crystalline and gas phases are given. The spectra are well resolved and by comparison of the crystal spectra to their gas phase counterparts, the groups of peaks are mainly found around the same positions.

However, shifts of the NH 3 bend from 1580 to 1622 cm -1 and of the C O stretching mode from 1919 to 1853 cm -1 are displayed, represented by the red arrows in figure 5.3. These shifts are due to the electrostatic interactions between the amino and carboxylate groups present in the crystal. .6 shows a slightly faster convergence for IR over VCD. The convergence of the spectrum from the crystal compared to the spectrum from the gas phase is faster for IR but slower for VCD. If only one molecule of alanine is considered for the solid, the convergence of the crystal is faster in both IR and VCD. The supramolecular interactions of the crystal may add instability to the magnetic dipole moment, thus disturbing the convergence.

Furthermore, the crystalline phase shows convergence with a shorter amount of sampling time than the gas phase. In the previous chapter, it was shown that a simulation time of 4 ns was needed for convergence in the gas phase. The crystalline phase can be considered converged at a sampling ratio of 0.5 so after the system has sampled 15 trajectories of 40 ps, thus only 600 ps of sampling time.

Assessing moments and spectra from one trajectory

As in chapter 4, a trajectory (positions and velocities) from classical molecular dynamics based on force fields was exploited to calculate the electric and magnetic dipole moments using DFT within the CPMD software [8]. The DFT dipole moments were computed by our collaborator S. Jähnigen.

The spectra from both methods are shown in figure 5.7. Both IR and VCD spectra from FFMD are very similar to their DFT counterpart. As for zwitterionic alanine in the gas phase, the force field gives rise to more activity in the 1410-1450 cm -1 range and a lower intensity at 1925 cm -1 in the IR spectrum. The signs of all peaks in the VCD spectrum agree between DFT and FFMD. generally even improved respect to the values found for the gas phase. This is probably due to the force field describing the crystalline phase with more accuracy than the gas phase. To describe the interactions within the crystal, the force field's treatment of electrostatics via multipoles centered on each atom and the inclusion of polarisation through induced dipole moments prove to be important especially for the good directionality of the hydrogen bond. µ or its time derivative ì µ, and the magnetic dipole moment ì m, as obtained for a specific FFMD trajectory at 300 K but calculated from the force field or using density functional theory.

Hydrated alanine

Amino acids are generally solvated in water thus it is important to examine their behaviour in the hydrated phase. The study of alanine is continued here with a solvent of water. This will permit the analysis of the method using classical molecular dynamics based on force fields to determine VCD spectroscopy in this more realistic phase.

Vibrational spectra obtained from force field molecular dynamics

A molecule of alanine in the zwitterionic form was soaked in a previously equilibrated periodic box of length L = 18.6 Å. This procedure gives a periodic box containing 211 water molecules and 1 alanine molecule, see fig.

Masking the water contribution

Water has a large electric dipole moment, this leads to an IR spectrum showing large intensities (∼20 000 M -1 cm -1 ) especially compared to the intensities in the gas phase spectrum (∼400 M -1 cm -1 ), as shown in figure 5.10. The IR spectrum is computed from the average of 10 trajectories of 40 ps. Thus, when calculating the IR spectrum of hydrated alanine, the water spectrum hides the contribution of alanine. Experimentally this is avoided by subtracting the spectrum of pure water from the spectrum of hydrated alanine. Such a heuristic approach would be tricky in simulation because the pure solvent and the hydrated molecule do not share exactly the same periodic box. Thus the spectrum of water alone will not correspond to the spectrum of alanine and water, minus the contribution of alanine.

However if the same molecular dynamics trajectory is exploited and no water molecules are considered when evaluating the dipole moments, it becomes possible to have the local spectrum of alanine influenced by the water molecules. This can be achieved through scaling the dipole moments, as explained in section 2.6.4 of chapter 2 with the following damping function:

P K J (t) = 1 + e (| ì R K J (t)|-R 0 )/D -1
. This scaling employs a cutoff coefficient R 0 . Within the ChirPy program, for a cutoff of 0 Å, only the solute is considered.

For the VCD spectrum of hydrated alanine, the contribution to the spectrum from the water molecules should be zero as this solvent is achiral. However, the possible effects of chirality transfer in the vicinity of the solute should not be excluded. In this case, the solute could transmit chiral properties, including a VCD signal, to the solvent.

Choosing a cutoff around the alanine molecule

The electric and magnetic dipole moments were computed within the periodic box of hydrated alanine. As before, the Fourier transform of the appropriate time correlation functions was performed for each trajectory, then the resulting spectra were averaged. Here 30 trajectories of 40 ps at 300 K were exploited. Again, the nearest periodic image is directly chosen between particles to make the VCD spectrum origin independent. A cutoff It is shown here that if no cutoff is taken, thus R c equals the length of the periodic box L, then the intensity of the water peak takes over the whole spectrum (Fig. 5.11a).

Depending on the cutoff taken and thus on the number of water molecules included, the influence of water on the alanine spectrum can be varied. This is particularly clear between 1600 and 1900 cm -1 in figure 5.11b). This is the case also for the VCD spectra, though the VCD signal of an achiral molecule is zero, thus there should be no water signal. One reason for this can be that sampling of 30 × 40 = 1 200 ps is still not long enough for the contributions of the water molecules The normalised errors were determined as a function of the sampling ratio S r = N T /30

and they are shown in figure 5.14. The error was computed by comparing the spectra from the average of 1 < N T < 30 trajectories to the converged spectrum. The normalised errors were determined from the whole frequency range, from 0 to 4 000 cm -1 . The high fluctuation of the OH stretching band at 3500 cm -1 is thus also taken into account. This explains why the normalised error (figure 5.14) shows a slower convergence than the spectra from the averages of an increasing number of trajectories (figure 5.13) that only show the 1000 to 2000 cm -1 region.

The convergence of the VCD spectrum of hydrated alanine is very much more arduous for IR and also for the other phases, especially when water molecules are accounted for. This large error arises from the large fluctuation of the peaks corresponding to the water molecules, which is why the normalised error with a vanishing cutoff is considerably lower than the error with finite cutoffs. This large error from the signal of water molecules confirms the hypothesis that the system has not been sampled enough for the water contributions to average out to zero. This could explain why a residual signal for water is present in the FFMD spectrum. In the future, spectra from a longer sampling time need to be computed.

Assessing dipole moments and spectra from one trajectory

Following the same lines as in the previous chapter, the electric and magnetic dipole moments are determined from the same MD trajectory with DFT and with the force field methods. of R c = 3.85 Å was chosen to consider the alanine and its first solvation shell only. The spectra are rather similar and the peaks of both VCD spectra have the same signs. Again, as for zwitterionic alanine in the gas phase, the force field gives rise to more activity in the 1410-1450 cm -1 range and a lower intensity at 1925 cm -1 in the IR spectrum.

For each point along a FFMD trajectory the electric dipole moment ì µ, its time derivative ì µ, and the magnetic dipole moment ì m, were determined from the force field and using DFT. The values of the moments from the force field are plotted along the y-axis and the values of the DFT moments are plotted along the x-axis. This yields figures 5.16 to 5.18, where r is the associated Pearson correlation coefficient.

Effective mode analyses for crystalline and hydrated alanine

The moments acquired though classical molecular dynamics based on the force field have thus been validated for crystalline and hydrated alanine. Here an effective modes analysis is performed to determine the vibrational modes causing the peaks. Two FFMD trajectories of 40 ps at 300 K were analysed to retrieve the positions and velocities. A first trajectory of the alanine crystal was exploited, then a trajectory of hydrated alanine. While the trajectories were performed on the alanine molecule with its full environment, the power spectra and the effective modes were determined by retrieving the molecular data from a unique alanine molecule. Therefore, the IR spectra and the trajectories themselves fully account for the effect of the environment, similarly to the fragment analysis reported at the end of the previous chapter. These positions and velocities were transformed into internal coordinates to initiate the effective modes analysis. The internal coordinates are detailed in table 5.1.

The two modes in the upper panels refer to the hydrated form of zwitterionic alanine (h.), while the modes in the lower panels refer to its crystal form (c.). Note that all modes are very similar to each other, involving some kinds of wagging and twisting of the CH 3 and NH 3 fragments, mainly coupled to the vibrations of the CH bond that is not involved in CH 3 . Only mode 3 in the crystal form appears slightly different, as the umbrella mode of NH 3 can clearly be identified.

ω1=1044 cm -1 (h.) ω4=1134 cm -1 (h.)

ω3=1169 cm -1 (c.) ω1=1046 cm -1 (c.) Concerning the appearance of the double peaks in the region 1400-1650 cm -1 of the IR spectrum of crystal alanine in comparison to the single peaks structure of hydrated alanine, the corresponding effective modes are reported in figure 5.22. In the figure, two modes were selected, one around 1450 cm -1 and one slightly above 1600 cm -1 .

ω6=1444 cm -1 (h.) ω4=1621 cm -1 (h.) ω11=1566 cm -1 (c.) ω5=1451 cm -1 (c.) The modes on the left panels and on the right panels of figure 5.22 are very similar: on the left, both modes involve the umbrella motion of NH 3 coupled to the (out-of-phase) umbrella motion in CH 3 ; on the right, both modes involve the HNH bending motion in NH 3 , coupled to an out-of-plane vibration of the third NH bond in the fragment. Therefore, it appears reasonable to assume that despite the fact the IR spectrum of the crystal form of alanine presents a better resolved structure than its hydrated form, the source of the absorption bands is to be found in modes that are fundamentally very similar. Presumably, the constrained vibrations in the crystal reduce the effect of anharmonicity, thus allowing the absorption bands to be better resolved in frequency.

Finally, in figure 5.23, the effective modes responsible for the appearance of the absorption bands at around 1350 cm -1 (left panels) and 1850 cm -1 (right panels) are analysed. The effective modes at lower frequency appear to be slightly different. As in water, an umbrella motion involving the fragment CH 3 cannot be found, but rather a combination to HCH bending and CH vibration. As previously, this slight difference can be ascribed to the different rigidities of the environment in which alanine is embedded.

ω5=1313 cm -1 (h.) ω13=1835 cm -1 (h.)

ω13=1854 cm -1 (c.) ω4=1169 cm -1 (c.) Figure 5.23: Effective modes at around 1350 cm -1 (left) and 1850 cm -1 (right) for the hydrated (h.) form of zwitterionic alanine (upper panels) and for the crystal (c.) form of zwitterionic alanine (lower panels).

The effective mode at higher frequency allows to very clearly identify the asymmetric stretching motion of COO in the crystal, which becomes instead a wagging motion in water and which is coupled to vibrations in the other fragments of the molecule. It can thus be concluded that in water in the region spectral region 1000-2000 cm -1 the pure antisymmetric stretching motion of COO appears to be suppressed.

Alanine solvated in dinitrogen

A second solvent for alanine was studied, dinitrogen. The idea of this study is to examine a solvent that does not involve hydrogen bonding with the solute and for which the results are expected to be closer to those in the gas phase, except for temperature effects. At 300 K, dinitrogen is naturally found in the vapor phase, and its influence on the vibrational spectra of alanine would be extremely limited. To impose a dense liquid, the simulations were performed at 77.35 K. This time a molecule of neutral alanine was solvated in a periodic box of 18.27 Å. The density of the box is 0.8171 g/cm 3 .

The low temperature prevents an extensive exploration of the potential energy surface during the MD simulation, as the alanine molecule stays mostly vibrating around its starting structure. Thus the lowest energy structures of families 2 and 3 of alanine found in the exploration in chapter 4 were used as starting points for FFMD simulations in the gas and solvated phases. Family 1 was excluded as it is not found in the force field simulations.

For each family, 10 different simulations were performed. The first simulation used the structure of its corresponding family and the other nine structures used the last structure of the previous simulation. Each simulation was composed of an initial 5 ps MD simulation in the canonical ensemble to equilibrate the system at 77.35 K. Then a 200 ps MD simulation, for the gas phase and a 40 ps MD simulation for the solvated phase, in the microcanonical ensemble were carried out to retrieve the electric and magnetic dipole moments. Thus for each family a spectrum from the average of 10 spectra was obtained.

Vibrational spectra obtained from force field molecular dynamics

To achieve an averaged spectrum with sufficient sampling at low temperature, thus including both families, the thermal probability Z α / i Z i , with Z corresponding to equation (4.12), was determined for each configuration of both families at 77.35 K, using the energies and harmonic frequencies from the force field. The value for each configuration was summed within each family so as to obtain values of the thermal probability for both families. The spectra from the average of 10 trajectories was then multiplied by the corresponding thermal probability to give the spectra at 77.35 K in figure 5.24.

First, the spectrum from the gas phase is compared from the FFMD simulations at 300 K and 77.35 K, then the spectrum from the gas phase at 77.35 K is compared to the spectrum of alanine solvated in N 2 at 77.35 K.

In figure 5.24a, the bands are present at similar locations for both temperatures. As expected, due to the higher temperature, the peaks for the alanine in the gas phase at 300 K, are significantly wider.

In figure 5.24b, at 1250 and 1700 cm -1 , there is a dramatic change between the VCD spectrum of alanine in the gas phase and the VCD spectrum of alanine solvated in N 2 at 77.35 K. The IR spectra do not shown significant change, only a slight widening of the peaks. This can be explained by the change of configuration of the structure of family 2 when solvated in N 2 compared to the structure shown in chapter 4. The NH 2 group rotates by 180 • and the hydrogen of the hydroxyl group turns towards the amino group. This change does not occur when the structure of family 2 is in the gas phase at 77.35 K. The structure of family 3 remains identical to the structure shown in chapter 4 in the gas phase and solvated in N 2 at 77.35 K. This confirms the sensitivity of the VCD spectum to conformational changes, especially compared to the IR spectra where the change is minimal.

Phenylcyclohexanediol in the crystalline phase and solvated in DMSO

A final system was studied, phenylcyclohexanediol. This molecule has a diol group which contains two hydroxyl groups. These diols can form hydrogen bonds, which is of particular interest. Combining with theoretical results, more information can be extracted from the experimental data. This molecule was examined in the crystalline phase and solvated in DMSO (see figure 5.25). The unit cell of the phenylcyclohexanediol crystal is composed of 8 molecules that form dimers. Each dimer has 2 hydrogen bonds between the 4 hydroxyl groups, represented by the black dotted lines. This unit cell is rather large (19.767 × 5.645 × 19.924 Å) making FFMD particularly valuable to address such a complex system.

The periodic box, for the crystalline phase, has the same size as the unit cell defined from crystallographic data shown in figure 5.25a. For the solvated phase, a cubic box of 18.0 Å was chosen. Phenylcyclohexanediol was solvated in 44 DMSO molecules and the density of the system was 1.0336 g/cm 3 ; it is depicted in 5.25b. Both phases were studied at 150 K, exploiting 20 trajectories of 80 ps to collect the electric and magnetic dipole moments.

Vibrational spectra obtained from force field molecular dynamics

The main difference between solvated and crystalline phenylcyclohexanediol is the interaction of its hydroxyl groups. When it is solvated in DMSO, its hydroxyl groups interact with DMSO though in the crystalline phase they interact with one another. This major difference influences the IR and VCD spectra of figure 5.26, where the crystal spectra are compared to their solvated phase counterparts.

For solvated phenylcyclohexanediol, as a first step, a vanishing cutoff (R 0 =0) was chosen, thus leading to ignore the DMSO contribution.

The interest of this system is the OH groups, thus the spectra is shown from 2800 to 4000 cm -1 where the bands of the OH stretch are present. The peaks at 3000 cm -1 correspond to the CH stretches, which are very similar for both phases in the IR spectra. In the VCD spectra, then are present at the same frequency but with a change of sign which is probably due to the change of environment of the carbon rings.

The peaks in between 3600 and 3800 cm -1 correspond to the OH stretching modes. The difference in both the IR and VCD spectra is very noticeable and expected due to the very different environments and configurations of the OH groups. The large peak at 3700 cm -1 in the spectrum of phenylcyclohexanediol solvated in DMSO could be composed of multiple peaks. Adjustment of the smoothing factor needs to be tested in the future. In addition, as no solvent is considered in the IR and VCD spectra of phenylcyclohexanediol solvated in DMSO, future work could consider the effect of the DMSO solvent. In the following, the spectra of isolated phenylcyclohexanediol are shown. In this system there is no inter-molecular hydrogen bonding as there is only one molecule. Figure 5.27a, shows the influence of the various terms of the magnetic dipole moment in the 800 to 2000 cm -1 range and figure 5.27b in the 2800 to 4000 cm -1 range. In figure 5.27a, the inclusion of all the terms containing the induced dipole moment influences the VCD spectrum the most between 1200 and 1400 cm -1 though without any peak inversion, contrary to the alanine spectrum discussed in chapter 4. In figure 5.27b, the influence of the terms containing the induced dipole moment is minimal, we can deduce that the induced dipole moments are small on the OH groups, compensation effects being quite limited as well.

Comparing these spectra to the crystalline and solvated phases, the OH peaks are found to be narrower and shifted to higher frequencies. The peaks corresponding to the CH stretches are slightly shifted to lower frequencies. Both of these shifts are expected due to the change of environment.

Statistical errors and convergence rates

The convergence is examined by comparing the spectrum from the average of 20 trajectories of 80 ps to the average of a lesser number of trajectories. First, the convergence is observed directly through the spectrum, in figure 5.28, then through the normalised error, in figure 5.29.

Concluding remarks

Throughout this chapter, five systems in condensed phases were studied: alanine crystal, hydrated alanine, alanine in dinitrogen, phenylcyclohexanediol in DMSO and phenylcyclohexanediol crystal. For these systems, the IR and VCD spectra were determined with classical molecular dynamics based on force fields. For alanine and phenylcyclohexanediol, the statistical errors and convergence rates were determined and they reveal the extended sampling that is needed, though the crystal shows the shortest sampling time required to reach convergence. This further supports the use of classical force fields to simulate VCD spectra for such systems which are barely accessible with methods that explicitly account for electronic structure. The moments were also determined from one trajectory, using DFT and the force field methods. The comparison of these results showed the excellent agreement of the FFMD results with DFT. This extends our validation of polarisable force fields for VCD spectroscopy from the gas phase to a variety of condensed phases that are more relevant in experiments.

The spectra of alanine in the hydrated and crystalline phase were finally decomposed into the contributions of effective modes. This analysis can be used in the future to assess the peak positions and adjust the force field parameters if necessary.

For the phenylcyclohexanediol, the study can be furthered with the effective modes analysis and with the study of the effect of the solvation by DMSO. A comparison of the dipole moments would confirm for a different system the classical molecular dynamics method to achieve VCD spectroscopy.

Chapter 6

Conclusions and perspectives

Vibrational circular dichroism and its general theory have been presented. Various methods to obtain VCD spectra were discussed, these methods include anharmonicity and temperature effects with various degrees. The method employed for the implementation of VCD spectroscopy in classical molecular dynamics based on force fields naturally accounts for such effects. Methods to access normal (harmonic) and effective (anharmonic) modes were discussed. From the effective modes, it is possible to assign spectral peaks from molecular dynamics trajectories to specific molecular vibrations.

Trans-1-amino-2-indanol was examined as a case study for VCD on a complex fluxional system. Because of this fluxionality, many structures coexist within tiny energy differences. These structures can be explored to a certain extent with static calculations though more conformers can be accessed through molecular dynamics, which crosses barriers and produces conformations not biased by intuition. In this case study, the effects of the solvent on the spectrum were scrutinised. In this system, the main contribution to the spectrum arises from the local solvent, at a distance less than 1 Å from the solute. This leads to the solvent influencing only indirectly the calculated spectra through structural effects. Ab-initio molecular dynamics simulations were shown to be more dependent on the starting point than simulations using classical molecular dynamics based on force fields. This latter MD can explore for longer periods of time as its computational cost is lower. As this system shows important solvent mobility, its solvation cannot be taken fully into account with small, finite-size clusters. Therefore VCD spectroscopy from FFMD presents itself as a complementary method to the static picture using the clusterin-a-liquid model and AIMD simulations.

During this thesis, appropriate expressions for the electric and magnetic dipole moments were implemented in classical molecular dynamics based on polarisable force fields to obtain subsequently the IR and VCD spectra from the Fourier transform of the corresponding time correlations functions. Alanine in the gas phase under the zwitterionic and neutral forms was used as a benchmark system. AMOEBA force field was parametrised for this amino acid, which turns out to be the smallest chiral amino acid. It was shown that the inclusion of the polarisation of the atoms is important to capture the VCD signal. Moreover, important sampling times are needed to reach a converged spectrum. These sampling times exceed the capabilities of AIMD. The electric and magnetic dipole moments computed from FFMD were compared successfully to DFT calculations at each time step. This comparison validates the implementation in the gas phase.

With the FFMD spectra validated, the IR spectrum was decomposed into effective modes to assign molecular vibrations to their respective peak. This information is not directly accessible from the MD simulation, thus an effective mode analysis had to be performed. This analysis is used as an opportunity to decompose the contributions from each functional group to the spectrum.

As the experimental interest lies in condensed phases, the implementation of VCD spectroscopy through classical molecular dynamics was extended to crystalline and solvated phases. Using this implementation, alanine was studied as a crystal, hydrated and solvated in dinitrogen. The crystal and hydrated alanine were simulated at 300 K and alanine in dinitrogen at 77.35 K. The inter-molecular interactions in the alanine crystal and their influence on the VCD spectrum were demonstrated via a shift of the bands corresponding to the NH 3 and CO stretching compared to their counterpart bands in the gas phase. These interactions hinder the convergence of the VCD spectrum of alanine in crystalline form. The IR spectrum of hydrated alanine presents a signal for the bending of the water solvent, as expected. However this is also the case for the VCD spectrum which is more unexpected. This could be due to insufficient sampling time, as shown through the normalised errors. This could also be caused by the water contributions to the spectrum that had not enough time to average out to zero. The force field and DFT moments in the condensed phases were also compared and shown to correlate to a very high degree. This contributed to validate our implementation of VCD spectroscopy with polarisable force fields within the framework of classical MD under diverse chemical environments and thermodynamical conditions. The studies of alanine in dinitrogen and of phenylcyclohexanediol in the crystalline phase and solvated in DMSO confirm the sensitivity of the VCD signal to the conformations of the system and their diversity.

The present work can be extended along various directions. A first direction would be a longer exploration of the hydrated alanine. This would assess the possible convergence of the system at longer sampling times. Moreover the effect, if any, of the periodic box size of the alanine crystal could be examined. For phenylcyclohexanediol, the electric and magnetic dipole moments need to be compared to dipole moments obtained via DFT, and the spectra examined through an effective mode analysis. The effective mode analysis of all the systems can be used to check the values of the frequencies against experimental results and to adjust the force field parameters, if necessary. Some modes are more susceptible than others to a change of sign in the VCD spectra. As the intensities of the VCD spectra are a function of the cross correlation of the electric and magnetic dipole moments, if the angle ξ(i) between the electric dipole moment and the magnetic dipole moment is smaller than 90 • , the intensity of the corresponding peak is positive. If the angle is larger than 90 • , the peak is negative [1]. This means that modes corresponding to dipole moments with an angle close to 90 • are more susceptible to a change in the sign. These are called non-robust modes and can be changed by small perturbations (experimentally, e.g. solvent effects, and computationally, e.g. level of calculations). If ξ(i) differs from 90 • by at least 30 • , the mode corresponds to a robust mode and will not change sign with ease [2]. An extension of the implementation could be to determine the angle between the two dipole moments for each molecule, to assess the robustness of the various modes. This will help with the comparison to other theoretical and experimental results as the analysis will be more precise because only robust modes will then be considered. A second extension of the code could be the determination of the electronic polarisability, to establish the Raman optical activity spectrum which arises from the difference in intensity of Raman scattered right-and left-circularly polarised light [3,4].

During the MD simulations, the nuclei were considered as classical particles. However this leads to a slight overestimation of the frequencies, as in the harmonic approximation. Though, solving the quantum problem, thus Schrödinger's equation for nuclei is excessively expensive for the systems treated in this thesis. For system of this size, semiclassical approaches can be used and MD simulations can be corrected to account for vibrational delocalisation. For example, the MD trajectories can be initialized based on the harmonic normal mode coordinates and associated momenta, then given an appropriate quantum mechanical thermal energy, with an amplitude chosen so the excess energy equals the (harmonic) zero-point energy [5]. To approximate quantum-mechanical correlation functions, the ring-polymer molecular dynamics (RPMD) method can also be used. In this method, each of the N classical particles in the system is represented by a harmonic ring-polymer of n monomers and n is pushed toward infinity [6] .

The main interest of using force field molecular dynamics is its capacity to treat large, explicitly solvated and fluxional systems. Thus, interesting perspectives lie in the determination of the vibrational signatures of polypeptides, fibrils and proteins. Various advantages of VCD have been shown for these systems. First, VCD facilitates the study of proteins as there is no interference from the water absorption [7]. Second, VCD has an unusual sensitivity to fibril growth and development in solution [8]. Fibrils with supramolecular chirality show highly enhanced VCD intensities and thus can be easily studied with this spectroscopic method [9]. Moreover, metal ions are found in intimate association with nucleic acids in their natural environment and VCD can be used to determine the effect of these ions on DNA structure. This is important as they can stabilize and destabilize biological structures [10].

Keywords: Vibrational circular dichroism, classical molecular dynamics, fluxional molecules

Abstract: Vibrational circular dichroism (VCD) is the weak difference in absorption for chiral molecules between right-and left-polarised light in the infrared (IR) range. It has become an increasingly popular spectroscopic technique since its discovery in the 1970s and has promising applications in pharmacology owing to its ability to determine absolute configurations of chiral molecules. The shape of VCD spectra is highly sensitive to minor changes in conformation and molecular interactions, which makes it a valuable probe of conformational isomerism and solvation. Structural assignment of experimental VCD spectra usually requires successful comparison with theoretical calculations, which also implies that the solvent is described with sufficient chemical accuracy. This manuscript examines how computational modelling based on polarisable force fields can be exploited to assist in the elucidation of VCD signals. This modelling is applied to flexible molecules for which conventional quantum chemical approaches relying on static approximations are necessarily limited. In particular, anharmonicities and finite temperature effects are difficult to account for, and the solvent is typically modelled in a rather rudimentary way through continuum descriptions or very few explicit solvent molecules. In the work presented here, the AMOEBA force field implemented in the Tinker software package was employed to model electric and magnetic dipole moments of various molecules of interest, for applications in the gas phase, as well as the condensed phases of solvated or crystalline systems. Sampling and real-time dynamical information were achieved by means of classical molecular dynamics simulations. Anharmonic IR and VCD spectra obtained from molecular dynamics trajectories were interpreted using an effective mode analysis in which the internal coordinates are linearly decomposed by projection onto suitable effective modes, the corresponding weights being obtained in such a way that the peak of each mode is as localized as possible. This effective mode analysis sheds light onto the various contributions of the different functional groups of the molecule to individual peaks in the vibrational spectrum. As a first illustrative example, the case of trans-1-amino-2-indanol solvated in dimethyl sulfide was examined by combining static explorations at the level of density-functional theory, short ab initio MD (AIMD) trajectories, and extended sampling of the potential energy surface using the force field. Comparison with existing measurements confirms the importance of flexibility for this molecule, the active role of the nearby solvent molecules, and their contribution to the VCD signals. The determination of VCD spectra from the polarisable force field was then undertaken and its performance was assessed by comparison with electronic structure calculations, namely AIMD. This comparison was performed on the very same MD trajectory, using the alanine amino acid in the gas phase as our benchmark system. The various contributions to the dipole moments arising from the induced dipole moment or its time dependence are discussed. Having validated the methodology, the statistical convergence of the IR and VCD spectra was quantified by defining dedicated error measures. These errors show in turn that convergence of the VCD spectra is generally slower relative to IR spectra and that convergence is achieved only after extended sampling, way more than the amount achievable by AIMD. The force field approach was then extended and applied to alanine in crystalline and solvated phases. Here again, comparison with AIMD simulations confirms that electric and magnetic dipole moments can be reliably evaluated using the polarisable force field. This paves the way towards determining the VCD spectrum of chemically complex and flexible systems that are inaccessible via methods that explicitly account for electronic structure.

  Caractérisation structurale des molécules chirales : développement d'outils de simulation moléculaire de spectres de dichroïsme circulaire vibrationnel Mots clés : Dichroïsme circulaire vibrationnel, dynamique moléculaire classique, molécules fluxionnelles Le dichroïsme circulaire vibrationnel (VCD) est la faible différence d'absorption des molécules chirales entre la lumière polarisée droite et gauche dans le domaine de l'infrarouge (IR). Cette technique spectroscopique est devenue de plus en plus populaire depuis sa découverte dans les années 1970 et a des applications prometteuses en pharmacologie en raison de sa capacité à déterminer les configurations absolues des molécules chirales. Ceci est démontré dans la figure suivante où chaque enantiomère donne un signal identique sur le spectre infrarouge et opposé sur le spectre VCD : IR spectrum VCD spectrum

Figure 1 :

 1 Figure 1 : Spectres IR et VCD de l'alpha-pinène. La figure a été adaptée de la référence [1].

Figure 2 :

 2 Figure 2 : Courbes de contour de l'occurence des distances O(trans-AI)-DMSO et N(trans-AI)-DMSO du complexe 1:1 du trans-1-amino-2-indanol à 300 K determiné à partir de simulations avec le champ de force

Figure 3 :

 3 Figure 3 : Spectres IR et VCD de la forme (a) zwitterionique et (b) neutre de l'alanine moyennés sur 40 trajectoires avec un champ de force à 300 K. La figure (a) examine l'addition des différentes contributions au moment dipolaire magnétique et la figure (b) montre les spectres produits avec le moment dipolaire électrique et avec sa dérivée.

Figure 4 :

 4 Figure 4 : Erreur normalisée des moyennes des spectres IR et VCD de l'alanine isolé de N T trajectoires avec un champ de force à 300 K, en fonction de la progression du rapport d'échantillonnage S r = N T /40.

Figure 5 :

 5 Figure 5 : Corrélation du moment dipolaire électrique ì µ, de sa derivée ì µ et du moment dipolaire magnétique ì m de l'alanine hydratée.

Figure 1 . 1 :

 11 Figure 1.1: One chiral molecule, two enantiomers.

Figure 1 . 2 :

 12 Figure 1.2: (left) The different smells of each enantiomer of carvone; (right) Chirality of the snail's shell.

Figure 1 . 3 :

 13 Figure 1.3: Visual representation of deoxyribonucleic acid (DNA). This figure has been reprinted from reference [6].

Figure 1 . 4 :

 14 Figure 1.4: Spectrum of electromagnetic waves (Adapted from an image of Inductiveload [CC BY-SA 3.0] [9] )

Figure 1 . 5 :

 15 Figure 1.5: Main mechanisms involved in IR and Raman spectroscopies.

Figure 1 . 7 :

 17 Figure 1.7: Chiral matter absorbing left circularly polarised light but transmitting right circularly polarised light. This defines circular dichroism.

Figure 1 . 8 :

 18 Figure 1.8: (Color) 2D projection images showing dichroism in the photoelectron angular distributions obtained from camphor enantiomers at hν =9.2 eV using left-and right-circularly polarised light (labeled lcp and rcp, respectively). Top row: the difference (rcp-lcp) for pure S-(-) and R-(+) enantiomers; Bottom row: the difference between molecular enantiomers [S-(-)-R-(+)] with light of a fixed helicity (either lcp or rcp). This figure has been reprinted from reference [19].

Figure 1 . 10 :

 110 Figure 1.10: IR (lower frames) and VCD (upper frames) spectra in the amide I mode of poly-L-lysine region for 3 conformations of the polypeptide. This figure has been adapted from reference [23].

Figure 1 .

 1 Figure 1.11: (4S,7R)-(-)-4-isopropylidene-7-methyl-4,5,6,7-tetrahydro-2(1)H-indazoles.This figure is adapted from reference[55].

Figure 1 .

 1 Figure 1.12: (a) Optimized conformers, relative energies and Boltzmann populations for pyrazole diastereomer 2a. (b) Comparison of VCD (upper frame) and observed IR (lower frame) spectra for (-)-diastereomer (right axes) with calculated spectra for diastereomer 2a conformers A and B. These figures have been reprinted from Tur et al. [55].

Figure 1 . 14 :

 114 Figure 1.14: Comparison of the simulated IR and VCD spectra in the gas phase and with PCM solvent model for CDCl 3 with the corresponding experimental spectra of Pd(TOLBINAP)Cl 2 in CDCl 3 , as reproduced from reference [56].

Figure 2 . 1 :

 21 Figure 2.1: Light polarisation, shown with the E-field in red and B-field in blue. (A) Vertically polarised light with polarisation vector ê = ẑ; (B) E-field of right-circularly polarised light (RCP) with the B-field not shown for clarity; (C) fields at the origin for LCP, facing backward along the propagation direction; (D) fields at the origin for RCP. This figure has been reprinted from reference [1].

Figure 2 . 2 :

 22 Figure 2.2: Energy level diagram of vibrational transitions.

Figure 2 . 3 :

 23 Figure 2.3: Jacob's ladder of density functional approximations. Th figure has been adapted from reference [25].

  .4. The experimental study was performed in different solutions depending on the mode under scrutiny, namely, a C 2 Cl 4 solution was used for the C-H and C-D stretching modes and a CS 2 solution otherwise.

Figure 2 . 4 :

 24 Figure 2.4: Structure of 2,3-trans-d 2 -oxirane.

Figure 2 . 5 :

 25 Figure 2.5: Comparison of calculated, Hartree Fock ( ) and DFT (•) rotational strengths (R calc ) of 2,3-trans-d 2 -oxirane to experimental values (R expt ). Experimental rotational strengths are from solution spectra. Arrows indicate C-H and C-D stretching modes. The line is of unit slope. Adapted from Cheeseman et al. [32].

Figure 2 . 6 :

 26 Figure 2.6: Comparison of the Morse and the harmonic oscillator potentials. This figure has been reprinted from reference [36].

Figure 2 . 7 :

 27 Figure 2.7: Crystal structure of the DHB dimer (left) and experimental and calculated VCD and IR spectra of methylbenzyl-amine-borane (right). This figure has been adapted from Merten et al. [50].

  then this timecorrelation function is called an autocorrelation function. If the variables represent two different quantities (C(τ) = A(t i ) • B(t j ) ) then the time-correlation function is called a cross-correlation function. A and B are dynamical variables [65, Chapter 5.2].

  the present work, IR absorption spectroscopy is represented with the time autocorrelation function of the time derivative of electric dipole moment: C IR (τ) = ì µ(0) • ì µ(τ) To represent VCD, the cross-correlation function C VCD (τ) with A = ì µ(0) and B = ì m(τ) is used, where ì m is the magnetic dipole moment.

Figure 2 . 8 :

 28 Figure 2.8: Schematic representation of (R)-propylene-oxide (RPO). This figure has been adapted from Scherrer et al. [8].

Figure 2 . 9 :

 29 Figure 2.9: Dynamical IR and VCD spectra of a (R)-propylene-oxide molecule solvated in water at 340 K (1 RPO molecule in 53 H 2 O molecules). The intensities are in units of km cm mol -1 . We show the molecular correlation and damping cutoffs of R 0 = 3.75 Å and R 0 = 4.0 Å, which includes the hydrogen bonded water molecules. The shaded areas

i)

  was compared to the average spectrum from all the simulations (I VCD av ). The individual spectra and the average spectrum are considered for different values of R 0 which are given in figure 2.10. Three dif- ferent systems were studied, neat RPO with 13 RPO and 1 propanal molecules and two solvated systems: one with 7 RPO molecules in 34 H 2 O molecules (low solvation) and one with 1 RPO molecule in 52 H 2 O molecules (high solvation).

Figure 2 . 10 :

 210 Figure 2.10: Cutoff dependence of the mean standard deviation of the VCD spectra. The magnitudes at the molecular level, i.e., at small distances, originate from the different amounts of sampling. This figure has been reprinted from Scherrer et al. [8].

Figure 2 . 11 :

 211 Figure 2.11: Schematic representation of the four key contributions to a molecular mechanics force field: bond stretching, angle bending and torsional terms and non-bonded interactions. This figure has been reprinted from reference [55].

Figure 2 . 13 :

 213 Figure 2.13: Lennard-Jones potential as a function of the particle distance. This figure has been reprinted from reference [96].

Figure 2 . 14 :

 214 Figure 2.14: Two equal point charges repel each other, and two opposite charges attract each other, with an electrostatic force F. This figure has been reprinted from reference [97].

Figure 2 . 15 :

 215 Figure 2.15: Schematic representation of the deformation due to the induced electric field E from the charges of other atoms.

Figure 2 . 16 :

 216 Figure 2.16: Qualitative description of approaches to modelling polarisation. (a) Fluctuating charges model where atomic charges are modified to reproduce the effect of polarisation; (b) Drude oscillator model where a point charge is attached to each atom by means of a spring, and moved, in order to reproduce polarisation effects; (c) induced dipole model where an explicit dipole term is added to each point charge. This figure has been reprinted from Illingworth et al. [106].

Figure 2 . 18 :

 218 Figure 2.18: Sketch illustrating the common origin problem within periodic boundaries.This figure has been reprinted from Jähnigen et al.[7].

Figure 2 . 19 :

 219 Figure 2.19: Structure of (1S,2S)-trans-1-amino-2 indanol and atom numbering.

Figure 2 . 20 :

 220 Figure 2.20: IR spectra of the trans-AI monomer in the gas phase from a FFMD trajectory at 300 K of 1 ns and from the average of 10 MD trajectories of 100 ps.

Figure 2 . 21 :

 221 Figure 2.21: IR spectra of the trans-AI monomer from FFMD simulations at 300 K with different time steps.

Figure 2 . 22 :

 222 Figure 2.22: Monitoring the CO and N-H N-methyl-acetamide stretch modes by DMD simulations. This figure has been reprinted from Thaunay et al. [119].

DMSO

Figure 3 . 1 :

 31 Figure 3.1: Structure of (1 S ,2 R)-(-)-cis-1-amino-2-indanol (cis-AI), of the solvent, dimethyl sulfoxide (DMSO) and of (1S,2S)-trans-1-amino-2-indanol (trans-AI).

  The 5 lowest configurations of (1 S ,2 R)-(-)-cis-1-amino-2-indanol (cis-AI) with 2 molecules of explicit DMSO solvent. The values in the top right corners are the Gibbs free energies in kcal•mol -1 . Exp Calculated (b) Top: Experimental VCD spectrum compared with the simulation resulting from Boltzmann-averaged contributions of the 1:2 complexes. Bottom: Individual VCD contributions.

Figure 3 . 2 : 3 . 2 .Figure 3 . 3 :

 323233 Figure 3.2: Structures and corresponding spectra of the 1:2 complex composed of 1 molecule of (1S,2R)-(-)-cis-1-amino-2-indanol (cis-AI) and 2 molecules of explicit DMSO solvent. This figure was adapted from Le Barbu-Debus et al. [4].

Figure 3 . 4 :

 34 Figure 3.4: (a) Four families of configurations of the (S,S)-trans-AI (b) Experimental VCD spectrum of (S,S)-trans-AI and simulated spectra for different configurations with the lowest-energy configuration at the bottom and the highest-energy configuration at the top. Adapted from Le Barbu-Debus et al. [12].At first sight, inspection of figure 3.4 indicates that the experimental spectrum can be accounted for by the contribution of a limited number of structures. The best match is obtained for mono OH -eqg -G -, which is the second most stable complex corresponding to the light green spectrum. All features except the positive band at ∼1180 cm -1

  Experimental VCD spectrum of trans-AI and the associated simulated spectra of the monomer and the 1:1 complex obtained from the Boltzmann-weighted average of 8 and 17 spectra, respectively.

Figure 3 . 6 :

 36 Figure 3.6: Distribution of distance between the oxygen atom of trans-AI and the oxygen/sulfur atom of DMSO, obtained from the AIMD trajectories at different temperatures (BLYP-D3).

Figure 3 . 7 :

 37 Figure 3.7: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances for the 1:1 complex, with comparison between the BLYP (left) and B3LYP (right) functionals in the AIMD simulations.

Figure 3 . 8 :

 38 Figure 3.8: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances for the 1:2 complex, using the BLYP functional in the AIMD simulations. The left plot represents DMSO 1 and the right plot DMSO 2.

Figure 3 . 9 :

 39 Figure 3.9: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances in the trans-AI:(DMSO) 1:1 and the trans-AI:(DMSO) 2 1:2 equatorial complexes at 150 K in FFMD simulations.

Figure 3 . 10 :

 310 Figure 3.10: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances in the trans-AI:(DMSO) 1:1 and the trans-AI:(DMSO) 2 1:2 axial complexes at 150 K in FFMD simulations.

Figure 3 . 11 :

 311 Figure 3.11: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances for each of the five DMSO molecules in the 1:5 equatorial complex at 150 K, together with their superposition. The contour plots are divided in four regions labelled with the corresponding static 1:1 complexes.

Figure 3 . 12 :Figure 3 . 13 :

 312313 Figure 3.12: Dihedral angle C 2 C 3 C 4 C 9 of the 5-membered ring of trans-AI during the FFMD simulations at 300 K with one, two and five DMSO molecules.

Figure 3 . 14 :

 314 Figure 3.14: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances in the 1:1 complex at 300 K in the FFMD simulations.

Figure 3 . 15 :

 315 Figure 3.15: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances for both DMSO molecules (left and right) in the 1:2 complex at 300 K in the FFMD simulations.

2 2 Figure 3 . 16 :

 22316 Figure 3.16: Distances evolution as a function of time in the trans-AI:(DMSO) 2 1:2 complex at 300 K, as obtained from FFMD simulations.

Figure 3 . 17 :

 317 Figure 3.17: Contour plots of occurrence of O(trans-AI)-DMSO and N(trans-AI)-DMSO distances for the five DMSO molecules in the trans-AI:(DMSO) 5 complex during the FFMD simulations at 300 K.

Figure 3 .

 3 Figure 3.18 provides the evolution with time of the O(trans-AI)-O(DMSO) distance and figure 3.19 the evolution of the N(trans-AI)-O(DSMO) distance, as well as the corresponding radial distribution functions (RDF).Two main conclusions can be drawn from these graphs: (i) at short distances, one DMSO molecule is always hydrogen-bonded to the OH and NH 2 groups, although this interaction is dynamical with frequent exchanges, leading to a residence time of ca. 70 ps and 100 ps, respectively; (ii) the other molecules (three or four, depending on the observed conformation) solvate the trans-AI solute and are also in its vicinity. These observations point to the mobility of the solvent at room temperature.

  Radial distribution function, and its integral, between the oxygen atom of trans-AI and the oxygens atoms of 5 DMSO molecules. Distance between the oxygen atom of trans-AI and the oxygen atoms of DMSO molecules.

Figure 3 . 18 :

 318 Figure 3.18: Radial distribution functions and distance evolution as a function of time in the trans-AI:(DMSO) 5 1:5 complex, as obtained from the FFMD simulations at 300 K. Here the O(trans-AI)-O(DMSO) hydrogen bonds are examined.

Figure 3 . 19 :

 319 Figure 3.19: Radial distribution functions and distance evolution as a function of time in the trans-AI:(DMSO) 5 1:5 complex at 300 K. Here the N(trans-AI)-O(DSMO) hydrogen bonds are examined.

  (a) Radial distribution functions between trans-AI and oxygen of the 56 DMSO molecules in a periodic box. (b) Integrated radial distribution functions between trans-AI and oxygen of the 56 DMSO molecules in a periodic box.

Figure 3 . 20 :

 320 Figure 3.20: Radial distribution functions and their integrations for trans-AI in a box of 56 DMSO molecules at 300 K. The RDF between the oxygen of trans-AI and the DMSO molecules (red) shows an intense peak at 2 Å while the RDF between the nitrogen of trans-AI and the DMSO molecules (blue) shows an intense peak at 3 Å.

Figure 3 .

 3 Figure 3.21: AIMD study of non-local effects on the IR and VCD spectrum of trans-AI in DMSO solution by means of 1:1 clusters at the B3LYP-D3 level of theory; (a) radially resolved spectra (top: IR, bottom: VCD) of the chiral solute with spatial regions marked as "local" and "solvent effect", respectively; (b) IR (top) and VCD (bottom) spectrum of trans-AI without consideration of the solvent's polarisation (grey lines) and the difference spectrum as a consequence of the solvent (red lines).

Figure 4 . 1 :

 41 Figure 4.1: Visual depiction of (stable) neutral (a) and zwitterionic (a) forms of alanine in the gas phase.

Figure 4 . 2 :

 42 Figure 4.2: IR and VCD spectra of (a) zwitterionic and (b) neutral alanine averaged over 40 FFMD trajectories at T=300 K. Figure (a) emphasizes the influence of the polarisation of the molecule and of the time varying induced dipole moments moving around fixed positions term added to the magnetic dipoe moment. Figure (b) shows the spectra using the position and the velocity forms of the electric dipole moment.

Figure 4 . 3 :

 43 Figure 4.3: IR and VCD spectra from the average of 2 to 40 FFMD trajectories at 300 K of (a,b) the zwitterionic and (c,d) the neutral forms of alanine.

Figure 4 . 5 :

 45 Figure 4.5: IR and VCD spectra from the averages of 5 different trajectories of zwitterionic (left) and neutral forms of alanine (right), as obtained from FFMD trajectories at 300 K.

Family 1 (

 1 1 configuration) Family 2 (2 configurations) Family 3 (4 configurations) RMSD = 0.0 RMSD = 1.2257 RMSD = 0.9258 ∆E FF = 0.00 ∆E FF = 0.61 ∆E FF = 0.48 ∆E DFT = 0.00 ∆E DFT = 1.06 ∆E DFT = 1.45 µ tot FF = 0.450 µ tot FF = 1.370 µ tot FF = 1.859 µ tot DFT = 1.298 µ tot DFT = 1.626 µ tot DFT = 2.144

Figure 4 . 6 :

 46 Figure 4.6: Comparison of harmonic and anharmonic spectra obtained from static DFT calculations, for the (a) family 1, (b) family 2 and (c) family 3. The harmonic frequencies are multiplied by a scaling factor of 0.98. No scaling was applied to anharmonic frequencies.

Figure 4 . 7 :

 47 Figure 4.7: IR and VCD spectra obtained by averaging anharmonic DFT spectra from the 3 families using weights from the MD and DFT data. Represented in grey is the FFMD spectrum from the average of 40 trajectories.

Figure 4 . 8 :

 48 Figure 4.8: Comparison of the IR and VCD spectra of (a) zwitterionic and (b) neutral forms of alanine using the dipole moments calculated with the force field and with DFT from the same FFMD trajectory at 300 K.

Figure 4 . 9 :

 49 Figure 4.9: Correlation of the projections of the electric dipole moment ìµ or its time derivative ì µ, and the magnetic dipole moment ì m, as obtained for a specific MD trajectory but calculated from the force field (FF) or using DFT for the zwitterionic form of alanine.

Figure 4 . 10 :

 410 Figure 4.10: Correlation of the projections of the electric dipole moment ì µ or its time derivative ì µ, and the magnetic dipole moment ì m, as obtained for a specific MD trajectory but calculated from the force field (FF) or using DFT for the neutral form of alanine

Figure 4 .

 4 Figure 4.11: (a) IR spectrum of isolated alanine in zwitterionic form, as obtained from FFMD trajectories at 300 K and (b) its decomposition in effective modes, where the different power spectra are shown in color. The power spectra have arbitrary intensities.

Figure 4 . 12 :

 412 Figure 4.12: Effective modes associated to the spectral region 1000-1350 cm -1 of the IR spectrum of zwitterionic alanine in the gas phase. The sticks represent alanine in zwitterionic form. The arrows illustrate the collective atomic displacements along the effective modes.

Figure 4 . 13 :

 413 Figure 4.13: Effective modes associated to the spectral region 1350-1500 cm -1 of the IR spectrum of zwitterionic alanine in the gas phase.

Figure 4 . 14 :

 414 Figure 4.14: Effective modes associated to the spectral region 1500 -1900 cm -1 of the IR spectrum of zwitterionic alanine in the gas phase.

  torsion H2-C1-H3-H4 torsion H2-C1-H4-H3

  ) and the umbrella mode (mode 2 in figure 4.17). identified by the analysis are very similar to those in CH 3 , namely they are a HNH bending mode (mode 1, depicted as arrows in the inset of figure 4.18), the umbrella mode (mode 2 in figure 4.18), and a combination of a HNH bending mode with an out-of-plane vibration of the remaining NH bond (mode 3 in figure 4.18).

Figure 4 . 18 :

 418 Figure 4.18: Black line: IR spectrum of isolated alanine in its zwitterionic form. Coloured lines: power spectra determined by the effective modes analysis, corresponding to the 3 modes of the NH 3 fragment. Insets: stick representation of the NH 3 fragment, along with the effective modes, depicted as arrows.

Figure 4 . 19 :

 419 Figure 4.19: Black line: IR spectrum of isolated alanine in its zwitterionic form. Coloured lines: power spectra determined by the effective modes analysis, corresponding to the 2 modes of the COO fragment. Insets: stick representation of the COO fragment, along with the effective modes, depicted as arrows.

Figure 5 . 1 :

 51 Figure 5.1: Visual depiction of (a) crystalline [2] and (b) hydrated alanine [3].

Figure 5 . 2 :

 52 Figure 5.2: Visual depiction of the alanine crystal in its periodic box.

Figure 5 . 3 :

 53 Figure 5.3: IR and VCD spectra averaged over 40 trajectories with a periodic box containing 16 alanine molecules, as obtained from FFMD simulations at 300 K. The VCD spectrum of alanine in the gas phase is multiplied by 3, for visibility. The red arrows highlight the shifts of peaks corresponding to the NH 3 bend and the C=O stretch, between the two phases.

Figure 5 . 5 :

 55 Figure 5.5: IR and VCD spectra of alanine crystal in a periodic box containing 16 molecules, as obtained from FFMD trajectories at 300 K. The averages from 2 to 29 trajectories of 40 ps are displayed. The number of averaged trajectories increases with a step of 3. Represented in grey is the average of 30 trajectories.

Figure 5 . 6 :

 56 Figure 5.6: Normalised error of the average in the IR and VCD spectra of crystalline alanine from N T FFMD trajectories of 40 ps at 300 K, compared to the average from 30 trajectories, as a function of the sampling ratio S r = N T /30. The corresponding data for the gas phase system are superimposed in this figure.

Figure 5

 5 Figure 5.6 shows a slightly faster convergence for IR over VCD. The convergence of the spectrum from the crystal compared to the spectrum from the gas phase is faster for IR but slower for VCD. If only one molecule of alanine is considered for the solid, the convergence of the crystal is faster in both IR and VCD. The supramolecular interactions of the crystal may add instability to the magnetic dipole moment, thus disturbing the convergence.

Figure 5 . 7 :

 57 Figure 5.7: Comparison of the IR and VCD spectra of the alanine crystal using the dipole moments calculated with the force field and with DFT from the same FFMD trajectory at 300 K.

Figure 5 . 8 :

 58 Figure 5.8: Correlation of the projections of the electric dipole moment ì µ or its time derivative ì µ, and the magnetic dipole moment ì m, as obtained for a specific FFMD

5 . 9 . 139 Figure 5 . 9 :

 5913959 Figure 5.9: Visual depiction of hydrated alanine in its periodic box.

Figure 5 . 10 :

 510 Figure 5.10: Comparison of the intensities of IR spectra of water and of alanine in the gas phase.

Figure 5 . 11 :

 511 Figure 5.11: IR and VCD spectra of hydrated alanine averaged over 30 trajectories with cutoff values at R c = 0.0 Å, 3.85 Å and 4.6 Å, as obtained from FFMD trajectories at 300 K.

Figure a) has

  Figure a) has been enlarged in to produce figure b).

Figure 5 . 13 :

 513 Figure 5.13: IR and VCD spectra of hydrated alanine from the average of 2 to 30 trajectories with cutoff of 3.85 Å, as obtained from FFMD trajectories at 300 K.

Figure 5 . 14 :

 514 Figure 5.14: Normalised error of the average in the IR and VCD spectra of hydrated alanine from N T FFMD trajectories of 40 ps at 300 K, compared to the average from 30 trajectories, as a function of the sampling ratio S r = N T /30.

Figure 5 . 15 :

 515 Figure 5.15: Comparison of the IR and VCD spectra of hydrated alanine using the dipole moments calculated with the force field and with DFT from the same FFMD trajectory at 300 K, processed with a cutoff of R c = 3.85 Å.

  Here three figures comparing the force field and DFT moments are shown. The first figure(5.16) represents the moments of the whole system of alanine and the 211 water molecules. The second figure (5.17) displays only the moments of the alanine, ignoring the moments of the water molecules. The third figure (5.18), depicts the opposite, showing only the moments from the 211 water molecules. The Pearson correlation coefficients for the dipole moments of the water molecules in figure 5.18 are very high, above 0.96. This possibly comes from water being a system that has been studied deeply and for which the force field parameters are very accurate, in particular the multipolar distribution. The values of r from the dipole moments of only alanine, shown in figure 5.17, are slightly lower. The values of r from the dipole moments of water are reflected in those of alanine and water, in figure 5.16, as there are considerably more molecules of water than of alanine. In all three figures the values of r are very high.

Figure 5 . 18 :

 518 Figure 5.18: Correlation of the electric dipole moments ì µ, the velocity form of the electric dipole moments ì µ and the magnetic dipole moments ì m of the water molecules, ignoring the dipole moments from the alanine molecule.

Figure 5 . 21 :

 521 Figure 5.21: Effective modes in the region 1000-1200 cm -1 for the hydrated (h.) form of zwitterionic alanine (upper panels) and for the crystal (c.) form of zwitterionic alanine (lower panels).

Figure 5 . 22 :

 522 Figure 5.22: Effective modes in the region 1400-1650 cm -1 for the hydrated (h.) form of zwitterionic alanine (upper panels) and for the crystal (c.) form of zwitterionic alanine (lower panels).

Figure 5 .

 5 Figure 5.24: (a) IR and VCD spectra of alanine in the neutral form in the gas phase at 77.35 K and 300 K. (b) IR and VCD spectra of alanine in the neutral form in the gas phase and in N 2 at 77.35 K.

Figure 5 . 25 :

 525 Figure 5.25: Visual depiction of phenylcyclohexanediol in the crystalline phase and solvated in DMSO.

Figure 5 . 26 :

 526 Figure 5.26: IR and VCD spectra averaged over 20 FFMD trajectories of 80 ps at 300 K of crystal phenylcyclohexanediol and solvated in DMSO. For visibility the spectra of the crystalline phase were multiplied by 0.1 for the IR spectrum and 0.02 for the VCD spectrum.

Figure 5 . 27 :

 527 Figure 5.27: IR and VCD spectra of isolated phenylcyclohexanediol averaged over 20 FFMD trajectories at T =300 K. Both figures emphasize the influence of the polarisation of the molecule and of the time varying induced dipole moments moving around fixed positions term added to the magnetic dipole moment.

Figure 5 . 28 :

 528 Figure 5.28: IR and VCD spectra averaged of phenylcyclohexanediol in crystalline form and solvated in DMSO. The averages from 2 to 20 FFMD trajectories of 80 ps at 300 K are displayed. The number of averaged trajectories increases with a step of 2. Represented in grey is the average of 20 trajectories.

Titre:

  Caractérisation structurale des molécules chirales : développement d'outils de simulation moléculaire de spectres de dichroïsme circulaire vibrationnel ....................................................................... Mots clés : Dichroïsme circulaire vibrationnel, dynamique moléculaire classique, molécules fluxionnelles Résumé:Le dichroïsme circulaire vibrationnel (VCD) est la faible différence d'absorption des molécules chirales entre la lumière polarisée droite et gauche dans le domaine de l'infrarouge (IR). Cette technique spectroscopique est devenue de plus en plus populaire depuis sa découverte dans les années 1970 et a des applications prometteuses en pharmacologie en raison de sa capacité à déterminer les configurations absolues des molécules chirales. La forme des spectres VCD est très sensible aux changements mineurs de conformation et aux interactions moléculaires, ce qui en fait une sonde de l'isomérie conformationnelle et de la solvatation. L'attribution des spectres VCD expérimentaux nécessite en général une comparaison avec les calculs théoriques, ce qui implique également que le solvant soit précisément décrit. Ce manuscrit décrit comment une modélisation basée sur des champs de force polarisables peut être exploitée pour interpréter les signaux VCD. Cette modélisation est appliquée à des molécules flexibles pour lesquelles les approches conventionnelles de chimie quantique reposant sur des approximations statiques sont limitées. En particulier, il est difficile de prendre en compte les anharmonicités et les effets de température finie, et le solvant est typiquement modélisé de manière simplifiée par un continuum ou par un nombre limité de molécules explicites. Dans le travail présenté ici, le champ de force AMOEBA mis en oeuvre dans le logiciel Tinker a été utilisé pour modéliser les moments dipolaires électriques et magnétiques de diverses molécules pour des applications en phase gazeuse, ainsi que dans des systèmes solvatés ou cristallins. L'échantillonnage et les informations dynamiques ont été obtenus par simulations de dynamique moléculaire classique. Les spectres IR et VCD anharmoniques obtenus à partir des trajectoires ont été interprétés à l'aide d'une analyse des modes effectifs dans laquelle les coordonnées internes sont décomposées linéairement par projection sur des modes effec-tifs appropriés, les poids correspondants étant obtenus de telle sorte que le pic de chaque mode soit le plus localisé possible. Cette analyse permet de mettre en évidence les contributions des différents groupes fonctionnels de la molécule aux pics individuels du spectre. Comme premier exemple, le cas du trans-1-amino-2indanol solvaté dans le diméthylsulfoxyde a été examiné en combinant des explorations statiques au niveau de la théorie de la fonctionnelle de la densité, des trajectoires courtes de MD ab initio (AIMD) et un large échantillonnage de la surface d'énergie potentielle à l'aide du champ de force. La comparaison avec les mesures existantes confirme l'importance de la flexibilité pour cette molécule, le rôle des molécules de solvant proches et leur contribution aux signaux VCD. La détermination des spectres VCD à partir du champ de force polarisable a ensuite été entreprise et ses performances ont été évaluées par comparaison avec des calculs de structure électronique. Cette comparaison a été effectuée sur la même trajectoire MD, en utilisant l'alanine en phase gazeuse comme système de référence. Les différentes contributions aux moments dipolaires sont discutées. Après avoir validé la méthodologie, la convergence statistique des spectres IR et VCD a été quantifiée via des mesures d'erreurs dédiées. Ces erreurs montrent à leur tour que la convergence des spectres VCD est plus lente que celle des spectres IR et quelle n'est atteinte qu'après un échantillonnage bien plus important que celui réalisable par AIMD. L'approche par champ de force polarisable a ensuite été appliquée à l'alanine en phases condensées. Là encore, la comparaison avec les simulations AIMD confirme que les moments dipolaires électriques et magnétiques peuvent être évalués de manière fiable. Ce travail ouvre la voie à la détermination du spectre VCD de systèmes chimiquement complexes inaccessibles via des méthodes qui tiennent compte explicitement de la structure électronique. Title: Structural characterisation of chiral molecules: tool development for molecular simulations of vibrational circular dichroism spectra ......................................................................................................

  is the static electric dipole moment of the molecule in its final state. The offdiagonal elements, ì µ fi and ì µ if , are the electric transition dipole moments, which enable coupling between the states.

Table 4 .

 4 6: List of CH 3 internal coordinates used for the specific analysis of the methyl group power spectrum.

Acknowledgements

The variations of the normalised error, shown as a function of the sampling ratio S r = N T /40, demonstrate that the sampling must be performed for a certain amount of time to reach convergence, see figure 4.4. When the sampling ratio reaches 0.5, hence when the system has been sampled for 20 trajectories of 200 ps which corresponds to 4 000 ps, the spectrum is considered to be converged.

To check whether there is an influence on the spectra depending on which trajectories are chosen, 5 trajectories were selected randomly between 40 trajectories for the zwitterionic and neutral forms of alanine. This was performed 3 times (sample 1-3). The average of 5 trajectories for each sample is shown in figure 4.5. Choosing different trajectories influences the intensity of the peaks especially for alanine in the neutral form. This confirms that the sampling needs to be performed over multiple trajectories to correctly access the spectral variations of the system. ω2=1150 cm -1 ω3=1200 cm -1 ω12=1706 cm -1 ω8=1486 cm -1 Figure 4.16: Effective modes associated to the spectral region 1000-1800 cm -1 of the IR spectrum of neutral alanine in the gas phase, highlighting to the major differences with the IR spectrum of zwitterionic alanine in gas phase.

Alanine_aminoacid_newParams

Frequency (cm

Alanine_gaz_plus_CH3-CH3

Frequency (cm ω1=1465 cm -1 ω2=1467 cm -1

CH3 in zwitterionic alanine Figure 4.17: Black line: IR spectrum of isolated alanine in its zwitterionic form, as obtained from FFMD trajectories at 300 K. Coloured lines: power spectra determined by the effective modes analysis, corresponding to the 3 modes of the CH 3 fragment. Insets:

stick representation of the CH 3 fragment, along with the effective modes, depicted as arrows.

IR analysis of the amino group

For the effective modes analysis of the amino group of zwitterionic alanine, the following internal coordinates are used. 4.7: List of NH 3 internal coordinates in zwitterionic alanine used for the specific analysis of the amino groups power spectrum.

The IR absorption bands of the NH 3 fragment are slightly blue-shifted with respect to those of CH 3 , and are mainly distributed in the region 1500-1700 cm -1 , as shown by the coloured bands in the lower panel of figure 4.18. Despite the blue shift, one might expect that the vibrations within the NH 3 fragment are coupled to those of CH 3 since the absorption bands partially overlap. Therefore, NH 3 contributes, similarly to CH 3 , to the large absorption region of zwitterionic alanine in gas phase in the region 1300-1600 cm -1 (black line in the upper panel of figure 4.18) as well as to the isolated band at 1700 cm -1 . The bands are not well localised in frequency as it was the case for the CH 3 fragment, which might be a manifestation of strong anharmonicity and coupling of NH 3 with other modes that have not been considered in this fragment analysis. The effective modes

Comparison to experimental results

The previous IR and VCD spectra of crystal alanine are compared to experimental data [1] in figure 5.4. The experimental data were obtained from ground crystals as nujol oil mull. These data were only available up to 1800 cm -1 . [1] and theoretical results from the force field MD approach at 300 K.

In the IR spectra, the peaks match between 1400-1800 cm -1 , below this range, the peaks are shifted. For the VCD spectra, a similar pattern of a positive then negative peak is found around 1400 cm -1 .

Statistical errors and convergence rates

The spectra of the alanine crystal shown above result from the average of 30 trajectories. To determine if the system has been sampled for a sufficient amount of time, the convergence is now examined. To show when the intensities of each peak stabilise at their converged values after a certain sampling time, the averages from an increasing number of trajectories (N T ) are plotted. Starting from the top of the figure 5.5, the spectra arise from the average from 2 trajectories (N T = 2). Then the number of trajectories increases by 3 stepwise. Each spectrum was compared to the reference spectrum obtained from the entire set of 30 trajectories.

The average of the spectra from 30 trajectories was considered converged and was taken as our reference to which other spectra with lesser sampling are compared. Thus the spectra from 1 < N T < 30 trajectories were compared to this averaged spectrum to determine their normalised error, with the method and equations shown in chapter 4.

This error is shown in figure 5.6 as a function of the sampling ratio S r = N T /30.

to average out to zero. The second reason can be the presence of chirality transfer from the alanine to the solvent, at least in the vicinity of the solute.

In the following, the comparison between the experimental and theoretical spectra is examined. The experimental spectrum [1] has been background subtracted with a water spectrum and the theoretical results correspond to the FFMD spectrum depicted in figure 5.11 with a cutoff value of R c = 0.0 Å, as the water is also subtracted in the experimental spectrum. These two results are compared in figure 5.12 within frequency range of 1000 to 1800 cm -1 .

IR VCD

Theoretical spectrum with Experimental spectrum Figure 5.12: IR and VCD spectra of hydrated alanine from experimental [1] and theoretical (FFMD) results at 300 K. The experimental results are background subtracted with a water spectrum and a cutoff value of R c = 0.0 Å is used for the theoretical spectra.

The IR spectra are quite similar especially above 1400 cm -1 . The low resolution of the theoretical spectrum is probably due to the smoothing performed during the Fourier transform of the correlation functions. This smoothing parameter will need to be adjusted in future computations. The main discrepancy is at 1050 cm -1 . The intense double peak in the experimental VCD spectrum is also represented in the theoretical results, similarly to the spectra of the alanine crystal. However there is a slight shift for the negative peak.

Statistical errors and convergence rates

The spectra of hydrated alanine shown above are the result from the average of 30 trajectories. The averages of an increasing number of trajectories (N T ) are plotted to examine the convergence and are compared to the full average in figure 5.13. Between successive spectra three new trajectories are added. A power spectrum with an arbitrary intensity and a vectorial depiction of each effective mode were acquired, they are compared to the FFMD IR spectrum in figure 5.19 and 5.20. The same set of internal coordinates is used as for the zwitterionic form of alanine in the gas phase. The goal of the analysis reported below is to comment on those similarities and differences based on the effective modes analysis. The main differences between the absorption spectra of hydrated and crystalline alanine are the following:

Alanine_solid

Frequency (cm

1. in the region 1000-1200 cm -1 , the two main bands are well-resolved in the crystal form, as well as the shoulder at around 1180 cm -1 , with the identification of three distinct effective modes bands, whereas, in the hydrated form, the two main bands are predicted to arise from the superposition of three modes, while a fourth mode is responsible for the shoulder, which seems, however, suppressed in water;

2. in the region 1400-1650 cm -1 , the lower intensity bands at the sides on the central most intense bands are both split into two contributions in the crystal form, while only one band appears to be present on each side of the intense band when alanine in solvated in water;

3. in the region around 1700 cm -1 , a large and not well-resolved band appears in the hydrated form of alanine, which is completely absent in the crystal form.

Additionally, if compared to the zwitterionic form of alanine in the gas phase, the band at around 1300 cm -1 is slightly blue-shifted in water and in the crystal, while, at the same time, the band at 1900 cm -1 is red-shifted when the molecule is not isolated.

For the decomposition of the hydrated spectrum, there is no effective mode for the peak at 1700-1800 cm -1 as this corresponds to the water asymmetric stretch, not accounted for in the analysis.

In order to analyse the points just discussed, in figure 5.21 some of the effective modes giving rise to the IR absorption bands in the region 1000-1200 cm -1 are reported. For the normalised error, in figure 5.29a, the two phases of phenylcyclohexanediol are compared, then, in figure 5.29b, the errors of phenylcyclohexanediol crystal are shown against those of the alanine crystal. The convergence is faster for the crystalline phase than for the solvated phase. This was also the case for crystalline alanine (see figures 5.6 and 5.14) and reflects the improved convergence that results from the presence of multiple molecules with strong interactions within the crystal, making the molecules less fluxional. When comparing the two crystals, again the alanine crystal shows lower IR errors, this time against phenylcyclohexanediol, but higher VCD errors.

The sampling time necessary for the spectra of the phenylcyclohexanediol crystal and of phenylcyclohexanediol solvated in DMSO to converge is 10 trajectories of 80 ps. This time is longer than for the alanine crystal and the hydrated alanine in the case of vanishing cutoff. This is due to the high flexiblity of phenylcyclohexanediol resulting in a large number of configurations that need to be all explored.

The dipole moments of the phenylcyclohexanediol systems remain to be compared to DFT dipole moments to confirm their accuracy. An effective modes analysis will also be performed to examine, with more detail, in the future, which band is assigned to each molecular vibration.