Weight-Averaged Anharmonic Vibrational Analysis of Hydration

Article pubs.acs.org/JPCB

Weight-Averaged Anharmonic Vibrational Analysis of Hydration Structures of Polyamide 6 Bo Thomsen,† Tomonori Kawakami,‡ Isamu Shigemoto,‡ Yuji Sugita,†,§,∥,⊥ and Kiyoshi Yagi*,†,§ †

Theoretical Molecular Science Laboratory, RIKEN, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan Advanced Materials Research Laboratories, Toray Industries, Inc., 2-1 Sonoyama 3-chome, Otsu, Shiga 520-0842, Japan § RIKEN iTHES, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan ∥ RIKEN Advanced Institute for Computational Science, 7-1-26 Minatojima-Minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan ⊥ RIKEN Quantitative Biology Center, 6-7-1 Minatojima-Minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan ‡

S Supporting Information *

ABSTRACT: Structures of polyamide 6 are investigated for different hydration levels using molecular dynamics (MD) simulations and quantum vibrational calculations. The MD simulations have shown that hydration leads to an increase in the diffusion coefficient, accompanied by a growth of water clusters in the polymer. The IR difference spectra upon hydration are calculated using a weight-averaged method incorporating anharmonicity of the potential energy surface. The predicted IR difference spectrum for the amide A band is in quantitative agreement with the experiment [Iwamoto, R.; Murase, H. J. Polym. Sci., Part B: Polym. Phys. 2003, 41, 1722−1729]. The proposed method, combined with experimental IR difference spectra, makes it feasible to elucidate the atomistic structure of hydrated polymer materials. Raman spectra, respectively, of the α and γ polymorphs of crystalline PA6. The frequency-scaled calculated spectra were in agreement with the experimental spectra recorded from 2000 to 400 cm−1. The agreement allowed for assignment of the experimentally recorded spectrum to combined motions of the crystalline material. Although these results are promising for the treatment of crystalline systems with periodic boundary conditions, their extension to armorphous hydrated materials and vibrational motions with a strong anharmonic character is currently unexplored. The harmonic approximation is generally not sufficiently accurate for vibrational motions with a high degree of anharmonicity, for example, hydrogen stretches. Anharmonic calculations require ab initio potential energy surfaces (PESs) and a solver of the vibrational Schrödinger equation (VSE) to compute highly accurate spectra.18−24 A system with N atoms contains 3N − 6(5) vibrational degrees of freedom, or modes. The generation of PES and subsequent solution of the VSE both scale exponentially with respect to the number of modes, making the full-dimensional calculation only tractable for systems of less than 10 atoms. An efficient approximation to circumvent the scaling problem is to expand the PES in terms of mode couplings. The accuracy of this approximation depends on the order of the couplings,

1. INTRODUCTION Understanding water−polymer interactions is essential to design and control functions of polymer materials.1−3 The hydration of polyamide 6 (PA6), also known as nylon-6 or polycaprolactam, has been considered by a mechanism proposed by Puffr and Šebenda.4 Their model contains three types of water molecules forming hydrogen bonds (H-bonds) to, (1) the CO group of PA6, (2) the N−H group of PA6, and (3) other water molecules without sharing any H-bonds with PA6. The third type of water molecule is denoted “free water”. When PA6 is poorly hydrated, the water molecules initially form H-bonds to CO, resulting in N−H without Hbond partners. These N−H donors subsequently form H-bonds to additional water molecules. Finally, “free water molecules” start to aggregate in a fully hydrated PA6. The polymer−water dynamics of PA6 has been probed using neutron and X-ray scattering,5−7 NMR,8−12 and IR spectroscopy.5,13,14 The Hbonds between PA6 and water molecules are directly probed using IR spectroscopy by monitoring the vibrational frequencies of amide groups (e.g., CO stretching, N−H bending, and N−H stretching modes), which are highly sensitive reporters of the H-bond formation. However, the interpretation of the spectra is nontrivial and thus requires extensive theoretical studies. Explicit quantum mechanical calculations of the harmonic vibrational spectrum for crystalline periodic systems are possible using the density functional theory (DFT).15 With this approach, Quarti et al.16 and Milani17 calculated the IR and © 2017 American Chemical Society

Received: January 12, 2017 Revised: March 15, 2017 Published: June 9, 2017 6050

DOI: 10.1021/acs.jpcb.7b00372 J. Phys. Chem. B 2017, 121, 6050−6063

Article

The Journal of Physical Chemistry B

2. METHODS 2.1. Anharmonic Vibrational Structure Calculations for a Chromophore. The bulk system is assumed to contain certain functional groups responsible for the IR absorption, denoted as chromophores. Here, we assume that other groups do not contribute to the absorption profile. The vibrational motion of the chromophore consisting of N atoms is described using a set of normal coordinates generated by constructing and diagonalizing a partial mass-weighted Hessian

which in turn depends on the coordinates used to represent the PES. One of the recent examples of using an efficient set of coordinates is the work by Bowman and co-workers, who reported a theoretical IR spectrum of bulk water25−28 and clathrate hydrates29,30 using normal modes localized to individual water molecules.31 To model systems with a larger diversity in composition, more drastic approximations are needed. Partial Hessian diagonalization32 is such an example, which has been extended to calculate the anharmonic frequencies of complex molecules.33 In the optimized coordinate vibrational self-consistent field (oc-VSCF) method,34 the coordinates are optimized by minimizing the vibrational ground-state energy with respect to a coordinate transformation matrix using a Taylor expanded potential. It has been shown that optimized coordinates increase the accuracy of VSE solvers.34−36 More recently, a flexible adaption of local coordinates of nuclei (FALCON) method 37 has been proposed, which is initiated by a set of partial Hessian diagonalizations for groups of atoms in a given molecule, and modes of the surrounding atoms are added iteratively until the Hessian is (block) diagonalized within a given threshold. Schemes for localizing the vibrational modes have shown promise for solving the VSE with high efficiency and accuracy for large systems. Unlike simple molecular clusters or solids, complex systems have an enormous number of structures (such as H-bond patterns) that contribute to the vibrational spectrum. Recently, we have proposed a weight-averaged approach to calculate the IR and Raman spectra of sphingomyelin (SM) bilayers at different phases.38,39 In this method, molecular dynamics (MD) simulations are first carried out to sample H-bond patterns of amide groups in SM and water molecules. Then, for each pattern, the MD trajectory is used to construct a cluster model, and the harmonic vibrational spectrum is computed using DFT. Finally, the spectra of all the patterns are summed with their statistical weights, which are determined from the MD trajectories. The calculated spectrum reproduced the experiment in the amide I region (CO stretching vibration) and provided firm evidence that SM molecules form H-bond clusters in a bilayer.39 In this study, the weight-averaged method is extended to the anharmonic theory. A set of local coordinates is created for the amide groups of PA6 and water molecules. These coordinates are then used to generate the PES and to solve VSE by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2)24 method. This method is, in particular, powerful for studying the N−H and O−H stretching modes of water molecules and PA6, respectively, and resolving the vibrational overtones, wherein anharmonicity plays an essential role. The following section describes the method for calculating the anharmonic frequency and IR absorption cross-section for the vibrational transitions. The spectral information is combined with structural information obtained from MD simulations to calculate the vibrational absorption coefficient of a bulk system. Section 3 gives the computational details of MD simulations and quantum vibrational calculations for hydrated PA6 systems. Section 4 contains the analysis of the PA6 MD simulation and the calculated IR difference spectra. A discussion on the current results and their predictions about the atomistic hydration structures of PA6 is given in Section 5. The concluding remarks and perspective applications of this method are discussed in Section 6.

(1)

LHL = Ω

where ∂ 2Eel 1 mimj ∂χi ∂χj

Hij =

(2)

Eel is the electronic energy, χi is the Cartesian coordinate, and mi is the mass associated with the ith atom in the chromophore. Equation 1 yields a set of 3N coordinates {Ql}, which forms the basis for calculating the vibrational properties. The coordinates are obtained as row vectors of L, with corresponding harmonic vibrational frequencies, ωl, from the diagonal elements of Ω. The set of normal coordinates can be further reduced to a set of F active coordinates. The VSE in atomic units for a chromophore is given in terms of {Ql} as ⎡ F ⎤ ∂2 ⎥Ψ = E Ψ Ĥ Ψn = ⎢ ∑ + V ( Q , ..., Q ) n n 1 F ⎥ n 2 ⎢⎣ ∂ Q ⎦ f f =1

(3)

where n denotes a state vector with integer elements nf. Note that the PES implicitly depends on the environment of a chromophore. The VSCF18,19 method is often used to generate a zeroth-order approximate solution to the VSE. A set of coupled one-mode equations is solved in the VSCF method until convergence with respect to energy is reached for all modes. The total VSCF wave function is a product of the onemode functions F

Φn =

F

M

∏ ϕn(f )(Q f ) = ∏ ∑ c (mf ()n ) χm(f()n ) f

f

f =1

f

(4)

f =1 m

Each of the one-mode wave functions is expanded as a linear ( f) combination of basis functions, χm(n . The one-mode equation f) for mode Qf for optimizing the free parameters, c(f) m(nf), is given as ⎡ 2 ⎤ ⎢ ∂ + V̅ (Q )⎥ϕ(f )(Q ) = ϵ(f )ϕ(f )(Q ) nf nf f ⎥ nf f f ⎢⎣ ∂Q 2f ⎦

(5)

ϵ(nff) is the one-mode energy, and ϕ(nff)(Qf) is the one-mode function for the state with vibrational quantum number nf. The mean field operator, V̅ (Qf), is an anharmonic potential in which the contribution from all other modes is integrated out F

F

V̅ (Q f ) = ⟨∏ ϕn(g )(Q g )|V (Q 1 , ..., Q F )| ∏ ϕn(g )(Q g )⟩ g

g≠f

g

g≠f

(6)

The VSCF method does not account for vibrational correlation corrections, thereby resulting in large errors for correlated states. VQDPT224 has been developed to incorporate the correlation effects. The method combines the contribution of 6051


Article


Figure 1. Chemical structure of the chains of PA6 used in the current study, each of the MD simulations contain four of these 52-mers of PA6 as well as a set amount of water molecules.

the static correlation in vibrational configuration interaction20 methods with the dynamic correlation of vibrational secondorder Møller−Plesset (VMP2).21 VQDPT2 divides the vibrational state space into two subspaces P = {p} ∧ Q = 1 − P = {q}

Γ(ν − νsi) =

⎛

∑ ⎜⎜

(7)

qϵQ

1

(0) (0) ⎝ Ep − Eq

⟨p|Ĥ |q⟩⟨q|Ĥ |p′⟩ 2

+

⎞ ⎟ − Eq(0) ⎟⎠ 1

E(0) p′

Es − E0 hc0

ct =

i

t

(13)

nt V

(14)

3. APPLICATION TO PA6 3.1. Molecular Dynamics. NAMD40 was used to run MD simulations of five PA6 membranes. Each contained 0, 32, 63, 126, and 252 water molecules, corresponding to 0.0, 2.3, 4.5, 8.6, and 15.9 wt % water in the membrane, respectively. All systems contain four PA6 chains consisting of 52 repeating units and capping groups at both terminals. See Figure 1 for the chemical structure. The CHARMM 36 lipid force field41 was used for the aliphatic carbon chain in PA6. The amide parameters of the force field were adapted from the CHARMM 2242 force field and the SM force field.43 The details of the PA6 force field are given in the Supporting Information (SI). The water molecules were modeled using the TIP3P parameters.44 The MD simulations were carried out using the smooth Particle-Mesh Ewald method.45,46 The time step for integration in all calculations was set to 1 fs. The nonbonding interactions were reduced to 0 between 10 and 12 Å, using a switching function. The neighbor list was updated every 10 fs using a pair list cutoff length of 14 Å. A Langevin thermostat with a coupling time set to 1 ps−1 was used to control the target temperature. The isotropic pressure was maintained by the Nosé− Hoover Langevin piston method, 47−49 with the oscillation period time and damping time set to 0.2 and 0.1

(10)

N

1 c ∑ ∑ σs Γ(ν − νsi) V i=1 s i

t

st

nt is the number of a cluster type, t, at a given time, which is obtained from an MD trajectory by grouping the chromophores to cluster types and counting their occurrences.

where ε0 is the permittivity of vacuum, x, y, and z correspond to the Cartesian directions, and μ̂ ζ is the dipole operator of the ζ component. 2.2. Weight-Averaged Vibrational Spectrum. Assuming that the system contains Nc chromophores, the absorption coefficient of the system is written as α (ν ) =

(12)

where ct is an ensemble average concentration of a cluster type, t, calculated as

(8)

(9)

2π 2νs ∑ |⟨s|μζ̂ |0⟩|2 3hε0c0 ζ = x , y , z

)

∑ ct∑ σs Γ(ν − νs ) t=1

Es and E0 are the energies of states s and 0, respectively, h is the Planck constant, and c0 is the speed of light in vacuum. The absorption cross-section between states s and 0 for an isotropic system is given by σs =

γ2 4

Nt

α (ν ) ≅

Ĥ (2) eff is a block diagonal matrix, where each block corresponds to a separate P-space. If a P-space only contains a single vibrational state, VQDPT2 is equivalent to VMP2. The VSE is solved by diagonalizing the Ĥ (2) eff matrix. The vibrational transition energy from the vibrational ground state, n = 0, denoted 0, to a vibrational excited state, n ≠ 0 denoted s, in wave numbers is given as

νs =

(

π (ν − νsi)2 +

γ is a broadening factor that can be chosen independently for each system and state. Note that the absorption coefficient, α(ν), is proportional to the absorbance, the quantity often reported in an experiment, within a constant factor. The chromophore and its environment with the same chemical composition (e.g., the chromophore and its first hydration shell) are denoted as a cluster. Assuming that all of the chromophores of a certain cluster type have vibrational transitions with the same wave number and absorption crosssection, eq 11 is rewritten as a sum over cluster types, Nt

The P-space consists of states with strong, static correlation internally. This space is chosen on the basis of a set of target states. The Q-space contains the states contributing dynamical correlation to the P-space states. The method for selecting the P and Q spaces is detailed in previous studies by one of the authors.24,36 The VQDPT2 Hamiltonian is given as follows (2) (Ĥ eff )pp ′ = ⟨p|Ĥ |p′⟩ +

γ 2

(11)

V is the volume of the system, the bracket represents an ensemble average, and Γ is a Lorentzian lineshape function 6052


Article

The Journal of Physical Chemistry B ps, respectively. Each system was equilibrated in the following procedure: First a 1000 step minimization was carried out. Then five sequential NVT calculations raising the temperature from 10 to 600 K over 500 ps were run, wherein the temperature was incremented in equidistant steps every 1 ps. Lastly, a 900 ps NPT dynamics simulation at 300 K and 1 atm was conducted. The production runs were done over 1 μs in the NPT ensemble at 1 atm and 300 K. The atomic coordinates were saved every 10 ps in the production run, resulting in 100 000 structural snapshots used for the analysis of each trajectory. The diffusion coefficient (D) of water was calculated by the Einstein formula (R t − R 0)2 = 6Dt

(15)

where the left-hand side corresponds to the mean squared displacement (MSD). t denotes time, R0 is the position of a water molecule at time zero, and the brackets indicate an average over all water molecules in the system. The position of a water molecule at a given time, Rt, is corrected for the motion of the center of mass (CM) of the heavy atoms (O, C, N) of the entire system R t = R(Wt ) − R(CM t )

(16)

Figure 2. Examples of the first hydration shell of an amide and water chromophore in the hydrated PA6 system. In the case of an amide, the cluster type highlighted by the dotted squares is 1A−A−1A, whereas the water chromophore showcases a 1A1W−W−2A cluster type.

where R(Wt) is the CM position of the water molecule at time t, and R(CMt) is the CM of all heavy atoms at time t. Averaging over every 20 ns of the MD trajectories and subsequently fitting the last 14 ns of the average to a linear function generate the diffusion coefficient. 3.2. Clustering the PA6 and Water System. The current study uses both a fine and a coarse clustering method of the MD trajectory data, both based on H-bonded structures of the chromophores. In this study, the chromophores are selected as the amide groups of PA6 and the water molecules, as these are expected to give the dominant contribution to the IR difference spectra in the region from 1400−4000 cm−1. The H-bond between the acceptor (A) and donor (D-H) is defined as R(A‐D) ≤ 3.2 Å and ∠(DHA) ≥ 130

3.2.2. Coarse Clustering. The coarse clustering method organizes the fine clusters with amide as a chromophore according to the following four types of H-bonds: V denotes a H-bond site lacking a partner, AA denotes an acceptor or donor atom of the amide participating in a H-bond with another amide, AWA denotes CO accepting a H-bond from water, and AWD denotes N−H donating a H-bond to water. Using the ensemble-averaged concentrations of the amide chromophores (eq 14), the probability in percent, pT, for each of the H-bond types, T, can be calculated as n c ∑ tT t × 100% = pT 3Camide (18) t ϵT

(17) 50

according to the definition of Desiraju and Steiner. Oxygen is defined to accept two H-bonds, and polar hydrogen is capable of donating to one oxygen acceptor. If more acceptors or donors are found, the ones with the shortest acceptor−donor distance are chosen. All clusters containing a capping group are discarded, because including these structures would overestimate the effect of the terminals. 3.2.1. Fine Clustering. Fine clustering is used to generate the clusters and concentrations for calculating the IR difference spectrum. Figure 2 describes the naming conventions and structure of the cluster types used in the fine clustering method. The cluster types are named by the functional groups donating H-bonds to the chromophore, the chromophore itself and the functional groups the chromophore donates H-bonds to. For example, a cluster type 1A−A−1A means that the chromophore accepts a H-bond from an amide group (1A−), and the chromophore is an amide (A), which donates a H-bond to another amide (−1A). Another cluster type 1A1W−W−2A means the chromophore accepts two H-bonds, one from water and another from amide (1A1W−), and the chromophore is water (W), which donates both of its OH bonds to two amide groups (−2A). Note that there is no distinction between which hydrogen atom of water is donating a H-bond, as there is no difference between the two hydrogens because of symmetry.

Camide is the sum of the ensemble-averaged concentrations of amide chromophores in the system. ntT is an integer denoting how many H-bonds of type T the chromophore t participates in. As an example, the time-averaged concentration of the 1A− A−1W cluster type would contribute, with ntT = 1, to the probabilities pAA, pV, and pAWD. 3.3. Electronic and Vibrational Structure Calculations. On the basis of the sum of their ensemble-averaged concentrations in the five trajectories, 40 dominant cluster types were chosen for the anharmonic calculation. N-methylacetamide (NMA), CH3NHCOCH3, is used to model the strands of PA6 in the electronic structure calculations. The B3LYP functional was used, with the 6-31++G** basis set for the amide group and water. The smaller 6-31+G* basis set was used for the methyl groups of NMA. The structures of each cluster type were optimized, and the Hessian matrix was calculated and diagonalized for the chromophore coordinates. The mode spaces of the chromophores were restricted to the water-bending mode and two stretching modes, and amide I− III and A modes for amide. The grid surfaces for the one-mode and two-mode PESs and dipole surfaces were constructed using 6053


Article

The Journal of Physical Chemistry B 10 and 100 single point calculations for each set of surfaces. All electronic structure calculations were done using Gaussian09.51 The VSCF calculations were done using a basis set of the 10 lowest-energy harmonic oscillators in each mode. VQDPT2 calculations were then conducted, targeting in the case of water the bending fundamental and overtone and the symmetric and asymmetric stretch. In the case of amides, the fundamental of amide A and the fundamentals, first overtones, and combination bands of amide I−III were targeted. The P and Q spaces used in VQDPT2 were constructed by considering three generations with k = 4 in the algorithm described in ref 36. The anharmonic vibrational frequencies and absorption cross-sections of the 40 cluster types are given in the SI. 3.4. IR Difference Spectra. The absorption coefficient for the PA6 membrane with different levels of hydration is calculated by eq 13. The harmonic spectrum is calculated on the basis of the frequencies arising from diagonalization of the local Hessian and intensities using the dipole derivatives from the electronic structure calculations. A uniform broadening factor of 40 and 80 cm−1 is used for the Lorentzian function in eq 12 in the amide I−II and amide A regions, respectively. These two factors are chosen so as to reproduce the overall shape of the experimental spectra. The IR difference spectra are calculated by subtracting the spectrum of the dry membrane and the spectrum of the membrane containing 8.6 wt % water unless otherwise noted. This level of hydration is selected because the absorption of water in PA6 has been reported to be 7.4 wt % in ref 8 and 4−10 wt % in references therein. The theoretical difference spectra are compared with the experiment by Iwamoto and Murase,14 who studied the change in absorbance as water evaporated from a fully hydrated 5 μm thick PA6 membrane. The experimental difference spectra were taken between the membrane 0, 1.5, 3.3, and 39 min after hydration, and one was taken 72 min after hydration. The spectrum 0 min after hydration is used for comparison with the theoretical spectra, unless otherwise noted. The data presented here was obtained through WebPlotDigitizer52 using a highresolution scan of the experimental IR difference spectra, given in Figures 2 and 5 of ref 14.

Figure 3. Snapshots of the unit cells of the PA6−water system at different water concentrations. Note that the difference in size of the unit cell in different water concentrations is an effect of adding more water into the membrane. The figures were generated using VMD.56

Table 1. Water Diffusion Coefficients and Densities of the PA6−Water Systems at Different Water Concentrationsa water concentration diffusion coefficient (10−10 cm2/s) density (g/cm3)

2.3 wt % 4.5 wt % 8.6 wt % 2.4

2.5

1.09

1.10

13 1.11

15.9 wt %

exp.8

31

1−10

1.11

1.08

Diffusion coefficients are found by fitting a straight line through the data presented in Figure 4 from 6 to 20 ns. The density stems from the mass divided by the average volume of the PA6−water systems. a

4. RESULTS 4.1. Hydrated Structures of PA6. Representative snapshots of the MD trajectories of the increasingly hydrated PA6 are shown in Figure 3. At low hydration levels, water molecules appear alone or in small clusters, all interacting strongly with PA6 chains through H-bonds. Increasing the water concentration causes the formation of long water chains connected by H-bonds. However, water molecules in such chains do retain their strong interaction with the PA6 chains. The calculated density of the PA6 membrane and water, shown in Table 1, is close to the experimental density of amorphous PA6. The density is insensitive to the increase in the concentration of water due to the swelling of the PA6 membrane making room for the additional water molecules. In the movie included in the SI, water molecules diffuse in jumps from one H-bound moiety to another. A jump is followed by a period of constrained rotational motions, retaining the H-bonds formed in the moiety, before another jump is made. Experimental studies have reported the diffusion coefficient of water in PA6 between 1.0 × 10−10 and 1.0 × 10−9 cm2/s at 25 °C.8 The calculated MSDs of the water molecules are given in Figure 4, and the resulting diffusion coefficients are given in Table 1. The results are close to the experimental

values. In the low hydration levels of 2.3 and 4.5 wt %, the diffusion coefficient remains unchanged around 2.5 × 10−10 cm2/s. Increasing the water concentration to 8.6 wt % increases the diffusion coefficient approximately five-fold to 1.3 × 10−9 cm2/s. The rise in the water diffusion coefficient continues when the water concentration is increased to 15.9 wt %, where the diffusion coefficient is found to be 3.1 × 10−9 cm2/s. The behavior of the diffusion constant with respect to hydration found here reflects the experimental findings,8 which reports an exponential rise in water diffusion as the water concentration in PA6 rises above 5 wt %. Using the coarse clustering method described in Section 3.2.2, the H-bonds of the amide groups in PA6 have been calculated for each water concentration in the current study. The results are given in Figure 5. In the dry PA6 membrane (0.0 wt %), half of the H-bond acceptor and donor sites lack a partner. The amount of vacant H-bond sites and H-bond sites occupied by other amides decreases when water is added to the membrane. Puffr and Šebenda4 suggested that water favors forming H-bonds to CO over N−H. This trend is reflected in the fact that pAWA occurs slightly more than twice as 6054


Article


Figure 4. MSD of water inside the PA6 membrane at various levels of hydration. The diffusion constant given in Table 1 is calculated by the linear fitting of a range between 6 and 20.

Figure 5. Probability of the amide group of PA6 accepting or donating a H-bond. pAA is amide accepting or donating a H-bond to another amide. pAWA represents an amide accepting a H-bond from water, whereas pAWD is an amide donating a H-bond to water. pV is the probability of an amide lacking a H-bond acceptor or donor. The error bars are the standard deviation in counts over the given trajectory.

frequently as pAWD in Figure 5, although the dynamical process of hydration remains a question, that is, whether CO or N−H is preferentially hydrated upon water absorption. The number of water molecules participating in a cluster of a given size was counted in each frame. The counts were converted to probabilities of water molecules participating in water clusters of a certain size. These probabilities were averaged over the MD trajectories and given in Figure 6. The water clusters with more than five (5+) water molecules were considered as one group in this process. For 2.3 and 4.5 wt %, the probabilities shift slightly favoring dimers and trimers more. Adding 8.6 wt % water to the membrane causes a rise in the probability of the clusters with 5+ water molecules at the expense of the probability of lone water. The clusters with 5+ water molecules constitute over 65% of the water in PA6 in the case of 15.9 wt% water in the membrane. Puffr and Šebendas’ model also introduces a type of water molecules, denoted free water, which only interacts with other water molecules in the

membrane. The percentages of water bound to the membrane were found using the time-averaged concentrations of water clusters. Approximately 93% of the water molecules are found to be H-bonded to PA6 chains for 2.3, 4.5, and 8.6 wt % water in PA6. The amount changes to 81% for 15.9 wt %, suggesting the rise in concentration of free water as suggested by Puffr and Šebendas’ model. The statistics suggest that the formations of chains of water molecules are important in facilitating faster diffusion of water through the PA6 membrane. These structures are predominantly introduced when there is more than 4.5 wt % in the model membrane. 4.2. Assessment of the Accuracy of Anharmonic Calculations Using N-Methyl-Acetamide (NMA). Harmonic and anharmonic vibrational frequencies of NMA, calculated by the scheme described in Section 2.1, are given in Table 2 along with experimental values.53 The harmonic results display a mean absolute deviation (MAD) of 15 cm−1 for the amide I− III bands and their overtones, with a maximum error of 40 6055


Article


Figure 6. Probability for a water molecule to be in a cluster constituted of X water molecules for the four different levels of hydration used in this study. Error bars are taken as the standard deviation of the water cluster counts in the given trajectory.

solvation experienced by the NMA molecule in a QM/MM box, compared with the inclusion of only water in the first hydration shell in this model. 4.3. IR Difference Spectrum of PA6 in the Amide A Region. The harmonic, anharmonic, and experimental14 difference spectra for the amide A region and the representative cluster types are shown in Figure 7. The harmonic result for this region is, as predicted in Section 4.2, in poor qualitative and quantitative agreement with that of the experiment. The anharmonic spectrum is in better qualitative agreement with that of the experiment, although the fourth peak of the anharmonic spectrum is negative, which is different from the positive valley of the experimental spectrum around 3300 cm−1. Zhou et al.13 have also reported the IR difference spectrum of this region, showing good agreement with the report by Iwamoto and Murase. The sharp peak at 3296.9 cm−1 is only found in Iwamoto and Murase’s study and might be due to the difference in sample processing. Both studies have independently assigned the experimental IR difference spectra. Zhou et al. assigned the difference spectrum to the hydrogen stretches of water molecules inside the PA6 membrane on the basis of the theoretical spectra calculated for water clusters of various sizes using a classical force field and the experimental IR spectrum of liquid water. Iwamoto and Murase assigned the spectral features above 3400 cm−1 to OH stretch of water and the peak at 3248 cm−1 to the amide A mode of amide groups forming H-bonds with water. This is based on the N−H Hbond to water being stronger than that to amide, thereby causing a red shift of the amide A peak. The shoulder of the 3248 cm−1 peak, centered at 3187 cm−1, was assigned to the combination band of amides I and II. As discussed in Section 4.2, we do still expect the theoretical frequencies to be shifted due to limitations in our model. The numbered theoretical peaks are, therefore, assigned as follows: Peaks 2 and 3 are assigned to the experimental peaks at 3465.5 and 3426.1 cm−1, respectively. Peak 4 is assigned to the valley between 3426.1 and 3296.9 cm−1. Peak 5 is assigned to the peak at 3248.1 cm−1. Peaks 5 and 6 are assigned to the experimental peaks at 3296.9 and 3248.1 cm−1, respectively. Lastly, peak 7 is assigned to the experimental peak at 3093.1 cm−1.

Table 2. Calculated Harmonic and Anharmonic Frequencies of the Amide I−III and A Bands in NMAa harmonic amide III fundamental amide II fundamental amide I fundamental amide A fundamental amide III overtone amide II overtone amide I overtone amide II + III combination band amide I + III combination band amide I + II combination band

anharmonic

experiment53

freq.

abs.

freq.

abs.

freq.

1251 1518 1740 3665 2502 3036 3480 2769

79.6 242 309 28.0

1240 1504 1726 3506 2474 3002 3437 2767

77.1 244 302 21.4 0.06 0.76 3.99 0.31

1258 1500 1728 3498 2504

2991

2968

0.40

2971

3258

3233

0.11

3440 2758

a

Harmonic overtones and combination bands are calculated as sums of the fundamental frequencies. The experimental data are from the resonance Raman spectrum of ref 53.

cm−1 for the amide I overtone. The anharmonic result exhibits an MAD of 8.9 cm−1 for the same motions with a maximum error of 30 cm−1 for the amide III overtone. The comparison for these modes slightly favors the anharmonic theory. The error in the amide A band from the harmonic theory is over 100 cm−1, whereas it is only a few cm−1 in the anharmonic case. This demonstrates the need for the anharmonic theory for the accurate description of the stretching motions of hydrogen. The theoretical absorption cross-section of combination bands and overtones are zero in the harmonic theory, whereas the anharmonic theory allows the description of these transitions. In a recent direct dynamics of NMA by Schwörer et al.,54 similar results were reported for NMA in vacuum using the B3LYP functional. However, the results of the solvated system are significantly different from the clusters describing the interactions between water and NMA, that is, the amide I peak of solvated NMA is calculated at 1628 cm−1, whereas the current study finds the amide I bands of clusters 2W−A−1W and 1W−A−1W to be 1664 and 1695 cm−1, respectively. The experimental reference55 for the solvated amide I band is 1625 cm−1. The difference might be explained by the increased 6056


Article


Figure 7. Comparison of the harmonic, anharmonic, and experimental14 IR difference spectra of wet and dry PA6 in the amide A region. The calculated spectra are broadened by a Lorentzian function with a width of 80 cm−1. The thick dashed lines indicate the contribution of cluster types with water as a chromophore, and the thin straight dashed lines denote ΔAbsorbance = 0. The representative structures given in (B)−(D) are the structures found to absorb most and closest to the peak locations found in the calculated spectra.

In the current study, there are 40 clusters contributing to the spectrum. It is, therefore, not possible to give an assignment directly to a certain cluster type. It is, however, possible to discern the general trends in the number of H-bonds and to some extent the molecules participating in the H-bonds with the chromophore. The two OH stretch modes for water are separated by around 100−200 cm−1. However, no direct assignment to symmetric or asymmetric stretches can be made on the basis of the coordinates of the partial Hessian diagonalization. The structures in Figure 7B−D are merely representative for the structures absorbing most and closest to the peak position. These structures are, however, generally representable for the cluster types that are found to absorb in the area with respect to the number of H-bonds and H-bond partners. Figure 7A shows that the region is dominated by water O−H stretch absorption as predicted by Zhou et al.13 The stretches of hydrogen lacking a H-bond partner are generally found at high frequencies around peak 1. As expected, the absorption strength of these peaks is very weak. When a water molecule donates two H-bonds, the peak center of its highest frequency can be found around peaks 2 and 3b, depending on the presence of a H-bond donor to water and the strength of H-

bonds formed. Peak 5b represents the low frequency O−H stretch of water molecules found to absorb around peaks 2 and 3b. On the other hand, peak 6b is the O−H stretch of water forming a strong H-bond by donation, either by accepting/ donating two H-bonds or by accepting/donating one bond. In the latter case, the O−H stretch absorbing around peak 6 is that of the hydrogen donating a H-bond. Iwamoto and Murase14 assigned the shoulder corresponding to peak 6 to a combination of amides I and II. However, such combination bands are found to have small absorption cross sections in any cluster type, and thus are not prominent in the calculated spectrum. Contributions to the difference spectrum due to shifts in the amide A band can still be found, although the region is dominated by water absorption. The absorption cross-section of the amide A bands of cluster types without a H-bond acceptor, which are found around 3500 cm−1, is found to be so small that these cluster types are hardly visible in the calculated difference spectrum. The amide A bands of clusters donating a H-bond to water are found to be blue shifted when compared with similar structures donating a H-bond to an amide group. The central negative peak 4 corresponds to the 1A−A−1A and 6057


Article


Figure 8. Comparison of the harmonic, anharmonic, and experimental14 IR difference spectra of wet and dry PA6 in the amide A region. The calculated spectra are broadened by a Lorentzian function with a width of 40 cm−1. The thick dashed lines indicate the contribution of cluster types with water as a chromophore, and the thin straight lines denote the ΔAbsorbance = 0 line. The representative structures given in (B)−(D) are the structures found to absorb most and closest to the peak locations found in the calculated spectra.

−A−1A cluster types. The change of 1A−A−1A to 1A−A−1W results in a blue shift of the amide A frequency to 3408 cm−1 (see Tables S6−S8). Peak 3a is therefore assigned to amide chromophores donating a H-bond to water. In contrast, peak 5a corresponds to amide clusters with CO accepting two Hbonds. Iwamoto and Murase14 and Zhou et al.13 have found a gradual blue shift around 3090 cm−1 in the IR spectrum with respect to hydration resulting in a positive peak in the IR difference spectrum. Iwamoto and Murase14 suggested that the blue shift was due to the replacement of amide−amide H-bonds with amide−water H-bonds, as the amide II band is expected to blue shift when N−H forms a H-bond to water. Two overtones are observed in this region for the anharmonic spectra, namely, the water bending and amide II overtones. The water bending overtones seem to contribute most in this area because the loss of amide concentration counteracts the rise in absorption due to the amide II overtone. Therefore, peak 7 is assigned to the water bending overtones in this study. 4.4. IR Difference Spectrum of the PA6 Amide I−II Region. The calculated and experimental14 IR difference

spectra for the amide I−II region are shown in Figure 8. The harmonic and anharmonic spectra are both in agreement about the overall shape of the spectrum. However, the shape of the theoretical and experimental spectra is quite different in the regions around the third and sixth peak in Figure 8A. This difference will be discussed in detail in Section 5. The anharmonic spectrum is assigned to amide I−II motions of PA6, except for the diffuse fifth peak, which is assigned to water bending motions. In the following discussion, peaks 3 and 6b will be considered the central peaks of the amide I and II band, respectively, and any calculated blue or red shift is with respect to this center. Iwamoto and Murase14 assigned all bands in this part of the IR difference spectrum to changes in amide I and II absorption profiles. The first bands (1700−1600 cm−1) are assigned to the amide I modes of CO lacking a H-bond donor and with water as a H-bond donor for the negative and positive peak, respectively. In the anharmonic calculation, four peaks (1−4) are found to correspond to the amide I motion. Loss of a Hbond partner leads to a blue shift of the amide I band in the anharmonic model. The exchange of H-bond donor from 6058


Article


Table 3. Assigned Structures for the Amide Absorptions and Their Anharmonic Calculated Band Centers in cm−1 Rounded to the Nearest 5a vibrational mode

H-bonds of amide chromophore

band centers (cm−1)

amide A

free NH N−H···OH2 N−H···OC (2A)b-N−H···OC free CO CO···(H−N/H−OH) CO···[2x(H−N/H-OH) + (H−N,H-OH)] N−H···OC N−H···OH2 free NH

∼3500 3390−3430 3310−3370 3310−3340 1710−1725 1685−1705 1655−1680 1555−1585 1535−1570 1505−1535

amide I

amide II

a

See Tables S6−S8 and Figures S7−S9 for the anharmonic frequencies of each cluster. bThe CO group of the amide chromophore accepts two Hbonds (i.e., the group (2A)-N−H···OC corresponds to the structures 2A−A−1A, 1A1W−A−1A and 2W−A−1A).

Figure 9. (Left) Experimental difference spectra at different times after hydration14 and (right) the anharmonic IR difference spectra of the amide A region of PA6 with different water concentrations in the membrane compared to those of the experimental reference at 0 min after hydration. All anharmonic spectra are calculated with a broadening factor of 80 cm−1.

5. DISCUSSION The model for water absorption in PA6 proposed by Puffr and Šebenda4 has been reproduced in the MD simulation. The water is initially found to form small clusters, whereas higher hydration levels lead to the formation of long water chains and increased amounts of free water in the membrane. The change in the hydration structure brings about a sudden increase in the diffusion coefficient of water, which is consistent with experimental observations.8 The calculation of the IR difference spectra with the harmonic theory did not produce sufficiently accurate results in the amide A region. This is both due to the strong anharmonic nature of the hydrogen stretching modes and the lack of overtone and combination bands in the harmonic theory. Anharmonic calculations suggest that both blue and red shifts of the amide A band can happen upon the hydration of PA6. The former is caused when the H-bond acceptor of the N−H bond is exchanged from amide to water (4 to 3a in Figure 7C), whereas the latter is due to the formation of a second H-bond

amide to water causes a minor shift in the peak center, where the direction depends on the nature of the cluster. CO gaining an additional H-bond red shifts the amide I band. The experimental amide II bands (1600−1500 cm−1) are assigned by Iwamoto and Murase to the N−H groups with water as the H-bond acceptor (positive band) and to N−H groups with no H-bond acceptor (negative band). Considering the 1A−A−1A cluster type band, peak 6b, as the peak center, it is found that replacing the amide−amide H-bond from the N− H group causes a red shift. The only exception is when more groups are found to donate a H-bond to the CO group, in which case the amide II band either remains the same or does a slight blue shift. It should be noted that a blue shift of the amide II band corresponds to a stronger H-bond, whereas a red shift corresponds to a weaker H-bond. The amide II band of clusters donating a H-bond to water is blue shifted compared to the clusters lacking a H-bond acceptor. 6059


Article


Figure 10. (Left) Experimental difference spectra at different times after hydration14 and (right) the anharmonic IR difference spectra of the amide I−II region of PA6 with different water concentrations in the membrane compared to those of the experimental reference at 0 min after hydration. All anharmonic spectra are calculated with a broadening factor of 40 cm−1.

has two H-bond donors (Table 3 and Figure S8). The peak frequencies of the cluster accepting a H-bond from one water or amide are found to overlap, thereby explaining the noisy nature of the shoulder around 1640 cm−1. Compared with the H-bond strength of N−H···X, that of CO···X does not show clear tendency, as shown in Figure S3. Nonetheless, the data do indicate that the second H-bond is weaker than the first one and that the strength of two H-bonds combined together is stronger than that of a single H-bond. In the current calculation, the difference spectrum around the amide II band consists of two single positive and negative patterns. Note that such a reverse pattern of amide A is expected, because a stronger H-bond results in a more restrained N−H bending motion. However, in the experimental spectrum, only one such positive/negative pattern is observed. The reason for the pattern observed in calculation is due to an overlap of the absorption of chromophores donating a H-bond to water and amide. Disentangling the pattern with respect to H-bond acceptors and donors of the central amide chromophore is difficult, because the H-bond partner of the CO group strongly influences the peak center. It may also be due to the mixing of the N−H bending and C−N stretching motion in normal coordinates of amide II (see Figure S5). The use of a higher level of calculations may bring the amide II frequencies of clusters donating a H-bond to amide and those donating a H-bond to water closer, thereby causing peaks 6b and 7 to cancel out and increase the intensity of peak 6a. PA6 is a semicrystalline material under standard conditions. Water can enter only into amorphous parts of the PA6 polymer material.3 Recently, the NMR experiment by Reuvers et al.9 has suggested that water enters only 32% of the PA6 membrane, whereas the rest of the membrane is left dry with up to 25% being crystalline PA6. The present calculation accounts for the amorphous part where the water can freely enter but neglects the effect of any dry or crystalline regions in PA6. Improving the resolution of the IR difference spectra in the vicinity of

to the CO group by water (5a in Figure 7C) (see also Table 3). Iwamoto and Murase14 suggested that the red shift of the amide A band was caused by replacing amide with water as the H-bond acceptor, assuming that the N−H···OC H-bond was stronger than the N−H···OH2 H-bond. However, the present calculation does not support this assumption. In Figure S2, we compare the H-bond strength of N−H···X H-bonds for all clusters considered in this study. It is found that the strength of H-bonds to amide is in a range of 6.5−9.2 kcal mol−1, whereas that to water is in the range of 4.5−6.5 kcal mol−1. This behavior is consistent with the predicted blue shift of the amide A band. Here, we instead suggest that the red shift is caused by H-bond donation of water to the CO group, which strengthens the H-bond between the amide groups. Note that the features of the amide A region are commonly observed for all concentrations of water, as shown in Figure 9. The calculated spectrum in the amide I−II region shown in Figure 8 deviates from the experiment around peaks 3 and 6. In the case of peak 3, the experimental spectrum is found to be very noisy, as noted in ref 14. The experimental difference spectrum at 39 min changes a positive shoulder at 1640 cm−1 to a negative delta peak around the region of the amide I absorption. Furthermore, the experimental spectrum contains a positive delta peak in this shoulder at other times. The reason for this behavior may be due to a strong absorption of the amide I band in the crystalline part of PA6. Unfortunately, this noisy part was hard to resolve when digitalizing the spectrum and therefore we refer to the original figure. The current study used a model containing several different H-bond patterns for the central amide chromophore of PA6, compared with the three-state model (free, H-bond to water, H-bond to amide) used to assign the spectrum by Iwamoto and Murase.14 We find by analyzing the shift patterns of the various H-bond conformations that the shifts can be reduced to a model in which the CO group either does not have a Hbond donor, water is a donor, amide is a donor, or the group 6060



■

ACKNOWLEDGMENTS B.T. is supported by the Special Postdoctoral Researchers Program at RIKEN. This research is partially supported by the “Molecular Systems”, “iTHES”, “Integrated Lipidology”, and “Dynamic Structural Biology” projects in RIKEN (to Y.S.), the Center of innovation Program from Japan Science and Technology Agency, JST, JSPS KAKENHI Grant Nos. JP26220807 and JP26119006 (to Y.S.), and JSPS KAKENHI Grant No. JP16H00857 (to K.Y.). We used computational resources provided by the HPCI System Research Project (Project ID: hp150023), the RIKEN Integrated Cluster of Clusters (RICC), and MEXT SPIRE Supercomputational Life Science (SCLS).

strong IR peaks, as the amide I band of PA6, is a challenge. One way to circumvent this problem is to perform an isotope substitution of the regions where water can enter. This was done in an IR experiment conducted by Murthy et al.5 They measured the isotope-substituted amide A and amide II bands from water-accessible amide groups by soaking the membrane in D2O and then recording the IR spectrum as the membrane dries. This type of measurement may lead to IR difference spectra with higher resolutions that are directly comparable to the calculated spectra for the amorphous regions (Figure 10).

6. CONCLUSIONS MD simulation of the atomistic hydration structure of PA6 predicts the formation of long chains of water molecules. These chains appear to be responsible for the increase in diffusion observed when a PA6 membrane contains more than 5 wt % water.8 The comparison of the anharmonic IR difference spectrum, which is based on the H-bond statistics of the MD simulation, with the experimental spectrum aids in the verification of the MD simulation. The model predicts a red shift of amide A and amid I bands in a series starting from no H-bonds to H-bonds to water, amide, and then the amide chromophores accepting two H-bonds at the CO group. The order is reversed in the case of amide II. The difference spectrum associated with the order of the H-bond strength has revealed the atomistic picture of hydrated PA6. However, the accuracy in the amide I and II region is found to be still insufficient, leading to a discrepancy in the lineshape of the spectrum. Further investigation using a higher level of electronic structure, clustering, and vibrational calculations will be the scope of future work. The method for calculating the IR difference spectrum presented here extends the method suggested in ref 38. The method is trivially extendable to other systems containing many similar chromophores, for example, other polymer and solvent systems. Furthermore, it is possible to extend to other types of vibrational spectroscopy, for example, Raman spectra. With careful consideration of the target system one can obtain a theoretical atomistic account of the changes in a bulk material, which is verifiable by comparison of the theoretical IR difference spectra with experiment.

■

■

ABBREVIATIONS CM, center of mass; DFT, density functional theory; H-bond, hydrogen bond; IR, infrared; MAD, mean absolute deviation; MD, molecular dynamics; MSD, mean squared displacement; NMA, N-methyl-acetamide; oc, optimized coordinate; PA6, polyamide 6; PES, potential energy surface; VMP2, vibrational second order Møller-Plesset; VQDPT2, second-order vibrational quasi-degenerate perturbation theory; VSCF, vibrational self-consistent field; VSE, vibrational Schrödinger equation; SM, sphingomyelin

■

REFERENCES

(1) Lee, K. P.; Arnot, T. C.; Mattia, D. A Review of Reverse Osmosis Membrane Materials for DesalinationDevelopment to Date and Future Potential. J. Membr. Sci. 2011, 370, 1−22. (2) Zhang, H. W.; Shen, P. K. Recent Development of Polymer Electrolyte Membranes for Fuel Cells. Chem. Rev. 2012, 112, 2780− 2832. (3) Murthy, N. S. Hydrogen Bonding, Mobility, and Structural Transitions in Aliphatic Polyamides. J. Polym. Sci., Part B: Polym. Phys. 2006, 44, 1763−1782. (4) Puffr, R.; Šebenda, J. On the Structure and Properties of Polyamides. XXVII. The Mechanism of Water Sorption in Polyamides. J. Polym. Sci., Part C: Polym. Symp. 1967, 16, 79−93. (5) Murthy, N. S.; Stamm, M.; Sibilia, J. P.; Krimm, S. Structural Changes Accompanying Hydration in Nylon 6. Macromolecules 1989, 22, 1261−1267. (6) Murthy, N. S.; Akkapeddi, M. K.; Orts, W. J. Analysis of Lamellar Structure in Semicrystalline Polymers by Studying the Absorption of Water and Ethylene Glycol in Nylons Using Small-Angle Neutron Scattering. Macromolecules 1998, 31, 142−152. (7) Grubb, D. T.; Murthy, N. S. Real-Time X-ray Study of Nylon-6 Fibers During Dehydration: Equatorial Small-Angle Scattering is Due to Surface Refraction. Macromolecules 2010, 43, 1016−1027. (8) Reuvers, N. J. W.; Huinink, H. P.; Fischer, H. R.; Adan, O. C. G. Quantitative Water Uptake Study in Thin Nylon-6 Films with NMR Imaging. Macromolecules 2012, 45, 1937−1945. (9) Reuvers, N. J. W.; Huinink, H. P.; Adan, O. C. G.; Garcia, S. J.; Mol, J. M. C. Water Uptake in Thin Nylon 6 Films as Measured by Electrochemical Impedance Spectroscopy and Magnetic Resonance Imaging. Electrochim. Acta 2013, 94, 219−228. (10) Reuvers, N. J. W.; Huinink, H. P.; Fischer, H. R.; Adan, O. C. G. The Influence of Ions on Water Transport in Nylon 6 Films. Polymer 2013, 54, 5419−5428. (11) Litvinov, V. M. Diffusivity of Water Molecules in Amorphous Phase of Nylon-6 Fibers. Macromolecules 2015, 48, 4748−4753. (12) Reuvers, N. J. W.; Huinink, H. P.; Adan, O. C. G. Plasticization Lags Behind Water Migration in Nylon-6: An NMR Imaging and Relaxation Study. Polymer 2015, 63, 127−133. (13) Zhou, S. M.; Tashiro, K.; Ii, T. Moisture Effect on Structure and Mechanical Property of Nylon 6 as Studied by the Time-Resolved and Simultaneous Measurements of FT-IR and Dynamic Viscoelasticity

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b00372. Movie depicting the diffusion of water for 2.3, 4.5, 8.6, and 15.9 wt % (AVI) Information on the force field used to describe the PA6 polymer, and the absorption cross-sections and frequencies of the cluster types used in this calculation (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +81-48-462-1407. ORCID

Bo Thomsen: 0000-0002-5662-4440 Yuji Sugita: 0000-0001-9738-9216 Kiyoshi Yagi: 0000-0003-1120-9355 Notes

The authors declare no competing financial interest. 6061


Article

The Journal of Physical Chemistry B under the Controlled Humidity at Constant Scanning Rate. Polym. J. 2001, 33, 344−355. (14) Iwamoto, R.; Murase, H. Infrared Spectroscopic Study of the Interactions of Nylon-6 with Water. J. Polym. Sci., Part B: Polym. Phys. 2003, 41, 1722−1729. (15) Dovesi, R.; Orlando, R.; Erba, A.; Zicovich-Wilson, C. M.; Civalleri, B.; Casassa, S.; Maschio, L.; Ferrabone, M.; De La Pierre, M.; D’Arco, P.; et al. CRYSTAL14: A Program for the Ab Initio Investigation of Crystalline Solids. Int. J. Quantum Chem. 2014, 114, 1287−1317. (16) Quarti, C.; Milani, A.; Civalleri, B.; Orlando, R.; Castiglioni, C. Ab Initio Calculation of the Crystalline Structure and IR Spectrum of Polymers: Nylon 6 Polymorphs. J. Phys. Chem. B 2012, 116, 8299− 8311. (17) Milani, A. Unpolarized and Polarized Raman Spectroscopy of Nylon-6 Polymorphs: A Quantum Chemical Approach. J. Phys. Chem. B 2015, 119, 3868−3874. (18) Bowman, J. M. Self-Consistent Field Energies and Wavefunctions for Coupled Oscillators. J. Chem. Phys. 1978, 68, 608. (19) Gerber, R. B.; Ratner, M. A. A Semiclassical Self-Consistent Field (SC SCF) Approximation for Eigenvalues of Coupled-Vibration Systems. Chem. Phys. Lett. 1979, 68, 195−198. (20) Christoffel, K. M.; Bowman, J. M. Investigations of SelfConsistent Field, SCF CI and Virtual Stateconfiguration Interaction Vibrational Energies for a Model Three-Mode System. Chem. Phys. Lett. 1982, 85, 220−224. (21) Norris, L. S.; Ratner, M. A.; Roitberg, A. E.; Gerber, R. B. Møller−Plesset Perturbation Theory Applied to Vibrational Problems. J. Chem. Phys. 1996, 105, 11261. (22) Beck, M. H.; Jäckle, A.; Worth, G. A.; Meyer, H.-D. The Multiconfiguration Time-Dependent Hartree (MCTDH) Method: a Highly Efficient Algorithm for Propagating Wavepackets. Phys. Rep. 2000, 324, 1−105. (23) Christiansen, O. Vibrational Coupled Cluster Theory. J. Chem. Phys. 2004, 120, 2149. (24) Yagi, K.; Hirata, S.; Hirao, K. Vibrational Quasi-Degenerate Perturbation Theory: Applications to Fermi Resonance in CO2, H2CO, and C6H6. Phys. Chem. Chem. Phys. 2008, 10, 1781−1788. (25) Bowman, J. M.; Wang, Y.; Liu, H.; Mancini, J. S. Ab Initio Quantum Approaches to the IR Spectroscopy of Water and Hydrates. J. Phys. Chem. Lett. 2015, 6, 366−73. (26) Liu, H.; Wang, Y.; Bowman, J. M. Quantum Calculations of the IR Spectrum of Liquid Water Using Ab Initio and Model Potential and Dipole Moment Surfaces and Comparison with Experiment. J. Chem. Phys. 2015, 142, No. 194502. (27) Liu, H. C.; Wang, Y.; Bowman, J. M. Quantum Local Monomer IR Spectrum of Liquid D2O at 300 K from 0 to 4000 cm‑1 Is in NearQuantitative Agreement with Experiment. J. Phys. Chem. B 2016, 120, 2824−2828. (28) Liu, H.; Wang, Y.; Bowman, J. M. Transferable Ab Initio Dipole Moment for Water: Three Applications to Bulk Water. J. Phys. Chem. B 2016, 120, 1735−1742. (29) Homayoon, Z.; Conte, R.; Qu, C.; Bowman, J. M. FullDimensional, High-Level Ab Initio Potential Energy Surfaces for H2(H2O) and H2(H2O)2 With Application to Hydrogen Clathrate Hydrates. J. Chem. Phys. 2015, 143, No. 084302. (30) Qu, C.; Bowman, J. M. Ab Initio, Embedded Local-Monomer Calculations of Methane Vibrational Energies in Clathrate Hydrates. J. Phys. Chem. C 2016, 120, 3167−3175. (31) Wang, Y.; Bowman, J. M. Towards an Ab Initio Flexible Potential for Water, and Post-Harmonic Quantum Vibrational Analysis of Water Clusters. Chem. Phys. Lett. 2010, 491, 1−10. (32) Head, J. D. Computation of Vibrational Frequencies for Adsorbates on Surfaces. Int. J. Quantum Chem. 1997, 65, 827−838. (33) Hanson-Heine, M. W.; George, M. W.; Besley, N. A. Rapid Anharmonic Vibrational Corrections Derived From Partial Hessian Analysis. J. Chem. Phys. 2012, 136, No. 224102.

(34) Yagi, K.; Keçeli, M.; Hirata, S. Optimized Coordinates for Anharmonic Vibrational Structure Theories. J. Chem. Phys. 2012, 137, No. 204118. (35) Thomsen, B.; Yagi, K.; Christiansen, O. Optimized Coordinates in Vibrational Coupled Cluster Calculations. J. Chem. Phys. 2014, 140, No. 154102. (36) Yagi, K.; Otaki, H. Vibrational Quasi-Degenerate Perturbation Theory with Optimized Coordinates: Applications to Ethylene and Trans-1,3-Butadiene. J. Chem. Phys. 2014, 140, No. 084113. (37) König, C.; Hansen, M. B.; Godtliebsen, I. H.; Christiansen, O. FALCON: A Method for Flexible Adaptation of Local Coordinates of Nuclei. J. Chem. Phys. 2016, 144, No. 074108. (38) Yagi, K.; Li, P.-C.; Shirota, K.; Kobayashi, T.; Sugita, Y. A Weight Averaged Approach for Predicting Amide Vibrational Bands of a Sphingomyelin Bilayer. Phys. Chem. Chem. Phys. 2015, 17, 29113− 29123. (39) Shirota, K.; Yagi, K.; Inaba, T.; Li, P.-C.; Murata, M.; Sugita, Y.; Kobayashi, T. Detection of Sphingomyelin Clusters by Raman Spectroscopy. Biophys. J. 2016, 111, 999. (40) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (41) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D., Jr.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Llipids: Validation on Six Lipid Types. J. Phys. Chem. B 2010, 114, 7830−7843. (42) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586−3616. (43) Venable, R. M.; Sodt, A. J.; Rogaski, B.; Rui, H.; Hatcher, E.; MacKerell, A. D., Jr.; Pastor, R. W.; Klauda, J. B. CHARMM All-Atom Additive Force Field for Sphingomyelin: Elucidation of Hydrogen Bonding and of Positive Curvature. Biophys. J. 2014, 107, 134−145. (44) Price, D. J.; Brooks, C. L., III A Modified TIP3P Water Potential for Simulation with Ewald Summation. J. Chem. Phys. 2004, 121, 10096. (45) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089. (46) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577. (47) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511. (48) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (49) Feller, S. E.; Zhang, Y. H.; Pastor, R. W.; Brooks, B. R. Constant Pressure Molecular Dynamics Simulation: The Langevin Piston Method. J. Chem. Phys. 1995, 103, 4613. (50) Desiraju, G. R.; Steiner, T. The Weak Hydrogen Bond: In Structural Chemistry and Biology; Oxford University Press, 2001. (51) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (52) Rohatgi, A. WebPlotDigitizer, version 3.9; Zenodo, 2015. (53) Mayne, L. C.; Hudson, B. Resonance Raman Spectroscopy of NMethylacetamide: Overtones and Combinations of the CarbonNitrogen Stretch (Amide II′) and Effect of Solvation on the Carbon-Oxygen Double-Bond Stretch (Amide I) Intensity. J. Phys. Chem. 1991, 95, 2962. (54) Schwörer, M.; Wichmann, C.; Tavan, P. A polarizable QM/MM approach to the molecular dynamics of amide groups solvated in water. J. Chem. Phys. 2016, 144, No. 114504. (55) Kubelka, J.; Keiderling, T. A. Ab initio calculation of amide carbonyl stretch vibrational frequencies in solution with modified basis 6062


Article

The Journal of Physical Chemistry B sets. 1. N-methyl acetamide. J. Phys. Chem. A 2001, 105, 10922− 10928. (56) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38.

6063


Weight-Averaged Anharmonic Vibrational Analysis of Hydration

Recommend Documents