We present results from a computational framework integrating genetic algorithm and molecular dynamics simulations to systematically design isotope engineered graphene structures for reduced thermal conductivity. In addition to the effect of mass disorder, our results reveal the importance of atomic distribution on thermal conductivity for the same isotopic concentration. Distinct groups of isotope-substituted graphene sheets are identified based on the atomic composition and distribution. Our results show that in structures with equiatomic compositions, the enhanced scattering by lattice vibrations results in lower thermal conductivities due to the absence of isotopic clusters.

The rise of graphene since its discovery in 2004 and subsequent breakthroughs in 2-dimensional (2D) analogous nanostructures (such as silicene, phosphorene, etc.) have captured the attention of the scientific community over the last decade.1–3 The unique structural and transport properties of the 2D carbon nanomaterial, especially those associated with charge and heat transfer, have been addressed extensively via a suite of experimental, theoretical, and computational techniques. Likewise, the effect of defects in these materials is of notable interest and importance, primarily because the presence of mass disorder and perturbations to electronic structure can impede transport mechanisms.4–9 Several reports have conclusively shown that the remarkably high thermal conductivity of graphene composed purely of 12C is reduced by a factor of two for a nanostructure containing 50% 13C isotopes.9–11 Hence, the occurrences of defects in graphene, and in any nanomaterial for that matter, is of concern due to their variety (point defects and impurities could range from vacancies, isotopes, dopants, and functional groups) and distribution (random, clustered, superlattice type).12–17 A controlled synthesis process to create these materials is not sufficient to eradicate the presence of defects that could be naturally occurring, as is the case for isotopically impure nanomaterials. On the other hand, such strong deviations in material characteristics due to simply a mass disorder paves the way for creating defect engineered structures that possess designer targeted properties. One such example is the potential application of graphene for analog circuits and sensors in high-frequency communications requiring low level of flicker noise, which decreases with increasing defect concentration in the nanomaterial.18,19

The key design parameters for these nanomaterials are the type and fraction of defects, while the inconspicuous and often ignored factor is their distribution. Results of molecular dynamics (MD) simulations of thermal conductivity (λ) for 750 samples of graphene nanosheets containing 50% 14C isotope substitution distributed randomly are presented in Fig. 1. We find that while the predictions reflect a normal distribution with majority of the isotopically impure structures converging to a mean thermal conductivity, there are potential distributions of isotopes that result in lower and higher λ values. This observation illustrates the need to consider distribution in the defect engineered nanomaterial design. We do acknowledge that achieving deterministic designer choices for both defect concentration and distribution through advanced synthesis methods is still an experimental challenge, addressed sparsely in the literature. Although computationally doable, the sheer number of potential cases that need to be simulated to exhaustively inspect every possible combination of defect fraction and distribution makes such a brute force approach unfeasible, especially for system sizes ranging from few hundreds to thousands of atoms. Hence, a systematic procedure is required to narrow the search space and provide an efficient framework for designing defect engineered materials. Here, we employ a genetic algorithm (GA) technique integrated with MD simulations to recursively explore generations of nanostructure geometries and offer recommendations about suitable defect distributions as a design parameter.20 We choose the well-studied system of 50% isotope substituted graphene amongst other defect fractions, and seek to obtain 14C distributions that provide lowest thermal conductivities.

FIG. 1.

Thermal conductivity distribution of graphene sheets containing 50% 14C isotopes. Molecular dynamics predictions of thermal conductivity of randomly isotope substituted graphene structures shows that even for the same defect concentration, the transport coefficient can vary across a significant range, demonstrating the importance of atomic distributions in addition to mass fraction. A pictorial (not to scale) representation of the isotope substituted carbon nanomaterial is shown as inset.

FIG. 1.

Thermal conductivity distribution of graphene sheets containing 50% 14C isotopes. Molecular dynamics predictions of thermal conductivity of randomly isotope substituted graphene structures shows that even for the same defect concentration, the transport coefficient can vary across a significant range, demonstrating the importance of atomic distributions in addition to mass fraction. A pictorial (not to scale) representation of the isotope substituted carbon nanomaterial is shown as inset.

Close modal

Classical MD simulations are performed on a 50 nm × 4.3 nm rectangular graphene sheet containing 8160 carbon (12C) atoms with periodic boundary conditions applied in all directions. The simulation cell is sufficiently extended in the out-of-plane direction to ensure that atoms do not interact with their periodic images. All simulations are executed using the LAMMPS program.21 50% of the 8160 12C atoms are randomly selected and substituted with the heavier 14C isotopes. C-C interactions are modeled with the optimized Tersoff potential that is known to predict the transport phenomena and graphene thermal properties with reasonable accuracy.22 Each structure is energy minimized, initialized at 300 K and 0.1 MPa, and subjected to thermal equilibration under the isothermal-isobaric (NPT) ensemble for 7 nanoseconds (ns), under the canonical (NVT) ensemble for 7 nanoseconds (ns), and finally under the microcanonical (NVE) ensemble for 2 ns. Nóse-Hoover thermostat and barostat with coupling times of 0.1 picosecond (ps) are employed for the temperature and pressure constraints, respectively, and a timestep of 0.001 ps is used throughout. λ is computed with time averaging of the steady state temperature distributions using the Reverse Non-Equilibrium MD method (RNEMD) as explained in detail elsewhere.4,23 The simulations are performed under ergodic hypothesis for every structure examined during the search through the populations of isotope substituted graphene sheets across several generations of the GA scheme.

Our intent is to both explore the space of distributions, while looking out for specific configurations that minimize λ. We employ an in-house GA framework that can simultaneously deploy several molecular simulations across multiple computing clusters in an automated fashion.24 The GA operates by evolving an initially random population of candidate configurations through a process of selection, crossover, and mutation. Choosing a reasonable population size ensures that the space of distributions is well explored. Each configuration is represented as a “chromosome,” which essentially maps the atomic distributions into a binary code. Chromosomes with good fitness are more likely to be selected to pass on their genetic code, while those with poor fitness tend to be washed out of the population. Selection is a mechanism for choosing pairs of chromosomes in the population that will crossover segments of their genetic material, creating two “children” for the next generation. Those not selected for crossover still have a chance of being mutated. Mutation randomly alters bits of the chromosome, which can kick candidate configurations out of local optima into unexplored areas of the search space. The process of selection, crossover/mutation, and evaluation repeat for many generations with the overall population fitness improving, thus providing insight into which patterns of distributions result in lowered λ. Upon several successive generational searches, the GA narrows down the suitable nanostructure chemistries much faster than a brute force exhaustive trial approach.

Each test graphene sheet in the population is represented as a chromosome with 1 representing 12C and 2 for 14C. The nanostructure is thus represented as a 1-dimensional array with a size of 8160, which is one chromosome in the population. The GA method creates a population of trial graphene sheets based on a mersenne twister random number generator.25λ of each structure is computed via MD. Nanostructures with lowest thermal conductivities (best candidates) are preferred since they are more “fit” to transfer their genetic information, although some other structures are also allowed to transmit isotope distribution to ensure that the algorithm is not restricted to local minima landscape.

A tournament selection method is employed where for a participant size k (chosen to be k = 4), the scheme selects k random population members and determines the best candidate that becomes one of the parents of an offspring for the following generation. We cross the two chosen parents by a simple cutting point defined between 0 and 8160; all atoms before the cutting point are taken from parent A, and the rest are taken from parent B. Each structure undergoes mutation where each atom has a small probability of being flipped between 12C and 14C. This stage is an important feature of the GA because it prevents the search from being limited to exploring local minima of the λ predictions and facilitates the examination of structures that while not the best candidates currently, could subsequently become one with potentially minor changes. Fig. 2 shows a schematic representation of the computational procedure. We use a population size of 40 for each generation. Each MD simulation is deployed on 128 processors on TACC Stampede (https://www.tacc.utexas.edu/stampede/). Each generation consisting of 40 MD simulations is simultaneously deployed on TACC Stampede. This implementation is done in a fault-tolerant and automatic fashion using our in-house toolkit FT-AdaGio (Fault Tolerant Adaptive Grid Allocator) that is capable of automatically submitting each simulation in a generation to the computing cluster, while respecting queue limits and queue space availability. Furthermore, the fault-tolerant implementation ensures that progress of each simulation is tracked (via file output) with resubmission of any simulation in case of hardware or communication failures. Once all simulations within that generation are completed, the framework packages and returns all relevant data to the local computer where results are analyzed, before the next generation is deployed. On an average, a single population takes about 8 h of wall time to execute.

FIG. 2.

Genetic algorithm process flowchart. (Left) The overall GA scheme starts with an initial population of materials that is evaluated with our in-house fault tolerant adaptive grid allocator (FT-AdaGiO) and undergoes an evolutionary cycle (right). Members of the population are chosen for crossover based on their fitness to produce additional candidate solutions. Finally, each additional member is subject to random mutations to introduce another genetic material. This process of evaluation and evolution is repeated until a desired fitness level is reached.

FIG. 2.

Genetic algorithm process flowchart. (Left) The overall GA scheme starts with an initial population of materials that is evaluated with our in-house fault tolerant adaptive grid allocator (FT-AdaGiO) and undergoes an evolutionary cycle (right). Members of the population are chosen for crossover based on their fitness to produce additional candidate solutions. Finally, each additional member is subject to random mutations to introduce another genetic material. This process of evaluation and evolution is repeated until a desired fitness level is reached.

Close modal

Over 21 successive generations of search and interrogation of the isotope substituted graphene structures, the GA-MD scheme identifies the candidates that provide the lowest λ. To easily visualize the diversity of the configurations explored, we perform principal component analysis (PCA) on the set of configurations. PCA is a method of linear dimension reduction that optimally explains the variance of every observation in the dataset. Once PCA is performed, every observation from the original dataset is represented as a linear combination of an additional set of basis (“principal component” or PC). Thus, each configuration (represented as an 8160 length vector) can now be equivalently represented in a much lower dimensional space (using a few PCs), enabling visual inspection and reasoning. We choose to visualize the configurations in 3D space by using the first three PCs of each configuration. Fig. 3 illustrates all the data generated and recorded in terms of the first three PCs. These three components explain 33.5% of the total variance in the material space. We identify four general areas, termed as regions I through IV in Fig. 3. Structures containing in excess of 50% isotopes of one kind, 14C and 12C, are located in regions I and III, respectively. The other two regions II and IV contain graphene sheets with exactly 0.5 fraction of each type of the two isotopes. The PCA and results from MD calculations suggest that graphene structures with lower λ occur in the regions where PC-1 assumes negative values. When the 2D structure of maximum fitness (corresponding to the minimum λ) is selected from each region, the isotope substituted graphene sheet from region IV exhibits the lowest λ. Similar to reports in the literature,4,8 our computations predict that equiatomic ratios of isotope substitutions provide relatively lower λ values, as in regions II and IV. Thus, we corroborate that both concentration and distribution of isotopes affects λ of mass disordered graphene. While previous thermal conductivity measurements of isotope substituted graphene have noted the same response to mass disorder,26 they are limited in controlling and determining the precise distribution of isotopes that, as we have shown above, does contribute to the heat transfer in these materials.

FIG. 3.

Representation of the material space from the first three principal components. Four general regions based on thermal conductivity predictions are identified through visual inspection of this isotope substituted graphene materials landscape. The color index indicates the fitness score of an isotopic substituted graphene sheet material, which is proportional to the thermal conductivity.

FIG. 3.

Representation of the material space from the first three principal components. Four general regions based on thermal conductivity predictions are identified through visual inspection of this isotope substituted graphene materials landscape. The color index indicates the fitness score of an isotopic substituted graphene sheet material, which is proportional to the thermal conductivity.

Close modal

The distributions of the 14C for representative structures from regions I and II are compared to understand the effects of isotope arrangement on λ. We define a normalized 14C–14C distance parameter d* = dbond/d14C-14C, where dbond is a mean bond distance (0.14 nm) between 12C-12C atoms in a graphene lattice, and d14C-14C is the average distance between any two 14C atoms in an isotope substituted structure. Accordingly, for a structure in region I, d* will assume a higher average value than for representative graphene sheets in other regions. Fig. 4 presents the distribution of d* along the length of the graphene by computing the distances across 100 equally spaced bins. We find that for the structure with a high fraction of 14C isotopes, the distribution reveals several notable peaks and valleys unlike that for the structure in region II comprising 50% 14C. Peaks are indicative of the proximity of 14C isotopes, while a valley represents the opposite. In the presence of a higher concentration of isotopes, several locations in the graphene sheet consist of 14C clusters (peaks) or the lack of isotopes which implies 12C aggregation (valleys). Atoms in these locations predominantly interact with like atoms, reducing the likelihood of energy scattering by lattice vibrations (phonons). On the other hand, for structures with exactly 50% isotope composition, the absence of significant peaks and valleys suggests a uniform 14C distribution throughout increasing the possibility for interactions between unlike isotopes followed by enhanced phonon scattering due to mass disorder. Thus, the nonexistence of clusters of similar isotopes in 50% substituted graphene structures results in extremely reduced λ relative to structures with higher fraction of 12C or 14C.

FIG. 4.

Distribution of normalized distance between two 14C atoms along the length of the graphene sheet. Representative low thermal conductivity structures containing 50% and greater than 50% mass fraction of 14C in the isotope substituted graphene are examined for the distribution of isotopes in them. A higher concentration of isotopes promotes clustering of like atoms that reduces energy scattering. Equiatomic composition is devoid of such atomic aggregation, and the enhanced scattering rate strongly reduces the thermal conductivity relative to other structures.

FIG. 4.

Distribution of normalized distance between two 14C atoms along the length of the graphene sheet. Representative low thermal conductivity structures containing 50% and greater than 50% mass fraction of 14C in the isotope substituted graphene are examined for the distribution of isotopes in them. A higher concentration of isotopes promotes clustering of like atoms that reduces energy scattering. Equiatomic composition is devoid of such atomic aggregation, and the enhanced scattering rate strongly reduces the thermal conductivity relative to other structures.

Close modal

In summary, we describe results from a computational framework integrating genetic algorithm-based exploration with molecular dynamics simulations to optimize isotope substituted graphene sheets for reduced thermal conductivities. A fault tolerant adaptive scheme is employed to seamlessly examine hundreds of 2D structures by large scale computations over tens of generations. A principal component analysis of the thermal conductivity predictions reveals the occurrence of four discrete regions of structures according to the fraction and distribution of isotopes present in them. A detailed investigation of the isotope arrangements in representative graphene sheets from two of the regions shows the existence of clusters of like atoms in structures with higher concentration of isotopes. On the contrary, the absence of such atomic aggregation in graphene sheets with equiatomic composition of the two isotopes leads to enhanced scattering of energy by lattice vibrations and significantly lower thermal conductivities.

This material is based upon the work supported by the National Science Foundation (NSF) under Grant No. CMMI-1404938. We thank NSF and TACC for computing resources via XSEDE CTS110007. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors' and do not necessarily reflect the views of NSF.

1.
A. A.
Balandin
,
Nat. Mater.
10
(
8
),
569
581
(
2011
).
2.
A. A.
Balandin
,
S.
Ghosh
,
W. Z.
Bao
,
I.
Calizo
,
D.
Teweldebrhan
,
F.
Miao
, and
C. N.
Lau
,
Nano Lett.
8
(
3
),
902
907
(
2008
).
3.
A. K.
Geim
and
K. S.
Novoselov
,
Nat. Mater.
6
(
3
),
183
191
(
2007
).
4.
G.
Balasubramanian
,
I. K.
Puri
,
M. C.
Bohm
, and
F.
Leroy
,
Nanoscale
3
(
9
),
3714
3720
(
2011
).
5.
F.
Banhart
,
J.
Kotakoski
, and
A. V.
Krasheninnikov
,
ACS Nano
5
(
1
),
26
41
(
2011
).
6.
S.
Broderick
,
U.
Ray
,
S.
Srinivasan
,
K.
Rajan
, and
G.
Balasubramanian
,
Appl. Phys. Lett.
104
(
24
),
243110
(
2014
).
7.
F.
Hao
,
D. N.
Fang
, and
Z. P.
Xu
,
Appl. Phys. Lett.
99
(
4
),
041901
(
2011
).
8.
F.
Leroy
,
J.
Schulte
,
G.
Balasubramanian
, and
M. C.
Bohm
,
J. Chem. Phys.
140
(
14
),
144704
(
2014
).
9.
H.
Malekpour
,
P.
Ramnani
,
S.
Srinivasan
,
G.
Balasubramanian
,
D. L.
Nika
,
A.
Mulchandani
,
R. K.
Lake
, and
A. A.
Balandin
,
Nanoscale
8
(
30
),
14608
14616
(
2016
).
10.
U.
Ray
and
G.
Balasubramanian
,
Chem. Phys. Lett.
599
,
154
158
(
2014
).
11.
S.
Srinivasan
,
U.
Ray
, and
G.
Balasubramanian
,
Chem. Phys. Lett.
650
,
88
93
(
2016
).
12.
A.
Bagri
,
R.
Grantab
,
N. V.
Medhekar
, and
V. B.
Shenoy
,
J. Phys. Chem. C
114
(
28
),
12053
12061
(
2010
).
13.
Z. Y.
Ong
and
E.
Pop
,
Phys. Rev. B
84
(
7
),
075471
(
2011
).
14.
M. T.
Pettes
,
I. S.
Jo
,
Z.
Yao
, and
L.
Shi
,
Nano Lett.
11
(
3
),
1195
1200
(
2011
).
15.
K.
Saito
and
A.
Dhar
,
Phys. Rev. Lett.
104
(
4
),
040601
(
2010
).
16.
J. H.
Seol
,
I.
Jo
,
A. L.
Moore
,
L.
Lindsay
,
Z. H.
Aitken
,
M. T.
Pettes
,
X. S.
Li
,
Z.
Yao
,
R.
Huang
,
D.
Broido
,
N.
Mingo
,
R. S.
Ruoff
, and
L.
Shi
,
Science
328
(
5975
),
213
216
(
2010
).
17.
A. Y.
Serov
,
Z. Y.
Ong
, and
E.
Pop
,
Appl. Phys. Lett.
102
(
3
),
033104
(
2013
).
18.
A. A.
Balandin
,
Nat. Nanotechnol.
8
(
8
),
549
555
(
2013
).
19.
M. Z.
Hossain
,
S.
Rumyantsev
,
M. S.
Shur
, and
A. A.
Balandin
,
Appl. Phys. Lett.
102
(
15
),
153512
(
2013
).
20.
A. E.
Eiben
and
J.
Smith
,
Nature
521
(
7553
),
476
482
(
2015
).
21.
S.
Plimpton
,
J. Comput. Phys.
117
(
1
),
1
19
(
1995
).
22.
L.
Lindsay
and
D. A.
Broido
,
Phys. Rev. B
81
(
20
),
205441
(
2010
).
23.
F.
Mueller-Plathe
,
J. Chem. Phys.
106
(
14
),
6082
6085
(
1997
).
24.
D.
Stoecklein
,
M.
Davies
,
N.
Wubshet
,
J.
Le
, and
B.
Ganapathysubramanian
,
ASME. J. Fluids Eng.
139
(3),
031402
(
2017
).
25.
M.
Matsumoto
and
T.
Nishimura
,
ACM Trans. Model. Comput. Simul.
8
(
1
),
3
30
(
1998
).
26.
S.
Chen
,
Q.
Wu
,
C.
Mishra
,
J.
Kang
,
H.
Zhang
,
K.
Cho
,
W.
Cai
,
A. A.
Balandin
, and
R. S.
Ruoff
,
Nat. Mater.
11
,
203
207
(
2012
).