A particularly attractive method to predict the dielectric properties of materials is density functional theory (DFT). While this method is very popular, its large computational requirements allow practical treatments of unit cells with just a small number of atoms in an ordered array, i.e., in a crystalline morphology. By comparing DFT and Molecular Dynamics (MD) simulations on the same ordered arrays of functional polyolefins, we confirm that both methodologies yield identical estimates for the dipole moments and hence the ionic component of the dielectric storage modulus. Additionally, MD simulations of more realistic semi-crystalline morphologies yield estimates for this polar contribution that are in good agreement with the limited experiments in this field. However, these predictions are up to 10 times larger than those for pure crystalline simulations. Here, we show that the constraints provided by the surrounding chains significantly impede dipolar relaxations in the crystalline regions, whereas amorphous chains must sample all configurations to attain their fully isotropic spatial distributions. These results, which suggest that the amorphous phase is the dominant player in the context, argue strongly that the proper polymer morphology needs to be modeled to ensure accurate estimates of the ionic component of the dielectric constant.

Polymeric materials are finding increased use in capacitor applications, and there has been significant motivation to predict the dielectric properties of new classes of polymers quantitatively without having to synthesize them. The addition of a small number of polar hydroxyl (–OH) groups is known to drastically increase the static relative permittivity (dielectric constant) of a nonpolar polymer, such as isotactic-polypropylene, in its semi-crystalline state.1–3 Our recent analysis of this system used Molecular Dynamics (MD) simulations coupled with Density Functional Theory (DFT) calculations to show that hydrogen bonds across spatially proximal –OH groups, combined with the presence of water, results in this higher dielectric constant, while at the same time not deleteriously affecting the dielectric loss.3 The critical role of hydrogen bonding and polarity has encouraged us to explore various polar side groups in our current work. A central assumption in our past work is that the morphology considered in the simulations is not critical to the prediction of the dielectric constant. This conclusion follows since the MD results (which deal with a semi-crystalline morphology, with both crystalline and amorphous components) are in good (qualitative) agreement with the DFT results, which only deals with the crystalline component.

Here, we focus on critically examining this central assumption that underpins much of the current understanding and predictive ability of dielectric properties. In particular, we predict the dielectric constant of polyethylene (PE) chains functionalized with a fixed amount (4.2 mol. %) of different polar groups in a series of DFT and MD simulations: specifically –CH3, –NH2, –NO2, –COOH, –OH, and –SH. We determined the crystal geometry for the functionalized PE using DFT. The crystal geometry obtained from DFT allowed us to calculate the dielectric constant using density functional perturbation theory (DFPT). Thereafter, we performed two flavors of MD simulations—one in the ground state crystal structure simulated by the DFT and the other in the self-assembled semi-crystalline state spontaneously formed by these polymers when they are cooled from the melt to a temperature below their equilibrium melting points. The ionic dielectric constants found in our first set of simulations are in quantitative agreement with the DFT calculations for the functionalized PE. This close agreement validates both the MD and the DFT protocols. However, the MD simulations using a semi-crystalline state yield ionic contributions to the dielectric constants that are up to 10 times larger than those obtained from crystalline simulations. An equilibrium statistical thermodynamical argument suggests that this reduced dielectric constant in the DFT calculations is a direct consequence of the fact that dipole moment fluctuations are significantly reduced in the crystalline state, while they are maximized in the amorphous state. Our results, therefore, emphasize the crucial role of polymer morphology in accurately predicting the dielectric constants of this commercially important class of materials.

We determine the relaxed crystal geometries of the functionalized PE using density functional theory (DFT). We started from the stable PE crystal structure with two chains stacked next to each other in a unit cell. In one of the two chains, every 4th –CH2 unit has a H atom replaced by a side group X, giving rise to the PE–X crystal (Figure 1). Note that the 25 mol. % functionalization simulated here is much higher than typical experimental systems with 0.8–8 mol. % functional groups. DFT,4,5 as implemented in the Vienna ab initio software package (VASP),6 is applied on the starting polymer geometry. The DFT relaxation is performed using the rPW86 functional in which the DFT-DF2 vdW correction is applied7 to capture the van der Waals interactions in the polymer correctly.8 The projector-augmented wave (PAW) pseudopotentials9 were applied, along with a tight energy convergence criterion of 10−8 eV and an energy cut-off of 600 eV. The relaxed polymer crystal structure arrangement obtained is used as an input for a subsequent density functional perturbation theory (DFPT)10,11 calculation. DFPT provides us with the dielectric constant tensor that includes the electronic component12 as well as the ionic (lattice) component.13 The dielectric constant values were obtained by determining the trace of the respective dielectric tensors. As expected, the calculated electronic contribution is found to be similar for all the polar groups. However, the ionic dielectric constant (εionic) is found to be different for the various functionalities (Figure 2, Table I).

FIG. 1.

PE–X crystal structure with two chains stacked next to each other in a unit cell used for DFT.

FIG. 1.

PE–X crystal structure with two chains stacked next to each other in a unit cell used for DFT.

Close modal
FIG. 2.

The ionic dielectric constant of PE–X systems calculated using DFT represented by the red histogram. The blue histogram represents the ionic dielectric constant calculated using MD for the crystal structures obtained from DFT. The green histogram represents the ionic dielectric constant calculated from MD for semi-crystalline systems.

FIG. 2.

The ionic dielectric constant of PE–X systems calculated using DFT represented by the red histogram. The blue histogram represents the ionic dielectric constant calculated using MD for the crystal structures obtained from DFT. The green histogram represents the ionic dielectric constant calculated from MD for semi-crystalline systems.

Close modal
TABLE I.

The dielectric constants calculated for different systems.

DFTelectronic DFTionic MDDFT MD
–CH3  2.45  0.02  0.02 ± 0.01  0.02 ± 0.01 
–OH  2.47  0.17  0.19 ± 0.01  1.02 ± 0.02 
–NH2  2.49  0.07  0.05 ± 0.01  0.98 ± 0.01 
–NO2  2.39  0.26  0.36 ± 0.02  0.79 ± 0.03 
–COOH  2.43  0.26  0.45 ± 0.03  0.48 ± 0.04 
–SH  2.55  0.07  0.05 ± 0.02  0.69 ± 0.02 
DFTelectronic DFTionic MDDFT MD
–CH3  2.45  0.02  0.02 ± 0.01  0.02 ± 0.01 
–OH  2.47  0.17  0.19 ± 0.01  1.02 ± 0.02 
–NH2  2.49  0.07  0.05 ± 0.01  0.98 ± 0.01 
–NO2  2.39  0.26  0.36 ± 0.02  0.79 ± 0.03 
–COOH  2.43  0.26  0.45 ± 0.03  0.48 ± 0.04 
–SH  2.55  0.07  0.05 ± 0.02  0.69 ± 0.02 

We use the DFT generated structures as the initial configuration for an MD simulation using the all atom OPLS-AA force field.14 To speed up the simulations, we exploit general-purpose graphical processing units (GPGPUs) to accelerate the van der Waals and long-range Coulombic calculations,15,16 as implemented in LAMMPS.17 Two infinite chains of 500 carbons each were prepared using the geometry obtained from DFT. The PE–X systems obtained were then allowed to equilibrate in the canonical (NVT) ensemble for 5 ns. The systems were then run for 50 ns to calculate the ionic dielectric constant (εionic) using ε ionic = 4 π M 2 M 2 3 V k B T , where M is the dipole moment of the system at a time step, V is the volume, T is the temperature, and kB is the Boltzmann’s constant. The ionic dielectric constants obtained from DFT and MD simulations for these purely crystalline samples agree with each other (Table I, Figure 2). However, the dielectric constant observed in experiments for PE–OH is found to be significantly larger than these predicted values.1 We predict an ionic contribution of ≈0.2, which makes the total dielectric constant ≈2.65 while the experiments yield values of more than 3.5. Similar results are also found for the PE–NH2 functional group (Figure 3).

FIG. 3.

Dielectric constant of PP–NH2 with 1.0, 1.4, 2.2, and 3.0 mol. % –NH2.

FIG. 3.

Dielectric constant of PP–NH2 with 1.0, 1.4, 2.2, and 3.0 mol. % –NH2.

Close modal

We postulate that this discrepancy arises because these DFT-inspired simulations do not consider the amorphous phase and the crystal-amorphous regions. To confirm this assumption, we studied the PE functionalized with 25 mol. % –OH groups further. The system is heated at 1000 K for 10 ns under NPT conditions, cooled to 500 K at a rate of 50 K/ns, and then allowed to equilibrate at 500 K for another 5 ns. Under these conditions, the system is molten. Subsequently, this melt is cooled to 300 K at a rate of 2 K/ns. At 300 K, the volume of the system is again equilibrated for 20 ns. The crystallinity of the final structure is then measured by calculating the local order parameter (LOP).18–23 Here, each carbon is assigned a bond orientation unit vector, which is calculated by connecting the midpoints of its two adjacent backbone bonds. The bond order parameter between the ith and the jth atom is given by A = 3 cos 2 ϕ 1 2 = 3 b i b j 2 1 2 where b i and b j are unit orientation vectors of the respective atoms. The order parameter for a carbon site is determined by averaging the bond order parameters for all carbons within a radius of 0.7 nm. Any carbon with a LOP of more than 0.80 is considered crystalline. The crystallinity of this system was found to be 21% ± 3%, which is significantly lower than the crystallinity for previously studied PE–OH systems (≈60%-70%) which admittedly had lower levels of functionality (maximum of 8.2 mol. %).3 The system is then run for 50 ns under NVT conditions to calculate the dielectric constant. The calculated εionic = 1.96 ± 0.03, which is approximately 10 times higher than the result from the fully crystalline structures, validates our hypothesis that the morphology of the materials is crucial for the accurate prediction of the dielectric constant of this class of polymers.

We now simulate semi-crystalline systems that are more relevant to experiment— with 4.2 mol. % functional groups. To generate PE–X we create a single polyethylene chain with 1000 backbone carbon atoms. The chain is first equilibrated at 500 K and 1 bar, where it is in a molten, amorphous state. We then randomly replace 4.2 mol. % H-atoms of the CH2 with the appropriate polar group to generate the amorphous PE–X. The PE–X systems were continuously cooled from 500 K to 300 K at a rate of 2.5 K/ns under isobaric conditions. The system’s volume is then allowed to stabilize at 300 K for 20 ns. Thereafter, it is equilibrated in the canonical (NVT) ensemble for 10 ns, and the net dipole moment is averaged for a further 50 ns.

Morphology analysis of the systems at 300 K suggests that the polar groups are predominantly localized in non-crystalline regions. More specifically, we quantify the distribution of polar groups by measuring the degree of crystallinity using the LOP discussed above. A carbon with a LOP less than 0.25 is considered amorphous; more than 0.80 is considered crystalline; everything in between is considered interphase. We observe for all the PE–X systems, the polar groups prefer to stay in non-crystalline regions (Figure 4), which agrees with Flory’s notion that these “defect” groups are rejected from the crystal.24 

FIG. 4.

The main figure shows the distribution of polar groups in crystalline, interphase, and amorphous regions. The inset shows the crystallinity in the system.

FIG. 4.

The main figure shows the distribution of polar groups in crystalline, interphase, and amorphous regions. The inset shows the crystallinity in the system.

Close modal

The inset of Figure 4 shows that the crystallinity using LOP of all the PE–X is comparable (≈60%-80%) except for PE–COOH. We compare the crystallinity of these PE–X systems with pure PE by crystallizing PE using the same protocols. The crystallinity of the PE system is calculated to be 79% ± 2% using LOP. This appears to not compare favorably with the crystallinity values obtained from differential scanning calorimetry, namely 48% and 46% for pure PE and 1.3 mol. % PE–OH, respectively.3 Clearly, a major source of discrepancy is that the theory and the experiments determine crystallinity in different ways. To obtain an experimentally relevant crystallinity, we compare the ratio of the heat of fusion as obtained from the simulations, ΔHf, with that of a perfect PE crystal (68.4 cal/g).25 To obtain the heat of fusion, we calculated the enthalpy difference between the semicrystalline sample and a corresponding fully amorphous sample at the same temperature (300 K). For our simulated PE system, ΔHf = 536.26 ± 24.63 kcal/mol, which is equivalent to 38.30 ± 1.74 cal/g. Thus, the crystallinity for the simulated PE system is found to be 56% ± 3%, which is within 20% of the experimental value.

Since hydrogen bonding plays a critical role, we analyze the systems for cluster formations. To quantify any pairing of the functional groups, we measure the distance between the polar groups. If the distance between them is less than 3.3 Å, they are considered to be a part of a cluster. The average number of clusters is then reported in Table II. PE–CH3 and PE–NO2, where the functional group does not hydrogen bond, are found to have low clustering. However, PE–NH2, PE–COOH, and PE–OH form clusters of the polar group.

TABLE II.

Percentage of groups in a cluster.

–CH3 –COOH –NH2 –NO2 –OH –SH
14% ± 5%  76% ± 2%  28% ± 3%  7% ± 5%  57% ± 2%  11% ± 4% 
–CH3 –COOH –NH2 –NO2 –OH –SH
14% ± 5%  76% ± 2%  28% ± 3%  7% ± 5%  57% ± 2%  11% ± 4% 

The εionic calculated from these semi-crystalline samples is shown in Figure 2. The dielectric constant from these semi-crystalline samples is found to be significantly larger than those estimated by DFT on the purely crystalline systems at a higher polar concentration in all the systems, except in the case of –COOH. The most significant difference is seen in PE–NH2 system where the εionic is calculated to be 20 times larger than for the semi-crystalline system as compared to the purely crystalline system. To understand the discrepancy in the results, we synthesized PP–NH226 and measured its dielectric constant. Figure 3 shows dielectric constants of several PP–NH2 containing 0.4, 1, 1.4, 2.2, and 3.0 mol. % –NH2, respectively, over a temperature range from −20 to 100 °C and a frequency range from 100 to 1 MHz. We find that our prediction of the dielectric constant using the semi-crystalline structure is much closer to the experimentally measured values. Note that our simulations do not contain water. As shown in our previous work,3 “bound water” can lead to an additional increase in the dielectric constant of the material. Regardless, we conclude that the morphology is crucial for an accurate quantitative prediction of the dielectric constant for polar polymers relevant to capacitor applications.

In order to understand the variation of dielectric constant with these different functional groups we are inspired by its definition for amorphous systems, M = 0 : ε ionic M 2 V . We choose to calculate the microscopic analog of this quantity, and thus we enumerate the dipole moment (μ) of each polar group (CH2–X) from the trajectory of the simulations. The volume of each polar group (CH2–X) is calculated by first assuming all atoms to be hard spheres. Using the known centers of mass, we find the minimum distance of the atoms of each element to the other atoms. The radius of the element is initially defined as half of the distance between the element and the closest atom. The assumed radius allows us to predict the volume occupied by the atoms within the simulation box, and subsequently the void volume of the system. We then change the radius of each element iteratively to minimize the void volume of the simulation box. Thus, the radius of the element in the polar group is estimated, and subsequently, the volume of the polar group. Figure 5 shows a high correlation between the calculated ionic dielectric constant with the ratio of the square of the dipole moment of each polar group relative to the calculated molar volume. One exception is the case of –NO2, which shows some discrepancy from the linear trend. Note also that we use the quantities relevant to a –COOH dimer due to its known propensity to dimerize. Therefore, we can predict the dielectric constant by calculating the “molecular” dipole moment and the “density” of these functional groups, an idea that also explains the results of Dong et al.27 

FIG. 5.

Comparison of the ionic dielectric constant to the square dipole moment by volume for various polar groups.

FIG. 5.

Comparison of the ionic dielectric constant to the square dipole moment by volume for various polar groups.

Close modal

We now provide a mathematical understanding of the role of morphology in the prediction of the dielectric constant. We go back to the definition: ε ionic = 4 π M 2 M 2 3 V k B T , and note that the significant quantities are M 2 and M , where M = i μ i , i.e., the sum of the individual dipoles in the system. For simplicity let us consider one chain, and make the simplified assumption that local dipoles are aligned along the bond direction. In this situation, the net dipole of the chain is equivalent, within a multiplicative constant, to the end-to-end vector of the chain. For a purely Gaussian chain, 〈M〉 = 0 and 〈M2〉 is equivalent to the mean square end-to-end distance of the chain. The statistics of a single chain under the action of a stretching force (as derived by Rubinstein and Colby28), clearly shows that the variance, M 2 M 2 , assumes its maximum value for a Gaussian chain (or alternatively in an amorphous phase) and decreases monotonically to zero at full extension (i.e., in a completely ordered crystalline domain). By analogy, for small enough variations about the amorphous state the εionic will be unaffected. Now consider a large enough perturbation where all the chains are stretched fully along the z direction but retain complete freedom along the x and y directions. In this case, it is easy to show that the εionic will be 2/3 of its value in the fully amorphous phase. The fact that the crystalline DFT results are up to an order of magnitude smaller must stem from the limited x and y mobility (“libration”) of the functional groups in the crystal.

Our results clearly show that simulating the right morphology, especially one where the dipoles in a sample can relax in a manner that mimics experimental reality, is necessary to quantitatively predict the dipole moment of the material. Thus, simulating a purely crystalline phase, where all the dipoles are strongly constrained, will not provide a good representation of the real, semi-crystalline polymeric material of relevance to experiment. In recent works, it has been postulated that large increases in the ionic dielectric constant can be achieved by using defective crystals with increased free volume.27,29 While this last idea is indeed appealing, we note that such fringed micelle type models for polymer crystals were the focus of debate over 50 years ago. Unequivocal experimental evidence suggests that nature avoids such structures and predominantly favors the formation of a densely packed crystal surrounded by amorphous phases with lower densities. At this time, we believe that the existence of a polymer in an amorphous phase, either by itself or as part of a semi-crystalline entity, maximizes the dielectric constant for this class of materials.

The authors acknowledge financial support for this work through a Multidisciplinary University Research Initiative (MURI) Grant from the Office of Naval Research under Contract No. N00014-10-1-0944. The simulations presented in this paper were performed on the Yeti HPC cluster system maintained by the Columbia University Information Technology Center.

1.
T. C. M.
Chung
, “
Functionalization of polypropylene with high dielectric properties: Applications in electric energy storage
,”
Green Sustainable Chem.
02
,
29
37
(
2012
).
2.
X.
Yuan
,
Y.
Matsuyama
, and
T. C. M.
Chung
, “
Synthesis of functionalized isotactic polypropylene dielectrics for electric energy storage applications
,”
Macromolecules
43
,
4011
4015
(
2010
).
3.
M.
Misra
,
M.
Agarwal
,
D. W.
Sinkovits
,
S. K.
Kumar
,
C.
Wang
,
G.
Pilania
,
R.
Ramprasad
,
R. A.
Weiss
,
X.
Yuan
, and
T. C. M.
Chung
, “
Enhanced polymeric dielectrics through incorporation of hydroxyl groups
,”
Macromolecules
47
,
1122
1129
(
2014
).
4.
P.
Hohenberg
and
W.
Kohn
, “
Inhomogeneous electron gas
,”
Phys. Rev. B
7
,
1912
1919
(
1973
).
5.
W.
Kohn
and
L. J.
Sham
, “
Self-consistent equations including exchange and correlation effects
,”
Phys. Rev.
140
,
1133
1138
(
1965
).
6.
G.
Kresse
and
J.
Hafner
, “
Ab initio molecular dynamics for liquid metals
,”
Phys. Rev. B
47
,
558
561
(
1993
).
7.
J.
Klimeš
,
D. R.
Bowler
, and
A.
Michaelides
, “
Chemical accuracy for the van der Waals density functional
,”
J. Phys.: Condens. Matter
22
,
022201
(
2010
).
8.
C. S.
Liu
,
G.
Pilania
,
C.
Wang
, and
R.
Ramprasad
, “
How critical are the van der Waals interactions in polymer crystals?
,”
J. Phys. Chem. A
116
,
9347
9352
(
2012
).
9.
P.
Blöchl
, “
Projector augmented-wave method
,”
Phys. Rev. B
50
,
17953
17979
(
1994
).
10.
S.
Baroni
,
S.
De Gironcoli
,
A.
Dal Corso
, and
P.
Giannozzi
, “
Phonons and related crystal properties from density-functional perturbation theory
,”
Rev. Mod. Phys.
73
,
515
562
(
2001
).
11.
G.
Pilania
and
R.
Ramprasad
, “
Dielectric permittivity of ultrathin PbTiO3 nanowires from first principles
,”
J. Mater. Sci.
47
,
7580
7586
(
2012
).
12.
F.
Bernardini
,
V.
Fiorentini
, and
D.
Vanderbilt
, “
Spontaneous polarization and piezoelectric constants of III-V nitrides
,”
Phys. Rev. B
56
,
10024
10027
(
1997
).
13.
X.
Zhao
and
D.
Vanderbilt
, “
First-principles study of structural, vibrational, and lattice dielectric properties of hafnium oxide
,”
Phys. Rev. B
65
,
233106
(
2002
).
14.
W. L.
Jorgensen
,
D. S.
Maxwell
, and
J.
Tirado-Rives
, “
Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids
,”
J. Am. Chem. Soc.
118
,
11225
11236
(
1996
).
15.
J. E.
Stone
,
J. C.
Phillips
,
P. L.
Freddolino
,
D. J.
Hardy
,
L. G.
Trabuco
, and
K.
Schulten
, “
Accelerating molecular modeling applications with graphics processors
,”
J. Comput. Chem.
28
,
2618
2640
(
2007
).
16.
J. A.
Anderson
,
C. D.
Lorenz
, and
A.
Travesset
, “
General purpose molecular dynamics simulations fully implemented on graphics processing units
,”
J. Comput. Phys.
227
,
5342
5359
(
2008
).
17.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
18.
R.
Haggenmueller
,
J. E.
Fischer
,
K. I.
Winey
,
V.
Pennsyl
,
R. V.
December
,
V.
Re
,
M.
Recei
, and
V.
February
, “
Single wall carbon nanotube/polyethylene nanocomposites: Nucleating and templating polyethylene crystallites
,”
Macromolecules
39
,
2964
2971
(
2006
).
19.
J.
Ramos
,
J.
Mart
, and
S.
Sanmart
, “
Following the crystallization process of polyethylene single chain by molecular dynamics: The role of lateral chain defects
,”
Macromol. Symp.
312
,
97
107
(
2012
).
20.
M. S.
Lavine
,
N.
Waheed
, and
G. C.
Rutledge
, “
Molecular dynamics simulation of orientation and crystallization of polyethylene during uniaxial extension
,”
Polymer
44
,
1771
1779
(
2003
).
21.
S.
Fujiwara
and
T.
Sato
, “
Molecular dynamics simulations of structural formation of a single polymer chain: Bond-orientational order and conformational defects
,”
J. Chem. Phys.
107
,
613
(
2008
).
22.
W.
Simiao
,
Y. U.
Xiang
,
K.
Bin
, and
Y.
Xiaozhen
, “
Molecular dynamics study of the nucleation behavior and lamellar mergence in polyethylene globule crystallization
,”
Sci. China Chem.
56
,
195
202
(
2013
).
23.
V.
Tomer
,
G.
Polizos
,
C. A.
Randall
, and
E.
Manias
, “
Polyethylene nanocomposite dielectrics: Implications of nanofiller orientation on high field properties and energy storage
,”
J. Appl. Phys.
109
,
074113
(
2011
).
24.
P. J.
Flory
,
Principles of Polymer Chemistry
,
Baker Lectures 1948
(
Cornell University Press
,
1953
).
25.
B.
Wunderlich
and
C.
Cormier
, “
Heat of fusion of polyethylene
,”
J. Polym. Sci., Part A-2
5
,
987
988
(
1967
).
26.
M.
Zhang
,
X.
Yuan
,
L.
Wang
,
T. C. M.
Chung
,
T.
Huang
, and
W.
Degroot
, “
Synthesis and characterization of well-controlled isotactic polypropylene ionomers containing ammonium ion groups
,”
Macromolecules
47
,
571
581
(
2014
).
27.
R.
Dong
,
V.
Ranjan
,
M.
Buongiorno Nardelli
, and
J.
Bernholc
, “
Atomistic simulations of aromatic polyurea and polyamide for capacitive energy storage
,”
Phys. Rev. B
92
,
024203
(
2015
).
28.
M.
Rubinstein
and
R.
Colby
,
Polymer Physics
(
OUP Oxford
,
2003
).
29.
Y.
Thakur
,
R.
Dong
,
M.
Lin
,
S.
Wu
,
Z.
Cheng
,
Y.
Hou
,
J.
Bernholc
, and
Q.
Zhang
, “
Optimizing nanostructure to achieve high dielectric response with low loss in strongly dipolar polymers
,”
Nano Energy
16
,
227
234
(
2015
).