Introduction

The solubility parameter, δ, concept provides a numerical estimate of the degree of interaction between materials, and can be a good indication of solubility, particularly for polymers [1–8]. As Bicerano emphasized in his book [8], the utilizations of polymers in many technological and industrial applications are critically dependent on their δ. The δ is often used in industry to predict the compatibility [1], permeation [1, 2], and swelling [2], bulk and solution properties [9–12] of polymers. In sensor applications δ can provide insights to predict the swelling changes of polymers in the presence of volatile chemical compounds [13, 14]. Furthermore, it is useful to optimize processing conditions, as well as select the suitable solvent or the optimum solvents combinations in coatings industry [2]. The δ for a pure compound is defined from Hildebrand-Scatchard solution theory as [3–7]

$$ \delta = \sqrt {{{E_{{CED}}}}}. $$
(1)

The δ is defined as the square root of the E CED. The E CED is the amount of energy needed to completely remove unit volume of molecules from their neighbors to infinite separation (an ideal gas), which is equal to the molar energy of vaporization (ΔH vap–RT) divided by molar volume (V) [4–6]. Equation 1 can then be converted to

$$ \delta = \sqrt {{\frac{{\Delta {H_{{vap}}} - RT}}{V}}}, $$
(2)

where ΔH is the molar enthalpy of vaporization, R is the gas constant, with all physical quantities refering to the same temperature T. The division of the Hildebrand parameter into three component Hansen parameters (dispersion force, polar force, and hydrogen bonding force) can be expressed in terms of the individual solubility parameters δ d, δ p and δ h [6, 8]

$$ \delta = \sqrt {{\delta_d^2 + \delta_h^2 + \delta_p^2}} $$
(3)

where δ d , δ h , and δ p are respectively representing Hansen parameter contribution from dispersion, hydrogen bond, and polar interactions [6, 8].

The δ of polymers is a thermophysical property [12–14], and it shows high dependency to the temperature, T [2, 11–16]. In fact, the T g of a polymer is closely related to its δ [8, 12, 15, 16]. By considering the T effects on the fractional free volume (v f) and the enthalpy (H) change, Eqs. 4 and 5 have been proposed to describe the relationship between δ and T respectively above and below T g for the polymers exhibiting a glass-rubber transition.

For a temperature above Tg, T ≥ T g, the polymer is at rubbery state [15, 16]

$$ \delta {(T)^2} = \left[ {{H_C}\left( {{T_g}} \right)\exp \left( { - 4{\alpha_l}\Delta T} \right) + R{T^2}\left( {\partial \ln {v_f}/\partial T} \right)} \right]/{\nu_l}(T)M, $$
(4)

where H c(T g) is a molar composite enthalpy at T g consisting of a conformational and an intermolecular interaction component, α l is the volume thermal expansion coefficient above T g, M is the molar mass, v l (T) is the specific volume above T g, and ΔT equals the temperature difference with reference to \( {T_{\text{g}}}\left( {\Delta T = T--{T_{\text{g}}}} \right). \) The v equals the reciprocal of its mass density (ρ); molar volume equals the molar mass divided by the mass density (ρ) (V = M/ρ). Thus, the molar volume can be expressed as V = vM. The molar mass of a polymer is a constant.

For a temperature below T g, T ≤ T g, a polymer at glassy state [15]

$$ \delta {(T)^2} = \left[ {{H_C}\left( {{T_g}} \right)\exp \left( { - 4{\alpha_l}\Delta T} \right) - M\int_T^{{{T_g}}} {\Delta {C_p}} (T)/dT} \right]/{\nu_s}(T)M, $$
(5)

where \( \Delta {C_{\text{p}}}(T) = {C_{\text{pl}}}(T)--{C_{\text{ps}}}(T) \), with C pl(T) and C ps(T) being the specific heat capacities at a constant pressure (p) of the polymer in the liquid and glassy states at T respectively, and v s (T) is the specific volume below T g.

Although theoretical works [8, 12, 15, 16] have greatly contributed to the understanding of the relationship between δ and T, a number of fundamental questions remain uncovered, for example, more experimental data on δ at different T are required to support Eqs. 4 and 5. However, determination of δ at different T from experimental methods is time consuming and labor intensive. Moreover, Eqs. 4 and 5 are too complex to be applied practically. The δ can be considered to be linear proportional to T separately in a temperature above and below T g [15, 16]. Thus, Eqs. 4 and 5 can be are respectively simplified to two linear functions: For T ≤ T g, when the polymer is at glassy state (solid phase)

$$ \delta (T) = {m_s}\left( {T - {T_g}} \right) + {\delta_g}, $$
(6)

where δ g is the δ at T g and m s is the thermal coefficient of δ at glassy state (solid phase); For T ≥ T g, when the polymer is at rubbery state

$$ \delta (T) = {m_l}\left( {T - {T_g}} \right) + {\delta_g}, $$
(7)

where m l is the thermal coefficient of δ at rubbery state (liquid phase).

Equations 6 and 7 show simple mathematical relations between δ and T. Given the m s (or m l ) and T g, δ at a given temperature can be calculated easily from the equations. The two equations and the material constants (m s, m l and T g,) is very important for the physical properties of the polymer. They can provide a fundamental understanding of temperature effects on the bulk and solution properties [1–8, 13, 14, 17] and also help in selecting or designing materials for specific applications. Moreover, they can be utilized in setting protocol for processes related to polymer–solvent interactions [1, 2, 15].

Literature [18–20] shows the possibility in calculating the δ of polymers via MD techniques. Additionally, simulations on bulk systems of moderate size (1000–5000 atoms) using recent systematically-derived “class II” forcefields such as “COMPASS” are capable of making predictions of δ with an accuracy that compares with experiment [21–25]. The forcefield describes approximately the potential energy hypersurface on which the atomic nuclei move. An accurate selection of forcefields is a key in enabling accurate predictions of the molecular interactions and material properties [17–25].

COMPASS is a powerful forcefield that supports atomistic simulations of condensed phase materials [23–25]. It has been parameterized and validated using condensed-phase properties in addition to various ab initio and empirical data for molecules in isolation [23, 24]. The forcefield enables accurate and simultaneous prediction of structural, conformational, vibrational, and thermophysical properties, that exist for organic molecules, inorganic small molecules, and polymers in isolation and in condensed phases, and under a wide range of conditions of temperature and pressure [23–25].

In this investigation, our objective is to model the temperature dependence characteristics of δ for emeraldine based polyaniline (EB–PANI) by using the COMPASS forcefield based MD techniques. The decrease of δ with increasing T for the polymer is plotted by running MD simulations at temperatures ranging from 250 up to 650 K and taking readings for δ at various temperatures. Due to the lack of published δ data at temperature other than 298 K, in order to verify the accuracy of our method, two kind of analysis has been are performed: (i) the δ value at 298 K for EB–PANI is predicted and compared with the literature reported data; (ii) the T g of the polymer is determined from the δ–T curve and compared with the experimental value. The temperature effects on the physical quantities (v, U coh and E CED) used in the definition of δ are also studied. The applications of the temperature dependence of δ on optimizing processing conditions will be addressed and discussed later in the discussion section.

Computational methodology

The use of atomistic methods to furnish a first principle calculation of E CED, is a desirable tool [18–25]. In atomistic simulations, the cohesive energy, U coh, is defined as the increase in energy per mole of a material if all intermolecular forces are eliminated [18, 19, 25]. The E CED corresponds to the cohesive energy per unit volume. It is a measure of the intermolecular forces within a system, and is estimated via the non-bonded van der Waals and electrostatic (includes hydrogen bond) interactions [18–20, 25, 26]. Equation 2 can be written in the form of:

$$ {E_{{CED}}} = \left( {{U_{{vdw}}} + {U_Q}} \right)/\nu M, $$
(8)

where U vdw and U Q are respectively van der Waals and electrostatic energy [8, 14, 18–20, 25]. Hence, the δ can be formally expressed in terms of

$$ \delta = \sqrt {{\delta_{{vdw}}^2 + \delta_Q^2}}, $$
(9)

where δ vdw and δ Q are respectively representing the contributions from van der Waals forces and electrostatic interactions and

$$ {\delta_Q} = \sqrt {{{E_Q}}} = \sqrt {{\frac{{{U_Q}}}{{\nu M}}}} $$
(10)
$$ {\delta_{{vdw}}} = \sqrt {{{E_{{vdw}}}}} = \sqrt {{\frac{{{U_{{vdw}}}}}{{\nu M}}}}, $$
(11)

where E Q is the electrostatic energy density and E vdw is the van der Waals energy density. Counting hydrogen bonding energy to electrostatic energy results in

$$ {\delta_Q} = \sqrt {{\delta_h^2 + \delta_p^2}}. $$
(12)

In the general case of a gas, liquid, or solid at a constant pressure, p, the volume thermal expansion coefficient, α, is given by

$$ \alpha = \frac{1}{{\nu M}}{\left( {\frac{{\partial \nu }}{{\partial T}}} \right)_p}. $$
(13)

Since the thermal expansion gradients of α for above and below the T g are different, and it can be used to predict the T g of a polymer [25–28]. Applying Eqs. 8–13 to a polymer, the thermal properties of solubility parameter and its related cohesive properties can be investigated and established by MD simulation.

Molecular modeling details

The polymer considered in this work is EB−PANI which consists of equal numbers of reduced (diamine) and oxidized (diimine) units: [(C6H4N)2(C6H4NH)2]n. It is regarded as one of the most useful conducting polymers due to its ease of synthesis, environmental stability and its semi-conductive characteristic upon protonic acid doping (Fig. 1) [25, 29–31]. Structurally, EB–PANI contains benzene rings and amine entity. The cohesive energy contributes from all types of interactions including the polar, the nonpolar (or dispersive), and the specific ones (or hydrogen bonding). The large amounts of benzene rings in the EB–PANI structure enhance the intermolecular interaction among the polymer chains [32–34]. These characteristics inherently make the physical properties of EB–PANI sensitive to temperature. EB−PANI is a good candidate for studying the thermal characteristics of δ and its related cohesive properties.

Fig. 1
figure 1

Non-oxidizing protonic acid (e.g., HCl, H2CO3) doping of EB–PANI: 2n is the number of protonic acid molecules, H+ is the proton, A− is an anion, e.g., chloride

Full size image

The molecular simulations are performed by the Forcite Plus module of the commercial software Materials Studio® 5.0 using COMPASS forcefield. Simulation approaches include MM and MD with canonical ensemble (NVT) and isothermal-isobaric ensemble (NPT). All MD simulations are run with a 1.0 fs time step. In all MD simulations, temperature and pressure were controlled by the Berendsen’s method using a half-life for decay to the target temperature of 0.1 ps (decay constant) and 0.1 ps for the pressure scaling constant. The non-bonded electrostatic and van der Waals forces are controlled by “atom based” summation and a cutoff value of 15.5 Å is used.

Building the amorphous polymer system

The EB−PANI polymer molecule is constructed using the polymer builder based on its stereoisomerism (tacticity) and sequence isomerism (connectivity). The tacticity of the repeat units (Fig. 2a) in the constructed polymer is isotactic. The connectivity of the monomer units (Fig. 2a) is head-to-tail. After geometry optimization, polymer chains are assembled into a three-dimensional (3D) unit cell subject to periodic boundary conditions by an “Amorphous Cell” module. The “Amorphous Cell” module provides a comprehensive set of tools to construct 3D periodic structures of molecular liquids and polymeric systems. The module builds molecules in a cell in a Monte Carlo fashion, by minimizing close contacts between atoms, whilst ensuring a realistic distribution of torsion angles for any given forcefield [35]. Using this algorithm, an amorphous polymer structure can be built with realistic conformations, while minimizing the number of close contacts [25, 35].

Fig. 2
figure 2

a Polymer monomer of EB–PANI with COMPASS forcefield assigned atomic charges. b Schematic illustration of a cubic unit cell for amorphous EB–PANI. Colors: carbon–gray, hydrogen–white, oxygen–red, nitrogen–blue, and hydrogen bond–red dotted line

Full size image

In order to suppress surface effects while keeping the number of particles in the model to a reasonable size, 3D periodic boundary condition (PBC), in which the system is considered to be surrounded on all sides by replicas of itself (Fig. 2b), has been employed. The length of the cubic unit cell is about 34 Å with 3688 atoms (4 polymer chains) (Fig. 2b). Each polymer chain in the amorphous model has 20 monomers and the polymer density is 1.23 g/cm3 [25]. The atomic charges are forcefield assigned [23–25] and hydrogen bonding interactions have been monitored as illustrated in Fig. 2b.

Geometry optimization and equilibration of polymer system

To create a final structure with realistic density and low-potential-energy characteristics, polymer system is subjected to geometry optimization and an equilibration cycles. The Forcite geometry optimization is employed to refine the geometry of polymer structure until it satisfies a certain criteria with the settings listed in Table 1. Equilibration of polymer system is conducted by an equilibrating protocol which consists of 13 intermediate MD stages (NVT annealings and NPT equilibration) as described in Table 1 [25]. The purpose of an annealing strategy in equilibration is to find the different local minimum energy structures by using the higher temperature periods of each annealing cycle to overcome energy barriers.

Table 1 Settings for geometry optimization and equilibration process of polymer system
Full size table

Glass-rubber transition of polymer

Once the cell structures have been optimized and equilibrated, NPT–MD simulations over a range of temperature above and below the T g are employed to model the glass-rubber transition of polymer. Experimentally, the T g value of EB−PANI range over 460 K interval [(473 K) [33] or (523 K) [34]]. A series of NPT–MD simulations over a temperature range from 250 up to 650 K with an interval of 25 K (each of 500 ps) are carried out for the prediction of T g . The simulation is continued, with the input of the final state from the previous simulation trajectory. With the simulation trajectory from the NPT–MD simulations, the specific volume, v, and cohesive energy, U coh, at various temperatures is calculated. The δ and its related physical quantities (E CED, E vdw, E Q), can be calculated by running the Forcite cohesive energy density calculation on the series of structures generated. The three components, van der Waals energy (δ vdw), hydrogen bonding energy (δ h), and polar interaction (δ p) contribute to the total solubility parameter, δ, are reported with their associated standard deviation.

Results and discussion

Validation of model in predicting solubility parameter

In order to validate the accuracy, the δ of EB−PANI at 298 K have been predicted and compared with the empirical data. Shacklette and Han determined the total solubility parameters of EB−PANI at 298 K by using empirical measurement and group additive based on Eq. 3. The value of each component is δ d = 17.40 (J/cm3)1/2, δ h = 10.70 (J/cm3)1/2, and δ p = 8.10 (J/cm3)1/2, which is equivalent to a total solubility parameter of δ = 22.2 (J/cm3)1/2 [36–39]. The solubility parameters of EB−PANI at 298 K are estimated in the current MD study as follows: δ vdw = 17.94 (J/cm3)1/2, δ Q = 14.07 (J/cm3)1/2, and the total solubility parameter of δ = 22.80 (J/cm3)1/2, which are in reasonable agreement with the reported values as listed in Table 2) [36–39]. The agreement between simulation results and empirical data indicates that it is possible to predict the δ at various temperatures by MD.

Table 2 Comparing the MD simulation predicted solubility parameters with experimental data (T = 298 K)
Full size table

Temperature dependence of specific volume and cohesive energy

Figure 3a shows that the specific volume as a function of temperature ranging from 250 to 650 K. There is a discontinuity at temperature 498 K, which represents a different thermal expansion in the rubbery and glassy states. Two best-fit lines have been drawn to represent the thermal expansion below and above T g. The T g is estimated as the point of intersection (approximately 498 K) between two fit lines. Again, it agrees well with the experimental data which is 473 K [33] and 523 K [34]. In addition to specific volume, the cohesive energy of polymer shows a relatively sudden change at the T g. To verify this argument, we plot the variation of cohesive energy against the temperature. As illustrated in Fig. 3b, an insignificant discontinuity is observed in the vicinity of T g. The T g as determined from the cohesive energy plot is 502 K, which is within the reported range from 473 K [33] to 523 K [34] in the literature. As also indicates from the other groups [27, 28], it demonstrates that T g can be estimated from change of physical entities along temperature obtained from MD studies.

Fig. 3
figure 3

Plots of specific volume and cohesive energy against temperature. There is a break in both plots of (a) specific volume vs. temperature and (b) cohesive energy vs. temperature. Experimental glass transition temperature of T g = 473 K [33] or 523 K [34]

Full size image

Temperature effects on the cohesive energy components related to solubility parameters

The thermal properties of total cohesive energy density, van der Waals energy, and electrostatic energy density used in defining δ, have been investigated. Figure 4a−c exhibit the various relationships between different physical entities and temperature (E CED vs. T, E vdw vs. T, and E Q vs. T). A noticeable break is detected from every plot. The corresponding temperature (493 K in Fig. 4a, 507 K in Fig. 4b, and 486 K in Fig. 4c) of these break are closed to the experimental T g data [33, 34]. These results indicate the cohesive energy density components versus temperature can also be used to determine T g. As formulated by Eq. 8, when an amorphous polymer undergoes a phase transition, the specific volume, cohesive energy, and cohesive energy density undergo a drastic change in the vicinity of the glass transition temperature, T g. The sudden jump at the glass transition temperature in these plots can be regarded as the evidence for Eqs. 4–7 from the previously theoretical works [8, 12, 15, 16].

Fig. 4
figure 4

Plots of cohesive energy density components related to solubility parameter against temperature. There is a discontinuous change in the vicinity of glass transition temperature, T g, in each plot of (a) cohesive energy density (b) van der Waals energy density (c) electrostatic energy density versus temperature. Experimental glass transition temperature of T g = 473 K [33] or 523 K [34]

Full size image

Temperature dependence of solubility parameters

Figure 5a shows the variation of the δ of EB−PANI as a function of T ranging from 250 to 650 K. An abrupt change is noticeable at a temperature of 495 K which indicates the occurrence of glass transition. The curve can be split into two linear sections in below and above the breaking point (T g). In each linear section, the solubility parameter decreases with increasing temperature. The corresponding temperature (495 K) of the breaking point agrees well with the T g value ranges (486 K–507 K) predicted from Fig. 3a–b and Fig. 4a–c, and it is in reasonable agreement with experiment T g ranges from 473 K [33] to 523 K [34]. Since the T g is an important characteristic for amorphous polymer, we believe the accuracy in getting T g from the calculated δ can justify the validation in estimating the δ with MD method.

Fig. 5
figure 5

a Plots of solubility parameter against temperature. There is a break at T = 495 K in the plot, which indicates the occurrence of glass-rubber transition. Experimental glass transition temperature of T g = 473 K [33] or 523 K [34]. b An example to determine the suitable processing temperature condition (T ptc) by using the δ–T curve. The EB–PANI has the best miscibility (or best solubility) with the target polymer (or solvent) at T ptc

Full size image

In applying Eqs. 6 and 7 to the current study, the thermal coefficient of solubility parameter at glassy state (m s) or liquid state (m l) and the solubility parameter at glass transition temperature (δ g) are calculated as:For T ≤ T g, the polymer at glassy state

$$ \delta (T){ } = { } - 0.00{45}\left( {T--{T_{\text{g}}}} \right){ } + { 22}.0{7}, $$
(14)

where the thermal coefficient of solubility parameter at solid phase m s = −0.0045 ((J/cm3)1/2  T −1) and the solubility parameter at T g , δ g = 22.07((J/cm3)1/2.

For T ≥ T g, the polymer at glassy state

$$ \delta (T){ } = { } - 0.00{66 }\left( {T--{T_{\text{g}}}} \right){ } + { 22}.0{7}, $$
(15)

where the thermal coefficient of solubility parameter at liquid phase m l = −0.0066 ((J/cm3)1/2  T −1). Equation 14 and Eq. 15 are useful for EB−PANI in the technical and industrial applications. For instance, it can be used to select a second polymer in formulating a polymer composite; incorporating a second component into conducting polymer film is one of the most important methods to develop new sensors [13, 14]. In comparison with modification of molecular structure of conducting polymers, the advantage of this technique is that it can avoid complicated chemical syntheses. Using Eq. 14, we can determine the miscibility of two polymers and the suitable processing temperature condition (T ptc) in fabrication as illustrated in Fig. 5b. In another scenario, the equations can help in the optimization of processing conditions, as well as selecting appropriate solvents or predicting the suitable temperature (see Fig. 5b) for combinations in coatings industry [2].

Conclusions

The molecular modeling methodology, which is capable of investigating the temperature dependence characteristics of δ of amorphous polymer, has been formulated. The temperature effects on δ are studied by running MD simulations at various temperatures. From the δ–T curve, the glassy phase and the rubbery phase are clearly seen. The MD predicted results over the significant range (from 250 up to 650 K) of temperatures are conveniently fitted to the two linear equations (Eqs. 5 and 6) of the form. Comparisons with experimental data confirm the validity of our method. The δ, v, U coh, and E CED show a similar temperature dependence characteristics and a drastic change around the T g. Thus, the temperature dependence of δ and its related cohesive properties (v, U coh, and E CED) may be utilized in determining T g of a polymer which has the glass-rubber transition. Furthermore, the applications of the temperature dependence characteristics of δ have been addressed and discussed.