Phenolic Polymer Solvation in Water and Ethylene Glycol, II: Ab Initio

Article pubs.acs.org/JPCB

Phenolic Polymer Solvation in Water and Ethylene Glycol, II: Ab Initio Computations Charles W. Bauschlicher, Jr.,*,† Eric W. Bucholz,† Justin B. Haskins,‡ Joshua D. Monk,‡ and John W. Lawson† †

Thermal Protection Materials Branch and ‡AMA, Inc., Thermal Protection Materials Branch, NASA Ames Research Center, Moffett Field, California 94035, United States ABSTRACT: Ab initio techniques are used to study the interaction of ethylene glycol and water with a phenolic polymer. The water bonds more strongly with the phenolic OH than with the ring. The phenolic OH groups can form hydrogen bonds between themselves. For more than one water molecule, there is a competition between water−water and water−phenolic interactions. Ethylene glycol shows the same effects as those of water, but the potential energy surface is further complicated by CH2−phenolic interactions, different conformers of ethylene glycol, and two OH groups on each molecule. Thus, the ethylene glycol−phenolic potential is more complicated than the water−phenolic potential. The results of the ab initio calculations are compared to those obtained using a force field. These calibration studies show that the water system is easier to describe than the ethylene glycol system. The calibration studies confirm the reliability of force fields used in our companion molecular dynamics study of a phenolic polymer in water and ethylene solutions.

I. INTRODUCTION Phenolic resins are used in a wide variety of applications requiring lightweight materials with thermal stability, strength, and fire retardancy. These resins are formed by curing lowmolecular-weight oligomers after dissolution in a suitable solvent. The original oligomers have a distribution of polymer lengths and are classified as either novalac or resole depending on the number of formaldehydes per phenol unit. Although many factors affect the quality of the final product, the choice of solvent is one of the most critical. Understanding solvent interactions with the phenolic polymer precursors will help guide the selection, and potentially the design, of solvents for improved resin curing. This work is Part II of a two-part computational study on the fundamentals of phenolic polymer−solvent interactions. The two considered solvents, ethylene glycol (EGL) and water (H2O), are commonly used for phenolic processing and are known to provide differing postcure resin morphologies. Water is the simplest example of a polar solvent, and ethylene glycol is the next most complicated with two sites for hydrogen bonding. In Part I, we performed1 molecular dynamics (MD) studies to evaluate a range of properties important for polymer solvation. Importantly, these properties were found to depend critically on subtle, nonbonded, and hydrogen-bonded intermolecular interactions. One open question for such simulations is whether the utilized interatomic force field adequately captured the complex, controlling interactions for such systems. In this article, we perform detailed ab initio calculations based on density functional theory (DFT) as well as higher levels of theory on simple model systems relevant for phenolic−solvent interactions. In particular, we report on results of studies on the © 2017 American Chemical Society

fundamental interactions between the solvents, water and ethylene glycol, interacting with models for the phenolic polymer. In addition to obtaining insight into the solvent− phenolic interactions, these studies are used to validate the force field utilized in the MD simulations performed in Part I.

II. METHODS Phenolic polymers are differentiated on the basis of the number of CH2OH side groups, with those with few side groups being denoted novalac, whereas those with a large number of side groups being denoted resole. The novalac and resole forms are differentiated because they are cured differently. However, the final cured resins are essentially identical. Our three models for the phenolic polymer are shown in Figure 1. The top two molecules are used to study the novalac form, whereas the bottom molecule is used to study the resole form. We use “N” and “R” to denote the novalac and resole forms, respectively. The number of C6 segments in the model is given after the letter denoting the form; thus, N2 signifies the two-segment novalac model. As one goal of this study is to determine the accuracy of our force field and of DFT methods to determine the interactions between water and ethylene glycol interacting with a phenolic polymer, we need to establish accurate results to be used for comparison. For the conformations of ethylene glycol and its dimer, the calculations are sufficiently small that the coupled cluster singles and doubles approach,2 including the effect of Received: January 11, 2017 Revised: March 13, 2017 Published: March 14, 2017 2852

DOI: 10.1021/acs.jpcb.7b00327 J. Phys. Chem. B 2017, 121, 2852−2863

Article

The Journal of Physical Chemistry B

and co-workers. We also add the empirical D2, D3, or D3(BJ) dispersion terms of Grimme.16−18 As noted above, electron correlation was also included by the MP2 and CCSD(T) approaches. The DFT calculations are performed using the 6-31+G* basis set of Pople and co-workers19 for oxygen and the 6-31G** basis set for the remaining atoms, which we denote 6-31(+)G**. The augmented correlation-consistent polarized triple zeta (aug-ccpVTZ) basis set of Dunning and co-workers20,21 is also used in the DFT calculations. The MP2 calculations are performed using the aug-cc-pVTZ and the quadruple and quintuple zeta analogues, aug-cc-pVQZ and aug-cc-pV5Z, respectively. The MP2 calculations are extrapolated to the complete basis set limit (CBS) using the X−3 approach of Helgaker et al.22 The basis set superposition error (BSSE) is computed using the counterpoise correction approach of Boys and Bernardi.23 Finally, we note that the aug-cc-pV double zeta (DZ) basis set is used in some MP2 calculations when the conformations of ethylene glycol were studied. Unless otherwise noted, all of the geometric parameters have been fully optimized. For the DFT techniques, the harmonic frequencies were computed to ensure that the structures correspond to minima. The reported binding energies do not include zero-point energy as we are unable to easily obtain vibrational frequencies for the force field results. The populations of the 10 conformers of ethylene glycol in the gas phase are determined using the computed relative free energy (ΔG) values at the DFT level. The populations for the liquid phase are determined using quantum molecular dynamics (QMD). These are Γ-point only calculations with an energy cutoff of 400 eV and a time step of 0.5 fs. The liquid phase simulations contained 50 ethylene glycol molecules and were carried out for 55 ps. The QMD approach was also applied to a single ethylene glycol molecule with a 70 ps simulation time. This gas phase QMD simulation was found to be in good agreement with the DFT ΔG approach. The Gaussian 09 program24 was used for all of the calculations except those involving the force field, which were preformed with a modified version of LAMMPS,25 and the QMD calculations that were performed using VASP.26−28

Figure 1. Molecules used as models of a phenolic polymer. The top two figures are the two- and three-segment models of the novalac polymer, whereas the bottom figure is the two-segment model used to study the resole polymer.

connected triples determined using perturbation theory,3 CCSD(T), can be used. This level of theory has been previously4 applied to the water and water dimer. However, for the models of a phenolic polymer, such calculations would be prohibitively expensive and therefore lower levels of theory are required. For the water dimer, previous calculations5 have shown that the basis set limit binding energy using second-order Møller−Plesset perturbation theory (MP2) is 5.0 ± 0.1 kcal/mol, which is in excellent agreement with the best estimate4 (5.0 kcal/mol) obtained from coupled cluster calculations including single, double, and triple excitations and a perturbative estimate of the quadruples. Therefore, we use the MP2 approach for as many of our calibration calculations as possible, but even this level becomes difficult for the largest systems of interest. Therefore, we attempt to establish which DFT method agrees best with our MP2 calculations and use it to test our force field on the largest systems. Our force field comes from several different sources. For the phenolic chains, the bonded and nonbonded interactions were computed using the all-atom optimized potentials for liquid simulation (OPLS-AA) force field developed by Jorgensen et al.6 The bond, angle, and dihedral parameters that were not provided by Jorgensen et al. were taken from Cornell et al.7 or from the GROMACS MD package.8 For ethylene glycol, the bonded and nonbonded interactions were taken from the OPLS-AA scaling electrostatic interaction (OPLS-AA-SEI) force field,9 whereas for H2O, the TIP3P/ew water model10 was used. For simplicity, we refer to the potentials as OPLS even though this represents only one component of the potential. DFT calculations were performed using a variety of functionals. These include the following: the B3LYP11 hybrid12 functional; the exchange and correlation function of Perdew, Burke, and Ernzerhof13 (PBE, which is denoted PBEPBE in Gaussian); and the M0614 and M06L15 functionals of Truhlar

III. RESULTS AND DISCUSSION For the larger systems, the number of conformations becomes very large and it would be extremely difficult to determine all possible structures. However, our goal is not to map out entire potential energy surfaces, but rather sample several different structures to understand the importance of the different types of bonds; for example, to compare structures with different number of phenolic−phenolic, phenolic−solvent, and solvent−solvent hydrogen bonds and to determine if the different methods can describe them with the same accuracy. In fact, we do not report all of the conformers that we found, but rather report a sampling of those showing different types of interactions. III.I. Water Dimer and Ethylene Dimer. The water and ethylene glycol dimers are shown in Figure 2. The water dimer has one hydrogen bond, whereas the ethylene glycol dimer has four H−O interactions, but because of the small HOH angle, none are classical hydrogen bonds. The binding energies are summarized in Table 1. The highest level of theory used to optimize the geometry is the MP2/aug-cc-pVTZ approach, and then, we extrapolate the MP2 binding energy to the CBS limit by performing calculations using the aug-cc-pVQZ and aug-ccpV5Z basis sets at the MP2/aug-cc-pVTZ geometry. The CBS 2853


Article


reliable. The CCSD(T)/aug-cc-pVTZ values are very similar to the MP/aug-cc-pVTZ values, showing that higher levels of correlation have only a small effect, as expected for both the H2O and ethylene glycol dimers. Using the MP2/CBS(4,5) values for comparison, we are able to establish the reliability of the more approximate methods. First, we consider the water dimer results. The DFT binding energies in the 6-31(+)G** basis set are too large. The larger aug-cc-pVTZ basis set reduces the binding energies because it significantly reduces the BSSE. For this basis set, the PBE functional has the smallest error but the errors for the B3LYP, M06, and M06L approaches are also quite small. For the larger basis set, we also try including the Grimme dispersion correction for the B3LYP, M06, and PBE functionals. The addition of the empirical correction leads to too large a binding energy. Overall, the best approach is B3LYP+D3(BJ), where the error is 0.23 kcal/mol, which is very similar to that obtained by the four best functionals (PBE, B3LYP, M06, and M06L) without any empirical corrections. Optimization of the geometry in the large basis set can be quite time consuming; therefore, we also performed the B3LYP/augcc-pVTZ and M06/aug-cc-pVTZ calculations using the B3LYP/ 6-31(+)G** and M06/6-31(+)G** geometries, respectively, and we find that the binding energies differ by only 0.01 kcal/mol from the analogous calculations using the aug-cc-pVTZ

Figure 2. Water and ethylene glycol dimers. For water, the hydrogen bond is shown. For ethylene glycol, the O−H distance is only slightly longer than that for water, 1.96 vs 1.92 Å, but the HOH angles are less than 90°; therefore, they are not consistent with classical hydrogen bonds, and the O−H interactions are noted with a smaller line.

values are reported using the aug-cc-pVTZ and aug-cc-pVQZ basis sets, which is denoted CBS(3,4), and the aug-cc-pVQZ and aug-cc-pV5Z basis sets, which is denoted CBS(4,5). These two extrapolated values are in good agreement with each other for both dimers. CBS(4,5) should be more reliable as it uses the two largest basis sets. For water, we note that the MP2/aug-cc-pV5Z value is very similar to the CBS value, suggesting that the best value is very reliable. Our best value is in excellent agreement with that reported by Feyereisen et al.5 Using the MP2/aug-ccpVTZ unscaled frequencies for a zero-point correction yields a best estimate of 2.87 kcal/mol for D0, which is in good agreement with the experimental value29 of 3.15 ± 0.03 kcal/mol. For ethylene glycol, there is a larger extrapolation effect, but it is sufficiently small that the MP2/CBS(4,5) value should be

Table 1. Summary of Water and Ethylene Glycol Dimer Binding Energies (in kcal/mol) as a Function of Level of Theorya b

c

binding energyd water

method

b

B3LYP/6-31(+)G** PBE/6-31(+)G** M06/6-31(+)G** B3LYP/aug-cc-pVTZ B3LYP/aug-cc-pVTZ B3LYP+D2/aug-cc-pVTZ B3LYP+D3/aug-cc-pVTZ B3LYP+D3(BJ)/aug-cc-pVTZ M06/aug-cc-pVTZ M06/aug-cc-pVTZ M06+D3/aug-cc-pVTZ M06L/aug-cc-pVTZ PBE/aug-cc-pVTZ PBE+D2/aug-cc-pVTZ PBE+D3/aug-cc-pVTZ PBE+D3(BJ)/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ CCSD(T)/aug-cc-pVTZ MP2/aug-cc-pVQZ MP2/aug-cc-pV5Z MP2/CBS(3,4) MP2/CBS(4,5) OPLS

geometry

c

opt opt opt opt B3LYP/6-31(+)G** opt opt opt opt M06/6-31(+)G** opt opt opt opt opt opt B3LYP/6-31(+)G** PBE/6-31(+)G** M06/6-31(+)G** M06/aug-cc-pVTZ opt MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ MP2/aug-cc-pVTZ opt

ethylene glycol

binding

error

binding

error

6.04 6.69 6.10 4.57 4.56 5.41 5.30 5.20 4.79 4.78 5.41 4.75 5.10 5.74 5.44 5.52 5.18 5.12 5.19(4.71) 5.18(4.70) 5.18 5.22 5.09 5.03 5.03 4.97 7.18e

1.06 1.72 1.12 −0.40 −0.41 0.44 0.33 0.23 −0.18 −0.19 0.44 −0.23 0.12 0.77 0.57 0.55 0.21 0.15 0.22(−0.27) 0.21(−0.27) 0.21 0.25 0.12 0.06 0.06 0.00 2.20

10.89 14.05 15.07 9.24 9.16 16.85 14.84 14.37 13.87 13.78 14.88 14.23 12.14 17.90 20.37 15.60 12.34 14.56 14.57(12.43) 14.73(12.59) 14.88 14.96 14.13 13.80 13.59 13.45 12.14

−2.55 0.60 1.62 −4.21 −4.29 3.41 1.40 0.93 0.43 0.34 1.43 0.78 −1.31 4.45 6.92 2.16 −1.11 1.12 1.12(−1.02) 1.28(−0.85) 1.43 1.51 0.68 0.35 0.14 0.0 −1.31

a

The error is the difference between a given level of theory and our best value (MP2/CBS(4,5)). Zero-point energy is not included. bThe method used to determine the binding energy. cThe method used to determine the geometry. “opt” indicates that the geometry is optimized at the same level of theory as the energy is evaluated. dThe values in parentheses are corrected for BSSE. eThe binding energy is 6.86 kcal/mol if the H2O geometry is fixed at that of the monomer. 2854


Article


suggest that the ethylene glycol dimer is more difficult to describe at the DFT level than the water dimer. The aug-cc-pVQZ and aug-cc-pV5Z calculations are very time consuming, and using these basis sets for the phenolic calculation is not currently practical. Therefore, a CBS extrapolation for the phenolic systems is also not practical. One alternative might be correcting the MP2/aug-cc-pVTZ results for the BSSE, and we report the BSSE corrected values in Table 1. The BSSE correction is too large, yielding binding energies that are too small. However, as we have observed on other occasions, the BSSE corrected value is too small by about as much as the uncorrected value is too large. Thus, if only the MP2/aug-ccpVTZ calculation is possible, using half the BSSE correction should yield the best binding energy. III.II. Phenolic−Water Interactions. We studied up to three waters interacting with our models for the phenolic polymer. The optimized structures are shown in Figures 3 and 4, and the

geometry. This very small energy difference supports using the smaller basis set for the optimization and represents a large saving in computational effort, with only a small loss of accuracy. Several MP2/aug-cc-pVTZ results are given, differing in the choice of the geometry. The B3LYP/6-31(+)G** geometry and M06 geometry from either basis set yield MP2 binding energies that are in excellent agreement with those obtained using the MP2-optimized geometry. Using the PBE/6-31(+)G** geometry yields a slightly smaller MP2 binding energy. As the DFT geometries yield MP2 binding energies that are similar to those obtained using the MP2-optimized geometry, their use represents a big saving in computer time with only a small loss in accuracy. The MP2 binding energy decreases as the basis set is improved because of the reduction in BSSE. Finally, we note that the error in the OPLS (TIP3P/ew) binding energy is 2.20 kcal/ mol, which is larger than that for virtually all of the ab initio calculations. Using the original TIP3P model would reduce the error in the dimer to about 1.5 kcal/mol but would yield a poorer description of bulk water. The ethylene glycol dimer results show some differences with those obtained for the water dimer. First, we note that the B3LYP/6-31(+)G** value is too small, and improving the basis set further reduces the binding energy. The PBE/6-31(+)G** value is 0.60 kcal/mol too large, whereas the PBE/aug-cc-pVTZ value is −1.31 kcal/mol too small. M06/aug-cc-pVTZ is the most accurate DFT approach, with an error of only 0.43 kcal/mol. For the water dimer, M06/aug-cc-pVTZ had a small error, being only slightly larger than the PBE/aug-cc-pVTZ result. Using the 631(+)G** geometries yields B3LYP/aug-cc-pVTZ and M06/ aug-cc-pVTZ binding energies that are very similar to those obtained when the geometry is optimized in the larger basis set. Adding the empirical dispersion yields binding energies that are too large. For B3LYP, all three of the dispersion corrections improve the binding energy, with D3(BJ) having the smallest error. For both the M06 and PBE functionals, the addition of the dispersion increases the disagreement with our best value. Thus, as found for the water dimer, the addition of the empirical dispersion does not seem to offer any improvement over using the M06 functional in the large basis set. All of the MP2/aug-cc-pVTZ binding energies using the DFT geometries differ somewhat from the MP2/aug-cc-pVTZ optimized results. This difference varies from 2.54 kcal/mol, when the B3LYP/6-31(+)G** geometries are used, to 0.15 kcal/ mol, when the M06/aug-cc-pVTZ geometries are used. We also note that the difference is about 0.3 kcal/mol for the PBE/631(+)G** and M06/6-31(+)G** geometries. As expected, the geometric differences mostly arise for the dimer, with the monomer geometry being quite similar for all DFT methods. The M06/aug-cc-pVTZ approach gives the best DFT description of the ethylene glycol dimer. The second best binding energy is at PBE/6-31(+)G**, which is clearly benefiting from some cancellation of errors, as improving the basis set significantly reduces the binding energy. For the water dimer, PBE/aug-cc-pVTZ had the smallest error, but M06/aug-ccpVTZ and B3LYP/aug-cc-pVTZ had only marginally larger errors. Thus, there is a larger variation in the DFT binding energies for ethylene glycol than that for water. Given the larger variation in the DFT binding energies, it is probably not too surprising to find a larger variation in the MP2 binding energies depending on the DFT approach used to obtain the geometry. The larger spread in the accuracy of the DFT binding energies and larger effect of DFT geometries on the MP2 binding energies

Figure 3. Water bonded to the two-segment phenolic model.

energetics is summarized in Table 2. For the two-segment model, it is possible to perform MP2 calculations, including optimization of the geometry at the MP2 level. Although it is possible to optimize the geometry, the extrapolation to the CBS limit would be very expensive; therefore, we determine our best value as the MP2 binding energy corrected with half of the BSSE correction, 2855


Article


which, for the dimers, was in very good agreement with our CBS value. We found three local minima for one water bonded to the N2 polymer. N2H2O A has the water hydrogen bonded to the phenolic OH group, whereas N2H2O B and N2H2O C have the water interacting with the ring part of the polymer. For N2H2O B, the water is inside the “V” created by the two rings, whereas N2H2O C is on the outside of the V. For N2H2O B, we were able to find two minima for the M06/6-31(+)G** and PBE/631(+)G** approaches. We use the more stable of the two as this has a structure more similar to the MP2- and B3LYP-optimized geometries. Starting MP2 or B3LYP from the M06 geometry of the less stable minima did not result in a second minima, but rather these levels of theory collapsed back to the more stable structure. The binding energies for B and C are significantly smaller than those for A. We should note that these small binding energies made it difficult to find the water−ring minima, as many attempts to locate them resulted in finding the water−OH minima. The N2(H2O)2 molecule has two waters bonded to the polymer and a hydrogen bond between the two waters. N2(H2O)3 A has a chain of three water molecules hydrogenbonded to each other and the end waters hydrogen-bonded to the polymer. N2(H2O)3 B has a two-water system like that in N2(H2O)2, with the third water hydrogen-bonded to the polymer. For the resole model, we consider two structures for a single water, one with water bonded to the polymer through the water O (R2H2O A) and the other where it is bonded through the H (R2H2O B). For the three-segment novalac model, we have two single-water cases, which differ in the orientation of the water. We have three models for the two-water case: (1) a single

Figure 4. Water bonded to the three-segment phenolic model.

Table 2. Summary of Water−Phenolic Binding Energy (in kcal/mol) as a Function of Level of Theorya reaction

method OPLS

geometry

OPLS

H2O + N2 → N2H2O A H2O + N2 → N2H2O B H2O + N2 → N2H2O C H2O + N2H2O → N2(H2O)2 2H2O + N2 → N2(H2O)2 H2O + N2(H2O)2 → N2(H2O)3 A H2O + N2(H2O)2 → N2(H2O)3 B 3H2O + N2 → N2(H2O)3 A 3H2O + N2 → N2(H2O)3 B H2O + N3 → N3H2O A H2O + N3 → N3H2O B H2O + N3H2O → N3(H2O)2 A H2O + N3H2O → N3(H2O)2 B H2O + N3H2O → N3(H2O)2 C 2H2O + N3 → N3(H2O)2 A 2H2O + N3 → N3(H2O)2 B 2H2O + N3 → N3(H2O)2 C H2O + R2 → R2H2O A H2O + R2 → R2H2O B

9.21 4.51

B3LYP

B3LYP/TZ B3LYP

PBE

M06

PBE

M06/TZ

MP2-BSSE

BESTb

B3LYP

PBE

M06

MP2

M06

MP2

M06

11.54 7.07 4.91 14.90 26.44 11.45

9.49 5.08 3.77 13.09 22.58 9.89

9.12 6.24 4.62 15.49 24.61 8.80

9.12 6.44 4.71 15.23 24.36

10.61 7.20 5.05 14.32 24.93 10.76

10.71 7.28 5.28 14.43 25.14

9.49 5.69 4.12 12.89 22.37 9.55

9.58 5.69 4.24 12.89 22.47

10.05(10.15) 6.45(6.49) 4.58(4.76) 13.60(13.66) 23.65(23.81) 10.15

3.12

5.40

4.17

5.43

5.59

4.49

5.04

34.58 28.87 16.00 15.17 4.51

28.77 23.33 13.43 12.85 3.07

37.89 31.85 16.62 16.21 5.44

32.46 26.75 14.69 14.43 4.16

33.42 30.05

35.69 30.51

31.92 26.87

33.81 28.69

7.22

6.61

5.48

6.64

5.66

11.90

10.67

9.07

10.51

8.86

19.31 21.34 26.02 6.98 8.04

20.51 22.60 26.66 6.19 7.36

16.50 18.91 22.50 4.35 5.65

22.06 23.25 27.13 8.19 8.66

18.85 20.35 23.55 6.65 6.96

14.36 23.56 11.53

8.69 4.49 3.76 15.62 24.31 10.27

7.31 1.19 2.08 12.89 20.20 8.57

5.35

4.56

35.10 28.91 14.12 14.12 5.20

c

9.81 5.76 4.64 17.71 27.51

M06

MP2/TZ

a The DFT calculations use the 6-31(+)G** basis set, unless otherwise noted. TZ denotes the aug-cc-pVTZ basis set. Zero-point energy is not included. bMP2 with half the BSSE correction added at the M06 geometry. The values in parentheses use the MP2-optimized geometry. cCollapses to the solution where the water is interacting with the phenolic OH.

2856


Article

The Journal of Physical Chemistry B water bridges two polymer OH groups, with the second water bonded to the central polymer segment, N3(H2O)2 A, (2) a single water bridges two polymer OH groups with the second water bonded to the first water molecule, N3(H2O)2 B, and (3) where there is a two-water bridge between two polymer OH groups, N3(H2O)2 C. Thus, our models show a variety in the bonding and therefore should serve as a good test of the accuracy of different levels of theory. We first consider the N2 model with one and two water molecules where we can perform MP2 geometry optimization. For these species, we find that the MP2 binding energy using the M06 geometry is in best agreement with the MP2-optimized result. The B3LYP and PBE geometries yield somewhat larger errors. Given the similar MP2 binding energies using the M06 and MP2 geometries, it is not surprising that these two geometries yield similar BSSE corrections, and hence best estimates that agree to about 0.1 kcal/mol per H2O−N2 bond. On this basis, we believe that MP2 at the M06 geometry corrected for half the BSSE yields our best estimate for larger systems where an MP2 optimization would be very expensive and/or the CBS extrapolation is not possible. Comparing lower levels of theory to our best estimate shows that for the N2(H2O)n species we find it difficult to differentiate between the B3LYP, PBE, and M06 functionals in the 631(+)G** basis set. On the basis of one water, PBE would be the functional of choice, but B3LYP performs much better for the two-water case. The B3LYP and PBE binding energies can be either larger or smaller than the best estimate, whereas the M06 results are consistently larger. Expanding the basis reduces the binding energies, as expected. Now, both B3LYP/aug-cc-pVTZ and M06/aug-cc-pVTZ binding energies are smaller than the best estimate, with M06/aug-cc-pVTZ being about as reliable as the B3LYP/6-31(+)G** approach. For the three-segment cases, we are unable to perform the MP2 calculations, but on the basis of the two-segment results, we will use the M06/aug-cc-pVTZ result for comparison. The OPLS results agree well for both first and total binding energies for the N2(H2O)2, N2(H2O)3, and N3(H2O)2 species and correctly find the three-water chain (N2(H2O)2 A) to be more than 5 kcal/mol more stable than the two-water chain plus one additional water (N2(H2O)2 B). It also correctly orders the three N3(H2O)2 species, finding the twowater bridge (C) to be the most favorable. The OPLS also agrees with B3LYP and M06 for the H versus O bonding for the resole model. Overall, the OPLS results are in good agreement with our best values, MP2 when practical and M06/aug-cc-pVTZ in other cases. The only significant difference is for N2H2O C, where the OPLS does not find the weak minima; however, we caution that this difference might be a result of just missing a shallow minimum. III.III. Phenolic−Ethylene Glycol Interactions. The bonding of the ethylene glycol to the polymer is much more complicated than for the water. This arises because the ethylene glycol has two OH groups and two CH2 groups that can interact with the polymer. In addition, the ethylene glycol can have several conformations, which further increases the number of possible interactions. The binding energies for different conformations are discussed below. These multiple bonding options mean that there are numerous local minima on the potential surface, and these minima can differ from method to method. We illustrate this for the resole R2 model in Figure 5. The first optimization used the M06 functional, and this gave the M06 A structure shown in the figure. Starting from this structure, the molecule was reoptimized at the B3LYP level using the

Figure 5. Local minima for the R2eg B species for different levels of theory. The relative energetics is given for each level of theory in kcal/ mol.

program defaults and (1) the starting second derivatives from M06 and (2) with the B3LYP second derivatives computed every step; these optimizations yielded the B3LYP A and B3LYP B structures, respectively. Three OPLS optimizations were performed, starting from M06 A, B3LYP A, and B3LYP B structures, which gave OPLS A, B, and C geometries, respectively. Finally, M06 was reoptimized starting from B3LYP A, which gave the M06 B structure. The small differences in the relative energy show how flat the potential energy surfaces are. It is also interesting to note changes in geometry. OPLS A has the ethylene glycol on the “outside” of the “V” created by the two C6 units even though it was started from M06 A that has the ethylene glycol inside the V, whereas OLPS-AA B has the ethylene glycol inside the V, even though it was started from B3LYP A, which has the ethylene glycol above the phenolic CH2 bridging group. Clearly, the potential is very flat with numerous minima, which make it difficult to compare the different levels of theory. To make the best comparison possible, we exchanged optimized structures between the different methods and reoptimized to find the best local minima. This complexity is not restricted to the resole model. In Figure 6, we show 10 local minima that we found for one ethylene glycol with a two-segment N2 model for the novalac form. There are probably many others that we did not find. We should note that not all minima are found at all levels of theory. This can arise because the minima do not exist at some levels of theory or that they are sufficiently shallow that we did not obtain them in a 2857


Article


Figure 6. Ethylene glycol bonded to the two-segment phenolic models.

addition, we consider fewer structures for the three-segment model and a maximum of two ethylene glycol molecules; see Figure 7. The binding energies are summarized in Table 3. The first thing to note is that the MP2 binding energy computed using the M06 geometry yields best agreement with the MP2-optimized result for N2eg C. The MP2 binding energy using the B3LYP geometry is consistently smaller than that obtained using the M06 geometry, as found for the water system. The best estimate (computed as MP2 with half of the BSSE correction) from the M06 geometry is only 0.17 kcal/mol smaller than the MP2optimized result corrected with half of the BSSE. Thus, the MP2 results using the M06/6-31(+)G** geometry can be used to compute our best estimate, with only slightly larger errors than those for the water case.

given optimization. For the N3eg2 system, starting from the same geometry, M06 and OPLS gave different structures. However, when the optimized structures were exchanged, the OPLS found the M06 geometry and M06 found the OPLS geometry. That is, both minima existed on both potential energy surfaces, even though the first optimization yielded different minima. We should note that the numerous local minima are not a result of numerical problems. For example, structures N2eg B and N2eg C are similar in energy at the M06/6-31(+)G** level, but both are still local minima with all positive frequencies even if the optimization thresholds are tightened. That is, B and C are both local minima at the M06/6-31(+)G** level of theory and are not a computational artifact. The ethylene glycol is larger than water, which means that the calculations are longer, and as a result, we optimized only the MP2 geometry for 1 of the 10 structures that we found. In 2858


Article


31(+)G**, and they are consistently too small. The fact that the larger basis set yields binding energies that are too small might indicate that a dispersion correction might improve the binding energies in the larger basis set. However, the dimer studies suggest that the dispersion would make the dimer interaction worse. Clearly, it might be difficult to obtain a completely balanced treatment between the solvent−solvent and solvent− solute interactions. The OPLS gives a similar order of the 10 structures as the M06; the biggest difference is that M06/6-31(+)G** finds a minimum for B, whereas OPLS does not. However, the inverse is true for structure I, where it is a minimum in OPLS, but not in M06. Whereas the relative energies are in good agreement with our best estimates (MP2 or M06/6-31(+)G** when MP2 is not possible), the OPLS total binding energies are smaller than those of the best estimate. For the N3 species, a comparison of the OPLS with the best estimate (M06/6-31(+)G**) shows reasonable agreement. The OPLS correctly finds that the first binding energy of ethylene glycol to N3 model is larger than that of the second. OPLS has the N3eg2 B and C structures close in energy and the A structure less stable. As in other cases, there are small differences between the OPLS and M06, for example, the OPLS has C slightly more stable than B, whereas the reverse is true for M06. For the R2 case, as noted above, M06 and OPLS find somewhat different structures for the B form. However, the OPLS, B3LYP, and M06 find the A and B structures to be close in energy, with the energy difference being smaller for the OPLS than that for the DFT methods. Up to this point, we have focused on the most stable conformation of the ethylene glycol for each system considered, but the existence of many low-lying ethylene glycol conformations can further complicate the polymer−ethylene glycol potential energy surfaces. The first step in investigating this complexity is the study of ethylene glycol monomer

Figure 7. Ethylene glycol bonded to the three-segment phenolic model.

Comparing our best estimate to the DFT results shows that the M06 functional in the 6-31(+)G** basis set is by far the best. Expanding the basis reduced the binding energy and the M06/ aug-cc-pVTZ results are actually inferior to the M06/631(+)G** results. Clearly, the smaller basis set is benefiting from a cancellation of errors. The B3LYP/6-31(+)G** results are significantly too small, and expanding the basis makes the B3LYP binding energies even smaller. The PBE/6-31(+)G** results are closer to our best estimate than those of B3LYP/6-

Table 3. Summary of Ethylene Glycol−Phenolic Binding Energy (in kcal/mol) as a Function of Level of Theorya reaction

method OPLS

geometry

OPLS

eg + N2 → N2eg A eg + N2 → N2eg B eg + N2 → N2eg C eg + N2 → N2eg D eg + N2 → N2eg E eg + N2 → N2eg F eg + N2 → N2eg G eg + N2 → N2eg H eg + N2 → N2eg I eg + N2 → N2eg J eg + N3 → N3eg eg + N3eg → N3eg2 A eg + N3eg → N3eg2 B eg + N3eg → N3eg2 C 2eg + N3 → N3eg2 A 2eg + N3 → N3eg2 B 2eg + N3 → N3eg2 C eg + R2 → R2eg A eg + R2 → R2eg B

8.66 Cc 13.64 C 8.79 5.92 12.03 6.07 8.18 12.85 15.87 9.19 13.54 13.57 25.05 29.41 29.43 9.80 9.73

B3LYP

B3LYP/TZ B3LYP

9.09 11.06 13.27 8.26 6.85 3.13 10.75 4.36 2.33 10.87 15.90 6.77 10.11 11.12 22.67 26.01 27.02 6.15 7.32

7.26 8.02 11.08 6.46 5.44 1.99 8.91 2.37 −0.18 8.00 13.68 4.87 7.93 9.59 18.56 21.61 23.27 4.79 5.34

PBE

M06

PBE 10.53 13.97 15.91

5.69 4.40 13.28

M06/TZ M06

12.48 17.09 17.98 C 8.24 4.70 12.80 7.02 J 15.98 19.81 12.38 13.21 13.04 32.20 33.03 32.86 10.28 12.37

MP2/TZ B3LYP

10.54 14.60 16.08 6.92 3.67 11.31 5.64 13.67 17.82 11.02 11.73 12.11 28.84 29.55 29.93 8.74 10.32

11.44 17.37 17.81 12.21 9.10 4.99 13.56 7.14 7.52

PBE

M06

18.28

13.15 18.43 18.81

MP2-BSSE MP2

M06

19.11

11.03 15.37 16.11

BESTb

MP2

M06

16.15

12.09 16.90 17.46(17.63)

9.33 5.60 13.99 8.58

7.80 4.08 11.68 6.72

8.56 4.84 12.84 7.65

16.67

13.96

15.32

a The DFT calculations use the 6-31(+)G** basis set unless otherwise noted. TZ denotes the aug-cc-pVTZ basis set. Zero-point energy is not included. bThe MP2 corrected with half of the BSSE correction at the M06 geometry. The value in parentheses uses the MP2 geometry. cIndicates that the optimization resulted in the structure C geometry.

2859


Article


is a transition state at both the M06 and MP2/aug-cc-pVTZ levels of theory. We find a stationary structure at the OPLS level. As we find this stationary point using the “fire” optimization method, which includes temperature, we believe that this is a minimum and not a saddle point. M06 has minima for conformer 3 using the larger aug-cc-pVTZ bias set but not the smaller 631(+)G** basis set, and the OPLS does not find the conformer 9 structure. We should note that Boussessi et al.32 did not find a stable conformer 8 and their CCSD(T)-F12/AVTZ relative energies for the other nine conformers agree with our CCSD(T)/aug-cc-pVTZ results to within 0.05 eV. The maximum difference between the CCSD(T) and our MP2 results in the aug-cc-pVTZ basis set is only 0.21 kcal/mol, which suggests that high-order correlation effects are small. The CBS(3,4) and CBS(4,5) results are in good agreement and differ only slightly from the aug-cc-pVTZ results, which suggests that basis set effects beyond aug-cc-pVTZ are small. This suggests that the CCSD(T) results are accurate to about 0.05 kcal/mol. A comparison of the more approximate methods with CCSD(T) shows that DFT and MP2/aug-cc-pVDZ are in good agreement with the most accurate results and even M06/6-31(+)G** and OPLS are in good agreement with the best results. Cramer and Truhlar30 and Howard et al.31 computed the relative populations for the 10 conformers at 298 K in the gas phase using their ΔG values, including a correction for the degeneracy. We have performed analogous calculations using our ΔG values computed from the geometries and unscaled frequencies computed at the B3LYP/aug-cc-pVTZ and PBE/ aug-cc-pVTZ levels; see Table 5. Because of the very small energy separations between the conformations, we tightened the thresholds for the geometry optimization and integral calculation and used the “ultrafine” grid. As a result of these changes, we find that our B3LYP/aug-cc-pVTZ results are slightly different from those of Howard et al. who used the same approach. Our populations for conformers 1−3 are more similar to the best values of Cramer and Truhlar than to those of Howard et al., who also used the B3LYP/aug-cc-pVTZ approach. The population of conformer 2 is higher for PBE than for B3LYP, and this can be traced to the smaller separation for PBE, 0.05 versus 0.34 kcal/ mol. We also performed the PBE calculations at 413 K and compare them to QMD calculations performed with VASP using the same PBE functional. The results are similar but not identical. Part of this difference probably arises from the harmonic

conformations. This has been investigated previously; see, for example, Cramer and Truhlar,30 Howard et al.31 (and references therein), and Boussessi et al.32 After accounting for symmetry, there are 10 unique conformations of the ethylene glycol monomer. We first explored the possible conformation at the MP2/aug-cc-pVDZ geometry. The 10 minima found at this level of theory are shown in Figure 8, and their relative stability is

Figure 8. Two views of the 10 conformations of ethylene glycol.

shown in Table 4. Starting from the optimized MP2/aug-ccpVDZ geometries, the 10 conformations were reoptimized at other levels of theory. Using the MP2/aug-cc-pVTZ geometry, the MP2 calculations were repeated using the aug-cc-pVQZ and aug-cc-pV5Z basis sets, and these MP2 binding energies were extrapolated to the basis set limit. Also, using the MP2/aug-ccpVTZ geometry, the CCSD(T)/aug-cc-pVTZ energies were computed. For MP2/aug-pVDZ, B3LYP/aug-cc-pVTZ, and PBE/aug-ccpVTZ, all 10 of the conformers are minima, whereas at other levels of theory, one or more of the conformations are a transition state. The most notable difference arises for conformer 8, which

Table 4. Summary of Conformations of Ethylene Glycol (in kcal/mol) as a Function of Level of Theorya structureb

PBE

B3LYP

CCSD(T)

basis setc

6-31(+)G**

M06 aVTZ

aVTZ

aVTZ

aVTZd

aVDZ

aVTZ

aVQZd

aV5Zd

CBS(3,4)

CBS(4,5)

OPLS

1 2 3 4 5 6 7 8 9 10

0.00 0.71 2.73 3.14 3.33 3.63

0.00 0.34 0.99 2.67 2.82 2.88 2.99

0.00 0.45 0.99 2.60 2.79 2.82 3.01

0.00 0.46 0.99 2.61 2.80 2.84 3.03

0.00 0.47 1.02 2.59 2.80 2.84 3.03

0.00 0.46 0.99 2.62 2.82 2.85 3.04

f

0.00 0.50 1.25 2.48 2.90 2.50 3.67 2.77

2.96 3.30

0.00 0.43 1.02 2.58 2.78 2.81 2.99 3.19 3.12 3.55

0.00 0.43 0.94 2.61 2.78 2.80 2.99

3.41 4.12

0.00 0.34 0.81 2.35 2.43 2.40 2.58 2.70 2.76 3.14

0.00 0.32 0.80 2.56 2.64 2.59 2.78

f

0.00 0.05 0.71 2.84 2.69 2.47 2.57 2.49 2.88 3.05

3.13 3.60

3.14 3.61

3.15 3.61

3.15 3.62

g

e

3.11 3.49

MP2

f

3.10 3.57

4.19

a

Zero-point energy is not included. bSee Figure 1. cThe basis sets are described using a short-hand notation, where aVXZ indicates the aug-cc-pVXZ basis settlement. dNot optimized; the aug-cc-pVTZ geometry is used. eHas an imaginary frequency; if displaced along the vector corresponding to the imaginary mode, the structure collapsed to structure 1. fHas an imaginary frequency; if displaced along the vector corresponding to the imaginary mode, the structure collapsed to structure 2. gCollapsed to structure 2 during the optimization. 2860


Article


geometry, the OPLS optimization was started from the M06optimized structures and M06 was started from the OPLS geometries. The results are summarized in Table 6. The binding energies are given with respect to both the ethylene glycol with the lowest energy (conformer 1) and the conformer used in the starting geometry. Using the binding energies with respect to the most stable ethylene glycol conformation, it is easy to see that for the M06 approach, conformations 3 and 9 converge to the most stable structure. For OPLS, conformations 2, 3, and 9 all converge to the most stable structure regardless of the starting structures. We also find different final structures depending on the choice of the starting geometry. This work again shows how complex the ethylene glycol−phenolic surface is. We find multiple minima depending on the orientation of the ethylene glycol relative to the phenolic and even multiple minima depending on the conformation of the ethylene glycol. Given the extreme complexity of this surface, it is difficult to even guarantee that the methods are fully sampling all of the possible minima; however, for those structures, we have found that the OPLS yields a reasonable description. III.IV. First Solvation Shell Structure. As noted above, the solvent−phenolic OH interaction is much larger than the solvent−ring interaction. The N3 structure shows that the phenolic OH groups interact with themselves, thus reducing the potential for solvent−phenolic interactions. As the phenolic ring area is much larger than the OH area and the phenolic OH groups interact with themselves, the question arises as to the size of the total solvent−phenolic interaction and what part of it arises from the ring versus that from OH. To address this question, we take the geometry from MD simulations (described in our companion paper (10.1021/acs.jpcb.7b00326)) on a ninesegment phenolic model surrounded by 7900 water or 1700 ethylene glycol molecules. The geometry of the phenolic model and the first shell of 141 waters is shown in Figure 9. We evaluated the energy of the structure with the first shell of waters using the M06/6-31(+)G** approach. We also evaluated the energy of just the phenolic molecule and just cluster of waters. No geometric relaxation is allowed in any of the three systems. This yields a binding energy of 54.7 kcal/mol. This procedure was repeated using only the five water molecules closest to the phenolic OH groups. This yields a binding energy of only 13.2

Table 5. Populations of the Ethylene Glycol Conformer gas phase

liquid phase

conformation

B3LYPa

PBE

PBE

QMD

QMD

QMD

temp 1 2 3 4 5 6 7 8 9 10

298 54.7 27.7 13.1 0.5 1.3 0.6 0.5 0.6 0.6 0.4

298 46.6 37.2 12.7 0.2 0.8 0.5 0.4 0.9 0.4 0.4

413 40.1 31.5 16.4 0.7 2.7 1.5 1.3 2.7 1.6 1.6

413 47.3 26.9 14.4 0.0 0.8 0.1 0.2 3.8 2.1 4.6

298 21.1 25.9 14.7 0.5 3.4 3.0 0.3 11.1 6.3 13.8

413 25.9 17.1 12.5 2.5 7.7 3.9 2.3 7.4 8.6 12.0

a

The B3LYP and PBE approaches were computed using the aug-ccpVTZ basis set using Gaussian 09, whereas the QMD calculations were performed using a plane wave basis set using VASP.

approximation used in the ΔG formulation; however, small differences in energies between the PBE/aug-cc-pVTZ and PBE/ plane wave formulations can also contribute to the differences in the populations. Using the QMD approach, the liquid phase was simulated at 298 and 413 K, and this is given in the table. The first thing to note is that the population is more spread out in the liquid phase than in the gas phase. As the population is more spread out in the ethylene glycol liquid, one can expect that a phenolic molecule in liquid ethylene glycol is more likely to be found interacting with different conformers than the gas phase ethylene glycol results would suggest. Given the presence of multiple conformers in the liquid, it is important to investigate the interaction of different conformers with the polymer. Therefore, we consider the effect of the ethylene glycol conformation on the binding energy of ethylene glycol to the two-segment model of the novalac form, that is, the N2 model. We start from the N2eg C model, which has the largest binding energy. The existing ethylene glycol is replaced with the other nine ethylene glycol conformations, and the new merged system is optimized. The M06/6-31(+)G** and OPLS models are used. In addition to starting from “replaced”

Table 6. Summary of Binding Energies of the Conformations of Ethylene Glycol (in kcal/mol) to the Two-Unit Phenolic Modela conformationb

M06/6-31(+)G**

geomc binding energy 1 2 3 4 5 6 7 8 9 10

replaced d

OPLS OPLS

replaced

M06

Δn

Δ1

Δn

Δ1

Δn

Δ1

Δn

Δ1

17.97 13.06 19.44 12.56 12.68 15.06 12.32 10.93 21.37 17.39

17.97 12.34 17.96 9.83 9.53 11.73 8.69 7.68 17.96 13.28

17.97 18.68 19.44 15.61 16.03 16.23 16.53 12.77 21.36 15.81

17.97 17.97 17.96 12.89 12.88 12.90 12.90 9.52 17.95 11.69

11.50 12.00 12.75 12.52 12.94 12.54 13.71 1.22 12.75 1.41

11.50 11.50 11.50 10.05 10.04 10.04 10.04 −1.55 11.50 −2.78

11.50 12.00 12.75 12.52 7.92 8.13 7.91 10.91 12.75 13.02

11.50 11.50 11.50 10.05 5.02 5.63 4.24 8.14 11.50 8.82

a Zero-point energy is not included. bIndicates the ethylene glycol conformation in the starting geometry. cIndicates the starting geometry: “Replaced” indicates that the ethylene glycol from N2eg C was replaced with the free ethylene glycol conformer, “M06” indicates that the M06optimized geometry was used as the starting geometry, and “OPLS” indicates that the the OPLS geometry was used for the starting geometry. dThe ethylene glycol binding energy: Δ1 means the binding energy is with respect to the most stable ethylene glycol conformer (1), whereas Δn means that the conformer n energy is used to compute the binding energy.

2861


Article


Figure 9. A nine-segment phenolic model surrounded by a shell of water molecules. The phenolic molecule has folded so the head and tail are adjacent.

kcal/mol. The analogous values for the OPLS force field are 93.9 and 23.4 kcal/mol. This difference in binding energy between the M06 and OPLS values is much larger than that for the smaller models discussed above, and this is probably mostly a result using the OPLS geometry for M06 calculations. We should note that attempts to optimize the geometry at either level resulted in dramatic changes as these structures are not stable without the surrounding water molecules. Despite these limitations, it is clear that while the water−phenolic OH interaction is much stronger than the water−ring interaction for a single water molecule, the sheer number of water−ring interactions means that they contribute more to the total binding energy than the stronger water−OH interactions. We performed equivalent calculations for ethylene glycol and the nine-segment model. The M06/6-31(+)G** value is 16.0 for the first shell of ethylene glycol molecules and 8.6 kcal/mol for the five nearest ethylene glycol molecules, whereas the analogous OPLS values are 145.2 and 47.0 kcal/mol. The difference between M06 and OPLS is much larger than that for water; thus, we focus on the OPLS values. As in the case of water, the ring binding is larger than the OH binding because of the much larger number of molecules interacting with the rings. Even these models with only the first shell of solvent molecules are very simple as the orientation of the solvent molecules around the phenolic polymer is determined by both the solvent−solvent and solvent−phenolic interactions. In this regard, we note that the geometric relaxation at the OPLS level of the system with only the first shell resulted in a 364 kcal/mol energy lowering for the water system and 550 kcal/mol lowering for the ethylene glycol system. Thus, any discussion of the solution is deferred to the companion paper (10.1021/acs.jpcb.7b00326); however, the current ab initio calculations give some insight into the difference between water and ethylene glycol, such as the ethylene glycol potential surface having many more local minima than those in water, and show the competition between solvent−phenolic OH and phenolic OH−phenolic OH interactions. These ab initio calculations also support the use of the OPLS-derived potential.

potential. This arises because ethylene glycol has several possible conformers and two OH and two CH2 groups that lead to many more possible interactions with the phenolic polymer than those of water. This surface is further complicated by the fact that geometry optimization starting from different geometries or with different parameters in the optimization can yield different results. In spite of these problems, we have been able to compare various levels of theory for selected systems of phenolic polymer models interacting with water or ethylene glycol. Although not complete, we have sampled a sufficient number of structures to have some confidence in our calibration study. Given the size and complexity of these systems, it is difficult to definitively compute a benchmark value for comparison, but we have established what we feel are reliable values for comparison, which are based on MP2 calculations on the smaller systems. However, we must admit that excluding the dimers our best estimates have some uncertainty. For H2O−phenolic interactions, we find that M06/aug-ccpVTZ is in best agreement with our best estimate, which is derived from MP/aug-cc-pVTZ corrected for 50% of the BSSE error. For the water dimer, PBE/aug-cc-pVTZ has the smallest error, but the error for the M06/aug-cc-pVTZ is only marginally larger. Thus, the M06/aug-cc-pVTZ approach seems to be an excellent compromise between the H2O−H2O and H2O− phenolic interactions. For the ethylene glycol systems, the errors are larger. M06/6-31(+)G** gives the best description of the ethylene glycol−phenolic interaction. Improving the basis set makes the agreement worse. The M06/aug-cc-pVTZ approach has the smallest error for the ethylene glycol dimer, whereas the M06/6-31(+)G** binding energy is significantly too large. Thus, it is more difficult to obtain a balanced description of the ethylene glycol system than that of the water system. Overall, the binding energies computed with the OPLS potential are in reasonable agreement with our best estimates and the DFT functional of choice. The relative energetics of different structures follows our best estimates, and even the total binding energies are reasonable. Given the size of the errors, we suspect that to go beyond the OPLS would require the generation of a polarizable force field.

IV. CONCLUSIONS We find the phenolic−ethylene glycol potential energy surface to be much more complicated than that of the phenolic−water 2862


Article


■

(18) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465. (19) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-Consistent Molecular Orbital Methods. 25. Supplementary Functions for Gaussian Basis Sets. J. Chem. Phys. 1984, 80, 3265−3269 and references therein. (20) Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (21) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron affinities of the First-Row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796−6806. (22) Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-Set Convergence of Correlated Calculations on Water. J. Chem. Phys. 1997, 106, 9639−9646. (23) Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (24) Frisch, M. J.; et al. Gaussian 09, revision C.01; Gaussian, Inc.: Pittsburgh, PA, 2009. (25) Plimpton, S. J. Fast Parallel Algorithms for Short-Range Molecular Dynamic. J. Comput. Phys. 1995, 117, 1−19. (26) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561; Ab Initio Molecular-Dynamics Simulation of the Liquid-MetalAmorphous-Semiconductor Transition in Germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251−14269. (27) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (28) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (29) Rocher-Casterline, B. E.; Ch’ng, L. C.; A Mollner, A. K.; Reisler, H. Communication: Determination of the Bond Dissociation Energy (D0) of the Water Dimer, (H2O)2, by Velocity Map Imaging. J. Chem. Phys. 2011, 134, No. 211101. (30) Cramer, C. J.; Truhlar, D. G. Quantum Chemical Conformational Analysis of 1,2-Ethanediol: Correlation and Solvation Effects on the Tendency To Form Internal Hydrogen Bonds in the Gas Phase and in Aqueous Solution. J. Am. Chem. Soc. 1994, 116, 3892−3900. (31) Howard, D. L.; Jørgensen, P.; Kjaergaard, H. G. Weak Intramolecular Interactions in Ethylene Glycol Identified by Vapor Phase OH-Stretching Overtone Spectroscopy. J. Am. Chem. Soc. 2005, 127, 17096−17103. (32) Boussessi, R.; Senent, J. L.; Jaïdane, N. Weak Intramolecular Interaction Effects on the Torsional Spectra of Ethylene Glycol, an Astrophysical Species. J. Chem. Phys. 2016, 144, No. 164110.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Charles W. Bauschlicher Jr.: 0000-0003-2052-332X Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Bucholz, E. W.; Haskins, J. B.; Monk, J. D.; Bauschlicher, C. W.; Lawson, J. W. Phenolic Polymer Solvation in Water and Ethylene Glycol, I: Molecular Dynamics Simulations. J. Phys. Chem. B 2017, DOI: 10.1021/acs.jpcb.7b00326. (2) Bartlett, R. J. Many-Body Perturbation Theory and Coupled Cluster Theory for Electron Correlation in Molecules. Annu. Rev. Phys. Chem. 1981, 32, 359−401. (3) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A Fifth-order Perturbation Comparison of Electron Correlation Theories. Chem. Phys. Lett. 1989, 157, 479−483. (4) Lane, J. R. CCSDTQ Optimized Geometry of Water Dimer. J. Chem. Theory Comput. 2013, 9, 316−323. (5) Feyereisen, M. W.; Feller, D.; Dixon, D. A. Hydrogen Bond Energy of the Water Dimer. J. Phys. Chem. 1996, 100, 2993−2997. (6) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (7) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (8) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (9) Kony, D.; Damm, W.; Stoll, S.; van Gunsteren, W. F. An Improved OPLS-AA Force Field for Carbohydrates. J. Comput. Chem. 2002, 23, 1416−1429. (10) Price, D. J.; Brooks, C. L., III A Modified TIP3P Water Potential for Simulation with Ewald Summation. J. Chem. Phys. 2004, 121, 10096−10103. (11) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (12) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (13) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. Errata: 78, 1396 (1997). (14) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215− 241. (15) Zhao, Y.; Truhlar, D. G. A New Local Density Functional for Main-Group Thermochemistry, Transition Metal Bonding, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Phys. 2006, 125, No. 194101. (16) Grimme, S. Semiempirical GGA-type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (17) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parameterization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, No. 154104. 2863


Phenolic Polymer Solvation in Water and Ethylene Glycol, II: Ab Initio

Recommend Documents