Fragment Molecular Orbital Based Parametrization Procedure for

Article Cite This: J. Phys. Chem. B XXXX, XXX, XXX−XXX

pubs.acs.org/JPCB

Fragment Molecular Orbital Based Parametrization Procedure for Mesoscopic Structure Prediction of Polymeric Materials Koji Okuwaki,† Yuji Mochizuki,*,†,‡ Hideo Doi,† and Taku Ozawa§ †

Department of Chemistry and Research Center for Smart Molecules, Faculty of Science, Rikkyo University, 3-34-1 Nishi-ikebukuro, Toshima-ku, Tokyo 171-8501, Japan ‡ Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan § JSOL Corporation, 2-5-24 Harumi, Chuo-ku, Tokyo 104-0053, Japan ABSTRACT: In the analyses of miscibility behaviors of macromolecules and polymers, dissipative particle dynamics (DPD) simulations are generally performed. In these simulations, the so-called χ parameters describing the effective interactions among particles are crucial. It has been known that such parameters can be obtained within the classical or empirical force field frameworks. However, there is a potential problem that charge transfer and polarization occasionally occur. Additionally, satisfactory reference parameters are not available for some cases. Therefore, we developed a new procedure to evaluate the set of parameters by using the ab initio fragment molecular orbital (FMO) method which can provide the set of interaction energies among segments as polymer units. Moreover, we evaluated the anisotropy of molecules by using the FMO-based effective interaction parameters for three standard binary mixture systems (hexane−nitrobenzene, polyisobutylene−diisobutyl ketone, and polyisoprene−polystyrene). The calculated values showed good agreement with the experimental values with about 10% errors. parameters are closely connected to the so-called χ parameters involving the interactions in the Flory−Huggins theory.6,13,14 Unfortunately, it has been well-known that the evaluation of reliable χ parameter is not an easy task. There are mainly two representative routes of χ parameter predictions. The first route depends on the solubility parameter (SP) values. The SP value method was devised by Hildebrand,15 and various empirical estimation models such as atomic group contribution models,16,17 Bicerano,18 and so on were developed. Note that the SP values could be obtained from molecular simulations for single molecules as well.19 The second route is based on the evaluation of interactions between heterogeneous molecules, due to the contact energy between segments20 and the difference in cohesive energies of the aggregated model.21 Hereafter, we focus on the contact energy method proposed by Fan et al.20 By Fan’s method,20 it is possible to directly determine the interactions between the contacting particles in a molecular level approach. The temperature-dependent χ parameters are evaluated by a Monte Carlo (MC) procedure (as shown later), where the pairwise contact energies are calculated for all possible combinations of basic repeating units of polymer (or segments). Note that an amount of configurations for each pair should be computed for each pair. In ref 20, the classical force

1. INTRODUCTION Polymeric materials, such as synthetic plastics and rubbers, have been used and distributed widely. In recent years, better performances of the polymeric materials have been required toward the sophistication of products. Structures of the mesoscale region (which is intermediate between nanoscale and macroscale) are well-known to have a major effect on the properties of the final products of material. Therefore, the prediction and control of mesoscale structures is an important issue. To predict mesoscopic properties and structures, several coarse-grained simulation techniques such as dynamic mean field theory1,2 and dissipative particle dynamics (DPD)3−8 are applied. In these methods, polymer chains are generally mimicked as coarse-grained bead−spring models. The coarsegrained simulation can handle long time scales of large scale molecular systems, relative to atomic simulations such as molecular dynamics. To date, a number of molecular systems have been studied with DPD simulations, some of which are addressed as below. Arai studied structural analysis of a telechelic polymer solution,9 whereas Yamamoto et al. simulated the mesoscopic structures of the polyelectrolyte membrane10 as well as the spontaneous vesicle formation of amphiphilic molecules.11 Shillcock et al. investigated the equilibrium structure and lateral stress distribution of amphiphilic bilayers.12 In DPD simulations, the parameters describing effective interactions among particles or segments are crucial. Such © XXXX American Chemical Society

Received: August 24, 2017 Revised: December 5, 2017

A

DOI: 10.1021/acs.jpcb.7b08461 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B field (FF) parameters were used to estimate the crucial pairwise contact energies. The use of FF could, however, be problematic when charge transfer and polarization are substantial. In other words, the interactions in real systems would not be sufficiently modeled by FF parameters in certain cases. In order to solve the above-mentioned problem, it is desirable to calculate the pairwise interaction energies by quantum mechanical methods such as ab initio molecular orbital (MO) calculations so that both polarization and charge transfer can be incorporated. In 2015, Sepehr et al. evaluated the interaction parameters of DPD for Nafion as a demanding system for classical FF modeling, based on MO calculations.22 Although the usual MO-based evaluation of χ parameters is attractive, its applicability would be limited due to the increased computational costs when the molecular sizes of segment pairs become large for coarser modeling. In this paper, we report a new approach to calculate the set of χ parameters, by using the fragment molecular orbital (FMO) method.23 The FMO method has been known as an efficient MO framework applicable to large scale molecules (e.g., proteins and molecular clusters) with chemical reliability, through highly parallel executions. The list of interaction energies among fragments is directly obtained once the FMO calculation completes. Additionally, detailed analyses on the nature of interactions can be performed when necessary. These features of FMO are convenient for the present subject, as discussed later. The rest of this paper is configured as follows. In section 2, the χ parameter evaluations are described in detail. In section 3 the scheme of test calculations is briefed. In section 4, the results of the test calculation are presented and discussed.

calculations, care should be taken for potential limitations due to the usage of FF parameters, as mentioned already. For further description of our procedure, the stepwise procedures of ref 20 are summarized as follows. (i) Preparation of Proper Structures for Segments. The fundamental structures of segments as repeating units in polymer components are constructed. These structures may be optimized using molecular mechanics. (ii) Configuration Setup and Calculation of Interaction Energies among Segments. For each segment pair, the interaction energies between segments are calculated for a large number (say a few thousands to ten thousands) of configurations. Methods of generating a configuration are described in the following sequence. Figure 1 shows the generated configurations of

2. THEORETICAL METHOD OF χ PARAMETER EVALUATION 2.1. Conventional Parameter Calculation Method. In the Flory−Huggins lattice theory, the free energy change (ΔG) for a binary system (consisting of two components) is expressed as13,14

segments, and Figure 2 illustrates the flow of the parameter calculation from configuration generation: refer to these two figures when needed, hereafter.

φ φ ΔG = 1 ln φ1 + 2 ln φ2 + χφ1φ2 RT x1 x2

Figure 1. Generated configurations of segment (see text). Blue and red spheres represent “segment 1” and “segment 2”, respectively.

(1)

where φi and xi are the volume fraction and the chain length (i = 1, 2 for the two components), respectively. The first and the second terms on the right side describe the entropy change, and the third term corresponds to the enthalpy change. The χ parameter is defined as

χ=

Z ΔE12 RT

(2)

(3)

Figure 2. Flow of the parameter calculation from generated configurations of segments. Blue and red spheres represent “segment 1” and “segment 2”, respectively. The average interaction energy per pair Eij is calculated by using the Metropolis MC method for the generated configurations. The symbols of “circle” and “cross” indicate acceptance and rejection in the MC sampling. The coordination number Zij is multiplied to the resulted energy.

Eij in eq 3 is the average interaction energy between the segments i and j in the two components, and ΔE12 represents the energy gain per segment due to the mixing. These relations imply a scale-down of problem from mesoscale to nanoscale (segment). Fan et al. proposed the procedure to calculate Z and ΔE12, based on MC simulations with classical FF.20 Although the χ values could be computed only through theoretical

(ii-1) The geometric centers of both segments are placed near the origin of the Cartesian coordinate frame. (ii-2) The particular orientation and translation direction of segment 2 (for segment 1 fixed) are determined by randomly choosing the three Euler angles and the vector from the origin to the surface of a unit sphere.

where Z is the coordination number of the model lattice and the contact energy ΔE12 is given by the following equation: ΔE12 = E12 −

1 (E11 + E22) 2

B


Article

The Journal of Physical Chemistry B

(MP2),32−34 can be incorporated in an additive manner. Equation 5 is then rewritten to define the pairwise interaction energies ΔẼ IJ as

(ii-3) Segment 2 is moved following the orientation and vector determined at step ii-2, until the van der Waals surfaces of the molecules come in contact with each other. (ii-4) The acceptance of the above-generated configuration is judged with the Metropolis MC algorithm.24 uvv′ is the probability for the transition from configuration v to configuration v′ at temperature T, just described as ⎛ ΔE ⎞ uvv = exp⎜ − vv ′ ⎟ ′ ⎝ RT ⎠

E=

I>J

∑ EIJ − (Nf I>J

(4)

− 2) ∑ EI I

I

(6)

EI′ in eq 5 is a fictitious monomer energy excluding the ESP contributions. Note that ΔẼ IJ has been called the pair interaction energy (PIE)35 or inter-fragment interaction energy (IFIE),36,37 and these values are useful quantities in analyzing the nature of fragment interactions in a given target system.23,26,27 For more detailed analyses, the technique of pair interaction energy decomposition analysis (PIEDA)35,38 is available, where a certain pair energy is expressed as the sum of contributions from electrostatic interaction (ES), exchange− repulsion (EX), charge transfer with some other mixed effect (CT + mix), and (attractive) dispersion interaction (DI): the correlated calculation just as MP2 is required to evaluate DI. Because the contact energies between segments (eqs 2 and 3) are evaluated with the list of pair interaction energies by the FMO calculations, the quantum mechanical effects of both polarization and charge delocalization can be taken into account, unlike the classical FF framework employed in ref 20. When the segment sizes grow, such larger segments may be divided into appropriate fragments and the associated interaction energies are computed through proper summation. This ability of FMO23,26,27 provides an applicability to large segments straightforwardly. Additionally, PIEDA35,38 is usable to characterize the specific segment interactions when needed. (ii) Reevaluation of Distance between Nonbonded Atoms. Fan et al.20 used the set of van der Waals radii of atoms to generate the configurations of a target segment pair. However, there is a presumable problem that closer contacts are necessary to describe certain quantum mechanical interactions such as hydrogen bonding. In our method, the distance parameters between nonbonded atoms in each segment pair are reevaluated when needed, through a series of MO-based geometry optimizations. (iii) Evaluation of the Surrounding Environment with Anisotropy. In ref 20, the coordination number Z in the numerator of eq 2 was calculated from a spatial filling model, and Z was multiplied to the contact energy averaged by the MC-based sampling. Namely, isotropy was assumed for the surrounding environment in the Flory−Huggins sense.13,14 However, such an assumption might not always stand, depending on both the molecular shapes and the mutual orientations of segment pairs in actual cases. For example, if a specific portion in a certain segment pair ij provides a strong stabilization, this would lead to an overestimation in the average interaction energy Eij through biased samplings. To overcome this issue, the space around a center segment is split into six regions from its center of mass (see Figure 3), and then the anisotropy factors are introduced (for each pair). The sampling weight of each region l (=1−6), dl, is calculated by the following formula from the results of the MC sampling:

where ΔEvv′ in the numerator is the energy difference between sequential configurations v and v′. The average value of the energies of accepted configurations by the MC selection is defined as the interaction energy of each segment pair (E11, E12, and E22). (iii) Calculation of Coordination Number Z. As shown in eq 2, the χ parameter depends on not only interaction energy Eij but also on the coordination number Z corresponding to the number of segments adjacent to the segment of interest. Namely, the value of Z should have a significant effect on the final value of χ as well. In ref 20, the number of segments which can be placed around the given center segment is calculated by repeating the procedures ii-2 and ii-3. For a binary system with tentative component labels 1 and 2, four different combinations are possible due to the following situations. The center segment can be segment 1 or 2 and the surrounding segment also can be segment 1 or 2: or there are four possible relations, symbolically described as 1−1, 1−2, 2−1, and 2−2. The average of the four combinations is considered in evaluating Z, because all segments are assumed to be spheres of an identical size in the Flory−Huggins theory.13,14 2.2. Proposed Method. In our proposed method, the pair energies Eij in eq 3 are evaluated by the FMO method23 originally developed by Kitaura et al. in 1999.25 This is the fundamental appealing point over Fan’s method.20 The details of our method are given as follows. (i) FMO-Based Contact Energy. The typical FMO calculation is briefed as follows (see ref 23 as well as refs 26−28 when necessary). A certain target molecular system is divided into appropriate fragments called “monomers”. The Hartree−Fock (HF) procedure is used to determine the MO set of each fragment under the environmental electrostatic potential (ESP) at the monomer stage.29−31 The grand iteration for fragment monomers is iterated until the self-consistent-charge (SCC) condition is satisfied. The intrafragment processing is parallelized, and the task of monomer list is processed in parallel as well (or two-layer parallelism). Once the monomer SCC iteration completes, a series of HF calculations of dimers (pairs of monomers) are then performed under the ESP of monomer SCC, with a similar parallelism with the case of monomers. The polarizations and charge transfers (or charge delocalizations) are taken into account at the monomer stage and the dimer stage, respectively. The total electronic energy E is obtained by the following simple formula: E=

∑ ΔEIJ̃ + ∑ EI′

dl =

(5)

where EI is the energy of monomer I, Nf is the number of monomers, and EIJ is the energy of dimer IJ. Note that the energy corrections due to the electron correlation through, for example, the second-order Møller−Plesset perturbation theory

Al A max

(7)

where Al is the number of acceptations in the region l and Amax is the number of acceptances in the region with the most acceptances. The number of trials is roughly 10 000 times. C


Article


set as 2000 for each segment pair (1−1, 2−2, and 1−2). The list of pair interaction energies (IFIE36,37) is obtained for the generated configuration, by a series of FMO calculations with the ABINIT-MP program.34 Then, the MC sampling (at a given temperature) is used to evaluate the average contact energy with the regional weights of eq 9. This procedure is applied to each segment combination: segments 1−1, segments 1−2, segments 2−1, and segments 2−2 (see also Figure 2). Subsequently, the χ parameter was calculated using the following equation (which is slightly different from eq 2): χ= Figure 3. Division of space around a center segment for scaling. Six regions are defined according to a cube. Blue and red spheres represent “segment 1” and “segment 2”, respectively.

(10)

where Zij is the coordination number of segment j around segment i, and Sij is a scaling factor (which associates with anisotropy) of the combination of segment j around segment i. The numerator of the right-hand side of eq 10 corresponds to the energy gain due to mixing of segments, because (E12Z12S12 + E21Z21S21) is the energy when the center segment has contact with maximum allowable different kinds of segments while (E11Z11S11 + E22Z22S22) is the energy when the center segment has contact with only segments of its kind. Note here that E12 and E21 as well as S12 and S21 are the same values as each other in the present method. If the central segment is a unit of polymer chain, an excluded volume of connecting portion as 2 (as suggested in ref 20) should be subtracted from the corresponding Z. (v) Evaluation of the Phase Transition Upper Critical Temperature. In the Flory−Huggins theory,13,14 the upper critical solution temperature of a binary mixture (Tc) can be used as a reliability indicator of the temperature dependence of estimated χ paramerters. The theoretical value of the critical point χc is calculated by the following equation:13

Difference in the numbers of acceptations may correspond to the difference in free energy as below: ⎛ A ⎞ −ln⎜ l ⎟ ≅ ΔG ⎝ A max ⎠

(8) 13,14

However, in the Flory−Huggins theory, it was assumed that the segment molecules are isotropic and ΔG is zero. Therefore, the following S is adopted as a parameter representing the anisotropy of the system: S=

∑ l

dl N

(E12Z12S12 + E21Z 21S21) − (E11Z11S11 + E22Z 22S22) RT

(9)

where N is the number of regions. If a given system is completely isotropic in an ideal case, dl is the same number in all the regions, and thus S is 1. Conversely, if the system has some anisotropy, S becomes smaller than 1. The S values are evaluated 100 times with different divisions of the space. Subsequently, the results are averaged, and the averaged value is used as a scaling factor for the next step. (iv) Generation of Conf iguration and Calculation of Eij and Z. Figure 4 sketches the workflow of our proposal for the evaluation of χ parameters. The generation of configuration and also the calculation of Z are performed (recall section 2.1, iii) using the J-OCTA program,39 where the contact radius refinements for nonbonded atoms are made with MO calculations. By considering the relative computational cost of FMO to FF, the number of generated configurations is typically

( χ = c

1 na

+ 2

1 nb

2

)

(11)

where na and nb (degrees of polymerization) are the quotients of the polymer molecular weight (Mw) and the segment molecular weight: this value is set to 1 for the case of solvent. A series of χ values may be sequentially evaluated (by a typical step size of 1 K) for the temperature range near χc. The

Figure 4. Computational flow of FMO-based evaluation of χ parameter (see text). D


Article

The Journal of Physical Chemistry B temperature when the calculated value of χ coincides with χc corresponds to the theoretical Tc value.

Table 1. Small Molecular Models Used To Reevaluate Distances between Nonbonding Atoms

3. SCHEME OF TEST CALCULATION To check the usability and reliability of our method, the sets of χ parameters were calculated for three standard binary systems (which were tested in ref 20) as (i) hexane−nitrobenzene (solvent−solvent) (ii) polyisobutylene−diisobutyl ketone (polymer−solvent) (iii) polyisoprene−polystyrene (polymer−polymer) Figure 5 illustrates the fundamental molecular structures of these systems.

atom type

definition

molecule

c3 c2 ca hc ha no o on

sp3 carbon sp2 carbon aromatic carbon hydrogen bonded to aliphatic carbon hydrogen bonded to aromatic carbon nitrogen of the nitro group oxygen with one atom oxygen of the nitro group

C2H6 C2H4 C6H6 C2H6 C6H6 CH3NO2 CH3CHO CH3NO2

Table 2. Reevaluated Distances between Nonbonding Atoms molecule 1

molecule 2

distance (Å)

c3

c3 c2 ca hc ha no o on c2 ca hc ha ca hc ha no

3.65 3.57 3.63 2.90 2.88 3.38 3.42 3.42 3.57 3.51 3.00 3.00 3.50 2.79 2.79 3.28

c2

Ca

molecule 1 molecule 2 ca hc

ha

no

o on

o on hc ha no o on ha no o on no o on o on

distance (Å) 3.42 3.42 2.37 2.40 3.03 3.20 2.86 2.58 3.29 2.90 3.17 3.64 3.12 3.12 3.30 3.30

segments (polyisobutylene, polyisoprene, and polystyrene), two monomers were combined as a single segment consisting of two fragments. Interaction energies between these polymer segments were calculated by summing up the interactions between the fragments belonging to different segments. For example, if one segment consisted of fragments 1 and 2 and the other consisted of fragments 3 and 4, the interaction between the two segments ΔẼ seg was calculated by the following equation:

Figure 5. Representative molecular structures of segment pairs used for test calculations. (i) Hexane−nitrobenzene (solvent−solvent) system. (ii) Polyisobutylene−diisobutyl ketone (polymer−solvent) system. (iii) Polyisoprene−polystyrene (polymer−polymer) system.

̃ = ΔE13 ̃ + ΔE14 ̃ + ΔE23 ̃ + ΔE24 ̃ ΔEseg

(12)

where ΔẼ IJ was the interaction energy between fragments I and J (refer to eq 6). Each solvent molecule (hexane, nitrobenzene, and diisobutyl ketone) was treated as a single segment as well as a single fragment. Once the FMO calculations completed, the lists of average interaction energies Eij (and Sij) were obtained. Note that all the necessary Zij values were separately computed (see section 2.2, iv, and also section 2.1, iii). Finally, the χ parameters were evaluated for the respective temperatures (with the given range), by the MC sampling with the prepared Zij and Sij factors. In some cases, Fan’s FF-based procedure20 was also used for comparison, where the general AMBER force field (GAFF43) and the Dreiding FF44 were adopted. A breakdown analysis in interaction energies was made with PIEDA35,38 for the hexane−nitrobenzene system.

For this testing, we reevaluated the distances between nonbonding atoms for the generation of segment pair configurations: recall the descriptions in section 2.2, iii. Several small molecules were employed in such reevaluations, as compiled in Table 1. The geometry optimizations were done at the dispersion-corrected B97D40 density functional level with the 6-31G(d) basis set,41 by using the Gaussian 09 program.42 The reevaluated distances are listed in Table 2. The differences from the original values20 should be reflected in the radial distribution functions, as will be shown in section 4. For a certain binary system, 2000 pair configurations for the respective combinations of 1−1, 2−2, and 1−2 labels were generated and subjected to the standard FMO-MP2/6-31G(d) calculations with ABINIT-MP program.27 For the polymer E


Article


4. RESULTS AND DISCUSSION 4.1. Hexane−Nitrobenzene (Solvent−Solvent) System. Prior to discussing the results of χ parameters, the effect of reevaluated distances is addressed with the first test case, or the hexane−nitrobenzene system45 (Tables 1 and 2). Figure 6

Figure 6. Radial distribution functions between nonbonding atoms of accepted configurations in the MC sampling for the nitrobenzene− nitrobenzene pair. Each peak is aligned to 1. Label of “Original” means the van der Waals contact based result by Fan’s original scheme,20 and “Present” corresponds to the present result with reevaluations through MO-based geometry optimizations.

Figure 7. Calculated χ values of hexane−nitrobenzene system as a function of temperature. Symbols of square, cross, and triangle correspond to the results of FMO with scaling (S), FMO without scaling, and GAFF,43 respectively. Dotted horizontal line indicates χc shown in Table 3. Curve with squares (FMO with scaling) can be best fitted with χ = −2.23 + (1.212 × 103)/T.

plots the radial distribution functions (RDFs) between nonbonding atoms of accepted configurations in the MC sampling for the nitrobenzene−nitrobenzene pair. In Figure 6 , RDF by our method is compared with that by Fan’s original scheme.20 As expected from the reevaluation, the RDF by the present method indicates that the accepted configurations are increased in short distances. It is also suggested that the sampling efficiency is effectively improved by the reevaluation of distance parameters. Table 3 shows the values of Eij, Sij, and Zij at 300 K as well as the Tc. It is found that the interaction energy of the hexane−

polymerization is 1 (refer to eq 11). Thus, the formal χc is 2.0, as given in Table 3. Our estimated Tc is 286 K, whereas the reported value of ref 20 was 268 K. Our value is in better agreement with the experimental value of 293 K.46 The PIEDA35,38 results are compiled in Table 4, where three characteristic configurations are taken: (A) the most stable configuration, (B) the most unstable configuration, and (C) the configuration which has the largest energy difference between the FMO-IFIE and FF (Dreiding44) values, as illustrated in Figure 8. Interestingly, the configurations with parallel alignment are stable all three pairs. For the hexane−hexane pair (label 1−1), the dispersion interaction provides the leading stabilization (refer to the DI column in configurations A and C of Table 4), as expected from the nonpolar nature. In contrast, the stabilization for the nitrobenzene−nitrobenzene pair (label 2−2) consists of the electrostatic (ES) attraction and charge transfer (CT + mix) as well as the dispersion interaction. Both the polarizable benzene ring with quadrupoles47 and the electronegative nitro group may be responsible for this result. The electrostatic repulsion between two nitro groups leads to the sizable destabilization in configuration B. For the hexane− nitrobenzene pair (label 1−2), configuration A has the stabilization of −7.1 kcal/mol from ES, CT + mix, and DI contributions whose values are smaller than the corresponding values of the nitrobenzene−nitrobenzene pair. Although there should be some overestimation in interaction energies evaluated at the FMO-MP2/6-31G(d) level, this value would be comparable to −6.7 kcal/mol obtained by the highly correlated coupled cluster calculations,48 suggesting a practicality of the present procedure. Finally, the energy difference between FMO and FF (Dreiding44) in configuration C is addressed (see again Table 4 and Figure 8). For the hexane−hexane pair, the magnitude of difference (0.55 kcal/mol) is considerably smaller than those of the two other pairs (1.86 and 2.09 kcal/mol) in which charge transfer contribution is substantial. In other words, the difference may be primarily attributable to the amount of

Table 3. Interaction Energies (Eij), Scaling Factors (Sij), and Coordination Numbers (Zij) for Hexane (1)−Nitrobenzene (2) System at 300 K, as Well as Calculated Critical Temperature (Tc) Tc (K) pair (i−j)

Eij (kcal/mol)

Sij

Zij

χc

exptl

FMO

1−1 1−2 2−1 2−2

−0.77

0.88

−1.78

0.83

2.0

293

286

−3.28

0.75

10.6 10.4 10.7 10.6

hexane pair (E11) is considerably smaller than that of the nitrobenzene−nitrobenzene pair (E22) and also that of the hexane−nitrobenzene pair (E12): the reasons for these differences will be addressed later. S11, S22, and S12 are 0.88, 0.75, and 0.83, respectively, indicating that the interaction between nitrobenzene molecules has some anisotropy. This situation may be attributed to the polarized nature due to the nitro group and the planar structure of a nitrobenzene unit. The calculated χ values of the hexane−nitrobenzene system are plotted against temperature in Figure 7, where a fitted curve is visible. For comparison, the results of GAFF43 and also of FMO without the scaling factor are included, indicating overestimations relative to the result of FMO with scaling. Since this system is of a solvent−solvent mixture, the degree of F


Article

The Journal of Physical Chemistry B Table 4. PIEDA Results for Hexane (1)−Nitrobenzene (2) Systema PIEDA pair (i−j) 1−1

1−2

2−2

A B C A B C A B C

IFIE

ES

EX

CT + mix

DI (MP2)

q (I → J)

Dreiding

diff between IFIE and Dreiding

−2.12 −0.19 −1.44 −4.12 −0.01 −2.99 −7.13 2.36 −4.68

−0.13 0.06 0.18 −0.71 0.16 −1.14 −2.09 2.84 −2.53

0.20 0.03 0.13 0.33 0.02 0.15 0.94 0.14 0.10

−0.38 −0.05 −0.29 −1.06 −0.05 −0.70 −1.20 −0.29 −0.81

−1.81 −0.23 −1.45 −2.68 −0.14 −1.31 −4.79 −0.32 −1.43

−0.0001 −0.0002 0.0001 −0.0076 0.0005 −0.0049 0.0013 0.0001 0.0070

−2.41 −0.29 −1.99 −3.09 −0.34 −1.13 −6.06 2.47 −2.59

0.29 0.10 0.55 1.04 0.33 1.86 1.07 0.10 2.09

a

ES, EX, CT + mix, and DI contributions to IFIE are shown in kcal/mol. q (Mulliken charge transfer) is given in e. Results by Dreiding FF44 are included as well. Labels A, B, and C are according to Figure 8.

300 K. These values mean relatively weak interactions and high isotropy in this system, relative to the above-discussed hexane− nitrobenzene system. The χ values against temperature are shown in Figure 9, and the Tc values depending on the

Figure 9. Calculated χ values of diisobutyl ketone−polyisobutylene system as a function of temperature. Symbols of square, cross, and triangle correspond to the results of FMO with scaling (S), FMO without scaling, and GAFF,43 respectively. Three dotted horizontal lines indicate χc shown in Table 6. Curve with squares (FMO with scaling) can be best fitted with χ = −0.39 + (3.138 × 102)/T.

Figure 8. Characteristic configurations of hexane−nitrobenzene system: (A) the most stable configuration; (B) the most unstable configuration; (C) the configuration which has the largest energy difference between the IFIE and FF values.

Table 6. Critical Temperatures of Diisobutyl Ketone (dib)− Polyisobutylene System

charge transfer (q), and this fact indicates the necessity of quantum mechanical treatment (just as FMO) in the evaluation of χ parameters. 4.2. Polyisobutylene−Diisobutyl Ketone (Polymer− Solvent) System. Table 5 shows the values of Eij, Sij, and Zij at

Tc (K)

Table 5. Interaction Energies (Eij), Scaling Factors (Sij), and Coordination Numbers (Zij) for Polyisobutylene (1)− Diisobutyl Ketone (2) System at 300 K pair (i−j)

Eij (kcal/mol)

Sij

Zij

1−1 1−2 2−1 2−2

−0.95 −1.13 −1.13 −1.38

0.92 0.92 0.92 0.88

7.6 6.9 9.7 9.1

Mw(dib)

χc

exptl

FMO

22 700 285 000 6 000 000

0.57 0.52 0.50

292 319 329

328 346 354

polymerization degrees (reflected in Mw) are listed in Table 6. The experimental Tc values in the three degrees are 292, 319, and 329 K,13 respectively, and they are compared with the present estimations of 328, 346, and 354 K in this order. The errors fall within about 10% range, indicating an acceptable reliability of our method. It is notable that the introduction of the scaling factor is effective. A high Z dependence was G


Article


our FMO estimations with scaling again show acceptable agreement (10% error) with the experimental values.49−51

observed in ref 20, and the best values of 295, 315, and 320 K were obtained with the Z value of 6.5: a similar situation was reported for the polyisoprene−polystyrene system. 4.3. Polyisoprene−Polystyrene (Polymer−Polymer) System. The calculated values of Eij, Sij, and Zij at 300 K are given in Table 7, and the χ values are plotted in Figure 10. The

5. SUMMARY In this paper, we reported a new method to estimate the χ parameters as an improvement with a quantum mechanical framework. Namely, the list of interaction energies between segments were evaluated by a series of FMO calculations.25−27 Some related improvements such as an introduction of anisotropy scaling were made as well. The present method was tested with the three simple systems of binary mixtures.20 The theoretical Tc values were in agreement with the experimental data within 10% errors for all systems, showing a practical reliability. Since our procedure could be free from FFs, several demanding systems (for which reasonable empirical parameters were not available) would be tractable. In fact, works to model the Nafion−water system as well as a lipid−water system52 had been carried out. Although several improvements are still necessary, we believe that the present FMO-based method of the χ parameter determination has promising potential with wide-range applicability to various mesoscopic systems under a nonempirical framework. Works to implement the present method as a workflow system (FCEWS, FMO-based Chi-parameter Evaluation System) are underway.

Table 7. Interaction Energies (Eij), Scaling Factors (Sij), and Coordination Numbers (Zij) for Polyisoprene (1)− Polystyrene (2) System at 300 K pair (i−j)

Eij (kcal/mol)

Sij

Zij

1−1 1−2 2−1 2−2

−1.16 −1.78 −1.78 −2.59

0.84 0.80 0.80 0.76

8.0 7.0 9.0 7.8

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Yuji Mochizuki: 0000-0002-7310-5183 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors would thank Dr. Kosuke Ohata (JSOL Corp.) for discussion and technical assistance and also Dr. Yuto Komeiji (AIST Japan) for fruitful comments on the manuscript. PIEDA35,38 was done with an MHIR local version of ABINIT-MP.27 This work was supported in part by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) as a social and scientific priority issue no. 6 (Accelerated Development of Innovative Clean Energy Systems) to be tackled by using post-K computer.

Figure 10. Calculated χ values of polyisoprene−polystyrene system as a function of temperature. Symbols of square, cross, and triangle correspond to the results of FMO with scaling (S), FMO without scaling, and GAFF,43 respectively. Three dotted horizontal lines indicate χc shown in Table 8. Curve with squares (FMO with scaling) can be best fitted with χ = 0.47 × 10−1 + 1.918 × 103/T2.

values of E12 and E22 are slightly larger in stabilization than those of the polyisobutylene−diisobutyl ketone system because of the existence of polarizable π electrons in the styrene unit. The calculated Tc values for the respective polymerization degrees are compiled in Table 8. One can see that Tc considerably varies depending on the Mw values. Fortunately,

■

Table 8. Critical Temperatures of Polyisoprene (pip)− Polystyrene (ps) System at 300 Ka Tc (K)

Mw pip

ps

χc

exptl

FMO

1000 2000 2700 2700

1000 2700 2100 2700

0.34 0.15 0.15 0.12

243 329 408 448

255 420 420 489

REFERENCES

(1) Fraaije, J. G. E. M.; van Vlimmeren, B. A. C.; Maurits, N. M.; Postma, M.; Evers, O. A.; Hoffmann, C.; Altevogt, P.; GoldbeckWood, G. The Dynamic Mean-Field Density Functional Method and Its Application to the Mesoscopic Dynamics of Quenched Block Copolymer Melts. J. Chem. Phys. 1997, 106, 4260−4269. (2) Hasegawa, R.; Doi, M. Dynamical Mean Field Calculation of Grafting Reaction of End-Functionalized Polymer. Macromolecules 1997, 30, 5490−5493. (3) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. Europhys. Lett. 1992, 19, 155−160. (4) Koelman, J. M. V. A.; Hoogerbrugge, P. J. Dynamic Simulations of Hard-Sphere Suspensions Under Steady Shear. EPL (Europhysics Lett. 1993, 21, 363. (5) Mackie, A. D.; Avalos, J. B.; Navas, V. Dissipative Particle Dynamics with Energy Conservation: Modelling of Heat Flow. Phys. Chem. Chem. Phys. 1999, 1, 2039−2049.

According to eq 11, χc has the same value (0.15) when Mw’s of pip− ps are 2000−2700 and 2700−2100. Therefore, Tc for FMO also has the same value. a

H


Article

The Journal of Physical Chemistry B (6) Groot, R. D.; Warren, P. B. Dissipative Particle Dynamics: Bridging the Gap between Atomistic and Mesoscopic Simulation. J. Chem. Phys. 1997, 107, 4423−4435. (7) Yamazaki, N.; Motoyama, M.; Nonomura, M.; Ohta, T. Morphology of Microphase Separated Domains in Rod−coil Copolymer Melts. J. Chem. Phys. 2004, 120, 3949−3956. (8) Wijmans, C. M.; Smit, B.; Groot, R. D. Phase Behavior of Monomeric Mixtures and Polymer Solutions with Soft Interaction Potentials. J. Chem. Phys. 2001, 114, 7644−7654. (9) Arai, N. Structural Analysis of Telechelic Polymer Solution Using Dissipative Particle Dynamics Simulations. Mol. Simul. 2015, 41, 996− 1001. (10) Yamamoto, S.; Hyodo, S. A Computer Simulation Study of the Mesoscopic Structure of the Polyelectrolyte Membrane Nafion. Polym. J. 2003, 35, 519−527. (11) Yamamoto, S.; Maruyama, Y.; Hyodo, S. Dissipative Particle Dynamics Study of Spontaneous Vesicle Formation of Amphiphilic Molecules. J. Chem. Phys. 2002, 116, 5842−5849. (12) Shillcock, J. C.; Lipowsky, R. Equilibrium Structure and Lateral Stress Distribution of Amphiphilic Bilayers from Dissipative Particle Dynamics Simulations. J. Chem. Phys. 2002, 117, 5048−5061. (13) Flory, P. J. Principles of Polymer Chemistry; Cornell University Press: Ithaca, NY, 1953. (14) Paul, D. R. Polymer Blends; Academic Press: New York, 1978. (15) Benesi, H. A.; Hildebrand, J. H. A Spectrophotometric Investigation of the Interaction of Iodine With Aromatic Hydrocarbons. J. Am. Chem. Soc. 1949, 71, 2703−2707. (16) Van Krevelen, D. W. Properties of Polymers: Their Estimation and Correlation with Chemical Structure; Elsevier: Amsterdam, 1976. (17) Coleman, M. M.; Serman, C. J.; Bhagwagar, D. E.; Painter, P. C. A Practical Guide to Polymer Miscibility. Polymer 1990, 31, 1187− 1203. (18) Bicerano, J. Prediction of Polymer Properties; CRC Press: Boca Raton, FL, 2002. (19) Prathab, B.; Subramanian, V.; Aminabhavi, T. M. Molecular Dynamics Simulations to Investigate Polymer−polymer and Polymer− metal Oxide Interactions. Polymer 2007, 48, 409−416. (20) Fan, C. F.; Olafson, B. D.; Blanco, M.; Hsu, S. L. Application of Molecular Simulation To Derive Phase Diagrams of Binary Mixtures. Macromolecules 1992, 25, 3667−3676. (21) Jawalkar, S. S.; Adoor, S. G.; Sairam, M.; Nadagouda, M. N.; Aminabhavi, T. M. Molecular Modeling on the Binary Blend Compatibility of Poly(vinyl Alcohol) and Poly(methyl Methacrylate): An Atomistic Simulation and Thermodynamic Approach. J. Phys. Chem. B 2005, 109, 15611−15620. (22) Sepehr, F.; Paddison, S. J. Dissipative Particle Dynamics Interaction Parameters from Ab Initio Calculations. Chem. Phys. Lett. 2016, 645, 20−26. (23) Fedorov, D. G.; Kitaura, K. The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems; CRC Press: Boca Raton, FL, 2009. (24) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087−1092. (25) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Fragment Molecular Orbital Method: An Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 1999, 313, 701−706. (26) Fedorov, D. G.; Nagata, T.; Kitaura, K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 2012, 14, 7562−7577. (27) Tanaka, S.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Fukuzawa, K. Electron-Correlated Fragment-Molecular-Orbital Calculations for Biomolecular and Nano Systems. Phys. Chem. Chem. Phys. 2014, 16, 10310−10344. (28) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 2012, 112, 632−672.

(29) Nakano, T.; Kaminuma, T.; Sato, T.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Application to Polypeptides. Chem. Phys. Lett. 2000, 318, 614−618. (30) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Fragment Molecular Orbital Method: Use of Approximate Electrostatic Potential. Chem. Phys. Lett. 2002, 351, 475− 480. (31) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. A New Hierarchical Parallelization Scheme: Generalized Distributed Data Interface (GDDI), and an Application to the Fragment Molecular Orbital Method (FMO). J. Comput. Chem. 2004, 25, 872−880. (32) Fedorov, D. G.; Kitaura, K. Second Order Møller-Plesset Perturbation Theory Based upon the Fragment Molecular Orbital Method. J. Chem. Phys. 2004, 121, 2483. (33) Mochizuki, Y.; Nakano, T.; Koikegami, S.; Tanimori, S.; Abe, Y.; Nagashima, U.; Kitaura, K. A Parallelized Integral-Direct SecondOrder Møller-Plesset Perturbation Theory Method with a Fragment Molecular Orbital Scheme. Theor. Chem. Acc. 2004, 112, 442−452. (34) Mochizuki, Y.; Koikegami, S.; Nakano, T.; Amari, S.; Kitaura, K. Large Scale MP2 Calculations with Fragment Molecular Orbital Scheme. Chem. Phys. Lett. 2004, 396, 473−479. (35) Fedorov, D. G.; Kitaura, K. Pair Interaction Energy Decomposition Analysis. J. Comput. Chem. 2007, 28, 222−237. (36) Amari, S.; Aizawa, M.; Zhang, J.; Fukuzawa, K.; Mochizuki, Y.; Iwasawa, Y.; Nakata, K.; Chuman, H.; Nakano, T. VISCANA: Visualized Cluster Analysis of Protein - Ligand Interaction Based on the Ab Initio Fragment Molecular Orbital Method for Virtual Ligand Screening. J. Chem. Inf. Model. 2006, 46, 221−230. (37) Mochizuki, Y.; Fukuzawa, K.; Kato, A.; Tanaka, S.; Kitaura, K.; Nakano, T. A Configuration Analysis for Fragment Interaction. Chem. Phys. Lett. 2005, 410, 247−253. (38) Tsukamoto, T.; Kato, K.; Kato, A.; Nakano, T.; Mochizuki, Y.; Fukuzawa, K. Implementation of Pair Interaction Energy Decomposition Analysis and Its Applications to Protein-Ligand Systems. J. Comput. Chem., Jpn. 2015, 14, 1−9. (39) See: http://www.j-octa.com/ (accessed Dec 1, 2017). (40) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (41) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (42) 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, revision B.01; Gaussian, Inc.: Wallingford, CT, 2002. (43) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (44) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. Dreiding: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897−8909. (45) Moore, W. J. Physical Chemistry, 5th ed.; Prentice-Hall: 1999. (46) Atkins, P.; de Paula, J. Atkins’ Physical Chemistry, 2nd ed.; Oxford University Press: 1982. (47) Battaglia, M. R.; Buckingham, A. D.; Williams, J. H. The Electric Quadrupole Moments of Benzene and Hexafluorobenzene. Chem. Phys. Lett. 1981, 78, 421−423. (48) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M. Intermolecular Interactions of Nitrobenzene-Benzene Complex and Nitrobenzene Dimer: Significant Stabilization of Slipped-Parallel Orientation by Dispersion Interaction. J. Chem. Phys. 2006, 125, 124304. (49) McIntyre, D.; Rounds, N.; Campos-Lopez, E. Thermodynamic and Structural Parameters in Tri-Block Polymers. Polym. Prepr. (Am. Chem. Soc., Div. Polym. Chem. 1969, 10, 531−537. (50) Koningsveld, R.; Kleintjens, L. A. Thermodynamics of Multicomponent Polymer Systems. Macromolecular Chemistry−8; I


Article

The Journal of Physical Chemistry B Elsevier: 1973; pp 197−230. DOI: 10.1016/B978-0-408-705165.50013-1. (51) Koningsveld, R.; Kleintjens, L. A.; Schoffeleers, H. M. Thermodynamic Aspects of Polymer Compatibility. Pure Appl. Chem. 1974, 39, 1−32. (52) Doi, H.; Okuwaki, K.; Mochizuki, Y.; Ozawa, T.; Yasuoka, K. Dissipative Particle Dynamics (DPD) Simulations with Fragment Molecular Orbital (FMO) Based Effective Parameters for 1-Palmitoyl2-Oleoyl Phosphatidyl Choline (POPC) Membrane. Chem. Phys. Lett. 2017, 684, 427−432. (53) Okuwaki, K.; Doi, H.; Mochizuki, Y. An automated framework to evaluate effective interaction parameters for dissipative particle dynamics simulations based on the fragment molecular orbital (FMO) method. J. Comput. Chem. Jpn. 2017, accepted.

■

NOTE ADDED IN PROOF A short paper about FCEWS has been available during the proof production see ref 53.

J


Fragment Molecular Orbital Based Parametrization Procedure for

Recommend Documents