Part II
New Development in Density-Functional Theory

Density-Functional Approach to Atomic and Electronic Structures and their
Dynamics in Nano-Structures
Atsushi Oshiyama1, Junichi Iwata1, Taisuke Ozaki2, Tsuyoshi Miyazaki3, and Eiji Tsuchida4
1
Department of Applied Physics, The University of Tokyo
2
Institute for Solid State Physics, The University of Tokyo
3
Computational Materials Science Unit, National Institute for Materials Science
4
Nanomaterials Research Institute, National Institute ofAdvanced Industrial Science and Technology
The main purpose of this project is to develop the density-functional approach which enables us
to perform more accurate and larger scale static and dynamical calculations on massively parallelarchitecture computers, to clarify underlying physics and chemistry in phenomena in nanomaterials and structures, and to predict new phenomena in nano-world. We have adopted a realspace scheme as a principal methodology to solve the Euler equation in the density-functional
theory (DFT). We have also incorporated Krylov subspace method or the density-matrixoptimization method to make the computation cost scaling linearly to the system size N. Thoroughthe computations using these calculation schemes on the K computer at Kobe, Japan and also onother multi-core parallel architecture computers, We have clarified and predicted various
properties in nano-materials and structures. In particular, ACM Gordon Bell Prize in 2011 wasawarded to our real-space DFT calculations for silicon nanowires consisting of more than 100,000
atoms. The obtained results in this project are disseminated as 98 papers in scientific journals andpresented in 31 invited talks of major international conferences. The following is a summary of
those achievements (Ref. number cited below corresponds to the number in the publication list).
1.

Development of Real-Space Density-Functional-Theory (RSDFT) Approach [1,6]
To perform accurate density-functional calculations efficiently on the multi-core massively
parallel architecture computers, we have adopted a real-space (RS) finite-difference scheme in
DFT combined with the first-principles pseudo-potential method. We have developed a highly
efficient computation code named RSDFT. In the real-space scheme, discrete grid points are
introduced in real space, and the Hamiltonian in the Kohn-Sham (KS) equation in the DFT isexpressed as a matrix represented at the grid-point space. Differential operators for kinetic energy
are replaced by finite-difference operators with sufficiently high orders. The Hamiltonian matrix
therefore becomes sparse and Fast-Fourier Transformation (FFT), which is heavy communication
burden among all the compute nodes in the conventional plane-wave basis-set scheme, isunnecessary in the Hamiltonian operation. Hence, the communication becomes essentially local
only among neighboring compute nodes so that the total communication cost is tremendously
reduced in RSDFT.
The basic algorithm to solve Kohn-Sham (KS) equation is almost the same as that used inconventional O(N3) scaling method. The main procedures consist of three parts: ConjugateGradient minimization, Gram-Schmidt ortho-normalization, and subspace diagonalization. The
latter two are the bottlenecks which have O(N3) scaling. Those procedures are basically linear
algebraic calculations. In high performance computing (HPC), it is important to perform thecomputation using the data on the cache, the nearest memory to the compute node, since to fetch
the data from far memories or the hard disk costs a lot of time. Then when a ratio of the numbersof computation operations to the number of data required for the computation is high in an
algebraic calculation, it can be performed using mainly the data on the cache and thus the HPC is
feasible. In the matrix times matrix calculation, this ratio becomes high with increasing the target

45

size, compared with the inner product of vectors and also with the matrix times vector product.
In the Gram-Schmidt ortho-normalization in RSDFT code, eigenvectors, {¥it i.e.}, Kohn-Shamorbitals, are gathered in several groups so that apparent inner products of the vectors are converted
to the matrix-times-matrix calculations. This conversion renders our ortho-normalization to beefficiently done: The computational performance has reached 80 - 90 % of the theoretical peak
performance by using Basic Linear Algebra Subprograms (BLAS) level 3 on the K computer in
Kobe, Japan.
In addition to the ground state properties, we have made RSDFT applicable for Car-Parrinello
molecular dynamics (CPMD) simulation. Several benchmark tests on the K computer show thatthe real-space scheme is also suitable for CPMD implementation. We have further tuned the realspace CPMD code on the K computer. We are now capable of performing CPMD simulations fora few-thousands atom systems with practical CPU cost.
Our nano-meter scale targets consist of tens of thousands of atoms. In studies of such largesize systems, electronic energy-band structures provide important information to understandoptical and electrical characteristics of the systems. However, the calculation of energy bands inthe nano-system often becomes a heavier task than the self-consistent field calculation since many
eigenvalue problems with a large number of k-point parameters need to be solved and thecalculations must be performed in double-complex precision. In collaboration with people in thefield of applied mathematics, we have implemented in RSDFT code the Sakurai-Sugiura (SS)method for eigenvalue problems in which we are capable of obtaining eigenvalues in a specified
energy region. Hence the SS method is applied to the computation of the eigenvalues near the
Fermi level, thus substantially reducing the computational costs. In the SS method, we have tosolve linear equations with multiple shifts many times. We have also introduced an efficientalgorithm to solve the linear equations with multiple shifts.
2.

Application to RSDFT SchemeElectron States and Electric Current Distribution in Si Nanowires (NWs) [6,15]: Current
miniaturization of field effect transistors (FETs) is suffered from the limitation of the performance
inherent to the planar structures. Nanowire structures, i.e., all-surrounding gate structures arepromising in the next-generation FETs [Int. Tech. Roadmap for Semiconductors:
http://www.itrs.net/].
We have performed RSDFT calculations for various cross-sectional shapes such as circles,
ellipsoids and dumbbells with several cross-sectional dimensions and crystallographic wire
directions. The obtained energy bands show benefits of [110]-SiNW with the ellipsoid shape [15].
Combining the obtained energy bands and Koh-Sham orbitals from the RSDFT calculations witha simple transport mode [K. Natori, IEEE Trans.
Electron Devices 55, 2877 (2008)], we have calculateed the electric current density underappropriate bias voltage. Figure 1 shows calculated
current density represented on the wire crosssections. This information is important since thecarrier scattering at the interface of the wire and theoxide insulator may be crucial in designing allsurrounding gate transistor. The RSDFT calculationsfor SiNW were awarded Gordon Bell Prize for its
high performance.
Silicene on Ag Substrates[17,18]: Graphene is a
Fig. 1: Current density distribution on the
hexagonally bonded carbon sheet. Its most intriguingcross section of silicon NW FET as aproperty is the linear energy dispersion near thefunction of gate voltage.

46

Fermi level (EF) at particular symmetry
points (Dirac cone) in the Brilloine zone
(BZ). This gives rise to unique phenomena such as the anomalous quantum
Hall effect and the unexpected magnetic
ordering. Moreover, the very high inplane mobility makes graphene a
promising material utilized for nextgeneration devices. However, fromviewpoints of the technology connectivity, Si, which is located just below C in
the periodic table, is a more attractive
element. There has been a lot ofexperimental efforts to synthesize Si
monolayer (silicene) mainly on Ag
substrates. In this project, we have
performed RSDFT calculations forsilicene on Ag and found that Dirac
Fig. 2: Calculated STM images (right parts) and top
electrons are absent in this system due toviews of atomic structures (left parts) of the silicene on
the mixing between π orbitals of silicene Ag(111) surfaces. The stable [(a), (c) and (e)] and theand Ag orbitals. Instead we have propose metastable [(b), (d) and (f)] structures with the 4 x 4periodicity [(a) and (b)], the 13 x 13 periodicity [(c)that the BN sheets or the hydrogen- and (d)], and the 2 x2 periodicity[(e) and (f)] are3
3
passivated Si (111) surface is a good shown. Each circle of the left parts depicts the position
candidate for the substrate which provide of the Si atom in the silicene layer: The positions of the
topmost Si and the bottommost Si are depicted by the
Dirac electrons.
largest(red) and the smallest (blue) circles; theFigure 2 shows the stable structures ofintermediate-size circles depict the Si atoms between
the monolayer silicene on Ag(111) sur- the top and the bottom. The large gray balls depict theface. Experimentally several structures positions of the substrate Ag atoms. The lateral unit cellswith super- periodicity are identified. Our
calculations reproduce all the experimenttally observed structures and also predict
new atomic structures. The total energiesof those structures are close to each otherwith the energy difference of a few tens ofmeV, awaiting the experimental confirmation. Calculated STM images are in excellent agreement with the experiments whenavailable (Fig. 2).
Fig. 3: Calculated energy bands (a) of the silicene on
Our calculations show that in the freestanding silicene there are Dirac cones at EF the hydrogen-covered Si(111) surface [inset of (a)]and contour plots of the mixed π+ (π+*) states (b)though the silicene plane is substantially [absolute values]. Enlarged energy bands near EF = 0
buckled (a half Å). On the contrary, we have are shown in another inset and the distributions of thefound that Dirac cones are absent in any mixed π+ (π+*) states at K point marked by thestable structure of the silicne on Ag. The square are shown on the plane perpendicular to thesurface. The blue and pink balls in the inset depict Sibuckling itself mixes σ and π orbitals but
and H atoms, respectively.
keeps the Dirac cone near EF. The mixing ofSi π plus σ orbiatls and Ag orbitals, however, produces bonding and antibonding states between
Si and Ag orbitals, thus removing linear-dispersion states from EF.

47

As an alternative for the Ag substrates, we consider hexagonal BN sheet and hydrogenterminates Si (111) surface. Figure 3 shows calculated energy bands for the stable silicene on Hterminated Si (111) surface. The calculations show that the distance between the silicene layer
and the substrate hydrogen is 2.44 Å and the binding energy is 140 meV per Si atom. The energy
bands unequivocally shows the existence of Dirac cone at EF. The character of the states isobviously Si π orbitals. Those calculations clarify that soft-landing substrates are good candidates
for preparing silicene with Dirac cone.
Moiré-Induced Electron Localization in Twisted Bilayer Graphene[25]: Bilayer graphene
is synthesized by either the exfoliation or
by heat treatments of SiC surfaces,
introducing a new degree of freedom,
i.e., the way of stacking of C atomic
layers. When each hexagonal sheet isstacked identically, i.e., an atom of eachlayer is on top of a partner atom of theadjacent layer, it is called AA stacking. In
graphite, each layer is rotated in the basal
plane from the adjacent layer by 60o. Fig. 4: Side and top views of geometry optimized TBLG
This is called AB stacking or Bernal with the twist angles θ of (a) 29.4o and (b) 8.26o. In (b),stacking. In bilayer graphene (BLG), two the corrugation ∆ is visible with the far and the near
distances, dfar and dnear.
layers may be stacked with a twist angle
θ being 0o < θ < 60o (θ=0o and θ = 60o correspond to the AA and AB stackings, respectively).
Twisted bilayer graphene (TBLG) is indeed observed experimentally. The layers with a tiny twist
angle θ generate a Moiré pattern with its period L = 3d / [2sin(θ / 2)] , where d is the C-C bondlength. Interlayer interaction in BLG is smaller than the intralayer interaction by an order ofmagnitude. Hence one may consider that the difference in the stacking plays only minor role inatomic and electronic structures of BLG. However, it is not true for twisted bilayer graphene(TBLG) with small twist angle θ, as is shown below.
We have performed RSDFT calculations for 59 TBLGs with various twisted
angles including the smallest θ=0.76o(22,708 atoms in a unit cell) using the
local-density approximation (LDA). We
have examined the validity of LDA byrepeating some of the calculations using
van der Waals (vdW)-DFT functional. The
results obtained are found to be essentially
identical to those by LDA. Figure 4 showsgeometry optimized TBLGs for two
representative twist angles θ. For θ=29.4o,
Fig. 5: Fermi velocity normalized to that of monolayer
graphene vF0, as a function of θ: the velocities in thethe two layers are almost flat [Fig. 4(a)]optimized
structures (vF); those for the flat structureswith the interlayer distance of 3.42 Å.
with the interlayer distances of 3.61 Å (vFfar) and 3.34 Å
When θ becomes small, we find Moiré (v near). Circles are experimental data [Luican et al.,Fpatterns in which the AA and AB stacking PRL, 106, 126802 (2011)]. The vertical dashed lines
regions as well as other stacking regions depict the position of the critical angle θe.
constitute a patchwork as in Fig. 4(b). The
period L of the patchwork pattern and the size of each of the AA and AB stacking regions increase48

as θ becomes close to 0o (or 60o). For small θ, the two layers are corrugated. The interlayer
distance is large (dfar) in the AA stacking region and small (dnear) in the AB stacking region.
The variation of the Fermi velocity vF as a function of θ is shown in Fig. 5. We have foundthat vF is almost the same as that in monolayer graphene for θ > 10o, and that it is dramaticallyreduced toward zero when θ < θe. When we adopt the twist angle at which the Fermi velocity
becomes half of the velocity in monolayer graphene, θe is about 3o. Although the value of θe has
ambiguity in a sense that the half-maximum setting is not uniquely adopted, the existence of thiscritical angle is evident. Analyses of the Kohn-Sham orbitals clarify that the reason for the
reduction of vF and then the emergence of the flat bands is ascribed to the electron localization
due to the Moiré pattern. For larger twist angle, the Kohn-Sham orbitals at EF are characterizedby π orbitals extending over the whole carbon sheet. However, when θ becomes as small as 1o,
the Moiré pattern becomes prominent and the Kohn-Sham orbitals at EF are localized in the AA
stacking region. The Moiré pattern which becomes prominent for small θ induces inhomogeneity.
The electron wave senses the inhomogeneity, and finally the Kohn-Sham orbitals are localized.
Moreover, the reduction in vF is not monotonic. As shown in the inset of Fig. 5, vF vanishes at θM=1.08o and then slightly increases for the smaller θ. There is a magic angle θM at which vF vanishes.
Magnetic Ordering with Non-Magnetic Elements on SiC Nanofacet[23]: Nanometer-scale
structures are now accessible through cutting-edge fabrication technique or a gift of selforganization phenomena. Common expectation that such nanostructures break through limitationsof current technology may rely on a fact that nano-scale shapes affect wave-functions of relevant
electron states. However, a way of designing nano-shapes which decisively affect the electronstates and therefore are intriguing in science and useful in technology is lacking. Our RSDFT
calculations here show that peculiar electron states which are localized near but extended along
step edges of the nano-facets of SiC exhibit either ferromagnetic (FM) or anti-ferromagnetic(AFM) spin polarization. Considering that SiC is an emerging material in the next-generation
environmental-friendly power electronics and that spin current is the new non-dissipative carrier
for the information in future, this finding based on our RSDFT scheme has great impact.

Fig. 6: Side view of the(11-2n) (n=12) nanofacet
on the 4H-SiC(0001) vicinal surface, which wehave found is most stable.
Green, yellow and purple
balls indicate the Si, C
and H atoms, respectively.

Figure 6 shows the obtained stable structure of the nano-facet on SiC (0001) surface. The
shape of the nano-facet is characterized by the magic facet angle ( ϕ = 15.3o ) and facet height withthe unit-cell dimension, being in excellent agreement with the nano-facet experimentally observed.
The energy bands of this nano-facet is shown in Figure 7. Near EF, we observe many states whichhave characters of dangling bonds (DBs) of surface Si atoms on the long terrace. The energy
bands more than 1.5 eV below EF is mainly ascribed to the valence bands of the bulk SiC, whereas
the bands at about 1 eV above EF are from the bulk conduction bands. Apart from those states, wehave found four distinct states in the energy gap, labeled as SE1, SE2, SE3 and SE4 in Fig. 7(a).

49

From the analyses of the Kohn-Sham (KS) orbitals of these states [Figs. 7(b)-(e)], we have foundthat they originate from the step-edge states: The KS orbitals of the SE1 and SE2 states are
localized near the step edges where the upper terraces have cubic symmetry, whereas the orbitalsof SE3 and SE4 states are near the step edges where the upper terraces have the hexagonal
symmetry They are localized near but extended along the step edges, exhibiting dispersion-freeflat bands.

Fig. 7: Energy bands along the ΓX parallel to the step edge of the stable (11-2n) nanofacet on 4HSiC (0001) surface obtained by our LDA calculations. The Fermi level is set to be zero. The SE1 andSE2, and also the SE3 and SE4 are almost degenerate. The Kohn-Sham orbitals of the SE1, SE2, SE3
and SE4 states at the Γ point are shown in (b), (c), (d), and (e), respectively, as isovalue surfaces. The
red and blue colors depict positive and negative values, respectively, of the Kohn-Sham orbitals.

This flat-band nature is indicative of the magnetic ordering if EF shifts downwards to SE states.
We have found that this is realized by the hydrogen passivation. Our detailed DFT calculations
have clarified that H passivation/adsorption on Si surface atoms and further edge Si atoms areenergetically favorable. Without hydrogen, the Si DBs are either doubly occupied or unoccupied.
The H adsorption annihilates this DB states and generates the low-energy lying Si-H bonding
states which are doubly occupied. Then the H adsorption shifts EF, depending on which sites are
H passivated. Our RSDFT calculations have indeed clarified that the H passivation is capable ofcontrolling the filling of the SE states and induce either FM or AFM chains along the step edgesof SiC nano-facets which are organized spontaneously.
3.

Development of Numerically Exact Low-Order Scaling Method
Many scientific problems can be cast to an eigenvalue problem, and it is known that thecomputational cost of solving the eigenvalue problem scales as O(N3), where N is the dimensionof the problem. The density functional theory (DFT) for molecules and solids is also cast to an
eigenvalue problem within semi-local approximations such as local density and generalized
gradient approximations, and thereby the computational cost must be O(N3). Thus, it is interesting
to question whether the computational cost of the DFT calculation can be reduced from O(N3) toa low-order scaling without introducing approximations. We have addressed the problem andfound that such a method can be actually developed by utilizing of locality of basis functions [3].
Our method directly computes selected elements of the density matrix by a contour integration ofGreen's function [T. Ozaki, PRB 75, 035123 (2007)], which allows us to evaluate only O(N)density matrix elements due to the locality of basis functions. The Green's function is evaluatedwith a nested dissection approach for sparse matrices as shown in Fig. 8. It turns out that thecomputational effort of the method scales as O(N(log2N)2), O(N2), O(N7/3) for one-, two-, and50

three-dimensional systems, respectively. Unlike O(N)methods developed so far the approach is a numericallyexact alternative to conventional O(N3) diagonalization
schemes in spite of the low-order scaling, and can beapplicable to not only insulating but also metallic systems
in a single framework. It is also noted that the well
separated data structure is suitable for the massively
parallel computation. We expect that the direction will befurthermore explored in the near future to extend theapplicability of DFT calculations for large-scale systems
without sacrificing the accuracy.
4.

Development of massively parallel O(N) method
One of ways to extend applicability of the density
functional theory (DFT) calculations towards large-scale
scale systems is to employ O(N) methods of which
computational cost is proportional to the number of atoms.
We have been developing such an O(N) method based on
Fig. 8: (a)-(c) Nested dissection of athe divide and conquer method and Krylov subspace sparse matrix, (d) binary tree
method [T. Ozaki, PRB 74, 245102 (2006)], and proven by structure for the interaction in (c).
applying the method to a wide variety of systems
that the O(N) method can be applied to not only
insulators but also metals. In the method, a largescale system is divided into a set of small clusters,
where each cluster centers on each atom, and thecalculation of each cluster is performed nearly
independently. After performing the calculations of
all the clusters, the Green’s function associated withall the central atoms are collected as a patch work,
which enables us to calculate density matrix
elements which are required in the calculation of
charge density. Thanks to the algorithm, the O(N)method is intrinsically suitable to parallelFig. 9: Atomic decomposition by acalculation. However, it is not so obvious how themodified bisection and inertial tensor
clusters should be distributed to MPI processes in methods.
the parallel calculation to minimize communication
volume and memory usage. It is important for increase of parallel efficiency to localize data thateach process accesses and to minimize the volume and the number of MPI communication
involved. We have addressed the issue and developed a spatial decomposition method based on amodified recursive bisection method and inertia tensor [29]. The method consists of three steps:(1) for a given number of MPI processes a binary tree structure is constructed by using a modified
recursive bisection method, (2) a principal axis for a set of atoms is determined by diagonalizing
an inertia tensor with load balancing weights, and these atoms are projected onto the principal
axis, (3) these atoms are grouped on the principal axis based on the binary tree structure generated
in the step (1). As shown in Fig. 9, performing the steps (2) and (3) according to the binary tree
structure generated in the step (1) allows us to spatially decompose atoms while keeping
localization of atoms. The dynamic load balancing is also possible by taking account of elapsed
time in previous molecular dynamics or geometry optimization step in constructing the inertial51tensor. As an example, how a supercell of diamond is decomposed is shown in Fig. 10, where thesupercell contains 16384 carbon atoms, 19 MPI processes are used, and atoms colored by thesame color are allocated to the same MPI process. It is found that the spatial decomposition isproperly performed even using 19 MPI processes being a prime number. We have implemented
the decomposition method in OpenMX, and measured the parallel efficiency of the O(N) methodon the K-computer. The parallel efficiency reaches at 68 % using 131072 cores on the K-computer,
where the benchmark system is the diamond superlattice consisting of 131072 carbon atoms. We
have been applying the newly developed methods to electrode-electrolyte interfaces [13, 28],
metal-insulator interfaces [20], and grain boundaries in permanent magnets.

Fig. 10: Atomic decomposition of adiamond superlattice where 19 MPI
processes are used.

5.

Application of Open-MX Code to Two-Dimensional Si Structures
Although it is believedthat two-dimensional honeyExpt.
Calculated(b)comb structures consisting ofsilicon atoms do not exist
experimentallydue
torelative instability of its
hybridized sp2 orbitals, a
recent experiment clearly
demonstrates that silicene,honeycomb structure ofsilicon atoms, can be
fabricated on ZrB2 (0001)
thin films [12]. We have
performed detailed studies Fig.11: omparison between angle-resolved photoemisisonon geometrical and electron- spectroscopy (ARPES) data and the theoretical band structure
ic structures of silicene on calculated by the unfolding method of silicene on ZrB2 (0001).
ZrB2 and a related two-dimensional structure by means of electronic structure calculations based on density functional
theories (DFT), guided by close collaboration with experiments. Theoretical chemical shift of Si2p states [12] and band structure calculations [26], as shown in Fig.11, strongly support theformation of silicene having a planar-like structure. The stability of the planar-like structure over
the regularly buckled structure can be understood by interaction between states of the silicene andsurface states consisting of the d-orbital of the top Zr atoms [19]. We have also proposed a possible
mechanism for the formation of the domain structure of silicene on ZrB2 [12, 27]. It is inferredthat the domain structure is induced by an instability of a phonon having a nearly zero frequency,
and is formed in such a way that the k-points having the zero frequency can be removed from thefirst Brillouin zone. The mechanism is verified by performing large-scale total energy calculations
as shown in Fig. 12. We have further explored a possible structure of multi-layer silicene, andfind that the MoS2 structure consisting of silicon atoms is stabilized with atoms in the inner layer
having a sixfold coordination, which results in cigar-shaped nematic orbitals originating from the52Si-sp2 orbitals [Gimbert et al., PRB 90, 165423 (2014)].

Fig. 12: Total energy ofsilicene on ZrB2 (0001) withthe domain structure.

6. Development of Order(N) Density-Matrix Optimization (CONQUEST) Approach
First-principles technique based on the density functional theory (DFT) is a powerful method
widely used to elucidate the physical properties of various materials at the atomic scale. However,the size of the systems modeled with the conventional first-principles methods is limited because
the computational cost increases as the cube of the number of atoms N in the simulation cell. To
overcome this problem, we have been developing a linear-scaling DFT code CONQUEST, whose
memory requirements and CPU time are both only proportional to N (O(N) method). The code ispowerful and efficient on massively parallel supercomputers like K computer in Kobe. We havedemonstrated that DFT calculations on million-atom systems are now possible using K computer.
(Fig. 13)
Fig. 13: Wall-clock time of O(N) DFT
simulations of silicon systems using the K
computer (1 CPU = 8 cores). The horizontal
axis shows the number of atoms, indicatingthat DFT simulations of million-atomsystems are possible. We can see that we cancalculate systems containing twice thenumber of atoms in the same time bydoubling the number of CPUs, meaning thatthe computational cost is O(N) and theklill l ffi i
f hdDevelopment of the efficient and accurate local orbitals used in CONQUEST: We haveintroduced multi-site local orbitals (support functions) into CONQUEST in order to reduce thenumber of support functions to the minimal basis while maintaining as much accuracy as possible.
The multi-site support functions are constructed by taking linear combinations of the Pseudo
Atomic Orbitals (PAOs) on both the target atom and its neighbor atoms within a cutoff. Since thecomputational cost is proportional to the cube of the number of support functions both in O(N)and conventional diagonalization methods, we can significantly reduce the computational timeby using the scheme of multi-site support functions.
To determine the linear-combination coefficients, we have used the local filter diagonalization(LFD) method proposed by Rayson et al. (Phys. Rev. B 80, 205104(2009)) and the optimizationof the coefficients to minimize the total energy. The multi-site support functions with large cutoffs

53

provide energies and forces with comparable accuracy with the original (primitive) PAO results.
When the coefficients are determined only with the LFD method, the multi-site support functions
with small cutoffs provide large deviations from the original results. However, the optimizationof the coefficients significantly improves the accuracy of the multi-site support functions. The
stability and convergence of the calculations with small cutoffs are also improved by the
optimization. Thus, the coefficient optimization enables us to use small cutoff for multi-site
support functions keeping reasonable accuracy, which will reduce the computational cost
significantly. We have demonstrated the accuracy and efficiency of the method for hydrated DNA
system (Fig. 14 and Table 1).

Fig. 14: Density of states of the hydrated DNA system by primitive and multi-site support
functions. Black, blue and red curves correspond to the results by primitive, multi-site (LFD) andmulti-site (LFD + opt) support functions, respectively. rLD is set to be equal to rMS. Deviations
from the DOS by the primitive support functions are shown for the occupied region.
Table 1: Computational times [sec.] for matrix construction, diagonalization and gradientcalculation with respect to the coefficients for the hydrated DNA system by multi-site support
functions (rLD-rMS) and primitive support functions.

(8.0-8.0)(11.0-11.0)(16.0-16.0)
Primitive
DZP

Matrix
Gradient
Diagonalization
construction
calculation
34.5
627.7
4.3
54.4
583.4
8.6
484.8
566.3
25.7
28.0

12317.4–

Stable and Efficient Linear-scaling First-principles Molecular Dynamics (FPMD): The
recent progress of linear-scaling or O(N) methods in the density functional theory (DFT) isremarkable. We expect that first-principles molecular dynamics (FPMD) simulations based on
DFT could treat more realistic and complex systems using the O(N) technique. However, very
few examples of O(N) FPMD simulations exist to date, and information on the accuracy andreliability of the simulations is very limited. We have shown that efficient and robust O(N) FPMD
simulations are now possible by the combination of the extended Lagrangian Born−Oppenheimer
molecular dynamics method, which was recently proposed by Niklasson (Phys. Rev. Lett. 2008,54

100, 123004), and the density matrix minimization method as an O(N) technique. Using ourlinear-scaling DFT code CONQUEST, we have investigated the reliable calculation conditions
for accurate O(N) FPMD simulations and have demonstrated that we are now able to do practicaland reliable self-consistent FPMD simulations of a very large silicon system containing 32768
atoms, using 512 CPUs of FX10 or 1024 CPUs of K computer.
7.

Application of CONQUEST to a Semiconductor Hetero-Structure
Using the ability of CONQUEST, we worked on the growth mechanism of Ge threedimensional (3D) islands, called hut cluster, which are formed when Ge is deposited on the Si(001) surface. The 3D islands are made of four {105} facets and they grow under some conditions
by increasing the length of thelonger side while keeping thewidth constant. To control themorphology of the nanostructured system, understanding
the mechanism of the growth isvery important, but we need totreat the systems containing at
Fig. 15: Adsorption sites (left) and energy map (right) of single
least a few tens of thousands of
Ge dimers adsorbed on large and small facets of Ge hut cluster.
atoms to model these systems.
We have performed DFT studyof adsorption of Ge dimers on the{105} facets of a Ge hut cluster inorder to identify stable nucleationsites. By using a linear scaling
technique, we were able toperform robust and efficient
structure optimization of large
systems (20,000 – 200,000 atoms).
We calculated the energy map ofsingle Ge dimers adsorbed on two
different facets of an elongated hut
cluster (Fig. 15). Our results show
Fig. 16: Energy change during facets completion: red & greythat ridges and edges between two circles – completion of a small facet starting from the basefacets are stable and that the upper and apex, respectively, and blue rectangles – completion of
positions are generally more stable a large facet starting from apex.
than
lower
positionsbyapproximately 0.2eV/atom. We have also investigated that the completion of {105} facets by
calculating the stability of ”imperfect” facets (Fig. 16). The results together with the energy mapof single Ge dimers on the facets suggest that a new facet layer grows from the apex to the bottom.
8.

New Algorithms in Static and Dynamical Density-Functional Calculations
An algorithm for increasing the stability limit of ab initio molecular dynamics
simulations[32]: The computational cost of ab initio molecular dynamics (AIMD) is still aserious obstacle, even on a supercomputer. If, however, the goal of the simulation is to obtainlow-energy conformations through simulated annealing, or to equilibrate the system prior to the
production run, the computational cost is minimized by using the largest possible time step. When
the Verlet method is used to integrate the equations of motion, the maximum size of the time step55

is given by hmax = T/π, where T is the period of the fastest oscillation in the system. In practice,
however, AIMD simulations often break down at ~ hmax /2 due to the strong anharmonicity of
the interatomic forces. We have recently shown that the stability limit of the Verlet method may
be significantly increased by setting an upper bound on the kinetic energy of each atom with onlya small loss of accuracy. The validity of our approach was demonstrated for molten LiF.
Accelerating electronic structure calculations with mixed precision algorithms: Matrix
diagonalization is often the computational bottleneck in large-scale electronic structure
calculations. We have recently demonstrated that a careful use of 32-bit floating point operations
in matrix diagonalization allows us to accelerate the electronic structure calculations by more than
30% on Intel/AMD CPUs, while keeping 64-bit accuracy. Moreover, most of the computationallyexpensive operations are performed by level-3 BLAS/LAPACK routines in our implementations,thus providing optimal performance on many platforms. This approach may also be extended tosimilar cases, e.g. density-functional perturbation theory and/or time-dependent density
functional theory. It will also be interesting to utilize the mixed precision arithmetic in conjunction
with Graphics Processing Units (GPUs). [Tsuchida & Coe, Comp. Phys. Commun. 183, 980 (2012)]
Accelerating ab initio molecular dynamics using generalized masses[10]: Mass tensor
molecular dynamics method was originally introduced by Bennett four decades ago for enhanced
sampling of the phase. The generalized Hamiltonian is given by
H (q, p) =1
∑ pi M ij−1 (q) p j + U (q)2 i, jwhere M is a positive definite, symmetric matrix of dimension 3Natom. If M is chosen carefully,
the efficiency of ab intio molecular dynamics can be significantly improved without changing
the results. Test calculations on liquid water consisting of 64 molecules showed a threefold
reduction in computational effort without making the fixed geometry approximation. We alsopresented a simple algorithm for estimating the optimal atomic masses using only the first
derivatives of the potential energy. This approach will be particularly useful for calculating
ensemble averages, e.g. radial distribution functions and free energies, in finite temperature
simulations. It would also be interesting to combine this approach with rare event methods, suchas hyperdynamics and metadynamics.
Collaborators
Kazuyuki Uchida (Kyoto Sangyo University), Zhi-Xin Guo (Xiangtan University), Keisuke
Sawada (Riken), Mauro Boero (Strasbourg University), Ken-ichi Koizumi (Institute for
Molecular Science), Shinnosuke Furuya (ARGO GRAPHICS Inc.), Yasuteru Shigeta (Universityof Tsukuba), S. Arapan (NIMS), A. Nakata (NIMS), D. R. Bowler (UCL), Yu-ichiro Matsushita(The University of Tokyo), Sukmin Jeong (Chonbuk University), Shinya Kyogoku (Toshiba
Corporation).
Publications and Invited Talks
Selected List of Publications (the number of papers published in the period of this project =98)
1.

2.
3.
4.

J.-I. Iwata, D. Takahashi, A. Oshiyama, B. Boku, K. Shiraishi, S. Okada and K. Yabana, ``A massivelyparallel electronic-structure calculations based on real-space density functional theory" J. Comp. Phys.
229, 2339 (2010).
Y. Gohda and A. Oshiyama, ``Stabilization mechanism of vacancies in group-III nitrides: Exchange
splitting and electron transfer" J. Phys. Soc. Jpn. 79, 083705 (2010).
T. Ozaki, “Efficient low-order scaling method for large-scale electronic structure calculations withlocalized basis functions”, Phys. Rev. B 82, 075131 (17 pages) (2010).
Y. Matsushita, K. Nakamura adn A. Oshiyama, ``Comparative study of hybrid functionals applied tostructural and electronic properties of semiconductors and insulators" Phys. Rev. B 84, 075205 (2011).

56

5.
6.

7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.

23.
24.
25.
26.

27.

28.
29.
30.
31.

S. Jeong and A. Oshiyama, ``Selective Alignment of Carbon Nanotubes on Supphire Surfaces: Bond
Formation between Nanotubes and Substrates" Phys. Rev. Lett. 107, 065501 (2011).
Y. Hasegawa, J.-I. Iwata, M. Tsuji, D. Takahashi, A. Oshiyama, K. Minami, T. Boku, F. Shoji, A. Uno,
M. Kurokawa, H. Inoue, I. Miyoshi, M. Yokokawa, ``First-principles calculations of electron states ofa silicon nanowire with 100,000 atoms on the K computer”, Proceedings of 2011 International
Conference for High Performance Computing, Networking, Storage and Analysis (SC2011)
M. Toyoda and T. Ozaki, “Exchange functional by a range-separated exchange hole”, Phys. Rev. A 83,
032515 (2011).
T. Ozaki and M. Toyoda, “Accurate finite element method for atomic calculations based on density
functional theory and Hartree-Fock method”, Comp. Phys. Comm. 182, 1245-1252 (2011).
M. P. Sena, T. Miyazaki and D. R. Bowler, “Linear Scaling Constrained Density Functional Theory in
CONQUEST”, Journal of Chemical Theory and Computation, vol. 7, p884-889 (2011)
E. Tsuchida, ``Ab initio mass tensor molecular dynamics”, J. Chem. Phys. 134，044112 (2011).
Y. Matsushita, S. Furuya and A. Oshiyama, ``Floating Electron States in Covalent Semiconductors"
Phys. Rev. Lett. 108, 246404 (2012).
A. Fleurence, R. Friedlein, T. Ozaki, H. Kawai, Y. Wang, and Y. Yamada-Takamura, "Experimental
evidence for epitaxial silicene on diboride thin films", Phys. Rev. Lett. 108, 245501 (5 pages) (2012).
T. Ohwaki, M. Otani, T. Ikeshoji, and T. Ozaki, "Large-scale first-principles molecular dynamics forelectrochemical systems with O(N) methods", J. Chem. Phys. 136, 134101 (9 pages) (2012).
D. R. Bowler and T. Miyazaki, “O(N) methods in electronic structure calculations” , Report on
Progress in Physics, vol. 75, 036503 (2012).
S. Kyogoku, J.-I. Iwata and A. Oshiyama, ``Relation between nanomorphology and energy bands of
Si nanowires" Phys. Rev. B 87, 165418 (2013).
K. Koizumi, M. Boero, Y. Shigeta, and A. Oshiyama, ``Atom-Scale Reaction Pathways and FreeEnergy Landscapes in Oxygen Plasma Etching of Graphene" J. Phys. Chem. Letters 4, 1592 (2013).
Z.-X. Guo, S. Furuya, J.-I. Iwata and A. Oshiyama, ``Absence of Dirac Electrons in Silicene on
Ag(111) Surfaces" J. Phys. Soc. Jpn. Letters 82, 063714 (2013).
Z.-X. Guo, S. Furuya, J.-I. Iwata and A. Oshiyama, ``Absence and Presence of Dirac Electrons in
Silicene on Substrates" Phys. Rev. B 87, 235435 (2013).
C.-C. Lee, A. Fleurence, R. Friedlein, Y. Yamada-Takamura, and T. Ozaki, "First-principles study oncompeting phases of silicene: Effect of substrate and strain", Phys. Rev. B 88, 165404 (2013).
H. Sawada, S. Taniguchi, K. Kawakami, and T. Ozaki, ``First-principles study of interface structureand energy of Fe/NbC", Modelling Simul. Mater. Sci. Eng. 21, 045012 (2013).
T. Otsuka and T. Miyazaki, ”A quantum chemistry study of Ds-Pa unnatural DNA base pair”, Int. J.
Quantum Chem. Vol. 113, p504-509 (2013).
T. Ohkubo, E. Tsuchida, M. Gobet, V. Sarou-Kanian, C. Bessada, Y. Iwadate, ``First-Principles
Molecular Dynamics Simulation and Conductivity Measurements of a Molten xLi2O-(1-x)B2O3
System， J. Phys. Chem. B 117，pp.5668-5674 (2013).
K. Sawada, J.-I. Iwata adn A. Oshiyama ``Magic Angle and Height Quantization in Nanofacets on
SiC(0001) Surfaces", Appl. Phys. Lett. 104, 051605 (2014).
Y. Matsushita and A. Oshiyama, ``Interstitial Channels that Control Band Gaps and Effective Masses
in Tetrahedrally Bonded Semiconductors", Phys. Rev. Lett. 112, 136403 (2014).
K. Uchida, S. Furuya, J.-I. Iwata and A. Oshiyama, ``Atomic corrugation and electron localization
due to Moiré patterns in twisted bilayer graphenes" Phys. Rev. B 90, 155451 (2014).
C.-C. Lee, A. Fleurence, Y. Yamada-Takamura, T. Ozaki, and R. Friedlein, "Band structure of silicene
on zirconium diboride (0001) thin-film surface: Convergence of experiment and calculations in the
one-Si-atom Brillouin zone", Phys. Rev. B 90, 075422 (2014).
C.-C. Lee, A. Fleurence, R. Friedlein, Y. Yamada-Takamura, and T. Ozaki, "Avoiding critical-point
phonon instabilities in two-dimensional materials: The origin of the stripe formation in epitaxial
silicene", Phys. Rev. B 90, 241402(R) (2014).
T. Ohwaki, M. Otani, and T. Ozaki, "A method of orbital analysis for large-scale first-principles
simulations", J. Chem. Phys. 140, 244105 (10 pages) (2014).
T.V.T. Duy and T. Ozaki, ``A three-dimensional domain decomposition method for large-scale DFT
electronic structure calculations", Comput. Phys. Commun. 185, 777-789 (2014).
A. Nakata, D. R. Bowler and T. Miyazaki, “Efficient Calculations with Multisite Local Orbitals in a
Large-Scale DFT Code CONQUEST”, J. Chem. Theory Comput., Vol. 10, 4813-4822 (2014).
M. Arita, D. R. Bowler and T. Miyazaki, “Stable and Efficient Linear Scaling First-Principles

57

Molecular Dynamics for 10000+ Atoms”, J. Chem. Theory Comput., Vol. 10, 5419-5425 (2014).
32. E. Tsuchida，``Stabilization of ab initio molecular dynamics simulations at large time steps”，JPS
Conference Proceedings 5，011019、2015.
33. A. Nakata, D. R. Bowler and T. Miyazaki, “Optimized multi-site local orbitals in the large-scale DFT
program CONQUEST”, Phys. Chem. Chem. Phys., published online (10.1039/C5CP00934K) (2015).

Selected List of Invited Talks in International Conferences(total number in the period of this project =31)
1.
2.

3.

4.
5.
6.
7.
8.

9.

10.

11.

12.
13.

14.

15.

A. Oshiyama, ``Materials Design through Computics: nanowires and Nanotubes" Inernational Focus
Workshop on Quantum Simulations and Design (September 27, 2011, Dresden, Germany)
T. Ozaki, "Low-order scaling methods for large-scale density functional calculations", Summer school
on Electronic Structure Analysis and Computation, Shanghai Jiao Tong Univ., Shanghai, China, June
6-10, 2011.
A. Oshiyama, ``Materials Design through Computics: Nanostructures of Silicon and Carbon" 10th Int.
Meeting on High Peformance Computing for Computational Science (VECPAR2012), July 17-20,
2012, Kobe, Japan).
A. Oshiyama, ``Real-Space Density-Functional Approach to Electronic Properties of Nanostructures"
Conference on Computaional Physics, (October 14-18, 2012, Kobe, Japan)
J.-I. Iwata, “First-principles electronic structure calculations with K computer”, Theory andApplications of Computational Chemistry (TACC2012), 2-7 Sep, 2012, Pavia, Italy
T. Ozaki, “Massively parallelized O(N) Krylov subspace method for large-scale DFT calculations”,
IUMRS-ISEM 2012, PACIFICO Yokohama, Japan, Sep. 23-28, 2012.
T. Ozaki, “Low-order Scaling Density Functional Methods Based on Quantum Nearsightedness”,
SIAM Conference on Parallel Processing for Scientific Computing, Savanna, USA, Feb. 15-17, 2012.
T. Miyazaki, “Atomic and electronic structures of BEDT-TTF and M(dmit)2 systems at highpressures: A first-principles study of organic molecular solids”, Int. Symp. of the Quantum Beam
Science Directorate, Japan Atomic Energy Agency, (Mito, Ibaraki, Japan), Sep. 23-27, 2012.
A. Oshiyama, ``Materials Design through Computics: Large-Scale Density-Functional Calculations
for Nanomaterials in the Real-Space Scheme" Int. Summer School on HPC Challenges in
Computational Sciences (New York City, USA, June 2013).
A. Oshiyama, ``Large-Scale Density-Functional Calculations in the Real-Space Scheme: Graphene
and Silicene" 7th Conf. Asian Consortium on Comnputational Materaisl Science (ACCMS7) (Nakhon
Ratchasima, Thailand, July 2013).
Z.-X. Guo, J.-I. Iwata and A. Oshiyama ``Large-Scale Calculations in the Real-Space Scheme:
Absence and Presence of Dirac Electrons in Silicene" 22th Int. Materials Research Congress (IMRC)
2013 (Cancun, Mexico, August, 2013)
E. Tsuchida，``Adaptive Finite Element Method for Ab Initio Molecular Dynamics Simulations”，
12th U.S. National Congress on Computational Mechanics，U.S.A., 2013/07/22.
A. Oshiyama, ``Large-Scale Real-Space Density-Functional Calculations: Moir¥'e-Induced Electron
Localization in Graphene and Floating States in SiC" 32nd International Conference on the Physics
of Semiconductors (Austin, Texas, USA, August 2014).
T. Ozaki, “Theoretical and experimental exploration of two-dimensional Si structures”, The 17th
Asian Workshop on First-Principles Electronic Structure Calculations, Yonsei University, Seoul,
Korea, Nov. 2-5, 2014.
T. Miyazaki, “Large-scale DFT simulations on nano structured materials using a linear-scaling
technique”, Collaborative Conference on 3D Research (CC3DMR) 2014, Songdo Convensia (Incheon,
South Korea), Jun. 23-27, 2014.

Prizes1.

2.

Y. Hasegawa, J.-I. Iwata, M. Tsuji, D. Takahashi, A. Oshiyama, K. Minami, T. Boku, F. Shoji, A. Uno,
M. Kurokawa, H. Inoue, I. Miyoshi, and M. Yokokawa, ``First-Principles Calculation of Electronic
States of a Silicon Nanowire with 100,000 Atoms on the K Computer” ACM Gordon Bell Prize (peakperformance) 2011
S. Arapan, ``JSST 2013 Outstanding Presentation Award” at the Int Conf on Simulation Technology,
JSST 201358

Studies on Structure Sampling and Non-Equilibrium Dynamics
Using Ab-Initio Molecular Dynamics Methods
Shinji Tsuneyuki1, Yoshihide Yoshimoto2,3, Jun Yamauchi4 and Minoru Otani5
Department of Physics, The University of Tokyo1,
Department of Applied Mathematics and Physics, Tottori University2,
Department of Information Science, The University of Tokyo3,
Department of Physics, Keio University4,
National Institute for Advanced Industrial Science and Technology5

1. Introduction
This group is mainly dedicated to the first-principles study on thermal physics and nonequilibrium dynamics in condensed matter, such as thermal conduction, solid-liquid phasetransition, electrochemical reactions at the electrode-electrolyte interface and so on. In thesimulation of these physical phenomena, where spaciotemporal fluctuation of the atomic structure
is so large, we need careful modeling and long-time/large-scale simulation of the atomic motion
to get statistically meaningful results. This is a hard requirement for first-principles simulations,firstly because the computational cost of the electronic structure calculation usually scales as N3
with N being the number of electrons, and secondly because it is essentially impossible toparallelize the dynamical simulation due to causality. On the other hand, for the simulation ofelectrochemical reactions, microscopic modeling method for the simulation is in itself
problematic due to applied voltage and non-conserving electron number. Comparison of thesimulation with spectroscopic experiments is also a difficult but indispensable task.
In this project, in order to overcome some of these difficulties, we took the following strategy.
Yoshimoto developed a highly parallelized first-principles electronic structure calculation
program package xTAPP and opened it to the public (Section 2). The program can be a powerfulbasis for the simulation of thermal physics and non-equilibrium dynamics. Tsuneyuki developed
an efficient scheme of modeling anharmonic interaction between atoms from first principles(Section 3). The anharmonic model potentials were then used for highly accurate calculation oflattice thermal conductivity. The program package ALAMODE was developed for all theprocedure and it was also opened to the public. Otani developed new useful methods forelectrochemistry by the extension of the effective screening medium method previously proposedby himself (Section 4), with which he realized simulation of electrochemical reaction underarbitrary and constant bias voltage. Finally, Yamauchi investigated the core-level X-ray
photoelectron spectroscopy from first principles, which enabled us to identify local defectstructures in semiconductors (Section 5).
Development of a plane wave basis first-principles electronic structure calculation package
xTAPP (Y. Yoshimoto)First-principles electronic structure calculation based on density functional theory and its local
density approximation is widely used in various research fields because the ratio of its accuracy
to the computational cost is excellent. The method is further classified by the used basis set.
Among them, the method with plane wave basis set using pseudo-potentials is commonly usedbecause it is good at atomic force filed calculations, can equally treat systems from sparse ones
to dense ones, and can control the computational accuracy due to the basis set easily.
2.

59As an implementation of the method, xTAPP was developed based on TAPP and it was opened
to the public at the following web sites.
http://xtapp.cp.is.s.u-tokyo.ac.jp
http://ma.cms-initiative.jp/ja/listapps/xtapp
The developed features of xTAPP in this project are as follows:












hybrid exchange correlation (XC) potentials with GPGPU computation
pseudo-potential library: all elements from H to Bi, except for radioactive ones (Tc, Pm)
and lanthanoids other than La, Ce, Nd, Gd
integration of maximally localized Wannier function code [1] and its application tospontaneous electronic polarization and projected DOS
surface under an electrostatic field by the effective screening medium method [2]total DOS, projected DOS, integrated charge distribution as a function of energy by using
the tetrahedron method
DFT-D[3] dispersion force
linkage with phonopy [4] : phonon related calculation
support for GRAPE-DR accelerator
TAPP input format is replaced with a modern one
user manual(22 page), reference manual for input files(33 page), tutorial and examples (53page), developers reference (7 page), Reference for formulations (72 page)xTAPP-test : examples of input files (93 patterns) and automated comparison program of
the test results

The package is opened to public under GNU General Public License version 3. This license is
a comfortable one in computer science, and it becomes easy for xTAPP itself to incorporate
foreign codes. The copyright of TAPP codes was cleaned and the code that is not suitable for the
public release was replaced. The consortium that maintains the TAPP code consented to makexTAPP open to the public. The replaced code includes the modernized code for sampled k points
with corresponding symmetry operations and it was utilized when the code to compute maximally
localized Wannier functions was integrated.
xTAPP implemented a fast computation of hybrid exchange correlation (XC) functionals byusing GPGPUs. The hybrid XC functionals (PBE, HSE, LC, and B3LYP etc.) are beyond local
density approximation and involve exact computation of exchange interactions. They improve,for example, description of water molecules. The obstacle for its wide use is the highcomputational cost of the exchange interaction. One solution of the issue is to use economical
computing devices such as GPGPUs.
The biggest issue in the use of GPGPUs is the relatively slow data transfer speed between main
memory and the device memory of a GPGPU.
The analysis of the structure of the computation from the computer science resolves the issue.
The computation is performed on each pairs of orbitals and its result is reduced. Namely, the orderof the arithmetic operations is O(N2) although the data size is O(N), where N is the number oforbitals. This means that the blocking of the computations is the appropriate approach: instead ofstep by step transferring each orbital pair to GPGPU and transferring back the computed potential,
we should transfer Nblk orbitals from each two series of orbitals to the device memory at a time,compute all the potentials for Nblk×Nblk pair block formed by the transferred orbitals and transfer
back the reduced potential. This improves the ratio of arithmetic operations to the transferred data
by Nblk times and thus we can hide the slow data transfer rate.

60

cut of wave number of the plane

3.6

4.0

4.8

5.0

5.4

6.4

FFT mesh

723

803

963

1003

1083

1283

Xeon X5690 [s]

378

549

994

1188

1583

2255

Radeon HD 6950 [s]

169

297

342

500

534

749

acceleration

2.23

1.84

2.91

2.37

2.96

3.01

wave basis [a.u.]

Table 1: 1 SCF time for Si 216 atom system

Based on this theoretical analysis, we adopted AMD Radeon 6950 (26,000 JPY) and OpenCL
language because of its cost benefit.
Table 1 shows the comparison of 1 SCF time of Si 216 atom system (Г point sampling) between
GPGPU and CPU(Intel Xeon X5690; 140,000 JPY). Nblk was 20. 8 GPGPU or 8 CPU were usedin parallel.
The higher the cut of the plane wave basis set becomes, the higher the obtained acceleration
ratio is in general. Up to 3 times acceleration was obtained. Because the cost of the Radeon
GPGPU was rather low, this is meaningful number. The depressions at FFT mesh of 803 and 1003
may be because the Radeon GPGPU is not good at computation of radix 5.
xTAPP only treats valence electrons using pseudo-potentials. Therefore a pseudo-potential
library is indispensable to make it adopted by a wide range of users. In this project libraries for(ultrasoft pseudo-potential / PBE type XC potential), (ultrasoft pseudo-potential / Ceperley-AlderPerdew-Zunger type XC potential), (norm conserving pseudo-potential / PBE type XC potential)
were prepared.
These libraries contain all elements from H to Bi, except for radioactive elements (Tc, Pm) andlanthanoids other than La, Ce, Nd, Gd.
The libraries were checked by developed test system for pseudo-potentials: the bulk band
diagrams to detect the ghost states, the reproducibility of lattice constants, bulk moduli, cohesiveenergies, vibration frequencies and equilibrium inter-atomic distance of diatomic molecules. The
result is included in the library.
The other extension of xTAPP includes effective screening medium method by M. Otani and
O. Sugino and computation code of maximally localized Wannier orbitals by K. Nakamura. To
integrate these functions, we have to check the soundness of existing functions. For this purpose,
we developed xTAPP-test, a test package with a specialized code to compare the results. xTAPPtest resolves the too many false-positive issue in the comparison between the results before andafter the integration of the target function.
Thus we can easily confirm the soundness ofD0 [10-6 m2/s] Ea [kJ/mol]
existing
functions.
xTAPP-test
is
PBE
3200
46
indispensable to keep the quality of xTAPP.
PBE0
70
31
The application of xTAPP includes themodeling of the low energy bands of dmit and
B3LYP
48
29
BEDT-TTF systems, organic semiconductor
D3LPY + DFT-D3
6.6
23
PTCDA on Si(001) surface (in cooperation
experimental [5]2.517with an experimental group) and modeling ofwater systems by a classical force fields using
Table 2: Self diffusion coefficients of water
thermodynamic downfolding method. The
up to 360 K
obtained force field for water was applied to the61self diffusion coefficients of water (Table 2), where diffusion constant D as a function oftemperature T is D = D0 exp(-Ea/kBT). To obtain good classical force fields both the hybrid XC
potential and dispersion forces are important.
The number of unique users registered at the xTAPP download site from April 2013 to March
2015 was more than 120.
3. Anharmonic lattice model and the lattice thermal conductivity from first principles (S.
Tsuneyuki)
Lattice thermal conductivity of crystals is governed by the anharmonicity of interatomic
potentials. It is an essential ingredient of the figure of merit for thermoelectric materials, where
lower thermal conductivity is preferred. On the other hand, higher thermal conductivity ispreferred in electronic devices for thermal stability. From the viewpoint of materials design forthese materials, we have developed a new scheme to calculate lattice thermal conductivity fromfirst principles.
So far three methods have been proposed for computational study of thermal conductivity: amethod based on Boltzmann’s transport equation and diagrammatic evaluation of lifetime ofphonons, molecular dynamics (MD) simulation with the Kubo-Greenwood formula and nonequilibrium MD simulation of heat transfer with temperature gradient. In any case, first-principles
simulation is highly difficult because of the large system size and long-time simulation or many
sampling points in k space integration necessary for reliable calculation. Thus we decided todevelop an accurate, efficient and versatile method for modeling anharmonic interatomic
potentials with first-principles MD based on DFT.
The model potential is just a Taylor expansion of the total energy of the system calculated by
DFT,
1
2

1
3!

Φ
,

Φ

⋯.

, ,

denotes the displacement of each atom in xyz direction from its potential minimum, Φharmonic force constant and Φ

⋯

is a

is called an anharmonic force constant. Practically the Taylor

expansion is terminated at a low order and the force constants are determined so that the model
potential minimizes the deviation of the force acting on each atom from that calculated from firstprinciples at some sampling structures under crystalline symmetry constraints [6]. In this study,
we adopt up-to the sixth-order force constants for the parameter fitting. For efficient sampling of
the potential energy hyper-surface, we use the locus of atoms obtained by first-principles MD at
a fixed temperature [7]: the temperature should be high enough to reflect anharmonic characterof the potential but low enough to avoid unphysical sampling. We find the third-order force
constants determined in this way are as accurate as those carefully obtained by numericaldifferentiation of the potential energy.
Once the force constants are determined, phonon dispersion, group velocity and the lifetime ofeach phonon mode are calculated, with which lattice thermal conductivity and its temperaturedependence is obtained within the relaxation-time approximation. We applied the scheme to Si,
Mg2Si, PbTe, Bi2Te3 and Ba8Ga16G30 and found that material dependence and temperature

62

Fig. 1 Temperature dependence of the lattice thermal conductivity of various materials.
Experimental data were taken from R. Kremer et al., Solid State Commun. 131, 499 (2004) (Si);J. J. Martin, J. Phys. Chem. Solids 33, 1139 (1972) (Mg2Si); I. Ravich et al., Semiconducting leadchalcogenides, Plenum Press (1970) (PbTe); C. B. Satterthwaite et al., Phys. Rev. 108, 1164 (1957)(Bi2Te3); M. Avila et al., Phys. Rev. B 74, 125109 (2006) (Ba8Ga16G30).

dependence of the thermal conductivity are quantitatively reproduced in a wide range from 100 to
104 W/mK as shown in Fig.1.
An apparent advantage of the use of the model potential instead of direct calculation of thermal
conductivity is that we can adopt system size large enough to include the effect of phonon
scattering by those with long wave length (typically, more than 104 atoms for MD or similar
number of sampling points in the reciprocal space for lifetime calculation), and MD simulation
long enough to match the life time of phonons (much more than 102 ps), although the modelingneeds only 102 atoms and a few ps MD simulation
for sampling.
Another advantage is that the system property
can be controlled by changing the configuration of
the model. For example, in case of Ba8Ga16G30, aclathrate containing rattlers inside, an effect of therattlers for its extremely low thermal conductivity
can be clarified by investigating the thermal
conductivity of an imaginary clathrate compound
without the rattlers. By comparing the simulations
with and without the rattlers, we found that the
rattlers cause tenfold reductions in the relaxation
time of phonons at low frequency by scattering
phonons (see Fig.2), while their impact on thegroup velocity is secondary [8].

63Fig. 2 Frequency dependence of
calculated life time of phonon modes for
Ba8Ga16G30 with (green) and without(purple) rattlers in the cage.

The program package ALAMODE was developed by T. Tadano for the present procedure andit was opened to the public.
http://github.com/ttadano/alamode; http://alamode.readthedocs.org
http://ma.cms-initiative.jp/ja/listapps/alamode/alamode

4. Extension of the effective screening medium method (M. Otani)
4-1 Introduction
An electrochemical reaction at the electrode-electrolyte interface is one of the most important
issues on performance and durability of energy harvesting devices such as a fuel cell, a secondary
battery, a photovoltaic cell and a photoelectrochemical cell. Although the first-principles
calculation technique is a reasonable choice to study the electrochemical reaction, the threedimensional periodic boundary condition, which is usually employed in a conventional firstprinciples calculation program, is not compatible with the electrochemical system. To make theprogram compatible with the system, we have developed the effective screening medium (ESM)
method [9] which allows us to simulate the electrochemical system by modifying the boundary
condition of the Poisson equation. Using the ESM method, we can calculate many
physical/chemical properties of the electrochemical systems in operando conditions and cananalyze the experimental result [10].
Although the ESM method enables us to simulate electrochemical system, there is a room toimprove the modeling of the system. For example, firstly the original ESM method needs a
vacuum region between a slab and the ESM to avoid the electronic wave function penetration into
the ESM region. It is difficult to define the applied bias potential between an electrode and theESM owing to the existence of the vacuum. Secondly conventional first-principles calculations
are done under fixed total number of electron (constant-N calculation), whereas theelectrochemical system is governed by an external potentiostat. Thus the calculation should bedone by the fixed Fermi energy (constant-μ calculation).
In this project we have been trying to develop the ESM method to solve above issues. In thefollowing subsections, we will describe how to remove the vacuum region and control the Fermi
energy during the first-principles calculations.
4-2 Smooth-ESM method
The reason why we need to introduce the vacuum region between the slab and the ESM is that
there is a discontinuity in the permittivity constant [see Fig.3 (a)]. When the electronic wave
function locates in the ESM region, we cannot solve the Poisson equation. To avoid this issue, we
introduce the following function form for the permittivity in the ESM region [see Fig.3 (b)]: (z)  e 2 a(zz ),1

where a denotes the screening length of the ESM region [11]. The value of the permittivity
smoothly changes from unity to infinity; we call this the smooth ESM. With this function form,
we can still solve the Poisson equation analytically and obtain the Green’s function in advance.
Thus the calculation cost remains comparable to that of the conventional first-principles
calculation.
Figure 3(c) shows a proof test calculation for a H atom approaching to the smooth ESM [12].
In the region far from the smooth ESM, the H atom does not feel any interaction from the ESM
as expected. Whereas it approaches to the ESM, there is an attractive interaction and it turns to a
repulsive interaction according as the H atom penetrates to the smooth ESM region. The originof the attractive interaction is the dipole-dipole interaction between the H atom and the ESM. The

64

ESM

ESM

(d)

(c)

Fig. 3 (a) Discontinuous permittivity constant change in the original ESM method and (b) smoothtransition of the permittivity constant in the smooth ESM method. The permittivity changes fromunity to infinity. (c) The total energy profile of an H atom approaching the ESM region. (d) Charge
redistribution of the H atom. The color of lines indicates the position of the H atom in (c).

magenta and blue color lines in Fig.3 (d) shows there is an induced dipole moment on the H atom.
The green and red color lines indicate that an electron on the H atom is delocalized when the Hatom penetrates to the ESM region. The delocalization of the charge leads instability to the Hatom; consequently there is a repulsive interaction between the atom and the ESM. Since theelectronic wave function can penetrate to the smooth ESM region, we do not need to introduce
additional vacuum region.
4-3 Constant-μ method
It is well known that, to realize the constant Fermi energy (constant-μ) calculation, we cannot
fix the Fermi energy during the self-consistent-calculation (SCF) loop. This intuitive strategy will
fail within a few SCF steps. To realize stable and robust calculation, we employed the extended
system technique which is usually used in the constant temperature and constant pressure
molecular dynamics simulation. We assume that the charge in the simulation cell is a dynamicalvariable; we call this the fictitious charge particle (FCP) [13]. This means that the system isconnected to a charge reservoir and can keep the Fermi energy constant by exchanging the FCP
with the reservoir. By solving the dynamics equation for the FCP, we can realize grand canonicalensemble with the following grand potential:  Etot   nFCP ,where Etot, μ and nFCP indicate the total energy and the Fermi energy, and the charge of the FCP.
By taking the derivatives of the grand potential, we can calculate the force acting on atoms and
FCP.
As a proof test calculation of the constant-μ method, we calculate an H atom diffusion on
Al(001) surface [see Fig.4 (a)] by combining the constant-μ method with the nudged elastic band
(NEB) method [12]. Figure 4 (b) shows the energy profile and the Fermi energy variation of the
H diffusion with the conventional constant-N calculation. The activation barrier for the diffusion
is about 0.3 eV. During the diffusion, the H atom is slightly polarized by the interaction between
the surface, thus the Fermi energy of the system varies as the diffusion. Figure 4 (c) shows theenergy profile and nFCP along with the diffusion path under constant-μ condition (three different
target Fermi energies (μext) are shown). Although overall energy profiles depending on the target
potential are consistent with the constant-N calculation, the activation barriers and nFCP slightly

65

(b)

(c)

constant−N
ext=−4.5 v
ext=−5.0 v
ext=−4.0 v

0.5
0.4

 (eV)

0.4
E (eV)(a)constant−N

0.5

0.3
0.2

0.33
0.32
0.31
4

0.3
0.2

0.1

0.1

0.0

0.0
0.050

nFCP (e/surf)

Fermi energy (eV)

−3.5
−4.0
−4.5

0.025
0.000
−0.025−5.0
AB
Diffusion path

C

A

B

C

Diffusion path

Fig. 4 (a) H atom diffusion on Al(001) surface. A and C sites are the most stable bridge site, B site is
the saddle point of the path. (b) The energy profile and the variation of the Fermi energy along thepath with constant-N calculation. (c) The energy profile and the variation of the excess charge of thesurface along the path with constant-μ calculation. The origin of the nFCP is the neutral system.

depend on the target potentials. This means that application of a bias potential affects the diffusion
path and activation energy.
4-4 Summary
In this project, we have developed the smooth-ESM and the constant-μ methods to improve themodeling of the electrochemical system and successfully shown simple proof test calculations
using the state of the art first principles calculation program. It is expected that we can simulate
more realistic electrode-electrolyte interface system and give a deep insight into theelectrochemical system experiments.
5. First-principles core-level X-ray photoelectron spectroscopy calculation of defects insemiconductors (J. Yamauchi)
For several decades, the behaviors ofdefects in semiconductors have been
investigated actively owing to industrial aswell as fundamental interest. However,atomic-scale phenomena are not yet
sufficiently understood. Although the mostfundamental information at this scale is theatomic structure of defects, it is usually very
difficult to determine the structure. This ismainly for two reasons. One is that defects
are located in the substrate and cannot bedirectly detected. The other is the weakFig. 5 X-ray photoelectron spectroscopy.
signals in measurements. The number of
impurity atoms is very low compared withthat of matrix atoms. Furthermore, the ratio of the number of defects of interest to the total numberof impurity atoms is also very small. Recently, to compensate for the weak signals, detailed Xray photoelectron spectroscopy (XPS) measurements (Fig. 5) have been performed using highintensity X-ray beams generated by large synchrotron radiation facilities such as SPring-8.

66

Therefore, the XPS has been expected to enable highly accurate measurements than before.
Since the XPS signal itself does not uniquely determine the defect structure, its identification
requires theoretical analyses for the XPS binding energies of the suggested models with highaccuracy. On the other hand, there have been few reliable first-principles core-level XPS
calculations for impurity defects in semiconductors, because the local potential boundary
condition of defect model systems has not yet been sufficiently evaluated. To obtain chemical
shifts in XPS binding energy with high accuracy, it is necessary to take a sufficiently large
supercell for a defect so as not to affect the local potential at the boundary. We report the firstprinciples core-level XPS calculations on the defects in silicon (Si) and silicon carbide (SiC)
considering the careful evaluation of the boundary condition. The tests for the supercell size
dependence revealed that the 512- and 1728-atom supercells for Si and SiC are required to achieve
the accuracy of 0.1 eV for the XPS binding energy. The XPS binding energies are calculated bythe SCF method. To calculate the XPS binding energy, special pseudo-potentials are necessary
to deal with the core state, because usual pseudopotentials ignore the core state for the saving of
the computational cost. There are two types of such potentials. One is the screened core holepseudopotential (SCH PP), which is generated using the atomic reference removing the electronof the core state. The other is the pseudopotential with core (PPwC), which deals with the corestate explicitly. We adopt the PPwC for As 3d and the SCH PPs for the other dopant elements.
Boron has exotic features and is the most important p-type dopant element for Si. Boron atom
is so small that it tends to make three-center bonding and prefers triangular shape. In nature,polyhedral units of boron clusters such as octahedral B6, icosahedral B12, and cubo-octahedral B12,
are available in elements and compounds. Thus, it is expected that the polyhedral B clusters exist
in Si crystals. We calculated the XPS binding energies for the defect models, which are suggestedby theories and/or experiments. Among them, we plotted the binding energies of B 1s near theexperimentally observed spectra in Fig. 6. In this figure, the vertical lines indicate the calculated
values for icosahedral B12 (ICO), substitutional B, <001> B-Si and octahedral B6 (B6), and thecurves are observed and deconvoluted experimental spectra. This figure clearly shows that theexperimental three peaks correspond to the ICO, substitutional B, and <001>B-Si, respectively.

Fig. 6 B1s core-level XPS binding energy ofB clusters in Si. Experimental anddeconvoluted data are from Mizushima andcoworkers. The vertical lines are calculated
binding energies. The substitutional B isadopted for the origin of the energy.

Fig. 7 As 3d core-level XPS binding energy ofdefects in Si. The energy origins are the neutral
and charged substitutional As for (a) and (b),respectively.

67

These results are consistent for the other experimental properties such as electric activity [14].
In Fig. 7, the calculated XPS binding energies of As 3d for As defects in Si are indicated in theneutral state (a) and the charged state. In this system, experimentally observed peak is -1.2 eV
relative to that of substitutional As. Within the neutral state only (a), where the neutral
substitutional As is set to the energy origin, we cannot find the defect suitable for the experimental
one. Including the charged states, where the energy origin is the binding energy for thesubstitutional As+, it is found that the vacancy related defects are good candidates for theexperimentally observed defects.
To investigate the dependence of the XPS binding energy on the substrate semiconductors, we
also performed the XPS calculation for the various dopants (B, Al, P, As) defects in SiC. It isfound that the relative range of XPS binding energy in SiC is about 4 eV, which is much larger
than that of 1 eV in Si. This is because, while the Si is the elemental semiconductor, the SiC is
the compound semiconductor and the deviation of the local potential is much larger than that in
Si. The SiC has many polymorphs. We calculated the defects in 3C- and 4H-SiC, and found thatthe dependence of the XPS binding energy on the polymorphs is not large.
References[1] K. Nakamura, H. Kagaku, 29, 432 (2008)[2] M. Otani and O. Sugino, Phys. Rev. B, 73, 115407 (2006)[3] S. Grimme, S. Ehrlich and L. Goerigk, J. Chem. Phys, 132, 154104 (2010)[4] Atsushi Togo, Fumiyasu Oba, and Isao Tanaka, Phys. Rev. B, 78, 134106 (2008)[5] H.Z. Weingärtner, Phys. Chem. 132, 129 (1982)[6] K. Esfarjani and H.T. Stokes, Phys. Rev. B77, 144112 (2008); Phys. Rev. B86, 019904 (2012).
[7] T. Tadano, Y. Gohda and S. Tsuneyuki, J. Phys. Condens. Matter 26, 225402 (2014).
[8] T. Tadano, Y. Gohda, and S. Tsuneyuki, Phys. Rev. Lett. 114, 095501 (2015).
[9] M. Otani and O. Sugino, Phys. Rev. B 73, 115407 (2006).
[10] M. Bagge-Hansen, et al., Adv. Mater. 27, 1512 (2015).
[11] I. Hamada, O. Sugino, N. Bonnet, and M. Otani, Phys. Rev. B 88, 155427 (2013).
[12] We implemented the smooth-ESM method and the constant-μ method into the PWscf
program in the Quantum ESPRESSO package: http://www.quantum-espresso.org/.
[13] N. Bonnet, T. Morishita, O. Sugino, and M. Otani, Phys. Rev. Lett. 109, 266101 (2012).
[14] J. Yamauchi, Y. Yoshimoto, and Y. Suwa, Appl. Phys. Lett. 99 191901 (2011).

6. Collaborators
Takashi Nakayama (Chiba University), Osamu Sugino (The University of Tokyo), Yoshitada
Morikawa (Osaka University), Kazuto Akagi (Tohoku University), Yoshitaka Tateyama (NIMS),
Yuji Suwa (Hitachi-Co), Yoshihiro Gohda (Tokyo Institute of Technology/The University ofTokyo)
7. Publications
 T. Suzuki, Y. Yoshimoto, K. Yagyu and H. Tochihara, “Adsorption of PTCDA on Si(001)-2 x
1 surface”, J. Chem. Phys. 142, Art. No. 101904 (2015).
 K. Nakamura, Y. Yoshimoto and M. Imada, “Ab initio two-dimensional multiband low-energy
models of EtMe3Sb[Pd(dmit)2]2 and κ-(BEDT-TTF)2Cu(NCS)2 with comparisons to single-

68band models”,

Phys. Rev. B86, Art. No. 205117 (2012).

 J. Yamauchi, Y. Yoshimoto and Y. Suwa, “Identification of boron clusters in silicon crystal by
B1s core-level x-ray photoelectron spectroscopy: A first-principles study”, Appl. Phys. Lett.
99, Art. No. 191901 (2011).
 Y. Yoshimoto et al., Software Package “xTAPP”;
http://xtapp.cp.is.s.u-tokyo.ac.jp
http://ma.cms-initiative.jp/ja/listapps/xtapp
 T. Tadano, Y. Gohda, and S. Tsuneyuki, “Impact of Rattlers on Thermal Conductivity of a
Thermoelectric Clathrate: A First-Principles Study”, Phys. Rev. Lett. 114, 095501 (2015).
 T. Tadano, Y. Gohda and S. Tsuneyuki, “Anharmonic force constants extracted from firstprinciples molecular dynamics: applications to heat transfer simulations”, J. Phys. Condens.
Matter 26, 225402 (2014).
 T. Tadano, Software Package “ALAMODE”;
http://github.com/ttadano/alamode
http://alamode.readthedocs.org
http://ma.cms-initiative.jp/ja/listapps/alamode/alamode M. Kawamura, Y. Gohda and S. Tsuneyuki, “Improved tetrahedron method for the Brillouinzone integration applicable to response functions”, Phys. Rev. B 89, 094515 (2014).
 I. Hamada, O. Sugino, N. Bonnet, and M. Otani, “Improved modeling of electrified interfaces
using the effective screening medium method”, Phys. Rev. B 88, 155427 (2013).
 N. Bonnet, T. Morishita, O. Sugino, and M. Otani, “First-principles molecular dynamics at aconstant electrode potential”, Phys. Rev. Lett. 109, 266101 (2012).
 J. Yamauchi and Y. Yoshimoto, “X-ray Photoelectron Spectroscopy for the Boron Impuritiesin Silicon: a First-principles Study”, AIP Conf. Proc. 1399, 89 (2011).
 J. Yamauchi, Y. Yoshimoto, and Y. Suwa, “A first-principles core-level XPS study on the boron
impurities in germanium crystal”, AIP Conf. Proc. 1566, 41 (2013).
 H. Kishi, M. Miyazawa, N. Matsushima, and J. Yamauchi, “First-principles core-level X-ray
photoelectron spectroscopy calculation on arsenic defects in silicon crystal”, AIP Conf. Proc.
1583, 226 (2014).
 N. Matsushima and J. Yamauchi, “First Principles XPS Calculation for the B defects in SiC”,
JPS Conf. Proc. 1 012027 (2014).
8. Invited Talks
 Yoshihide Yoshimoto, “Development and application of plane wave based first-principles
electronic structure calculation code "xTAPP"”, 4th International Workshop on Massively
Parallel Programming Now in Quantum Chemistry and Physics - Toward post-K computers,
Nov. 23 - 24, 2014, University of Tokyo.
 S. Tsuneyuki, “Toward Accurate Calculation of Electronic and Structural Properties ofMaterials from First Principles”, The 9th General Meeting of ACCMS-VO, OIST, Okinawa,Dec. 20-22, 2014.
 S. Tsuneyuki, “Atomistic Modeling of Materials Based on First-Principles Calculation”, The
10th NOBUGS conference, KEK, Tsukuba, Sep. 24-26, 2014.

69

 S. Tsuneyuki, “Computer Experiments for Materials Science”, International Symposium on
Recent Trend of Interdisciplinary Research of Physics, Earth and Space Science, Osaka Univ.,
Dec. 17-18, 2013.


S. Tsuneyuki, “Atomistic Modeling of Materials Based on First-Principles Calculation”, theSeventh General Meeting of ACCMS-VO, Tohoku Univ., Sendai, Nov. 23-25, 2012.



M. Otani, “Computational study of structure, vibrational property and electrochemical reaction ona biased metal-water interface”, International Symposium of Quantum systems in chemistry andPhysics, September 11-14, 2011, Kanazawa, Ishikawa, Japan。.



M. Otani, “First-principles molecular dynamics study of electrochemical reactions at electrodeelectrolyte interfaces”, ICMR Workshop on Charged Systems and Solid/Liquid Interfaces, July 611, 2014, Santa Barbara, CA, USA.



M. Otani, “First-principles molecular dynamics simulations on electrochemical reactions using
effective screening medium method”, 248th ACS National Meeting & Exposition, August 10-14,
2014, San Francisco, CA, USA.



M. Otani, “First-principles molecular dynamics study of electrochemical reactions at electrodeelectrolyte interfaces”, APS March meeting 2015, March 2-6, 2015, San Antonio, TX, USA.

70

Nanoscale Non-Equilibrium Electric Transport Dynamics Based on
Density Functional Theory
S. Watanabe1, K. Watanabe2, S. Souma3, T. Ono4
The University of Tokyo, 2Tokyo University of Science, 3Kobe University,
4
Osaka University (Present affiliation: University of Tsukuba)

1

1. Introduction
Nanoscale electric transport and related phenomena, in particular dynamical processes
involving realistic complex atomic structures, non-equilibrium electron/ion movements underapplied high local electric field, etc. are important in the light of both fundamental science andpractical applications. Our understanding on them, however, is still insufficient.
Our team has been conducting computational studies aiming at deeper understanding on
nanoscale electric transport phenomena, especially transport dynamics of various aspects. Our
research subjects include the analyses of interplay of electric transport with phonon/heat, ion
transport, spin, electronic excitation, etc., analyses on models more realistic than thoseexamined previously, and development of computational methodologies and codes for these
purposes. In this report, we describe the highlights of our research activities.
2. Overview
As mentioned above, our researches in this project include various topics. In this article, wewill briefly describe our achievement for six selected topics.
The first topic is the AC transport properties. From our study on the subterahertz AC responseof metallic carbon nanotubes using the nonequilibrium Green’s function (NEGF) method, we
have found many interesting features of the AC response, especially in the case of nanotubes
with an atomic vacancy. Combining the results for the nanotubes and single molecular bridges,
we have got a unified picture of the AC response of nanostructures.
The second one is the transient transport behavior. Since in this case we need much heavier
computation than in the cases of DC and AC transport behaviors, we have concentrated on very
simple models, such as a quantum dot with a single level between electrodes. On the other hand,
we have examined the effects of various factors such as the switching speed of the applied bias
voltage, electron-phonon interaction, and local heating.
The third one is the characteristics of graphene-based novel nanoscale devices. Firstly, we
have examined transport properties of graphene nanoribbons. More specifically, we haveexamined the effects of strain induced pseudo magnetic field on the locally strained graphene
channel field-effect-transistor operation quantitatively. Next, we have proposed a new type of
the pure spin current generation device that utilizes the quantum pumping effect and the edge
spin polarization in zigzag edged graphene nanoribbons.
The fourth one is the issues related to tantalum oxide based resistive switching devices. In
this device, the formation of conduction filaments is crucial in the device operation. So we haveexamined the features of the conduction filaments and ion movements.
The fifth one is the light-induced electronic excitations and dynamics of nanostructures.
Using sophisticated methodologies such as the Bethe-Salpeter equation and time-dependent
density functional theory (TDDFT), we have examined optical activation of dark excitons incarbon nanotubes, laser-assisted field emission from graphene nanoribbons and silicene
nanoribbons, excited-state nuclear forces and femtosecond laser reaction dynamics of a71

molecule, and optical properties of graphene nanoribbons and boron nitride nanoribbons.
Finally, a few methodologies and codes have been developed. The first one is a novel but
simple approach for the consideration of bias voltage within the Kohn-Sham formalism ofdensity functional theory. The second one is TDDFT methods to efficiently calculatesecond-order nonadiabatic couplings (NAC). For the computational code, we have developed
first-principles electronic-structure and transport calculation code based on real-space
finite-difference method. In the code, the problem of heavy computation due to the direct
inversion of large matrices to obtain the surface Green’s functions and self-energy terms ofelectrodes has been overcome by drastically reducing the matrix size using novel algorithm.
3. Research Results(Note: The numbers of cited references correspond to those in “5.[1](b) Selected publications.”)[1] AC Transport Properties
Considering high potential of metallic carbon nanotubes (M-CNTs) as interconnects in next
generation electronic devices, we have examined the subterahertz AC response of MCNTs using
the nonequilibrium Green’s function (NEGF) method. We have found several interesting
features of AC response, especially for MCNTs with an
atomic vacancy. First, using the NEGF method
combined with a simple tight-binding model, we havefound that the susceptance (exactly speaking, emittanceis examined instead of susceptance) does depend on theposition of the vacancy, which makes a strong contrast to
the position independence of the DC conductance (seeFig. 1).1) Next, using the NEGF method based on thedensity functional theory (DFT), we have found that thetwo states induced by the vacancy, the π-orbital-origin
vacancy state and the σ-orbital-origin dangling-bondstate, show different susceptance behaviors at firstsight.2) That is, besides capacitive peaks seen in both of Fig. 1. Emittance and DC conductancethe two states, satellite inductive peaks are clearly seen (inset) of a metallic carbon nanotube
adjacent to the main capacitive peak only around the with an atomic vacancy. Here, ddenotes the distance from the center ofdangling-bond state. Our analysis with a simple resonant the tube between the electrodescattering model reveals that the origin and the regions.1)
magnitude of these satellite inductive peaks can beunderstood by just one parameter, i.e., the lifetime of electrons at a defect state.
We have also examined AC transport properties of 1,4-Benzenedithiol molecule sandwiched
between semi-infinite one-dimensional gold electrodes using the NEGF based on DFT. In thecalculated susceptance spectrum, we have found two types of features; one corresponding toconductance peaks/dips and the other not related to those. We have found that the former can beunderstood similarly to the case of the above MCNTs, while the latter, which is seen at theenergies having small conductances, correlates with the DOS of the electrodes. This correlation
with the DOS is considered to be the same as the one seen in mesoscopic capacitors.
[2] Transient Transport Properties
As well as the AC transport behaviors, the transient current dynamics is also important inunderstanding the current behaviors under practical operations of novel nanoscale electronic
devices. We have examined the transient current dynamics using simple models and the NEGF
method.

72

First, we examined the transient current dynamics in a resonant tunneling model having asingle energy level with considering the switching-speed of the applied bias voltage.3) We foundthat with the decrease of the switching speed, the amplitude of the current overshoot decreases.
We also found that the crucial factor determining the relaxation time changes from the dwell
time of an electron in the resonant level to the voltage rise time.
Next, we have examined the effects of electron-phonon interaction and local heating on thetransient currents in a single-level quantum dot system weakly coupled to a localized vibration
degree of freedom using the NEGF method
under the wide-band-limit and lowest-orderexpansion approximations.4) We calculated
the time dependent phonon number using aphenomenological method employing thetime-dependent power transfer. We havefound that the dynamical behaviors in thissystem can be grouped into two responses: if
the energy corresponding to the applied bias
voltage eV is smaller than or equal to theFig. 2. Schematic illustration of the model structure
phonon energy Eph, first the system of graphene field effect transistor. The single layer
temperature T increases and then decreases; graphene is stretched along the armchair direction
alternatively, if eV is larger than Eph, T only over the gated channel region. The left and theincreases monotonically until a steady state right unstrained semi-infinite graphene regions are
assumed to be doped into n-type or p-type, while
is reached. The total electrical current is
the central strained graphene channel region issuppressed by phonon heating, while the applied by the gate voltage
heat transferring between dot andenvironment mitigates the effects of suchheating.
[3] Characteristics of Graphene Based
Novel Devices(a) Strained graphene field effect transistor
incorporating the pseudo magnetic field
effectRecently, the possibility of using lateralstrain to control the electrical conduction insingle layer graphene (SLG) has been
attracting great attention. Experimentally,such strain can be applied by bending or
stretching SLG on flexible substrate. From
the theoretical view point, it is important tonote that there are two distinct mechanisms
in the strain induced modulations ofelectronic properties in SLG. The first one isthat caused by the pseudo magnetic field at
the interface between unstrained and strained
regions of SLG. The second one is thatoriginated from the strain induced bandgap
opening. The former idea stimulated various
theoretical and experimental studies such as

Fig. 3．Gate voltage dependence of the drain current
density for two different doping configurations(top: p-i-p, and bottom: p-i-n) and three different
strain ratios.

73the quantum Hall effect and valley filtering device, and switching device. In the latter case,shear type strain is especially important since it allows SLG to possess the bandgap with theexperimentally achievable strain ratio. While the effect of such strain induced bandgap opening
in the strained graphene channel FET operation is similar to that in conventional semiconductor
MOSFET, the possible application of the former idea (pseudo magnetic field effect) to the FET
device operation has not been quantitatively explored theoretically so far. With such motivation,
we study numerically the performance of locally strained graphene based FET (Fig. 2) paying
special attention to the role of the pseudo magnetic field, on/off ratio of the drain current, andthe sub-threshold slope at the room temperature. By employing the non-equilibrium Green’s
function and the atomistic tight-binding methods, we have numerically studied the performanceof strained graphene based field effect transistors, where the strain around 8% is applied overthe channel region of graphene along the armchair direction to induce the pseudo vector
potential effect. As shown in Fig. 3, we have predicted that the gate voltage applied to thecentral strained graphene region can switch the drain current on and off with high on/off ratio
more than six orders of magnitude at the room temperature in spite of the absence of the
bandgap in the strained channel region. Steeper sub-threshold slope below 60 mV/decade is also
predicted to be achieved at the room temperature due to the mechanism similar to theband-to-band tunneling FET.5)(b) Graphene based pure spin current generator based on quantum pumping effectProactive use of the spin degree of
freedom in semiconductor is one of thepromising pathway for further developmentof the semiconductor technology. One of themost important subjects in such rapidly
developing field of spintronics is theefficient generation or injection of the spins Fig. 4. Illustration of ZGNR-based pure spininsemiconductors.
Especially,
the generator. Time-oscillating potentials are applied togeneration of pure spin current (i.e., four locally gated regions independently with eachspin-flow without charge current) is other.
important since the pure spin current does
not accompany the Joule heating and thus
enable us to transmit information at very
low energy consumption. So far various
types of devices have been proposed anddeveloped to generate the pure spin current,
including the spin Hall effect, the adiabatic
pumping in graphene, and the one using the
microwave induced magnetization dynamics.
We propose another type of the pure spincurrent generation device that utilizes thequantum pumping effect and the edge spin Fig. 5．Pumping phase dependence of the spinpolarization in zigzag edged graphene dependent pumped charges, where the pumping
nanoribbons (ZGNRs).6) ZGNR is known to phases (see Fig. 4) are set as φ1=φ2=0, φ3=φ, andshow the partly flat electronic band structure φ4= −φ, and the value of φ is varied.
near the Fermi energy, and thecorresponding wavefunction is localized near the edges of the ZGNR. Since such edge localizedstates give high density of states at the Fermi energy, more precise calculations taking into74

account the spin dependence of the electron-electron interaction predict the appearance of theedge spin polarization with opposite spins along the opposite edges. Moreover it has also been
predicted that the half metallic band structure can be induced by applying the transverse electricfield. Recent experiment has demonstrated the presence of edge spin polarization in ZGNR bymeans of the scanning tunneling spectroscopy. Our proposal is based on the idea to utilize suchedge spin polarization in ZGNR by means of the adiabatic quantum pumping effect to generate
the pure spin current. Quantum pumping effect is the idea to generate a dc current by applying
the local time periodic potentials to the conductor without applying the dc bias.
By employing the spin-density-functional tight-binding method and the adiabatic quantum
pumping theory, we have shown theoretically that the pure spin current can be generated inzigzag edged graphene nanoribbons (ZGNRs) through the adiabatic pumping by edge selective
pumping potentials (Fig. 4). The origin of such pure spin current is the spin splitting of the edge
localized states, which are oppositely spin polarized at opposite edges. In the proposed device,each edge of the ZGNR is covered by two independent time periodic local pumping gate
potentials with a definite phase difference, inducing the edge localized spin polarized current.
By making the pumping phase difference opposite between two edges, the total pumped charge
currents are canceled and then generates the pure edge spin current (Fig. 5). Such edge selective
time-dependent local potential assumed in the proposed device should be able to be applied bymeans of the scanning gate microscopy. The proposed method of generating pure spin current isuseful in the development of graphene spintronics.
[4] Conductive Filament Structures and Ion Movements in Tantalum-Oxide-Based
Resistive Switching Devices
Oxide-based resistive switching devices haveattracted much attention because of its potential
as next-generation memory cells. In particular,
amorphous-tantalum-oxide based devices show
good endurance and switching performances.
Their switching behavior is understood based on
the formation/rupture of conductive filaments(CFs), but, their atomistic structures andmicroscopic formation/rupture processes have notbeen well clarified yet.
We have examined the conduction paths inthe tantalum-oxide-based atomic switch using the
DFT calculations. In the Cu/amourphous-TaO2.5(a-TaO2.5)/Pt heterostructure, where Cu is
believed to be responsible for the CFs, Cu
single-atom chains have been found unstable. On
the other hand, Cu nanowires with widths larger
than the single atom have been found stable
during the molecular dynamics (MD) simulation
at 500 K for 15 ps, and can form conductive paths
Fig. 6. Atomic structures of Cu/a-TaO2.5/Pt(see Fig. 6).7) In the Pt/amorphous-TaOx/Pt
and local density of states near the Fermi
heterostructures, we have found the stronglevel. From top to bottom: the case with nocorrelation among Ta/O coordination number,
Cu wire, with a disconnected wire, with a
Ta–Ta/Ta–O bond lengths, and O concentration,
continued wire, and with a thicker continued
wire in TaO2.5.7)
and with the decrease of O concentration, Ta

75

atoms tend to merge together and finally form a continued Ta-rich region. This suggests that not
O vacancies, but the Ta–Ta bonding mainly contributes to the CF.8)
We have also examined the ion migration behaviors in TaOx using ab initio molecular
dynamics and nudged elastic band method. For the oxygen diffusion coefficient in TaOx, weshow that it strongly depends on x. We also show that at low O content (i.e., x = 1), the diffusion
coefficient of Ta becomes only ~0.5 order lower than that of O at room temperature. So the
diffusion of Ta cannot be ignored at low O content and at high temperature. For Cu at grainboundaries in TaOx, our calculations using surface models instead of grain boundaries suggestthat the presence of water strongly promotes its diffusion.
[5] Light-Induced Electronic Excitations and Dynamics of Nanostructures(a) Optical activation of dark excitons in carbon nanotubes
The optical properties of CNTs with hydrogen adatoms were investigated theoretically (Fig. 7).
We have shown that adsorbed hydrogen atoms make the triplet-exciton states optically active.
This effect is due to enhancement of the spin-orbit interaction as a result of distortion of sp2
orbitals so that they become more like sp3 orbitals. We
investigate the effects of spin-flipping on exciton statesby solving the Bethe-Salpeter equation under thetight-binding approximation. The absorption spectrum
was calculated using the solution of the Bethe-Salpeter
equation. It shows a satellite peak that is attributed totriplet excitons. Most importantly, we determined themicroscopic mechanism for the optical activation of
spin-triplet exciton states by adsorbed impurities, which Fig. 7. Photo absorption9)of H-adsorbed
offers a solution to the recent issue regarding the satellite (20,0) carbon nanotube.
peak observed in photo luminescence measurements.9)(b) Laser-assisted field emission from graphene nanoribbons and silicene nanoribbons
The energetics and dynamics of LAFE from nanoribbons were clarified by thetime-dependent density functional theory (TDDFT) simulations (Fig. 8). The emission
mechanism, which is common to the SiNR 10) and GNR 11), is found to be governed by threefactors: the electronic dipole matrix element, the orbital characteristics of the excited states, andthe energy levels of the excited states relative to the potential hump. However, the emission
currents from the SiNR are found to be much larger than from the GNR for the same laserparameters in spite of the larger work
function of the SiNR. This can be interpretedby the higher density of states, or morechannels for dipole transition in SiNRs thanthat of GNRs. We demonstrated the
dynamical correlation among the laser pulse,
the Kohn-Sham potential and emission
currents wiggling under the laser field in realFig. 8. Time evolution of emission currents
time and space on an atomic scale, andfrom SiNR at different points.10)
observed emission current being driven backand forth in the early stage of emission.
(c) Excited-state nuclear forces and femtosecond laser reaction dynamics of a molecule
We developed a simple and computationally efficient method to calculate excited-state76

nuclear forces on adiabatic potential-energy surfaces(APES) from linear-response TDDFT within areal-space framework.12) The Casida ansatz, which has
been validated for computing first-order nonadiabatic
couplings in previous studies, was applied to thecalculation of the excited-state forces. Our method isvalidated by the consistency of results in the lower
excited states, which reproduce well those obtained by
the numerical derivative of each APES. The
femtosecond laser reaction dynamics of biphenyl
molecule was investigated using TDDFT-MD
Fig. 9. Time evolution of torsional
simulation.13) The results (see Fig. 9) confirm that theangle of a biphenyl molecule.13)
probe pulse triggers a short-duration Coulomb explosion
and the kick pulse induces torsional motion of the two phenyl rings for a long duration by
selection of the laser parameters. The Coulomb explosion dynamics simulation confirms that the
F and Br atoms dissociate to the ion detector while maintaining their initial alignment withrespect to the phenyl rings, which is the fundamental basis for Coulomb explosion imaging ofmolecular torsion. Since the TDDFT-MD method takes account of the electric dipole-laserinteraction without assumptions, the simulation enables quantitative understanding oflaser-induced ultrafast molecular dynamics at a microscopic level.
(d) Optical properties of graphene nanoribbons and boron nitride nanoribbons
Dielectric function of BNRs and GNRs were calculated by the TDDFT.14) Dimensionality
dependence of dielectric functions has been clarified,indicating the signatures of the transition of plasmon excitaion
to particle-hole excitation in lower dimensions. We found an
apparent difference in low-energy absorption spectra between
BNRs and GNRs, reflecting their energy-band structures. We
carried out a systematic study on the quasiparticle band gapsof armchair boron nitride nanoribbons (ABNNRs) as afunction of the ribbon width, using a first-principles
many-body approach based on the GW approximation.15)
Because of the quasi-one-dimensional nature of ABNNRs, theenhanced Coulomb interaction effects and the reduced
screening greatly influence the quasiparticle band gaps, giving
rise to the GW correction in a range 2.5~3.0 eV (see Fig. 10).
Our GW calculations validate the qualitative correctness of the Fig. 10. Width-dependent GW
quasiparticle band gaps of
DFT study on BNNRs, together with the demonstration of the
H-terminated armchair boron
necessity of many-body approaches for quantitative nitride nanoribbons. 15)
understanding of quasiparticle band gaps.
[6] Developments of Methodologies and Codes and Their applications(a) Orbital-separation approach for considering bias voltage
We have developed a novel but simple approach, named orbital-separation approach, for theconsideration of bias voltage within the Kohn-Sham formalism of DFT.16) In the approach, theKohn-Sham orbitals around the Fermi level are separated into anode or cathode parts, anddifferent Fermi levels are applied in the determination of occupation numbers. This method canbe readily implemented in conventional first-principles codes.

77

We have applied this method to theelectronic charging of Au/vacuum/Au andAu/MgO/Au capacitors.16) Our results forvarious dielectric thicknesses have revealed
that the bulk optical and static dielectric
constants can be obtained accurately. We havealso performed finite-bias calculations on the
SrRuO3/BaTiO3/SrRuO3(SRO/BTO/SRO)
capacitors with 9 and 4 BTO unit cells. Figure
11 shows the calculated local inverse
permittivity profiles. In both systems, theinverse permittivity inside SRO converges tozero indicating a metallic response. On theother hand, a decrease in the permittivity isobserved at the SRO/BTO interfaces (deadlayer effect). When comparing the region
inside the BTO film, a clear difference is seen:Fig. 11. Inverse permittivity profiles of the SRO
for the 9-layer BTO, the permittivity is /BTO/SRO capacitors with (a) 9 BTO unit cells
positive as usual, but for the 4-layer BTO, theand (b) 4 BTO unit cells.
permittivity is found to be negative. Thisconfirms that ultrathin FE films with thickness of the order of a few unit cells exhibit negative
capacitance, at least if ferroelectric domain formation is suppressed. Since a layer of a few unit
cells is too thin to suppress tunneling current to a manageable amount in actual devices, we havebeen exploring the possibility of capacitance boost in the metal/PE/metal by sandwiching BTO
at the metal/PE interface.
(b) TDDFT methods for efficient calculation of nonadiabatic couplings (NAC)
We proposed a TDDFT method to efficiently calculate second-order NAC, using the Casida
ansatz adapted for the computation of first-order NAC. Calculations in the immediate vicinity ofJahn-Teller and Renner-Teller intersections showed that the results agreed well with theprediction from the Jahn-Teller/Renner-Teller models.17) We also presented extensive TDDFT
calculations on the quantization behaviors of first-order NAC in polyatomic systems and theangular NAC was found to show an oscillatory feature in the small-displacement limit (see Fig.
12), instead of the quantized value of 0.5 in the D3h symmetry.18) We further presented an
improved method for accurately calculating NAC from TDDFT, by the analytic evaluation of
the Kohn-Sham derivative matrix elements on the basis of density functional perturbation
theory.19) We elucidated the fact that the use of fractional
occupation, while being a natural extension for getting the
ground-state total energies or excitation energies, is not sostraightforward for getting NAC.20) Contrary to theconjecture that the Tamm-Dancoff (TDA) approximation
might cause the NAC results to deteriorate and violate thesum rule, the surprisingly good performance of the TDA
for NAC was revealed and attributed to that the TDA canpartially compensate for the local-density-approximation Fig. 12. Oscillating behavior oferror and give better excitation energies in the vicinity of angular nonadiabatic18) couplings inpolyatomic systems.
intersections of potential energy surfaces.21)78

(c) First-principles electronic-structure and transport calculation code based on real-space
finite-difference method
The nonequilibrium Green’s function (NEGF) approach has been used extensively inconnection with tight-binding models and first-principles methods employing localized basis
sets. However, the NEGF method has not been combined with the real-space formalism, whichis advantageous to massively parallel computer architectures and thus enables us to carry out
large scale transport simulations. This is because the direct inversion of Nx×Ny×Nz-dimensional
matrices to obtain the surface Green’s functions and self-energy terms of electrodes iscomputationally demanding due to the large numbers of grid points. Here, Nx, Ny, and Nz are
the numbers of grid points along the x-, y-, and z-directions in the unit cell of the electrodes andelectrons flow along the z-direction. The overbridging boundary-matching (OBM) method has
been employed with real-space approaches. In the formalism of the OBM method, the effect ofelectrodes is treated by the matching formula using Nx×Ny×Nf-dimensional ratio matrices,where Nf is the order of the finite-difference approximation for the kinetic energy operator in the
Kohn-Sham equation and is much smaller than Nz (Nf«Nz).
We have developed the real-space finite-difference (RSFD) NEGF scheme, in which thesurface Green’s functions and self-energy terms of the electrodes required in the NEGF method
can be obtained from the ratio matrices in the OBM method.22-26) The dimension of the matrices
for the surface Green’s functions and self-energy terms, which are directly inverted, issignificantly reduced from Nx×Ny×Nz to Nx×Ny×Nf in the present scheme while keeping the
rigorousness of the mathematical formulation. In order to demonstrate the applicability of the
RSFD NEGF method and the importance of the interpretation using scattering wave functions,the transport property of the C/BNNT where one carbon ring of (9,0) CNT is replaced with a
BN ring is examined.23,27,28) Inset of Fig. 13 shows the computational model. The central
finite-difference case, i.e., the case of Nf=1, is used. To determine the Kohn-Sham effective
potential, we use a conventional supercell under a periodic boundary condition in all directions
with a real-space grid spacing of ~ 0.24 Å; the dimensions of the supercell are Lx=13.34 Å,
Ly=13.34 Å, and Lz=4.32 Å for the electrodes and Lx=13.34 Å, Ly=13.34 Å, and Lz=8.64 Å forthe transition region, where Lx and Ly are the lateral lengths of the supercell in the x- andy-directions perpendicular to thenanotube axis, respectively, and Lz is
the length in the z-direction. The sizeof the matrices for the self-energy
terms is reduced from 75264 to 3136
using the present method.
The conductance spectrum of
the C/BNNT is shown in Fig. 13. We
have also calculated the conductancespectrum and channel transmissions
by the OBM method and found thatthese perfectly agree with the results
obtained by the present method
except for rounding error ofnumerical calculation. From Fig. 13, Fig. 13. Channel transmissions as functions of energy ofincident electrons. (Inset) Computational model where
we found that two channels actually
one BN zigzag ring is sandwiched between (9,0) CNT
contribute to the transport in theelectrodes. Large dark, small dark, and small light balls
are N, C, and B atoms, respectively. Reprinted withvicinity of the Fermi level and no
permission from Ref. 23).
significant peaks due to the resonant

79

tunneling through the edge states are observed in the conductance spectrum. The wave function
matching at the interface reveals that the rotational symmetry of wave functions with respect tothe tube axis plays an important role in electron transport. Since the states coming from andgoing to electrodes show threefold rotational symmetry, the states in the vicinity of the Fermi
level, whose wave function exhibits fivefold symmetry, do not contribute to the electrontransport through the BN ring. This result indicates that the real-space picture of the scatteringwave function, which is not necessary to be taken into account to calculate the transport
property in the NEGF formalism, helps us to interpret transport phenomena.
The spatial maps of the local density of states obtained by scanning tunneling spectroscopy(STS) can give us the images of standing waves, which provide important information about thedispersion relation of the electron scattering process at the potential barrier. By employing the
newly developed computational code in this project, we have carried out first-principles
calculation for obtaining the scattering potential of the defects on Ge(001) surfaces using acomputational model shown in Fig. 14.29) The Ge-Ge dimer at the center of the scattering region
is replaced by Ge-Si or Ge-Sn dimer as a defect. When the Si and Ge atoms are located at thelower and upper sides of the dimer, respectively, the dimer is referred as an SiL dimer. Other
dimers are named in a similar manner. By calculating the scattering wave functions using theOBM method, the standing waves in the spatial map of the local density of states are examined;
the waves correspond to the image of the differential conductance obtained by STS.
The charge density distribution of the scattering wave functions for the electrons
propagating from both the left and right electrodes is investigated with an energy of EF+0.55 eV.
To examine the period and phaseshift of the standing waves, thedensity of the standing waves is fittedby α(x)=Acos(2kxx+φ), where A, φ,and kx are the amplitude of thestanding waves, phase shift of thestanding waves, and Bloch vector inx direction, respectively. The fitting
is carried out using the density of thestanding waves above the lower
atoms of the dimers. The period ofFig. 14. Schematic image of computational model. Dark
the standing wave and its phase shiftspheres are Ge atoms and light spheres represent theare shown in Table I, in which they
dimer replaced by a Ge-Si or Ge-Sn dimer, which isagree with those obtained by theindicated by the arrow. Atoms are denoted by large andexperiment using STS. The scattering
small spheres according to the distance from the surface.
Reprinted from Ref. 29).
potential acts as a barrier when theelectronegativity of the upper atom
Table I. Amplitude and phase shift of the standing wavesof the dimer is larger than that of thearound the Ge-Si and Ge-Sn dimers.
lower atom, while it becomes a wellModelin the opposite case. The scattering
SiL
SiU
SnL
SnU
potential is related to the stabilization
Amplitude
3.908 1.692 7.169 1.031
of the π bands of the Ge(001) surface-7
3(10 e/Å /eV/spin)
due
to
the
differenceinPhase shift
0.221 -0.602 -0.650 0.142
electronegativity between Ge and the(πrad)
impurity.

80

4. Members(Affiliation at (the beginning, if later moved to anywhere else, and) the end of the period when the
person participated in our team is shown.)[1] Leader and principal members
Satoshi Watanabe (Professor, Graduate School of Engineering, The University of Tokyo)
Kazuyuki Watanabe (Professor, Faculty of Science, Tokyo University of Science)
Satofumi Souma (Associate professor, Graduate School of Engineering, Kobe University)
Tomoya Ono (Assistant professor, Graduate School of Engineering, Osaka University⇒ Associate professor, Center for Computational Sciences, University of Tsukuba)[2] Collaborative members
Akira Sakai (Professor, Graduate School of Engineering, Kyoto University)
Tomofumi Tada (Assistant professor, Graduate School of Engineering, The University of Tokyo⇒ Associate professor, Materials Research Center for Element Strategy, Tokyo Institute ofTechnology)
Takahiro Yamamoto (Assistant professor, Graduate School of Engineering, The University ofTokyo ⇒ Junior associate professor, School of Engineering, Tokyo University of Science)
Satoru Konabe (Postdoctoral fellow, Graduate School of Pure and Applied Science, Universityof Tsukuba)
Chunping Hu (Assistant professor, Faculty of Science, Tokyo University of Science)
Yasunobu Ando (Assistant professor, Graduate School of Engineering, The University ofTokyo)
Kanji Sasaoka (Assistant professor, Organization of Advanced Science and Technology, Kobe
University)
5. Publications and Invited Talks[1] Publications(a) Total number: 72 (English: 69, Japanese: 3)(b) Selected publications
1) D. Hirai, T. Yamamoto, and S. Watanabe, “Theoretical Analysis of AC Transport in
Carbon Nanotubes with a Single Atomic Vacancy: Sharp Contrast between DC and AC
Responses in Vacancy Position Dependence”, Appl. Phys. Exp. 4, 075103 (3 pages)(2011).
2) D. Hirai, T. Yamamoto, and S. Watanabe, “Anomalous satellite inductive peaks inalternating current response of defective carbon nanotubes”, J. Appl. Phys. 115, 174312(6 pages) (2014).
3) K. Sasaoka, T. Yamamoto, and S. Watanabe, “Wavelet Analysis of Quantum Transient
Transport in a Quantum Dot”, Appl. Phys. Lett. 102, 233107 (4 pages) (2013).
4) W. Liu, K. Sasaoka, T. Yamamoto, T. Tada, and S. Watanabe, “Inelastic transient
electrical currents and phonon heating in a single-level quantum dot system”, J. Appl.
Phys. 113, 123701 (9 pages) (2013).
5) S. Souma, M. Ueyama, and M. Ogawa, “Simulation based design of strained graphene
FET incorporating the pseudo magnetic field effect”, Appl. Phys. Lett., 104, 183103(2014).
6) S. Souma and M. Ogawa, “Pure spin current induced by adiabatic quantum pumping inzigzag-edged graphene nanoribbons”, Appl. Phys. Lett. 104, 213505 (2014).
7) B. Xiao, T. Gu, T. Tada and S. Watanabe, “Conduction paths in Cu/amorphous-Ta2O5/Pt
atomic switch: First-principles studies”, J. Appl. Phys. 115, 034503 (7 pages), 2014.

81

8)

B. Xiao and S. Watanabe, “Oxygen vacancy effects on an amorphous-TaOx-based
resistance switch: a first principles study”, Nanoscale 6, 10169-10178 (2014).
9) S. Konabe and K. Watanabe, “Mechanism for optical activation of dark spin-triplet
excitations in hydrogenated single-walled carbon nanotubes”, Phys. Rev. B 83, 145407(2011).
10) T. Higuchi, C. Hu, and K. Watanabe, “Energetics and Dynamics of Laser-Assisted Field
Emission from Silicene Nanoribbons: Time-Dependent First-Principles Study”, e-J. Surf.
Sci. Nanotech. 13, 115 (2015).
11) D. Utsugi, C. Hu, and K. Watanabe, “Laser-Driven Field Emission from Graphene
Nanoribbons: Time-Dependent Density-Functional Theory Simulations”, Appl. Phys.
Express 5, 105101 (2012).
12) J. Haruyama, T. Suzuki, C. Hu, and K. Watanabe, “Excited-state nuclear forces on
adiabatic potential-energy surfaces by time-dependent density-functional theory”, Phys.
Rev. A 85, 012516 (2012).
13) J. Haruyama, C. Hu, and K. Watanabe, “First-principles molecular-dynamics simulationof biphenyl under strong laser pulses by time-dependent density-functional theory”, Phys.
Rev. A 85, 062511 (2012).
14) S. Hagiwara, H. Goto, C. Hu, and K. Watanabe, “Optical Properties of Boron Nitride andGraphene Nanoribbons: A Time Dependent Density Functional Theory Simulation”, JPS
Conf. Proc. 1, 012072 (2014).
15) C. Hu, R. Ogura, N. Onoda, S. Konabe, and K. Watanabe, “Quasiparticle band gaps ofboron nitride nanoribbons”, Phys. Rev. B 85, 245420 (2012).
16) S. Kasamatsu, S. Watanabe, and S. Han, “Orbital-separation approach for considerationof finite electric bias within density-functional total-energy formalism”, Phys. Rev. B 84,
085120 (11 pages) (2011) (Editor’s suggestion).
17) C. Hu, O. Sugino, and K. Watanabe, “Second-order nonadiabatic couplings fromtime-dependent density functional theory: Evaluation in the immediate vicinity ofJahn-Teller/Renner-Teller intersections”, J. Chem. Phys. 135, 074101 (2011).
18) C. Hu, R. Komakura, Z. Li, and K. Watanabe, “TDDFT Study on quantization behaviors
of nonadiabatic couplings in polyatomic systems”, Int. J. Quant. Chem. 113, 263 (2013).
19) C. Hu, T. Tsukagoshi, O. Sugino, and K. Watanabe, “Nonadiabatic couplings fromtime-dependent density functional theory: Formulation by the Kohn-Sham derivative
matrix within density functional perturbation theory”, Phys. Rev. B 87, 035421 (2013).
20) C. Hu, O. Sugino, H. Hirai, and Y. Tateyama, “Reply to Comment on ‘Nonadiabatic
couplings from the Kohn-Sham derivative matrix: Formulation by time-dependent
density-functional theory and evaluation in the pseudopotential framework’ ”, Phys. Rev.
A 88, 056502 (2013).
21) C. Hu, O. Sugino and K. Watanabe, “Performance of Tamm-Dancoff approximation on
nonadiabatic couplings by time-dependent density functional theory”, J. Chem. Phys.
140, 054106 (2014).
22) T. Ono, S. Tsukamoto, Y. Egami, and Y. Fujimoto, “Real-space calculations for electrontransport properties of nanostructures”, J. Phys.: Condens. Matt. 23, 394203 (13 pages)(2011).
23) T. Ono, Y. Egami, and K. Hirose, “First-principles transport calculation method based onreal-space finite-difference nonequilibrium Green’s function scheme”, Phys. Rev. B 86,
195406 (13 pages) (2012).
24) T. Ono, Y. Fujimoto and S. Tsukamoto, “First-principles calculation methods for
obtaining scattering waves to investigate transport properties of nanostructures”, Quant.

82

Matt. 1, 4 (16 pages) (2012).
25) M. Heide and T. Ono, “Convergence of the Broyden density mixing method innoncollinear magnetic systems”, J. Phys. Soc. Jpn. 82, 114706 (10 pages) (2013).
26) T. Ono, “First-principles electronic-structure and transport calculation formalism using
real-space grid based method,” AIP Conference Proceedings 1558, 1524 (4 pages)(2013).
27) T. Kojima and T. Ono, “First-principles study on transport property of peapods,” Curr.
Appl. Phys. 12, S100 (5 pages) (2012).
28) H. D. Nguyen and T. Ono, “Doping Effect on Magnetism and Transport Property ofHeterojunction between Carbon and Boron Nitride Nanotubes”, J. Phys. Chem. C 117,
24115 (6 pages) (2013).
29) T. Ono, “First-principles calculation of scattering potentials of Si-Ge and Sn-Ge dimers
on Ge(001) surfaces”, Phys. Rev. B 87, 085311 (6 pages) (2013).
30) S. Souma, “Atomistic tight-binding simulation of spin-orbit coupled semiconductor
devices”, in “Spintronics in Nanoscale Devices” edited by E. R. Hedin and Y. S. Joe (Pan
Stanford Publishing), 101 (20 pages) (2013).
31) S. Souma, Y. Ohmi, and M. Ogawa, “Effect of lateral strain on gate induced control ofelectrical conduction in single layer graphene device”, J. Comput. Electron. 12, 170,(2013).
32) A. Mehdipour, K. Sasaoka, and M. Ogawa, and S. Souma, “Strain-induced modulationof anisotropic photoconductivity in graphene”, Jpn. J. Appl. Phys. 53, 115103 (2014).
[2] Invited talks(a) Total number: 35 (international: 25, domestic: 10)(b) Selected invited talks
1) S. Watanabe, T. Tada, S. Kasamatsu and T. K. Gu, “Nanoscale simulations on electronic
and atomic conduction in solid electrolyte/metal junction systems”, IUMRS 11th
International Conference in Asia (Qingdao, China, September 25-28, 2010).
2) S. Watanabe, T. Tada, S. Kasamatsu and T. K. Gu, “Ab Initio Based Simulations on
Electronic and Atomic Transport in Solid Electrolyte/Metal Junction Systems”, Materials
Research Society 2011 Spring Meeting (San Francisco, U. S. A., April 27, 2011).
3) T. Ono, “Spin-polarized current through graphene nanoflake”, The 6th Japan-Sweden
Workshop on Advanced Spectroscopy of Organic Materials for Electronic Applications(Kagaonsen, Japan, November 23-26, 2011).
4) C. Hu, O. Sugino, Y. Tateyama, and K. Watanabe, “Nonadiabatic couplings fromtime-dependent density functional theory: Progress and challenges”, ISSP-CMSI
international symposium on MAterial Simulation in Petaflops era (MASP2012)(Kashiwa, Japan, July 2, 2012).
5) S. Watanabe, B. Xiao, and S. Kasamatsu, “First-Principles Simulations on Electric andDielectric Properties of Metal-oxide Heterostructures for Novel Nanoscale Devices”,
IUMRS International Conference in Asia 2013 (IUMRS-ICA-2013) (Bangalore, India,
December 19, 2013).
6) Tomoya Ono, “First-principles electronic-structure and transport calculation formalism
using real-space grid based method” (Plenary talk), Computational and theoreticalanalysis of grid-based quantum many-body theory (Rhodes, Greece, September 21-27,
2013).
7) Tomoya Ono, “First-principles study on transport property of defects at surface and83interface”, The 16th Asian Workshop on First-principles Electronic Structure
Calculations (Beijing, China, October 28-30, 2013).
8) Tomoya Ono, “First-principles calculation on transport properties of carbon nanotubes
and fullerenes”, International Workshop of Computational Nano-Materials Design on
Green Energy, JSPS Core-to-Core Program Workshop (Osaka, Japan, June 1-3, 2014).
9) S. Watanabe and B. Xiao, “Density Functional Study on the Structures of Conductive
Filaments and Ion Migration Behaviors in Tantalum Oxide based Resistive Switching
Devices”, IEEE International Nanoelectronics Conference (INEC) 2014 (Sapporo, Japan,
July 30, 2014).
10) S. Watanabe and B. Xiao, “Density Functional Study concerning Tantalum Oxide based
Resistive Switching Devices”, 9th International Conference on Computational Physics(ICCP9) (Singapore, January 7-11, 2015).

84

New properties of materials probed by proton and muon, and their quantum
dynamics
H. Nakanishi 1, H. Goto 1, M. Geshi 1, M. Wilde 2,
W. Dino1*, H. Kasai1,2,3*, K. Fukutani2*, H. Kikuji1*
Osaka University 1, The University of Tokyo 2, National Institute of Technology, Akashi college3
We have been developing the first principles simulation methods and novel experimental
methods for the positive muon and proton states as hydrogen isotope nuclei states in thematerials, and investigating their behaviors and their induced new properties of materials. In
Chapter 1, we report our developing first principles simulation code “Naniwa”, for the positive
muon, proton and other small mass nuclei in the material environments, and its applications toreal material systems. In Chapter 2, we report a new algorithm for electron-electron correlationenergy calculations by superposition of nonorthogonal Slater determinants. In Chapter 3, wereport the physical properties of metallic hydrides under high pressure. In Chapter 4, we report
our experimental methods and their results for the catalytic reaction, diffusion, and desorption
dynamics of hydrogen in the near-surface region of palladium.

Ch.1 The first principles simulation code “Naniwa” and its applied results
1.1 Introduction
The first principles calculation based on the density functional theory is one of the mostsuccessful and powerful methods in the materials physics. And also the first principles
molecular dynamics methods have been used for various dynamical behaviors. In the case ofsmall mass atoms like hydrogen atoms, we should treat their nuclei (proton) in a quantummechanical manner, in addition to electrons. The quantum behaviors of the nuclei have crucial
roles to play in the materials related to today’s energy technology: fuel cell and lithium-ion
rechargeable battery, and also various catalytic reactions to produce hydrogen fuel, not only
from their academic interests. Positive muon is also one of our target quantum particles, because
it has been used as a magnetic probe particle in materials, -SR. The positive muon has thesame spin and same charge as a proton and behaves as the very light isotope of hydrogen which
has ninth mass of proton. It is important to analyze the quantum behavior of these particles in
the materials. In this chapter, we report our developing simulation codes “Naniwa” and its
application examples to real material systems.
1.1 Simulation code “Naniwa” and its algorithm.
We have been developing the quantum simulation code “Naniwa” for the small mass atom
nuclei, proton and positive muon on solid surfaces, in subsurfaces and bulk crystals. These
particle-case situations are different from that of electron case in the materials. We propose thefollowing suitable methods for the quantum state calculation for these particles in materials.
Here, we consider that the material are constructed from the nuclei and electrons. If weknow how these particles behave, we can know the almost material properties. The total
Hamiltonian of the electrons and nuclei in material is given byN 2  2 n  2  2(1)
   V (r,R )   H  

  2M 
n c .

I 1



I



  2m 
e

I

i 1



e



i

where MI and me are masses of the I-th nucleus and electron. Numbers of nucleus and electron
are N nc. and ne. The positions of the electrons and nuclei are indicated by(2)
r 1r ,2r ,3r ,, ri, rn  and R  R 1 , R 2 , R 3 ,, R I ,, R N ,enc.

* Co-Investigators (renkei-kenkyūsha).

85

where ri is the i-th electron position vector and R I is the I-th nucleus position vector. The
potential energy is given by the coulomb interaction between these particles;Nn N
n i 1I 1ZJ ZI e 2
 ZI e 2(3)e2V (r,R )   rR  r r   R R

n c .

e

i 1 I 1

i

n c .

e

i 1 j 1

IjI 1 J 1iJI

where Zl is atomic number of the I-th nucleus, and e is the elementary charge.
The total Schrödinger equation of electrons and nuclei is(4)
H  (r, R )  E  (r, R ) ,where E and  (r, R ) are the total eigen energy and the eigen wave function describing the stateof this material system. When we adopt the adiabatic approximation between the electrons andnuclei, the wave function, (r, R ) , can be written by,(5)r(R, ) n; R r( ) ; n (R )n ,and the total Schrödinger equation (4) can be divided into coupled two equations:(6)        V (r, R )  (r )  U (R )  (r ),     U (R )  (R )  E  (R )   2M 
 2m 
ne



N nc.

2

2
i

i 1



e

n ;R

n

n ;R



2

2
I

I 1

n



I

 ;n

 , n  ;nfor electron system and for nucleus system. Here,  n;R (r ) is the n-th eigen wave function ofelectron system with its eigen energy U n (R) in the case of the fixed nucleus positon R. By theuse of conventional ab initio electron states calculations or coming electron state calculation inCh.2 for various configurations of R, we can obtain the eigen energy U n (R) as a function of
R. And we solve the latter part of equation (6) and we can obtain the -th eigen function
 ;n (R ) of nucleus system in the case of electron system in its n-th state, and eigen energyE ;n , which is appeared as E in the eigen equation (4).
In the case of single small-mass-nucleus in the host materials, we described  ;n (R ) bylinear combination of the 3D-Gauss functions, which are located at mesh grid points   ,  ,   with its decay factors ( x ,  y ,  z ) :.
  x y z  exp   x x   2   y  y   2   z z    2
3



1

NG

 ( x, y, z )   C , 
 1



4



(7)

By use of the variation method, latter part of eq. (6) becomes matrix equation:
H  E S C  0, H  T  U .




 

 ,0

 

 

 ,

 

 

(8)

In the case of Gauss function based set, the overlap integral matrix elements, S  , and the
kinetic energy matrix elements, T , are analytically given by
y
 

2
   2   z    2 ,S  exp  x     
 




2









2

2





(9)




 2 
  S       2  1       2  1       2  1 .
T   
x
x


y
y


z
z


 2M  



 











(10)

When we use the discrete Fourier transform (DFT) of the potential energies, U n ( X j1 , Y j2 , Z j3 ) ,at sampling points ( X j , Y j , Z j ) , j1  1 N1 , j 2  1 N 2 , j3  1 N 3 :
1

2

3

~

 U n (k1 , k 2 , k 3 ) ,

(11)

DFT

U n ( X j1 , Y j2 , Z j3 )the potential energy matrix elements is also analytically given by
U    N1 1 / 2
S  

N 1 N 2 N 3 k1    N1 1 / 2

 N 2 1 / 2

k 2   N 2 1 / 2

 N 3 1 / 2~

k3   N 3 1 / 2

n

1

2

3



2

k 

k 
k 
1  2 k1 
1  2 k 2 1  2 k 3    i 1        i 2       i 3      . (12) L 
8 x  L y 
8 z  L z x x  Lz 
 Lx  y 

  U (k , k , k ) exp  8  L2

2

Hence the all matrix elements are given, and we can solve eq.(8) efficiently in Naniwa code [1].

86

1.2 Quantum states of hydrogen isotope adsorbed on metallic surface.
By use of our simulation code “Naniwa”, we investigate the hydrogen isotope states on themetallic surfaces. In the Table 1, we show the numerical results for the vibrational excitedenergies. Our simulation results show the good agreement with the experimental results. Theirdifferences are less than 10 %. The ratio of the excitation energy between deuterium andprotium is less than 2 . This reduced ratio is also shown in experiment. Even for suchfundamental physical quantity, the simple harmonic approximation around the minimum energy
point is failed, and such ab initio quantum simulation is necessary for the small mass atom
states.
NaniwaExperiment

Error

H on Pd(111)// 91 meV⊥ 114 meV//
96 meV [2]⊥ 124 meV [2]// 5.2%⊥ 8.1%

H on Cu(111)⊥ 135 meV⊥

129 meV [3]⊥ 4.7%

D on Cu(111)⊥ 104 meV⊥

96 meV [3]⊥ 8.3%

Isotope effect

1.29

1.34

3.4%

Table 1. Excitation energies of hydrogen isotope atoms adsorbed on metal surfaces for parallel (//) andperpendicular (⊥) vibrational modes.

In collaboration with the other group of (MEXT) Grant-in-Aid for Scientific Research on
Innovative Areas: “Frontier of Materials, Life and Particle Science Explored by Ultra Slow
Muon Microscope”, we investigate the positive muon and proton states on graphene [4].
Naniwa result for the hydrogen adsorbed energy, 0.58 eV, shows the good agreement with thatfrom TDS (Thermal Desorption Spectroscopy) experiments, 0.59 or 0.65 eV [5]. Only fromconventional electron states calculation, the calculated adsorbed energy is over estimated,because of lack of the zero-vibrational energy of proton. We found that an adsorbed state of
positive muon is also stable on graphene with lower adsorbed energy, 0.22eV. On graphene, apositive muon catches single electron on its 1s-orbital and becomes non-spin-polarized Mu(muonium), which will be observed by the coming ultraslow muon microscope.
In association with the experimental parts Ch.4, we investigate the hydrogen atom
absorption in hydrogen-covered Pd(110) (1 × 2) missing-row surface [6]. On this reconstructedsurface, the absorption of initially adsorbed hydrogen atoms is entirely assisted by the cominghydrogen molecule from outside, instead of the assistance of initially adsorbed H atoms. The
dissociation of hydrogen molecule is preferred to occur on H-vacant ridge site than on troughsite. Origin of such difference between these sites is clarified by tracing the position of thederived anti-bonding state in the transition states. We found that, for both cases, the dissociation
process is the rate-limiting step and hydrogen absorption is a non-activated process on thissurface.
Furthermore, we have been also developing the other version of Naniwa codes for scattering
problems, and investigate the dynamics of hydrogen isotopes (Mu, H, D, T) absorption into
Pd(111) subsurface [7]. We found the quantum correlation effects between hydrogen isotope
atom and surface lattice vibrations in adsorption process.
1.3 Quantum state of the hydrogen isotope impurity in the covalent bond crystals
After the achieve a measure of success in the case of metallic surface, we try to investigate
the quantum states of impurity hydrogen isotopes in covalent crystals, diamond structure crystalof silicon (Si) and -quartz (SiO2). In the case of Si, we found the impurity hydrogen state in
BC site as its ground state, and that in T site as its excited state, which correspond to those87

reported with the aid of the first principles path integral simulation by Miyake [8]. We noted that
correct results for these hydrogen impurity sites require the larger size of supercell in thesimulation, because of the impurity induced strong lattice stress. In this case, 3x3x3 supercell isneeded at the least. This simulation methods can be applied for revealing the hydrogen
behaviors in the metallic hydrides under high pressure appeared in Ch.3.
Nevertheless -SR experimental results show that a positive muon exist as a muonium in
SiO2 crystal, as a results of previous conventional ab initio study, the neutral muonium andneutral hydrogen atom impurity states are unstable in SiO2. We calculate the quantum states ofmuonium, hydrogen, deuterium and tritium atoms in -quartz by use of Naniwa code (Fig.1).
As a result of the 3x3x3 supercell calculation, which includes 81 silicon atoms and 162 oxygen
atoms, these neutral states are stable, whose ground state energies are -1.03eV, -1.15eV, -1.18eVand -1.20eV. The wave function is localized in SiO2 lattice cage. Calculated hyperfine constantis 4.47GHz, which are good agreement with experimental value [9]. In Fig.1, we found thesignificant quantum effects of Muonium as compared with other hydrogen isotope atoms, forexample, three degenerate ground states and larger vibrational excitation energy, which can beseen in energy gap from quantum number 3 to 4.

Fig.1. Eigen energies of a muonium (Mu), Hydrogen (H), Deuterium (D) and Tritium (T) atoms in -quartz
crystal. Upper insets show the wave functions of the 3 degenerate ground states of Mu. Lower right inset shows thewave function of the ground state of H. Cyan balls indicate silicon atoms. Red balls indicate oxygen atoms.

1.4 Conclusions
We have been developing the first principles simulation code “Naniwa” for positive muonand proton in various material environments. We have verified practical effectiveness of oursimulation methods for metallic surfaces, and for covalent bond crystals, against theexperimental given data. And we found some typical new quantum effects derived from small
mass of nuclei coupling with material environments. As for experimental methods by probe ofmuon being energetically developed in J-PARC (Japan Proton Accelerator Research Complex),
our simulation can provide the necessary corroborating information about muon site, muonstates (its motion, its spin, its valence electron and spin), their differences from proton(hydrogen), etc. in the various target materials. The quantum simulation must be powerful andnecessary tools to reveal the new properties of materials.
1.5 References[1] Japanese patent 4774523, US patent US8140467[2] H. Conrad, et al, J. Vac. Sci. Technol. A5, 452(1987).
[3] G. Lee, et al., Surf. Sci. 498 (2002) 229.
[4] A. D. Pant, et al., JPS Conf. Proc. 2, 010313(2014)[5] X. Zhao, et al., J. Chem. Phys. 124(2006)194704.

[6] A. A. B. Padama, et al., J. Chem. Phys. 140,
244707 (2014)[7] K. Shimizu et al.,, J. Phys. Soc. Jpn. 83, 013601(2014)[8] T. Miyake, et al., Phys. Rev. Lett. 81 (1998) 1873
[9] S. F. J. Cox, Rep. Prog. Phys. 72 (2009) 116501

88

Ch.2 A new algorithm for electron-electron correlation energy calculations bysuper- position of nonorthogonal Slater determinants
2.1 Introduction
In this study, a reliable and accessible algorithm for calculating electron-electron correlationenergies of few-particle systems with continuously controllable accuracy by the superposition ofnonorthogonal Slater determinants (SDs) [1] is presented. The key points of the proposed
procedure are updating of all the one-electron wavefunctions by employing correction functions
and optimizations on the basis of the variational principle without constraints on theorthogonality and normalizations. The accuracy and applicability of the present algorithm aredemonstrated through calculations of the ground- and excited-state energies for atoms andmolecules.
2.2 Optimization algorithm
An N-electron wavefunction Ψ(𝐫1 , 𝜎1 , 𝐫2 , 𝜎2 , ⋯ , 𝐫𝑁 , 𝜎𝑁 ) is expressed by a linear
combination of nonorthogonal SDs as follows:𝐿Ψ(𝐫1 , 𝜎1 , 𝐫2 , 𝜎2 , ⋯ , 𝐫𝑁 , 𝜎𝑁 ) = ∑ 𝐶𝐴 Φ 𝐴 (𝐫1 , 𝜎1 , 𝐫2 , 𝜎2 , ⋯ , 𝐫𝑁 , 𝜎𝑁 ) .

(1)

𝐴=1

Here, 𝐫𝑖 and 𝜎𝑖 represent the position and spin index of the ith electron, respectively. L is thenumber of SDs and Φ 𝐴 (𝐫1 , 𝜎1 , 𝐫2 , 𝜎2 , ⋯ , 𝐫𝑁 , 𝜎𝑁 ) is the Ath SD constructed by one-electron
wavefunctions(2)
𝜓𝑖𝐴 (𝐫, 𝜎) = 𝜙𝑖𝐴 (𝐫)𝛾𝑖 (𝜎)𝐴with 𝜙𝑖 (𝐫) and 𝛾𝑖 (𝜎) being nonorthogonal and unnormalized one-electron spatial
wavefunctions and spin orbital functions, respectively. The one-electron wavefunction 𝜙𝑖𝐴 (𝐫)is constructed as a linear combination of Gaussian basis functions. In the present procedure, theith one-electron wavefunction 𝜙𝑖𝐴 (𝐫) in the Ath SD is updated by
𝑁𝑐
𝐴(𝑛𝑒𝑤)(𝐫) = 𝐶𝐴 𝜙𝑖𝐴(𝑜𝑙𝑑) (𝐫) + ∑ 𝐶𝐿+𝜇 𝜉𝜇 (𝐫),
𝜙𝑖
𝜇=1

(3)

where 𝜉𝜇 (𝐫), 𝐶𝐿+𝜇 and Nc are a correction function, the expansion coefficient and the numberof correction functions, respectively. 𝜉𝜇 (𝐫) is also expressed as a linear combination ofGaussian basis functions. Thus, we have a linear combination of L+Nc SDs as the new
N-electron wavefunction. The coefficients 𝐶𝑖 (𝑖 = 1,2, … , 𝐿 + 𝑁𝑐 ) are obtained by solving thegeneralized eigenvalue equations obtained by adopting the variational principle with respect tothe total energy. Iterations of this updating procedure for all one-electron wavefunctions
belonging to all the SDs lead to the ground state with the superposition of nonorthogonal SDs.
This algorithm is the main result of this study. After the convergence of the ground state, an
excited state can be obtained by adding new SDs and updating the new one-electron
wavefunctions by the same procedure.
2.3 Applications to few-electron systems
The total energy for a carbon atom by the present algorithm is shown in Fig. 1. The ground-stateenergy is obtained using less than 100 nonorthogonal SDs with an error of 0.001%. The convergenceproperty is so smooth that the accuracy of the total energy is controllable by the number of employed SDs.
In contrast, the full configuration interaction (CI) method requires over 105 orthogonal SDs, andthis reduction in the numbers of SDs is a notable advantage of adopting nonorthogonal SDs [1].
Fig. 2 illustrates the total energy convergence of a Be atom for different values of 𝑁𝑐 . The
convergences are improved by multiplying the correction vectors, since the volume of the89search space for a one-electron wavefunction increases with 𝑁𝑐 . Figs. 3 and 4 show thenumbers of SDs to count 98% of the electron-electron correlation energy for few-electron
systems. Both sets of results illustrate that the number of nonorthogonal SDs to determine theground state by the present algorithm gradually increases with the system. Fig. 5 shows thecalculated potential energy curve of a HF molecule. Although for short bond lengths, the errors
obtained by coupled-cluster theory with singles and doubles (CCSD) and plus perturbative
triples (CCSD(T)) are a few milliHartrees, at longer bond lengths, the accuracy appears to
become depleted [2]. In contrast, the present calculation ensures practical exact ground states atall bond lengths, since no approximations are employed. In the case of 90 degree twisted
ethylene around the C-C bond axis, the same problem on counting the static correlation energy
arises. Fig. 6 indicates the noticeable advantage of the present algorithm. Non-negligible errors
in the total energy difference between the singlet and triplet states at 90 degree torsion angle
given by CCSD(T) and CAM-B3LYP are shown.

Fig. 1. Total energy convergence as a function of thenumber of SDs for a C atom with 3 correction
functions.

Fig. 2. Convergence performance of the total energy
for a Be atom with different numbers of correction
functions.

Fig. 3. Required number of SDs to determine theground state of atoms with 6-31G* basis set.

Fig. 4. Required number of SDs to determine theground state of small diatomic molecules.

Fig. 5. Potential energy curve of a HF molecule obtainedusing proposed algorithm with 6-31G** basis set.

90

Fig. 6. Total energy variation depending on thetorsional angle of ethylene.

2.4 Conclusion
A new calculation algorithm for searching the ground- and excited-state wavefunctions by thesuperposition of nonorthogonal SDs is proposed. The accuracy and applicability of the present
algorithm are demonstrated through calculations of the ground- and excited-state energies for
few-electron systems.
2.5 References[1] H. Fukutome, Prog. Theor. Phys. 80 (1988) 417.
[2] P. J. Knowles, et al.,, J. Chem. Phys. 133 (2010) 224106, K.Shimizu et. al.,J. Phys. Soc. Jpn. 83, 013601 (2014)

Ch. 3 Physical properties of metallic hydrides under high pressure
3.1 Introduction
Hydrogen at high pressure has attracted many physicists for long time, for examples,
metallization and superconductivity. In addition, Ashcroft’ suggested that metallic hydrides
might be metallic and realized the superconducting state at lower pressure than pure hydrogen[1,2]. Many researchers expected that the superconducting transition temperature (T c) ofmetallic hydrides might be higher than that of the previously found materials. Firstly, this wasinvestigated by means of first-principles calculations based on Density Functional Theory(DFT) using Allen-Dynes formula based on BCS theory. Many reports showed promising
results. However, experimental maximum Tc was 17 K of SiH4 [3]. The difference between
theoretical predictions from the DFT calculations and experimental results discourages theresearchers to study H-rich metals in view point of promising superconductors.
We focus on H-rich metals from different points of view. One is Metal-to-Nonmetal transition,and the other is effects of hydrogen for elastic properties of H-rich materials.
Basic calculations of all studies are performed using Quantum-Espresso code [4] to relax
atomic positions and cells and calculate phonon dispersions and elastic constants. The elastic
constants are calculated by means of my developed method [5]. To analyze accurate electronic
structures, we use Wien2k code [6].
3.2 Pressure induced Metal-to-Nonmetal transition of AlH3
AlH3 is a R3c structure at ambient pressure. With increasing pressure, it transforms to a
Pm3 n structure through an intermediate structure, a P 1 structure [7]. The electronic structureof the Pm3 n structure is semimetallic. Normally, band dispersion becomes wide and strongly
overlaps with increasing pressure. However, this semimetallic band structure changes tononmetallic ones and the band gap is open at 300GPa within the GGA result in spite of the factthat the band dispersion becomes wide. In order to obtain more accurate estimation of theMetal-to-Nonmetal transition pressure, we have performed higher accurate calculation using themethod of one-shot GW approximation [8-10]. From this result, we have estimated thetransition pressure of about 200GPa as shown in Fig.1. Using the leading-edge technique of ahigh pressure measurement this may be verifiable.
The reasons of pressure induced band gap opening are (1) the moderately large difference ofelectronegativity between aluminum and hydrogen and (2) the orthogonality between the 3s
states and 2s states of Al.. We have also examined several materials by replacing Al with B, Ga,
In, Tl, N, P, As, Sb, and Bi, and by replacing H with F assuming the same structure. However,Pressure induced band gap opening with the same mechanism did not found in those synthetic
materials.
The structural stability of higher pressure region than 164 GPa is confirmed by calculating91

phonon dispersion and enthalpy compared with the Pm3 n structure and many type of
optimized structures. This is the first finding of pressure induced band gap opening of hydrides
and the mechanism is different from Metal-to-Nonmetal transition of high pressure Li.

Fig.1: Band dispersion of AlH3 with Pm3n structure at 200 GPa. Black solid (Red dashed) lines indicate
GGA (GW approximation) results.

3.3 Chemical trend of elastic constants of metallic hydrides
It is difficult to identify positions and numbers of hydrogen in hydrides from experiments
except for neutron diffraction measurement. We therefore have calculated elastic constants for
hydrides with changing the number of hydrogen and have investigated the chemical trend andsome properties from the electronic structure analysis. The calculated elastic constants can becompared with experimentally observed ones by means of Brillouin scattering measurement.
We have calculated RH2 (R=Sc, Y, La, Ti, Zr, V, Nb, Ta) assuming the structure being the
CaF2 structure which is the typical structure of this chemical composition. In this structure theatom of R is at 4a site (0, 0, 0) and hydrogen does at 8c site. When hydrogens are added, theyoccupy the 4b site (1/2, 1/2, 1/2). When hydrogens are removed from RH2, they are done fromthe 8c site (1/4, 1/4, 1/4; 3/4, 3/4, 3/4). We have calculated the elastic constants for each numberof hydrogen from n=-4 to 4 with respect to RH2 (RH2+; -1 ≤  = n/4 ≤1).
(a)
(b)
(c)

(d)

(e)

(f)

Fig.2: The elastic constant C11-C12 for (a) YH2+(-1≤≤１;=n/4), (b) ZrH2+, and (c)NbH2+, and C44 for (d) YH2+,(e) ZrH2+, and (f) NbH2+.

The trend of the elastic constants depends on the group of metals. Figure 2 shows the elastic

92

constants, C11-C12 and C44, for YH2+, ZrH2+, and NbH2+ as representatives. Note that for n=-2
and n=2 their structures are always tetragonally distorted although the data is plotted in Fig. 2(a),(b), and(c). The origin of the tetragonal distortion is the crystal symmetry of RH2.5 and RH1.5. In
the positive region of C11-C12, the value of C44 is also positive. These show the existing regionof these metallic hydrides within the harmonic approximation, which is consistent with the Joint
Committee on Powder Diffraction Standards (JCPDS) data [11]. These results can be used toestimate the number of hydrogen and their positions for RH2+.
A general outline of the electronic structures of these materials can be understood that they
are quite similar with those having D4h symmetry. This means that these materials essentially
includes the electronic structure to distort tetragonal structure by changing the amount ofhydrogen. In fact, several RH2+ type materials are tetragonal structures [11].
3.3 Summary
We have found pressure induced Metal-to-Nonmetal transition for AlH3. This transition is thefirst finding for hydrides. The elastic constants have been calculated for several metallic
hydrides. If the elastic constants are measured from experiments, then the number of hydrogenand their positions can be estimated.
3.4 References
Universitaet Wien, Austria),2001.
[7] I. Goncharenko et al., Phys. Rev. Lett. 100, (2008)
045504.
[8] T. Kotani, M. van Schilfgaarde, Phys. Rev. B 81
(2010) 125117.
[9] M. van Schilfgaarde, T. Kotani, S.V. Faleev, Phys.
Rev. B 74 (2006) 245125.
[10] T. Kotani, M. van Schilfgaarde, S.V. Faleev,
Phys. Rev. B 76 (2007) 165106.
[11] http://www.icdd.com/[1] N.W. Ashcroft, Phys. Rev. Lett. 92 (2004)
187002.
[2] J. Feng et al., Phys. Rev. Lett. 96, (2006) 017006.
[3] M. I, Eremets, et al., Science 319 (2008)1506.
[4] P. Giannozzi, et al., J.Phys.:Condens.Matter 21,(2009) 395502.
[5] M. Geshi, Physica B 405, (2010) 517.
[6] P.Blaha, et al., WIEN2k, An Augmented plane
Wave þ Local Orbitals Program for Calculating
Crystal Properties Karlheinz Schwarz, Techn.

Ch.4. Investigation of the catalytic reaction, diffusion, and desorption
dynamics of hydrogen in the near-surface region of palladium
4.1 Introduction
In the final year of this project, we have continued our investigation into the hydrogen (H)transportation dynamics through palladium (Pd) surfaces, which are of crucial fundamental
importance for hydrogen storage and purification applications in hydrogen energy technology as
well as in industrial hydrogenation catalysis. Our previous research has provided detailed
insight into the near-surface hydride nucleation process and into the atomic-level hydrogen
absorption mechanism at the H2-exposed Pd(110) single crystal surface [1-4]. This established
our ability to systematically vary the spatial distribution of Pd-absorbed hydrogen underneaththe Pd(110) surface over a wide depth range. We also obtained evidence that the hydrogen
absorption kinetics are surface structure-sensitive. Based on this knowledge, we further pursued
our research following two specific objectives: 1.) To clarify the role of palladium-dissolvedhydrogen in the mechanism of Pd-catalyzed olefin hydrogenation, which has been subject ofcontroversial discussion in the literature for many years; and 2.) to explore the possibility tomanipulate the hydrogen transport between the interior and the surface of palladium throughdeliberate surface structural modification.
4.2 Clarification of the role of Pd-dissolved hydrogen in olefin hydrogenation catalysis on
palladium surfaces
Palladium is the most extensively applied industrial catalyst for the partial hydrogenation of93

unsaturated hydrocarbons (olefins) in petroleum refining and in organic synthesis. This reaction
adds two hydrogen atoms to the C=C double bond of the olefin adsorbed on the catalyst surface.
Surprisingly, this surface reaction is known to require hydrogen dissolved in the interior of thepalladium catalyst, but no consensus on the particular role of the Pd-absorbed hydrogen in the
catalytic mechanism has yet been reached in the literature. On one hand, it is possible that
Pd-absorbed H atoms become directly involved in the hydrogenation after they emerge from themetal interior onto the catalyst surface in an energetic state, because H in Pd bulk interstitialsites is ~0.3-0.4 eV less stable than in chemisorption sites on the surface [5]. This argument isanalog to the rather well-accepted explanation for the particular hydrogenation reactivity ofsubsurface hydrogen underneath alkyl-adsorbed Ni surfaces [6, 7]. On the other hand, one may
consider that – if present in sizeable populations (>50 %) - hydrogen in subsurface sites may
indirectly activate co-adsorbed reactants on the surface by altering the electronic structure of thecatalyst surface [8].
To elucidate this important mechanistic problem in hydrogenation catalysis on Pd, we here
applied our unique combination of hydrogen depth profiling via resonant 1H(15N,)12C nuclear
reaction analysis (NRA) and thermal desorption spectroscopy (TDS) [9] to reveal the hydrogen
distribution on and beneath a Pd(110) model catalyst surface under hydrogenation reaction
conditions for a typical conversion of cis-2-butene (C4H8) to butane (C4H10) and characterized
the reaction behavior of hydrogen/C4H8 co-adsorption systems on Pd(110) by means of TDS. A
full account of this study is published in Ref. [10].
Initially, we prepared different amounts and depth distributions of absorbed hydrogen (Habs)
in Pd(110) by controlling the H2 dosage and the Pd crystal temperature during the H2 exposure,then post-adsorbed a monolayer of C4H8 molecules at 115 K. The TDS traces of the C4H10hydrogenation product in Fig. 1 (a) show that the catalytic reaction proceeds efficiently between
160 and 250 K when absorbed hydrogen has been pre-charged into the Pd(110) surface, and thatthe product yield increases with the amount of pre-dosed H2. Notably, surface-adsorbed H on
Pd(110) (obtained after H2 exposure of 1 L = 1.33×10-4 Pa ∙ s) alone does not achieve any C4H10
production. This result is similar to observations at Pd nanoparticle catalysts [11] and forPd(111) [12, 13] and showcases the indispensable role of Pd-absorbed hydrogen for the catalytic
reactivity. Fig. 1 (b) summarizes the results of TDS reactivity measurements by showing thedependence of the butane yield on the amount of absorbed H (controlled by the H2 dosage) for
three different Pd crystal temperatures during the H2 exposure (Te = 115, 160, and 200 K).

Fig. 1: (A) Series of C4H10 TPD spectra recorded after sequential exposure to various amounts of H2 and to 1.25 L
C4H8 at 115 K. (B) Butane (C4H10) reaction yield as a function of the amount of Habs prepared at Te = 115 K (circles),
160 K (rectangles), and 200 K (diamonds). The solid lines denote the H 2 quantity desorbing in TPD measurementsfrom bare Pd(110) below 250 K. (B) is adapted from Ref. [10].

Evidently, the butane reaction yield depends linearly on the quantity of Pd-absorbedhydrogen at all Te’s and larger hydrogenation yields are obtained when Habs is prepared at lower
temperatures. The NRA H depth profiles of H2-exposed Pd(110) in Fig. 2 (a) reveal that this94

behavior relates to the initial depth distribution of the absorbed hydrogen, which tends to belocated closer at the surface at lower Te. The NRA H profiles in Fig. 2 (b) – (d), however,evidence that although the depth distributions of absorbed hydrogen prepared at the three Te’sare initially very different, the potentially large initial H concentration in the subsurface region
vanishes when the butene/hydrogen co-adsorption system on Pd(110) is stepwise annealed to
the reaction temperature of the catalytic hydrogenation (maximum at 218 K, Fig. 1 (a)).
Complementary TDS measurements detect only little H2 desorption during annealing to 218K,
hence the loss of the H in the near-surface region is due to diffusion into the deeper Pd bulk.
The remaining NRA profiles at 218 K (Fig. 2 (d)) reflect the signal of the saturated H
chemisorption layer on the Pd(110) surface at a coverage of 1.5 monolayers (ML).[4]Fig. 2: (a) Near-surface NRA H-depth profiles from Pd(110) exposed to H2 to prepare an equivalent of 10 ML ofabsorbed hydrogen at Te = 115 K (circles), 160 K (squares), and 200 K (diamonds). Triangles and a corresponding
Gaussian fit (solid line) indicate the resonance profile of surface hydrogen at saturation coverage (1.5 ML). (b)-(d):Evolution of the near-surface H distributions after co-adsorbing 1.25 L C4H8 at 115 K and briefly annealing to (b) 150
K, (c) 190 K, and (d) to the hydrogenation reaction temperature, TR = 218 K. The -yield of H in co-adsorbed C4H8
has been subtracted from the data in (b)-(d).

The combined data of Figs. 1 and 2 thus signify that the catalytic hydrogenation reaction canproceed efficiently although the H concentration in the Pd subsurface region is as small as 0.5
at. % (evaluated in depths larger than 5 nm, where surface hydrogen does no longer contributesignificant NRA signal). This clearly rules out the proposed indirect scenario that large H
concentrations in subsurface sites would influence the hydrogenation reactivity by activating
surface-adsorbed reactants through alterations of the surface electronic structure. On the other
hand, the solid lines in Figure 1 (b) denote the amount of absorbed H that desorbs from theclean Pd(110) surface in the temperature range of the hydrogenation reaction (up to 250 K). This
indicates that the butane reaction yield is linearly proportional to the number of absorbed H
atoms arriving at the Pd surface by diffusion from the Pd bulk at the reaction temperature.
Our results strongly support the interpretation that re-emergence of bulk-dissolved hydrogen
onto the H-saturated Pd surface triggers the catalysis by supplying reactive H species for thehydrogenation reaction, whereas the hydrogenation catalysis does clearly not depend on large
stationary H concentrations in Pd subsurface sites. We hence reach a clear decision in favor ofthe direct scenario for the role of Pd-dissolved H in Pd-catalyzed olefin hydrogenations, where
the Pd-absorbed hydrogen provides a reservoir of reactive hydrogen species that entertain thehydrogenation reaction after resurfacing from the Pd bulk.
4.3 Release control of Pd-dissolved hydrogen through surface structural modification
Our previous research into the atomic-level hydrogen absorption pathways at H2-exposed
Pd(110) with 15N NRA H depth profiling and thermal desorption spectroscopy withisotope-labeled surface hydrogen [4] has demonstrated that two H2 absorption pathways exist on
the Pd(110) surface: One is related to surface defects, whereas the other one involves regular95

terrace sites. In comparison to H2 absorption at Pd(100) [14, 15] and Pd(111) [16], where only
surface defects (such as step edges and terrace pits) appear to be active sites for hydrogen
ingress, the efficient absorption at regular terraces of Pd(110) seemed to suggest that hydrogen(H2) absorption is a structure-sensitive process, and that certain properties of the particular
H-induced pairing-row reconstruction of the Pd(110) surface enhance the hydrogen uptake
kinetics relative to the atomically more close-packed (111) and (100) surfaces. In order toexplore this hypothesis further, we have looked into the possibility to deliberately influence thehydrogen transport through palladium surfaces by structural modifications, exploiting the factthat the Pd(110) surface is susceptible to several CO (carbon monoxide) adsorption-induced(de-)reconstructions. While a detailed account of this investigation is published in Ref. [17], weonly provide a brief summary of the main results here.
The experimental study was performed by dosing 1000 L of H2 onto the clean Pd(110)
surface at either 120 or 160 K in order to H-saturate the surface chemisorption sites and to
populate the H absorption states [4] with a few MLs of hydrogen. Subsequently, various
amounts of CO (up to 3 L) were co-adsorbed at the same surface temperature and the H2desorption behavior was examined with TDS. The resulting structural phases of the Pd(110)
surface and the associated effects on the H2 thermal desorption spectra are summarized in
Figure 3. We verified that the H2 pre-dosage of Pd(110) has no influence on the CO-induced
structural modifications and that CO adsorbs on this surface with a constant sticking probabilityof 0.5 until it saturates after dosage of 2 L at a maximum coverage of 1.0 ML.[17]The bare low-temperature (160 K) H2-exposed Pd(110) surface initially exhibits a H-induced(2×1) pairing-row (PR) reconstruction. The low temperature region of the H2 TDS trace isdominated by two features (1 at 160 K and 3 at 200 K), which are due to desorption ofabsorbed H through surface defects and regular terrace sites, respectively [4]. Preferential
adsorption of small CO amounts (below 0.4 L or 0.2 ML) at defects eliminates the 1 desorption
channel at 160 K by blocking the exit sites for absorbed H, leaving now only one dominant
desorption feature of absorbed H from the (2×1) PR reconstructed terraces at 200 K (1).
A second step of CO-induced modification occurs in the dosage/coverage regime of 0.4 – 1.2
L (0.2 – 0.6 ML), which lifts the original H-induced PR reconstruction and restores a bulk
truncated (1×1) surface. This is accompanied by a second high temperature shift of the H2desorption feature from 200 K (1) to 270 K (2).

Fig. 3 (A) Ball models of the CO-induced surface (de-)reconstructions of Pd(110). White, gray and dark balls denotesurface, first and second subsurface Pd atoms, respectively. PR: (2×1) pairing-row; MR: (1×2) missing-row structures.
(B) Three steps of H2-TPD modification with increasing CO post-dosage. Pd(110) is sequentially exposed first to
1000 L H2 and various amounts of CO. Gas dosage is carried out (a) at 120 K, or, (b) and (c) at 160 K. Adapted from
Ref. [17].

A final stage of CO-induced structural modification takes place during the TDS heating ramp,
when the Pd(110) surface holding 0.3 – 1.0 L CO is heated above 340 K. This results in a (1×2)

96

missing-row (MR) reconstruction of the surface, which in turn leads to another substantial
high-temperature shift (by 105 K!) of the H2 desorption that now takes place in a dominant 3feature at 375 K. Remarkably, this final desorption temperature is higher than that of surfacechemisorbed hydrogen from the clean Pd(110) surface (~300 K), and because post-dosed CO
pushes chemisorbed H into the Pd interior [17], a distinction of surface and Pd-absorbed H as in
the H2 TDS trace from clean Pd(110) [4] is no longer possible after CO-modification.
A further noteworthy observation is that all H2 desorption features of the Pd(110) surface
described in Figure 3 occur while CO is still present on the surface (CO desorbs at even higher
temperatures than H2). With exception of the suppressed 160 K 1 feature, this clearly
distinguishes the present manipulations of the H2 desorption behavior as actually being due tosurface structural modification rather than as resulting from mere CO site blocking effects asreported for a few other surfaces [18-21]. We further revealed through NRA H depth profilingthat prior to desorption through the CO-modified surface above 200 K, all near-surface absorbedhydrogen, including H in first-layer subsurface sites, diffuses into the deeper Pd bulk [17].
In summary, we succeeded in demonstrating the possibility to deliberately raise thedesorption temperature of Pd-dissolved hydrogen in three distinct steps from 160 K of the clean
Pd(110) surface to 200 K, 270 K, and even up to 375 K through well-defined CO-induced
surface structural modifications. The individual H2 release modes (related to distinct surface
structural phases) are easily selectable by controlling the CO coverage. The largest achieved
temperature increase by remarkable 215 K shifts the H2 release from a cryogenic regime (160
K) to a range significantly above room temperature (375 K), which may open novel
opportunities for hydrogen storage technology of much reduced cooling energy cost.
4.4 Summary
Through direct visualization of the hydrogen distribution on and below a Pd single crystal
surface during catalytic conversion of butene (C4H8) to butane (C4H10) with NRA hydrogen
depth profiling, we have succeeded in finally clarifying the long-debated mechanistic role ofPd-dissolved hydrogen in Pd-catalyzed olefin hydrogenation. Pd-dissolved hydrogen is
identified as a resource for reactive hydrogen species after it emerges onto the H-saturated
surface by diffusion from the Pd bulk. Large stationary H concentrations in subsurface sites, on
the other hand, are not necessary to sustain the catalytic hydrogenation reactivity [10].
We further demonstrated the possibility to control the desorption dynamics of Pd-dissolvedhydrogen through well-defined modifications of the Pd(110) surface structure, achieving the so
far largest reported desorption temperature shift (215 K) of absorbed hydrogen from a cryogenic
regime (160 K) to above room temperature (375 K) [17]. The temperature at which Pd-absorbedhydrogen re-emerges at the Pd surface to either desorb as H2 or to induce a catalytichydrogenation reaction is a critical parameter that has the potential to improve process
conditions and catalytic selectivity.
In conclusion, we believe that the fundamental insight into the atomic-level dynamics ofhydrogen in catalytic reactions, diffusion, and desorption processes at palladium surfaces
obtained in this project will contribute significantly to the development of energy-efficientmaterials and process conditions in hydrogen storage applications and in industrialhydrogenation catalysis.
4.5

Referencesmaterials probed by proton and muon, and theirquantum dynamics”, Ch.4.

[1] The 2011 annual report on “New properties ofmaterials probed by proton and muon, and theirquantum dynamics”, Ch.4.

[3] The 2013 annual report on “New properties of[2] The 2012 annual report on “New properties ofmaterials probed by proton and muon, and their97

quantum dynamics”, Ch.4.

10853.

[4] S. Ohno, et al., J. Chem. Phys. 140 (2014)
134705.

[13] A. M. Doyle, et al., Angew. Chem. Int. Ed. 42(2003) 5240.

[5] Y. Fukai, The Metal Hydrogen System, Springer,
Berlin, New York, 2005.

[15] M. Wilde, et al., Surf. Sci. 482-485 (2001) 346.

[14] H. Okuyama, et al., Surf. Sci. 401 (1998) 344.

[6] A. D. Johnson, et al.,, Science 257 (1991) 223.

[16] G. E. Gdowski, et al., J. Vac. Sci. Technol. A 5(1987) 1103.

[7] S. T. Ceyer, Acc. Chem. Res. 34 (2001) 737.
[8] H. A. Aleksandrov, et al., Angew. Chem. Int. Ed.
53 (2014) 13371.

[17] S. Ohno, et al., J. Phys. Chem. C., submitted.

[9] M. Wilde, et al., Surf. Sci. Rep. 69 (2014) 196.

[19] G. Rupprechter, et al., Surf. Sci. 554 (2004) 43.

[10] S. Ohno, et al., J. Phys. Chem. Lett., submitted.

[20] M. D. Marcinkowski, et al., Nature Materials 12(2013) 523.

[18] C. Nyberg, et al., Surf. Sci. 256 (1991) 9.

[11] M. Wilde, et al., Angew. Chem. Int. Ed. 47(2008) 9289.

[21] S. Ogura, et al., J. Phys. Chem. C 117 (2013)
9366.

[12] M. Morkel, et al., J. Chem. Phys.119 (2003)

Ch.5 Refereed papers & Invited talks
5.1 Refereed papers
1. H. S. Kuncoro, M. Sakaue, H. Nakanishi, H. K. Dipojono, H. Kasai, J. Phys. Soc. Jpn. 80(2011) 024601.
2. M. K. Agusta, W. A. Diño, M. David, H. Nakanishi, H. Kasai, Surf. Sci. 605 (2011)1347.
3. Y. Kunisada, H. Nakanishi, W. A. Diño, H. Kasai J. Vac. Soc. Jpn. 55 (2012) 115.
4. A. Sarhan, O. Y. Ching, H. Nakanishi, K. Fukutani, H. Kasai, J. Appl. Phys. 113 (2013) 023703.
5. K. Shimizu, W. A. Dino, H. Kasai , J. Phys. Soc. Jpn. 83 (2014) 013601.
6. A. A. B. Padama, B. Chantaramolee, H. Nakanishi, H. Kasai, Inter. J. Hydrogen Energy 39(2014) 6598.
7. H. Goto and K. Hirose, J. Nanosci. Nanotechnol. 11 (2011) 2997.
8. A. Sasaki, M. Kojo, K. Hirose, H. Goto, J. Phys.: Condens. Matter 23 (2011) 434001.
9. A. Sasaki, K. Hirose and H. Goto, Curr. Appl. Phys., 12 (2012) S96-S99.
10. H. Goto, M. Kojo, A. Sasaki, K. Hirose, Nano. Res. Lett. (2013), 8:200
11. M. Geshi, T. Fukazawa, Physica B 411 (2013) 154.
12. S. Ohno, M. Wilde, K. Fukutani, J. Chem. Phys. 140 (2014) 134705.
13. M. Wilde, K. Fukutani, Surf. Sci. Rep. 69 (2014) 196-295.
Other referred papers: 53

5.2 Invited talks
1. H. Nakanishi, “Quantum states of the muon on solid surfaces and in subsurfaces”, JPS 2011 Autumn Meeting
Symposium: “Ultra slow muon spectroscopy: its innumerable possibilities”, Toyama University, Toyama, Japan,
September 23, 2011.
2. M. Geshi and T. Fukazawa, “Pressure Induced Band Gap Opening of AlH 3”, International Workshop ofComputational Nano-Materials Design on Green Energy (Joint Workshop of Interactive Materials Science
Cadet Program (IMSC), Osaka University and S-1 JSPS Core-to-Core Program (A) Advanced Research
Networks) Awaji Yumebutai International Conference Center, June 16-19, 2013
3. M. Wilde, “Subsurface hydride formation mechanisms at Pd(110)”, 9th International Conference on Diffusion
in Solids and Liquids (DSL-2013), Madrid, Spain, June 24-28, 2013.
4. H. Nakanishi, “Quantum states of impurity hydrogen isotopes in covalent crystals”, 10th International
Conference on Diffusion in Solids and Liquids (DSL-2014), Paris, France, June 23-27, 2014.
5. M. Wilde, “Atomic-level mechanism of hydrogen absorption at palladium surfaces”, Seminar at Chemical
Physics Group (Prof. Dr. H.-J. Freund), Fritz-Haber-Institut der Max-Planck-Gesellschaft, Berlin, Germany,
July 1, 2014.
6. H. Nakanishi, “First-principles calculations of probe quantum positive muon”, JPS 70th Annual Meeting（2015)
Symposium: Current status and perspective of first-principles calculation methods, Waseda University Tokyo,
Japan, March 21, 2015.
Other invited talks: 22

98

Reaction Analysis and Conformational Sampling forMultidimensional Large Scale Systems

1

T. Yamato, 2A. Mitsutake
1
Nagoya University, 2Keio University

AFM force-distance curve of membrane protein: A coarce-grained dynamics simulation
Regarding the structure and dynamics of single-biomolecules, atomic force microscopy (AFM)
measurements provide us detailed information we would not obtain otherwise. For example,
during the forced unfolding of a bacteriorhodopsin (bR) molecule with typicalseven-transmembrane architecture, we observe force-distance (F-D) curves with multiple peaks,
which have attracted a broad attention from many researchers both in experimental andtheoretical expertise. Unfortunately, however, the physicochemical mechanisms of theseprocesses have not been fully elucidated yet. To address the issue, we conducted a molecularsimulation study with a novel coarse-grained model.
From a theoretical viewpoint, the abovementioned peaks are thought to be direct reflectionsof energy barriers along the unfolding pathways of bR on its complex energy surface. Regarding
the energy landscape of membrane proteins, Popot and Engelman proposed a famous model
describing the protein folding process as consisting of two independent steps: (1) formation ofindividual helices and their insertion into the membrane and (2) the subsequent helix-helix
association within the membrane environment. Here, it is very likely that the interactions thatare responsible for the first step mainly characterize of the F-D curves with multiple peaks.
Based on the working hypothesis, we conducted computer simulations of bR with thecoarse-grained (C-G) model with rigid-sphere representation for peptide units. The surrounding
waters and lipid membrane were treated as implicit continuum environment. As a result, theexperimental peak positions of the F-D curves in the literatures were reproduced well by our
computer simulations, indicating that the essential features of the molecular interactions
involved in the forced unfolding processes of bR by the AFM apparatus were well captured in
our C-G model.
Figure 1.
We
havedevelopedanewcoarse-grained model of membraneproteins. Our simulation reproduced
experimental peak positions observedin the AFM force-distance curves ofbR. The overall system (left), F-D
curve (middle), and bacteriorhodopsin
molecule (right).

Energy exchange network of inter-residue interactions within a thermally fluctuating
protein molecule: A computational study
Molecular dynamics (MD) allows us to “see” how the protein behaves under physiologicalconditions. To study molecular mechanisms involved in protein functions, it is important to99

develop good tools to extract useful information from MD trajectories.
Here we propose our method to analyze the residue-residue interactions based on the energy
flow within a thermally fluctuating protein molecule. The “crosstalk” among interacting
residues was represented as a graph of energy exchange network (EEN). As an example, we
analyzed the properties of a small globular domain, PDZ3, which is well known to exhibit
single domain allostery. The removal of its C-terminal helix is known to decreases ligand
affinity by 21-fold without changing the overall protein structure; nevertheless, wedemonstrated that the helix constitutes an essential part of the EEN graph, indicating that theinter-residue interactions play crucial roles in the protein’s function.
Our computer program called “CURP” (CURrent calculatioin for Proteins), which was usedpreviously to analyze the atomic stress tensors in a small protein molecule, has been extended toinclude the energy flow analysis in this study.
N363

R312

T321

G329

I336

I338

V362

V315

G345

F337

I359

G356

G364

D366

P346

V365

D332

S371

R368

A375

E373

Q374

G324

A376

I377

A378

I327

K380

N381

A382

G333

T385

Ligand binding pocket
G322

H372

L379

S339

I341

E401

K355

Y397

E331

G335

E334

F400

E396

P394

E395

K393

R399

A402

S398
I314

E310

P311

S320

L349

R354

D357

R309

Q358

Q391

Y392

A390

S361

R313

D348

E352

L360

A347

S350

I316

T387

I389

R318

G319

Q384

G351

L353

Figure 2.
The structure of PZD3 domain (left) and its energy exchange network (EEN) (right).
Each node represents an amino acid residue, and interacting residue pairs with irEC greater than 0.015(0.008) {0.006}

are connected by thick red (thin blue) {thin green} edges. The

ligand-binding pocket was indicated by yellow box. Black rounded rectangles represent amino acidresidues located in the α3 helix.

CURP (CURrent calculation for Proteins)
The physical properties of condensed matter systems can be illustrated at the microscopic scaleby evaluating the continuum variables in atomistic simulations. Among condensed mattersystems, native protein molecules are particularly interesting in that they realize biological
functions via an energy relaxation process from the non-equilibrium state to the equilibrium
state in response to external perturbations such as light illumination and ligand binding. We cancharacterize native proteins in terms of continuum variables and provide biophysical grounds
for the molecular mechanisms of protein functions. The values of continuum variables vary100

from site to site in a protein molecule of interest, reflecting its inhomogeneity and anisotropy.
Therefore, we can expect to find functionally important sites by analyzing the 3D field ofcontinuum variables in a protein molecule. For instance, we analyzed electron-tunnelingcurrents in proteins and discovered a key residue for the electron transfer reaction in DNA
photolyases[1]. Other examples are the analysis of the energy-transfer pathway in photoactive
yellow protein[2,3] and the identification of the “epicenter” of the protein quake in photoactive
yellow protein via stress tensor analysis[4]. Note that the first quantity is related to the electron
flow, the second is related to the energy flow, and the third is related to the linear momentum
flow. The CURrent calculation for Proteins (CURP) program is written in Python and
FORTRAN, reads (1) the parameters of the force-field functions and molecular topology dataand (2) the atomic coordinates and velocities from the molecular dynamics trajectory, and thencalculates the flow of physical quantities. The AMBER format is supported for the current
version[5]. In the future, we plan to distribute an advanced version of the CURP program, andthe supported formats will be extended to include those of other molecular simulation packages.
References:[1] Y Miyazawa, H Nishioka, K Yura, T Yamato, Biophys J 94 (2008) 2194-203.
[2] T Ishikura, T Yamato, Chemical Physics Letters 432 (2006) 533-37.
[3] T Yamato, in: D.M. Leitner, J.E. Straub (Eds.), Proteins : energy, heat and signal flow,
Taylor and Francis, New York, 2009, p. 129-47.
[4] K Koike, K Kawaguchi, T Yamato, Phys Chem Chem Phys 10 (2008) 1400-5.
[5] T Ishikura, T Hatano, T Yamato, Chemical Physics Letters 539 (2012) 144-50.

Figure 3.
The Logo of the CURP program (left), calculation
flow (upper right), and an example of CURP
calculations (bottom right ) in combination with
conformational sampling by conventional MD
simulation.

Computational study on the roles of amino acid residues in the active site formation
mechanism of blue-light photoreceptors
The cryptochrome and photolyase families have different biological functions, despite sharing
high degrees of sequence identity and structural similarity between them. From the structural
point of view, there is no doubt that their enzymatic activities are particularly sensitive to thefine-tuning of their active site conformation.
It has been accepted that the active site methionine (hereafter denoted as M-site) and glutamic
acid (E-site) residues play important functional roles in the enzymatic activity of thecyclobutane-pyrimidine dimer (CPD) photolyase. To examine the functional roles of the M-site101and E-site residues, we performed in silico mutation at the M-site using MD simulation in asystematic manner, and focused on the hydrogen bonding between the E-site and the CPD
substrate. Fragment molecular orbital calculations with electron correlations demonstrated that
substitution of the M-site methionine with either alanine or glutamine always destabilized theinteraction energy between the E-site and the CPD by more than 12.0 kcal/mol, indicating thatthe methionine and glutamic acid residues cooperatively facilitate the enzymatic reaction in theactive site.
Figure 4.
The E-site glutamic acid and the
M-site methionine residue
cooperatively facilitate the hydrogenbonding formation between the CPD
and E-site, which is thought to beimportant for efficient electron transfer
reaction in the active site of CPD
photolyases.

Ligand migration in myoglobin: A combined study of computer simulation andtime-resolved x-ray crystallography
Ligand migration in myoglobin (Mb) is a fundamental biophysical process in which theinterplay between dynamics, thermal fluctuation, and ligand-protein interaction plays a crucialrole. Nevertheless, the details of the molecular mechanism involved have not been fully
understood. From structural point of view, it has been accepted that the hydrophobic cavities
called Xe pocket and Distal pocket (DP) play important roles in ligand migration in Mb.
However, no x-ray structure exhibits broad pathways for ligand migration movements
especially at the bottlenecks between different internal hydrophobic cavities, indicating that thethermal fluctuations around the average structure are essential for regulating the molecularfunction under the physiological conditions.
So far, a large number of studies have been devoted to illustrate ligand migration in Mb withexperimental techniques, computer simulation techniques and other theoretical approaches. By
using time-resolved X-ray crystallography at high resolution under the low-temperature condition,
Tomita and colleagues recently demonstrated the structural changes and the breathing motion
associated with the ligand migration movement of carbon monoxide in the wild-type sperm whale
myoglobin. Laser flash photolysis induces the structural change of the protein from thecarbon-monoxide-bound form (MbCO) to the carbon-monoxide- unbound form (Mb…CO). Under
cryogenic conditions, temperatures are low enough for high-resolution crystallography, and eachligand is supposed to reside inside each protein.
To illustrate the ligand migration mechanism, they proposed the self-opening model, in which theligand-protein interaction broadens the bottleneck between cavities and facilitates efficient ligand
movement. Here, we developed a new computational method: homogeneous ensemble displacement(HED), which effectively approximates the statistical ensemble of Mb…CO based on that of MbCO.
Using this method, we demonstrate the importance of the ligand-protein interaction for the ligand
migration. Although the structural changes were small, the presence of the ligand molecule in102the protein matrix significantly affected the 3D scalar field of the potential of mean force for theligand-protein interactions, in accord with the self-opening model proposed in the previous
x-ray study.

Figure 5.
Flash photolysis of myoglobin. Upon flash
photolysis, the CO molecule (red) is detachedfrom the heme group (blue). There areseveral internal cavities: four Xe-bindingsites (Xe1, Xe2, Xe3, Xe4), the distal pocket(DP) and the phantom site (Ph), indicated bydotted lines (green). In the Mb…CO form,
the ligand molecule migrates from one cavity
to another.

Figure 6.
PMF map for MbCO (left) and Mb…CO(right). (a,b) Xe1, (c,d) Xe2, (e,f) Xe3, (g,h) Ph,
and (i,j) DP and Xe4. Value of PMF, W, is
indicated by blue

, cyan, white, pinkand red. A 3Dinterpolation between grid points provided thescalar field of PMF for the molecular interior
space. Isoenergetic contour

in

Mb…CO (right) is indicated by broken lines in
MbCO (left). We indicate the protein interior inblack.

Atomic stress tensor analysis of proteins
We performed stress tensor analysis of polypeptides at the atomic level by molecular
dynamics simulations. Most of the approaches to stress tensor analysis of molecules in theliterature assume that the interatomic forces are pairwise; however commonly used force-field103

functions contain multi-body interactions. To address this issue, we derived new mathematicalformulae for atomic forces based on force-field functions, in which the total force acting oneach atom is expressed as a summation of pairwise interatomic forces. Then, a computer
program for stress tensor analysis was developed, and it successfully characterized themechanical features of molecules at the atomic level.
Figure 7.
We
derivedmathematicalexpressionsforpairwiseinteratomic forces based onforce-field functions containing
multi-body interactions. Then,
atomic stress tensors were
calculated based on Hardy’s
formula.

Normal Mode Analysis of Compressibility of Evolutionary Lysozymes
According to X-ray crystal structures, globular proteins have a precisely defined compact
structure, but there are many packing defects or cavities in the interior of a protein molecule
which cause internal motions or flexibility response to thermal or mechanical forces. Suchcavities are easily compressed by pressure and hence isothermal compressibility βT is an
important measure of protein dynamics or volume fluctuation. However, the direct measurementof βT is generally difficult for protein as seen in limited reports because the volume of protein
molecule is only slightly reduced by pressure and high pressure induces conformational change
or denaturation of protein. Therefore, most compressibility studies of protein have been carriedout with the adiabatic compressibility βs that can be more easily and accurately determined bysound velocity measurements. An important finding of these compressibility studies is that only
single amino acid substitution induces large changes in βs and that the mutant having larger βsvalue has higher enzyme activity. This means that the local structural changes by mutation
cause significant modification of overall protein dynamics to affect its function.
The isothermal compressibilities (βT) of chicken lysozyme and its six evolutionary mutants at
Thr40, Ile55, and Ser91 were calculated from their X-ray structures by means of a normal mode
analysis at 300K. The energy-minimization was carried out by a Newton’s method using the
program FEDER/2 and an energy function ECEPP. The βT value was calculated as a sum of thevolume fluctuation from individual normal modes within the potential energy of harmonic range.
The βT values of two extant lysozymes from chicken and quail (bobwhite) were 1.61 and 1.59
Mbar-1, respectively, but other ancestral mutants showed the larger βT values up to 2.08 Mbar-1.
These results suggest that the ancestral mutants were highly flexible and that such mutation
followed a non-neutral evolutionary pathway. No visible change was found in the tertiary
structure of all the mutants but most ancestral mutants induced a large increase in the total
cavity volume (Vcav). An observed positive linear correlation between βT and Vcav predicts thatcavity or atomic packing is an important determining factor of the volume fluctuation inmolecular evolution of protein.

104

Table 1 Compressibility and volume fluctuation of extant and ancestral lysozymes.
a) The values in parentheses represent the ratio (%) of <V> to total molecular volume.
b) Mean total cavity volume estimated with the probe radii of 1.05–1.20 Å.
c) Taken from BA Malcom et al. Nature 345 (1990) 86-88.

Molecular mechanism of long-range synergetic color tuning between multiple amino acidresidues in conger rhodopsin
The synergetic effects of multiple rhodopsin mutations on color tuning need to be completely
elucidated. Systematic genetic studies and spectroscopy have demonstrated an interesting
example of synergetic color tuning between two amino acid residues in conger rhodopsin’s
ancestral pigment (p501): —a double mutation at one nearby and one distant residue led to asignificant λmax blue shift of 13 nm, whereas neither of the single mutations at these two sites led
to meaningful shifts.
To analyze the molecular mechanisms of this synergetic color tuning, we performed
homology modeling, molecular simulations, and electronic state calculations. For the double
mutant (DM), N195A/A292S, in silico mutation analysis demonstrated conspicuous structuralchanges in the retinal chromophore, whereas that of the single mutant (SM), A292S, was almost
unchanged. Using statistical ensembles of QM/MM optimized structures, the intrinsic excitationenergy of retinal chromophore was evaluated forthe three visual pigments. As a result, the λmax
shift of DM from p501 was –8 nm, while that of
SM from p501 was +1 nm. Moreover, principal
component analysis successfully identified three
determinants for color tuning: bondlength
alternation, counterion displacement.
Figure 8. (Left) Protonated Schiff-base of retinal and amino acid residues at the sites 195 and 292
in rhodopsin. Although these two amino acid residues are separated far away, they cooperatively
work for the color tuning of rhodopsin. (Right) Average structures of the protonated Schiff-base ofretinal and the surrounding amino acid residues are shown. The thermal average structure wasobtained from MD simulations for p501 (red), SM(orange), and DM (green).

Relaxation mode analysis of proteinsIn this project, we have applied generalized-ensemble algorithms to the water confined toslit nanopores to investigate the liquid-solid phase transition of the system. We also have apply
relaxation mode analysis to the protein system to investigate the dynamics and kinetics. Here we
explain the results of relaxation mode analysis.
Biopolymers have flexible structures and several varieties of functions. The functions arederived from not only structures but also the dynamics of structural fluctuations. Therefore, the105

dynamic properties of structural fluctuations of biopolymers are important for understanding theinterrelationship between its movement and functions. Classical molecular dynamics (MD)simulations are a popular and powerful method for describing the structure, dynamics, andfunction of biomolecules in microscopic detail. Recent technological advances have allowed for
simulations on the timescales of the order of microsecond. Because longer MD simulations havebeen feasible in recent years, it is more important to develop an analysis method to extract goodorder parameters for identifying local-minimum-energy states automatically, and investigate
dynamics and kinetics of proteins.
Relaxation mode analysis (RMA) has been developed to investigate ``dynamic'' propertiesof spin systems and homopolymer systems. We applied RMA to biomolecular system in the
project. The equilibrium time correlation functions of the relaxation modes, Xp satisfy
X p (t)X q (0) = δ p,q e

where− λ ptX p (t)X q (0) denotes the equilibrium correlation of Xp at time t and Xq at time 0 and λp

is called the relaxation rate of Xp. The relaxation modes and rates are also given as left
eigenfunctions and eigenvalues of time-evolution operator of the master equation of the system,
respectively. For single linear polymers represented by the Rouse model, which has no excluded
volume interaction and no hydrodynamic interaction, the relaxation modes are identical to the
Rouse modes, which have been playing an important role in the theory of the polymer dynamics.
In the case of linear polymers, the relaxation phenomena have been studied systematically interms of the relaxation modes and rates.
In RMA, the relaxation modes and rates are estimated by solving a generalized eigenvalue
problemN∑ C (t + τ ) f
j=1

i, j

N

N

i=1

j=1 p,i

∑ ∑

0

p, j

=e

− λ ptN∑ C (t ) f
j=1

i, j

0

p, j

f Ci, j (t0 + τ ) fq, j = δ p,q ,Figure 8. (Left) Protonated Schiff-base of retinal and amino acid residues at the sites 195 and 292 inrhodopsin. Although these two amino acid residues are separated far away, they cooperatively workfor the color tuning of rhodopsin. (Right) Average structures of the protonated Schiff-base of retinal
and the surrounding amino acid residues are shown. The thermal average structure was obtained
from MD simulations for p501 (red), SM(orange), and DM (green).

, where Ci, j (t) denotes the equilibrium time-displaced correlation function of the position of
the i-th segment relative to the center of mass of the polymer. RMA is formulated as avariational problem equivalent to the eigenvalue problem of the time-evolution operator of thesystem, where a relaxation mode Xp is approximated by a trial function, which is constructed asa linear combination of relevant physical quantities that are time-evolved for t0/2. The parametert0 is introduced in order to reduce the relative weight of the faster modes contained in thephysical quantities. The slow relaxation modes and rates are obtained approximately. The
optimization of the normalized equilibrium autocorrelation function of the trial function
time-displaced by τ leads to the generalized eigenvalue problem. In practice, the timecorrelation matrices of structural fluctuations for two different times (t0 and t0+τ) arecalculated through simulations. Then, by solving a generalized eigenvalue problem for these106matrices, the relaxation rates and modes are obtained from the eigenvalues and eigenvectors,
respectively. Thus, we call this method the RMA method using a single evolution time t0/2.
We applied RMA to the Monte Carlo simulation for a heteropolymer, peptide, system witha small number of degrees of freedom and its effectiveness was demonstrated. In theapplications to "homo"polymer, polymer systems, the translational degrees of freedom areremoved from the conformations of the polymer sampled in a simulation. However, therotational degrees of freedom are remained because they are interested in the slow rotational
relaxation of the polymer. For "hetero"polymer, biomolecule systems, in contrast, both thetranslational and rotational degrees of freedom are removed from sampled conformations of abiomolecule because the structural fluctuations around its average conformation are interested.
We explained how to remove the contribution of the rotational degrees of freedom in RMA andapplied RMA to a five-residue-peptide, Met-enkephalin in gas phase because it is known tohave several metastable or stable conformations.
A long Monte Carlo simulation of Met-enkephalin in gas phase was performed. The results
were analyzed by the principal component analysis (PCA) and RMA. PCA extracts modes withlarge conformational fluctuation and is a static analysis method. We compared the results ofboth methods and showed the effectiveness of the relaxation mode analysis. From thesimulation, we observed three groups due to different hydrogen bond patterns of backbone
atoms, A and B, C and D, E and F. The clusters A and B in the same group have same hydrogenbonds between backbone atoms H2-5 and the global-energy-minimum state. The cluster C andD have H1-4 and are the second lowest-energy-minimum state and the cluster E and F haveH1-5 and are the third lowest-energy-minimum state. Difference between A and B, C and D, or
E and F is the arrangement of side chain of Tyr1 (see Fig.9). We obtained the different view ofthe free energy surfaces between PCA and RMA. In Fig.1 the free energy surfaces obtained by
RMA are shown. The directions of slow relaxation modes obtained by RMA correspond to slow
movements (transitions) between local-energy-minimum states. The free energy surface interms of the coefficients of the first, second, and third relaxation modes distinguished theclusters because they have long relaxation times due to their stabilities. The first relaxationmode, which is similar to the first principal component axis, corresponded to the transitions
between clusters A and C or clusters A and D. The second relaxation mode, which is notobserved in PCA, represented the movement of the side chain of Tyr1 or the transition between
clusters C and D. The third relaxation mode, which is similar to the second principal component
axis, corresponded to the transitions between clusters C and E or between clusters D and F. We
find that the transition between compact structures of H1-4 and H1-5 is faster than the
movement of the side
chain of Tyr1 or thetransition between
cluster C and D. From
the free energy
surfaces obtained by
RMA, we identify theclusters andinvestigate thetransitions between
these clusters while
Figure 9. The free energy surfaces along the first and second slowest
PCA gives only static
relaxation mode axes (a) and along the second and third slowest relaxation
information. RMA ismodes axes (b) obtained by RMA and six typical structures.
a powerful analysis

107method to investigate the dynamics and kinetics of protein systems.
Although RMA is a powerful method to extract slow modes, it requires a relatively highstatistical precision of time correlation matrices. Owing to this requirement, RMA needs a muchlong simulation than to occur many transitions between local-minimum-energy states, which
may not be a problem because very long simulations have been feasible in recent years.
Moreover, it cannot treat a large number of degrees of freedom directly. Therefore, thereduction of the degrees of freedom is necessary for the application of RMA to protein systems.
To apply RMA to a protein system by reducing its degrees of freedom, we proposed an analysis
method, which is referred to as the principal component relaxation mode analysis (PCRMA)
method. In this method, PCA is carried out first and then RMA is applied to a small number ofprincipal components with large fluctuations. This method can reduce the degrees of freedom
systematically. Even if we succeed in reducing the degrees of freedom by using PCA, it issometimes difficult to obtain reasonable relaxation modes and times by RMA using a single
evolution time t0/2. This occurs when we cannot make t0 so large that the fast modesin the selected principal components time-evolved for t0/2 become negligible, because the valueof t0 should be chosen so that the statistical precision of the correlation matrices C(t0) and C(t0 +τ) is maintained sufficiently high. Note that the statistical precision of C(t) becomes worse forlarger t. In such a case, the selected principal components that are time-evolved for t0/2 containfast relaxation modes that are not well described by a linear combination of the selectedprincipal components. To solve this problem, we also proposed the RMA method using multiple
evolution times, which can reduce the contribution of the fast modes efficiently. We used the
combination of the two proposed methods, which is the PCRMA method using multiple
evolution times.
We performed an all-atom molecular dynamics simulation of human lysozyme in aqueous
solution for 100 ns at 300 K. PCA was carried out on the mass-weighted coordinates for 1029
heavy atoms of the protein. We used 10 lowest PC modes for RMA. We reduced the degrees of
freedom from 3081 (3087-6) to 10 to apply
RMA. In Fig. 10, the three-dimensional
vectors in the first and second PC modes andin the first and second slowest PCRMA
modes are shown. Black lines represent the
average structure and arrows represent thedirections of the modes. For the first PC
mode, the residues with large fluctuations
are Pro103 and Gln104, which are located at
an exposed turn part. The first PC mode
corresponds to a local movement of the turnpart. On the other hand, for the first slowest
PCRMA mode, there is no local large
fluctuation near Pro103 and Gln104 andFigure 10. Three-dimensional vectors at Cα, C
the mode seems to be a collective mode.
β, and Cγ atoms for the first (a), second (c) PC
PCA sometimes extracts local movements
modes obtained by PCA and the first (b) andwith large and fast fluctuations located atsecond (d) slowest RM modes obtained byexposed parts (e.g., turn parts) as
PCRMA. The vectors are divided by the mass
lower-order PC. PCRMA succeeds inweight and multiplied by 8. The Cα atoms ofextracting a collective mode without local
Pro103 and Gln 104 are colored in yellow.
and large fluctuation as the slowest
relaxation mode.

108

PCRMA can reduce the degrees of freedom systematically. The combination method of
PCRMA and multiple evolution times can extract slow relaxation modes.
Finally, we applied PCA and RMA to the system, in which large conformational changes
are observed. PCA is difficult to extract order parameters to identify the local-minimum-energy
states for the system with large conformational changes because PCA extracts the modes withlarge structural fluctuation around an average structure. PC modes with large structuralfluctuations do not correspond to the transitions between the local-minimum-energy states. On
the other hands, RMA extracts slow relaxation modes. The local-minimum-energy states areusually stable and stay long time during a simulation. The order parameters with slow
relaxations may correspond to the direction between local-minimum-energy states. Slowrelaxation modes may correspond to good order parameters to identify local-minimum-energy
states.
We performed a simulation of chignolin at a transition temperature. Chignolin consists often amino acid sequences and is the test peptide to develop simulation methods.
There were many simulations for chignolin. From previous researches, it is found thatchignolin has the stable state (native structure) and a misfolded state. We consider that byincreasing temperature, many transitions between two states may occur in the order of hundred
ns simulation. The simulation at a transition temperature corresponds to a long time simulation.
We have performed a 750 ns molecular dynamics simulation of chignolin in aqueous
solution at 450K. Many transitions among the native, misfolded, and unfolded structures wereobserved. In Fig. 3 the free energy surfaces along the first and second lowest PC axes (a) andalong the first and second slowest RM axes (b) are shown. We obtained the different view of thefree energy surfaces of PCA and RMA. From the free energy surfaces of RMA, we identified
not only the stable and misfolded states but also intermediate state automatically. The native andmisfolded structures have hairpin-line structures. The stable and misfolded states have slightly
different hydrogen bond patterns of backbone atoms. The intermediate structures havecharacteristic structure; turn structure from Asp3 to Glu5. RMA neglects the conformationalchanges with fast relaxation of the terminal residues and extract the intermediate state while
PCA dose not identify the intermediate state. The stable and misfolded structures have the sameturn structure. From the results of RMA, it is found that chignolin folds to the native or
misfolded structures through the turn structure. From RMA, we clearly understand the folding
process and kinetics of the system. RMA is powerful and efficient to analyze the system withlarge conformational changes such as folding/unfolding and long simulations and investigate
dynamics and kinetics of the system.

10
5
0 2nd

Unfolded

15

Intermediate

10

Stable

5 2nd
0

-5
-10
-5

0

5 10 15
1st-5

Misfolded

-10

-5

0
1st

5

10

Figure 11. The free energy surfaces along the first and second lowest PC axes (a) and along thefirst and second slowest RM axes (b)109Research publications (Total # of papers and book chapters: 23)
T. Kaneko, J. Bai, K. Yasuoka, A. Mitsutake, and X.C. Zeng, “Liquid-solid and solid-solid
phase transition of monolayer water: High-density rhombic monolayer ice", Journal ofChemical Physics, 140, 184507 (2014)
T. Yamato and K. Yura, “Combinatorial approach of molecular dynamics simulations anddatabase analysis for the studies of protein structure and function”, SEIKAGAKU (The Japanese
Biochemical Society Journal), 85 (8) pp. 646-655, 2013.
T. Kaneko, J. Bai, K. Yasuoka, A. Mitsutake, and X.C. Zeng, “New Computational Approach to
Determine Liquid-Solid Phase Equilibria of Water Confined to Slit Nanopores", J. Chem. Theo.
Compt. 9 (2013) 3299
T. Tsuduki, A. Tomita, S. Koshihara, S. Adachi, T. Yamato, “Ligand migration in myoglobin:
A combined study of computer simulation and X-ray crystallography”, J. Chem. Phys., 136(16)165101 (9 pages) (2012)
Ishikura, T. Hatano, T. Yamato, “Atomic stress tensor analysis of proteins”, Chem. Phys. Lett.
139-140, pp. 144-150, (2012).
A. Mitsutake, H. Iijima, H. Takano,” Relaxation mode analysis of a peptide system: comparison
with principal component analysis” J. Chem. Phys.135, 164102(15 pages), 2011.
Invited talks (Total #: 22)
Ayori Mitsutake, “Effective Sampling Algorithms and Analysis Methods for Protein
Simulations’’, The 14th KIAS Conference on Protein Structure and Function, September 19,
2014, Korea
T. Yamato, “Role of reorganizable interactions between native residue contacts on the allostericcommunication in proteins: A computational study by the inter-residue energy conductivity
analysis”, Protein dynamics, Telluride research science workshop, Telluride, CO, Aug. 5-9,
2013.
A. Mitsutake, “Development of effective sampling algorithm and analysis method forbiomolecular simulations”, 4th Japan-Korea Seminar on Biomolecular Sciences Experiments andSimulations, Jan.2012, Nara
T. Yamato, “Ligand migration in myoglobin: A combined study of computer simulation andx-ray crystallography”, 4th France-Japan Joint Seminar, Imaging of spatiotempral hierarchies inliving cells -- an overview of dynamics from molecules to cells – , Hyogo, Japan, Jan. 6 --11,
2012.
T. Yamato, “Intramolecular Communication Chart of Proteins: Computational Analysis of
Energy Transfer Pathways”, Thermal Transport at the Nanoscale, Telluride Science ResearchCenter Workshop, Telluride, CO, USA, Jun. 21-25, 2010.
Awards (Total #: 2)
A. Mitsutake, Young scientist award of the Physical Society of Japan (2013).

110

Developments of first-principles multiscale simulation
for coupled dynamics of light-wave and electrons in dielectrics
Kazuhiro Yabana
Center for Computational Sciences, University of Tsukuba

1. Introduction
In this project, we have developed a first-principles computational approach to describeinteractions between light and matter based on time-dependent density functional theory.
Interactions between light and matter are usually described by two theories: Macroscopic
Maxwell equations describe propagation of a light wave in a medium, while quantum
calculation of susceptibilities such as dielectric function is one of central issues infirst-principles condensed matter calculations. Recent progresses of laser technologies,
however, require theories beyond them. Ultrashort laser pulses of even a few cycles are often
used, and pump-probe measurements are carried out in time domain. They require theories intime-domain, solving the time-dependent Schroedinger equation. Strong electromagnetic fields
of a laser pulse induce extremely nonlinear electron dynamics which no more allows
perturbative separation between macroscopic electromagnetic fields and microscopic electron
dynamics. As our solution to the problems, we have recently developed a multiscale theory forlight-matter interactions, describing electron dynamics in a solid using real-time time-dependent
density functional theory (TDDFT) and coupling it with macroscopic Maxwell equations in a
multiscale strategy. The theory provides a quite general and computationally feasible basis for anumber of problems, including ordinary macroscopic electromagnetism in a weak field limit
and being applicable to problems involving intense and ultrashort laser pulses.

2. Multiscale theory for coupled dynamics of electrons and light-wave
In our simulation, we combine two computations of different spatial scales. One is the
microscopic electron dynamics whose spatial scale is less than a nano-meter. The other is themacroscopic dynamics of light-wave whose spatial scale is about a micrometer. We start withthe explanation of the former.
2.1 Microscopic electron dynamics
In a microscopic spatial scale, we consider electron dynamics in a unit cell of a solid under aspatially-uniform, time-varying electric field. A core solver in our approach calculates
time-evolution of Bloch orbitals under the electric field:
2δE
1 ⎡ ! ! e ! ⎤ !
+
+
(
)
A
t
u nk − eφunk! + xc unk!
p
k
⎢
⎥
δn
∂t
2m ⎣
c
⎦
,
(1)
!
!
!
where the electric field is described using the vector potential A(t ) as E(t) = −(1 / c)dA (t ) / dt .

i"∂u nk!

=.
Before the irradiation of the electric field, the system is given by a ground state of static
Kohn-Sham equation. From the time-evolved Kohn-Sham orbitals, we may calculate electric111

Figure 1. Left panel shows the ground-state electron density of Si in [110] plane. The middle and right panels show
the change in the electron density from that in the ground state by the laser pulse. In the middle and right panels, thered color indicates the increase in the electron density while blue color indicates the decrease. Taken from [8].

current density flowing in the unit cell,

!
! !
*
!! !
−e
1 ' !* ! ! ! e !$ ! !
	
 
A
ψ
(
r,
t)
p
+
ψ
(
r,
t)
+
c.c.
j
(
r,
t)
=
(
,
)
,
I (t ) =
d
r
j
r
t
#
&
∑
)
,.
i
i
"
%
c
2m
Ω ∫Ω
(
+
i
(2)
!
Integrating the current over time, we may! obtain the polarization, P (t ) . !Thus our microscopic
solver provides us with the polarization P (t ) for a given electric field E(t) . We may regard
this microscopic calculation as a numerical constitutive relation. In fact, we may extract
information of dielectric function using this solver [10].
In Fig. 1, we show electron density in a typical time evolution calculation, an ultrashort laser
pulse irradiation on crystalline silicon. The laser pulse of 2.5 eV frequency is applied in [001]direction. The left panel shows the electron density in [110] plane in the ground state. The
middle panel shows a change of electron density from the ground state under irradiation of thelaser pulse. The red region indicates increase of electron density from the ground state, and theblue region indicates decrease. Electrons in the bond region are showing dipole oscillation. The
right panel shows electron density change after the laser pulse ends. The figure indicateselectrons in the bond region decrease by excitations.
2.2 Maxwell + TDDFT multiscale simulation
A propagation of light wave
is usually described by macroscopic! in! a! medium!
! ! Maxwell
we assume that B = H and there
equations for the four fields, D, E!, B, H!
! . In the! following,Eholds a local relation
betweenandasusual.
Namely, the macroscopic
D = E + 4πP! ! !
polarization
at point R , P R, t is determined by the electric field of the same point in the past,! !
E R, t ' , t ' < t . When the electric field is weak, the polarization and the electric field are relatedby the linear suceptibility χ αβ (t − t ') . When the polarization field is strong but still expressedin power series of electric field, we may use, for example, the third order susceptibility tensor,( 3)(t − t ' , t − t ' ' , t − t ' ' '). If the field is further strong, polarization may not be expressed as a
χ αβγδ
power series. We then need an approach in which we do not assume power series expansions.
For this problem, we have proposed a multiscale scheme, solving the coupled equation of

( )

( )

Figure 2. Coordinate system adopted
in multiscale simulation. At eachmacroscopic grid points shown left,
microscopic grid points shown rightare used.

112

Maxwell and time-dependent Kohn-Sham equations. The macroscopic equation to be solved isan ordinary wave equation and is given by

4π e
1 ∂ 2 AZ (t ) ∂ 2 AZ (t )−=−J Z (t )
2
2
2
c
c
∂t
∂Z
,

(3)

where we assume the propagation of the light wave is in Z-direction. We solve this equation in
macroscopic grid points, which is shown in the left side of Fig. 2. At each macroscopic grid
point, we prepare a microscopic grid points, shown in the right side of Fig. 2. The microscopic
grid points are used to solve the time-dependent Kohn-Sham equation (1).
As an example, we show a calculation of propagation of a laser pulse irradiating on a Si
surface [10]. In Fig. 3, we show snapshots of the time evolution. In solving Eqs. (2) and (1), themacroscopic coordinate Z is discretized into 256 grid points with the grid spacing of 250 au forthe medium region, Z > 0 . At each grid point, microscopic time-dependent Kohn-Shamequation (1) is solved.
The initial laser pulse at t = 0 has a peak intensity of 1×1011 W/cm 2 , the pulse length is
18 fs, and the average frequency is 1.55 eV. This is shown in the upper panel of the figure. The
middle panel shows the field when the center of the pulse has just reached the surface. In thebottle panel, one can see a transmitted and the reflected waves. In the right panels, electronic
excitation energy per atom is shown. The excitation is seen in the spatial region wheretransmitted wave exists. Although the laser frequency is far below the direct band-gap energy,
the excitation is seen after the transmitted pulse went away in the right-bottom panel. This iscaused by multiphoton excitations.

3. Physics applications
We have been applying the method to various problems. Among applications, laser damage
processes which are intimately related to laser machining are one of important areas that thepresent method will be useful. If a laser pulse of long pulse duration, picoseconds and more, isFigure 3. Snapshots of the vector potential divided by light speed, A / c < 0 (left panel), and of electronic
excitation energy per atom (right panels) at different times are shown as a function of macroscopic position. The
vacuum is at Z < 0 and the Si crystal is at Z > 0 . Top panels show initial starting field with a pulse on left
moving toward the Si surface. Middle panels show the field at a time when the middle of the pulse reaches thesurface. Lower panels show the field when reflected and transmitted pulses are well separated. The maximum
intensity of the incident laser pulse is set at 1×1011 W/cm 2 . Taken from [10].

113

used, avalanche mechanism is considered to dominate for electron excitations. In recent
researches, it has been realized that nonthermal, much cleaner processing is possible using
femtosecond laser pulses. In this case, the initial excitation process is dominated by one-electron
process such as multiphoton and tunneling excitations. Since the present multiscale simulation
may describe energy transfer from laser pulses to electrons in solids, we expect the present
simulation will be useful to understand basic mechanism of nonthermal laser processing.
We have made calculations of energy transfer from a laser pulse to a medium solving thecoupled Maxwell + TDDFT simulation. As the material to be investigated, we chose α-quartz
which is of practical importance. A few-cycle, intense laser pulse of 20 fs duration is irradiatednormally to the α-quartz. In Fig. 4, we show energy deposited to the α-quartz as a function of
the distance from the surface. The energy deposition to the medium is shown in a unit ofeV/atom. Calculations are achieved for various laser intensities around 1014W/cm2.
Experimentally, it has been known that ablation process starts to occur at around the intensity of
2x1014 W/cm2.
As shown in Fig. 4, the energy deposited to electrons in the medium increases as theintensity of the laser pulse increases. However, the amount of the energy deposition saturates ata certain value. Electrons inside the surface cannot be excited much even if we increase the laser
intensity. This is because the laser pulse is strongly reflected at the surface when the number of
excited electrons at the surface exceeds a critical density. As a measure to know what kind ofprocess will take place after the energy deposition to electrons, two horizontal lines, ablationand melting, are drawn. The line 'ablation' is drawn at the energy corresponding to the averagecohesive energy of α-quartz. It can be seen that the deposited energy exceeds the ablationthreshold at the surface when the intensity of the laser pulse exceeds 2.5x1014W/cm2. We
consider that this roughly coincides with measured threshold. At stronger intensities, the energy
deposition exceeds the ablation line to the depth of about 100-150 nm from the surface. Thisvalue is again consistent with the observed ablation depth. We thus consider that the present
approach succeeds to describe basic aspects of the initial stage of laser damage processes.

Figure 4. Energy deposited from a laser pulse of various intensities to α-quartz is shown as a function of distancefrom the surface. Two horizontal lines, ablation and melting, are shown to indicate damage thresholds.

114

4. Summary and future prospects
We consider the present approach solving Maxwell and time-dependent Kohn-Shamequations simultaneously is a quite general and comprehensive framework for light-matterinteractions. For weak electromagnetic fields, it may be regarded as a solver of Maxwell
equation with the accurate dielectric function given by microscopic TDDFT. As the intensity of
the light pulse increases, various nonlinear effects are automatically included such as second
and third order nonlinearities. It can even describe interactions of very strong laser pulses withsolids when electrons are strongly excited to form a plasma. Thus the present method may be
used for analyses of laser machining. The simulation requires, however, quite large computerresources. For the present scheme to be really useful, use of newest supercomputers underintimate collaborations with computer scientists is quite important.

Publications (Refereed journals and books)
1) M. Schultze, K. Ramasesha, C.D. Pemmaraju, S.A. Sato, D. Whitmore, A. Gandman, J.S.
Prell, L.J. Borja, D. Prendergast, K. Yabana, D.M. Neumark, S.R. Leone, Attosecond
band-gap dynamics in silicon, Science 346, 1348-1352 (2014).
2) S.A. Sato, Y. Shinohara, T. Otobe, K. Yabana, Dielectric response of laser-excited silicon at
finite electron temperature, Phys. Rev. B90, 174303 (8 pages) (2014).
3) G. Wachter, C. Lemell, J. Burgdoerfer, S.A. Sato, X.-M. Tong, K. Yabana, Ab Initio
Simulation of Electrical Currents Induced by Ultrafast Laser Excitation of DielectricMaterials, Phys. Rev. Lett. 113, 087401 (5 pages) (2014)4) M. Noda, K. Ishimura, K.Nobusada, K. Yabana, T. Boku, Massively-parallel electron
dynamics calculations in real-time and real-space: Toward applications to nanostructures ofmore than ten-nanometers in size, J. Comput. Phys. 265, 145-155 (2014).
5) S.A. Sato, K. Yabana, Efficient basis expansion for describing linear and nonlinear electron
dynamics in crystalline solids, Phys. Rev. B 89, 224305 (11 pages) (2014).
6) K.-M. Lee, C.M. Kim, S.A. Sato, T. Otobe, Y. Shinohara, K. Yabana, T.M. Jeong,First-principles simulation of the optical response of bulk and thin-film α-quartz irradiatedwith an ultrashort intense laser pulse, J. Appl. Phys. 115, 053519 (8 pages) (2014).
7) S.A. Sato, K. Yabana, Y. Shinohara, T. Otobe, G.F. Bertsch, Numerical pump-probe
experiments of laser-excited silicon in nonequilibrium phase, Phys. Rev. B 89, 064304 (8pages) (2014).
8) K. Yabana, Y. Shinohara, T. Otobe, J.-I. Iwata, G.F. Bertsch, First-Principles Calculations
for Laser Induced Electron Dynamics in Solids - Time-Dependent Density-Functional
Theory for laser matter interactions, Advances in Multi-Photon Processes and Spectroscopy,Vol. 21, pp. 209-244, Eds. S.H. Lin, A.A. Villaeys, Y. Fujimura, World Scientific (2014).
9) Y. Shinohara, S.A. Sato, K. Yabana, J.-I. Iwata, T. Otobe, G.F. Bertsch, Nonadiabatic
generation of coherent phonons, J. Chem. Phys. 137, 22A527 (2012).
10) K. Yabana, T. Sugiyama, Y. Shinohara, T. Otobe, G.F. Bertsch, Time-dependent density
functional theory for strong electromagnetic fields in crystalline solids, Phys. Rev. B85,
045134 (2012).
11) K.-M. Lee, K. Yabana, G.F. Bertsch Magnetic circular dichroism in real-timetime-dependent density functional theory, J. Chem. Phys. 134, 144106 (9 pages) (2011).

115

Conference talks (invited)
1) First-principles description of strong electromagnetic fields in solids, K. Yabana, invited
talk at focus session: Computer Simulation of Interaction of Electromagnetic fields andNanostructures, APS Meeting 2014, Colorado Convention Center, Denver, USA, Mar. 3-7,
2014.
2) Maxwell + TDDFT multiscale simulation for strong electromagnetic fields in solids, K.
Yabana, Quantum Dynamics Research Meeting, Vanderbilt University, Nashville, USA,
Mar. 10-12, 2014.
3) Multiscale description for strong electromagnetic fields in solids, K. Yabana, Gordin
Research Conference on TDDFT, Univ. New England, Biddeford, ME, Aug. 11-16, 2013.
4) Light-matter interactions in time-dependent density functional theory, K. Yabana, Advancesin time-dependent methods for quantum many-body systems, ECT*, Trento, Italy, Oct.
14-18, 2013.
5) Ab-initio theory for intense laser pulses in solids, K. Yabana, Korean Physical Society 2012
Spring Meeting, April 25-27, 2012, Daejeon, Korea.
6) Real-time electron dynamics in solids under strong electromagnetic fields, K. Yabana,
ISSP-CMSI International Workshop/Symposium on Material Simulation	  in Petaflops era(MASP2012), June 27,2012, Tokyo, Japan.
7) Time-Dependent Density Functional Theory for Electron Dynamics, K. Yabana,
BERTSCHFEST, UW, Seattle, Sept. 7-9, 2012.
8) Ab-initio description for the interaction of intense laser pulses with solids, K. Yabana,
Frontiers in Intense Laser-Matter Interaction Theory (FILMITh), Sept. 19-21, 2012, Max
Planck Institute fur Quantenoptik, Garching, Germany.
9) Multi-scale first-principles approach for strong electromagnetic fields in solids, K. Yabana,Int. Symp. on Computics: Quantum Simulation and Design,	  Oct. 11-13, 2012, Osaka,
Japan.
10) Time-dependent density functional theory for ultrafast electron dynamics in solid, K.
Yabana, Computational approaches to nuclear many-body problems and related quantum
systems, Feb. 12-16, 2013, RIKEN, Wako, Japan.
11) Real-time TDDFT Calculation in Molecules and Solids, K. Yabana, ISTCP-VII(International Symposium on Theoretical Chemistry and Physics), Waseda Univ. Sept. 2-8,
2011.
12) Real-Time TDDFT for Molecules and Solids, K. Yabana, Dynamics and Correlations in
Exotic Nuclei, YITP Kyoto Univ., Sept. 23, 2011.
13) Time-dependent density functional theory for femtosecond electron dynamics in dielectrics,
K. Yabana, The 14th Asian Workshop on First-Principles Electronic Structure Calculations,
Univ. Tokyo, Oct.31-Nov. 2.
14) Time-Dependent Density Functional Theory for Intense Laser Pulse Propagation in Solids,
K. Yabana, 3rd Int. CQSE Workshop on Atomic, Molecular and Ultrafast Science andTechnology, National Taiwan Univ., Jan. 7-8, 2012.
15) Computational Approach for Dynamics of Many-Fermion Systems - from Nuclear Physics
to Optical Science-, K. Yabana, Kyoto Univ. GCOE Symposium “Links among
Hierarchies”, Feb. 13-15, 2012.

116

Hybrid-density-functional theory for magnetic materials
Yoshihiro Gohda
Department of Materials Science and Engineering,Tokyo Institute of Technology, Yokohama 226-8502
1. Introduction
First-principles calculations based on density functional theory play an indispensable rolein materials science and condensed matter physics. While standard generalized gradientapproximation is highly successful, a better approximation is desired for some of materials
including magnetic ones with localized states. One of such approaches is hybrid density
functional theory. However, a problem of hybrid density functional theory is theadjustability of the mixing parameter. In this work, we tried to examine a scheme toremove such adjustability. As a result, we identified that the mixing parameter can be
determined using the generalized Koopmans condition, which should be satisfied by exact
single-particle states. We also found that this parameter can be highly dependent individual
electronic state, so that the reduction of the computational cost remains a problem to be
overcome. In addition, we performed first-principles calculations for various magnetic
materials rather in a standard manner: strain effect on the magnetic anisotropy in Y-Fe-B
compounds and development of a method for magnetic-anisotropy decomposition with low
computational costs; the magnetic order and magnetic anisotropy of Fe ultrathin films on
the Cu(001) surface; and gigantic magnetic anisotropy in the ε phase of Fe2O3.
2. Results
We first examined the generalized Koopmans theorem to apply the hybrid density
functional scheme. We derived first-principles mixing parameters for theexchange-correlation functional in the generalized gradient approximation and the exact
Hartree-Fock exchange. By applying our scheme for type-II antiferromagnetic NiO with adistorted rocksalt structure, we successfully reproduced its electronic structure as acharge-transfer insulator. In contrast with standard generalized gradient approximation and
HSE hybrid functional, our functional with optimized mixing parameter provide theposition of unoccupied 3d states that is consistent with experiments [Fig. 1]. On the other
hand, the type-II antiferromagnetic FeO exhibits strongly state-dependent mixing
parameters, even though the parameter determined by the states at the conduction band
minimum reproduces FeO as a Mott-Hubbard insulator. We also found that the mixing
parameters for the above compounds are larger than transition metal impurities in nitride
semiconductors due to larger locality of electronic states. Even though their charged states
are metallic, numerical calculation can be performed because wrong zero-density of statesof Hartree-Fock is caused by singular slope of density of states at the Fermi energy causing
practically no influences. Indeed, the optimized mixing parameter for bcc Fe is almost zero,117

which is consistent with the fact that the generalized gradient approximation is satisfactory
for this material. A strategy to determine the state-dependent mixing parameter more
precisely is to use a double-hybrid scheme where the correlation of density functional
theory is mixed with the MP2 correlation. However, this will significantly increase thecomputational cost, so that this problem remains to be solved. In addition to itinerant
electrons that can be appropriately described by the generalized gradient approximation,
localized electrons should be treated appropriately from first principles. For thedevelopment of a unified scheme, the knowledge obtained from the present results is
expected to be helpful. In addition, it should be pointed out that the present scheme todetermine the mixing parameters is applicable to determine the Hubbard U in the DFT+U
method, in particular for 3d transition metals where optimal U value significantly changes
depending on the chemical environment.
(b)

NiO

Density of states [arb. units]

Density of states [arb. units]

(a)

0

−10

−5

0

5

Energy relative to εVBM [eV]

10

FeO

0

−10

−5

0

5

Energy relative to εVBM [eV]10

Fig. 1: (a) Density of states for distorted-rocksalt NiO and (b) distorted-rocksalt FeO withthe type-II antiferromagnetic ordering. Solid lines are obtained with optimized mixing
parameter in hybrid density functional theory, whereas the dashed lines are results obtainedusing the generalized gradient approximation.
We also performed standard first-principles calculations for various magnetic materials,
where we focused mainly on the magnetocrystalline anisotropy. First, we calculated
Y-Fe-B compounds to examine the intersite couplings due to d electrons and the magnetic
anisotropy of Fe sublattice, because utilizing itinerant 3d electrons instead of localized 4f
electrons in exhibiting the magnetic anisotropy is an important new direction bothscientifically and technologically. In particular, the development of rare-earth-free magnets
will be made possible driven by fundamental understanding of the itinerant-electron
anisotropy. Since itinerant magnetic states are responsible to direct interactions amongmagnetic sites, d states are expected to be sensitive to lattice strain. In this study, we
examined strain effects on magnetic properties theoretically studied by first-principles
calculations for Y2Fe14B and YFe3, where Y is a prototypical f0 rare-earth element.
First-principles calculations on the basis of density functional theory within the generalized
gradient approximations are performed using the OpenMX code. The lattice constants a

118

and c of Y2Fe14B are changed from the equilibrium values independently, where we foundthe uniform compression enhances the perpendicular magnetic anisotropy. To clarify theorigin of this enhancement, we developed a method to decompose the magnetic-anisotropy
energy into contribution from each atomic site as well as from couplings among specificatomic orbitals, where the sum rule is satisfied by including off-site contributions in thesecond-order perturbation indirectly. As a result, we identified specific coupling among 3d
orbitals at the Fe j2 site to play an important role in enhancing the anisotropy, where boththe effect of the local density of states and the matrix elements among atomic orbitals are
essential [Fig. 2]. Our method enables us to study the anisotropy at microstructureinterfaces. In addition, our scheme can derive lattice spin models including the single ion
anisotropy, which can be used in micromagnetics and Monte Carlo simulations.

(a)

B

0.4
0.2
0
−0.2
−0.4

Y

LDOS (states/eV)

Fe

(b)

up

dx2-y2

down

†＊
dxy

0.4
0.2
0
−0.2
−0.4

＊
−4

−2

†

0

Energy (eV)2Fig. 2: (a) Crystal structure and spin-density distribution of Y2Fe14B. (b) Local density ofstates (LDOS) projected onto two 3d components of the Fe j2 site with (solid lines) andwithout (shades) compression. Pairs of daggers and asterisks indicate anisotropy coupling,respectively.
We also studied magnetic ultrathin films attracting attentions in fundamental surface
science. In contrast with the ground state of Fe bulk, the bcc structure with theferromagnetic order, distorted fcc Fe films have been grown onto Cu(001). In this system,
peculiar structural and magnetic properties of Fe have been reported. But its ground stateof magnetic structure is controversial. We performed first principle calculations in order tounderstand the electronic and magnetic structures of this system, where we have found thatthe ground state is not ferromagnetic but antiferromagnetic in the 4-monolayer (ML)system from the total energy analysis of allowable spin configurations. In the 5 to 7-ML
systems, an antiferromagnetic order is more stable than the ferromagnetic order, which isin good agreement with experimental results.
Furthermore, the origin of gigantic magnetic anisotropy in ε-type Fe2O3 was studied. We

119have found that the oxygen vacancies can play a significant role to exhibit magnetic
moments and large magnetic anisotropy.
3. Summary
We examined the dependence of the total energy of electronic systems on the number ofelectrons through the generalized Koopmans theorem. Deriving an exact condition to besatisfied for single-electron energy levels, we determined the mixing parameter appearing
in hybrid density functional theory. We applied this scheme for various magnetic materials
including distorted-rocksalt NiO as a charge transfer insulator and distorted-rocksalt FeO
as a Mott-Hubbard insulator.
We also performed first-principles calculations to investigate the magnetic anisotropy.
Strain effects for Y2Fe14B are clarified, where we developed a convenient and reliable
scheme to calculate single-ion anisotropy constants. This scheme is a promising tool to
combine first-principles calculations and spin lattice models for micromagnetics
simulations based on Landau-Lifshitz-Gilbert equations and Monte Carlo simulations.
4. Participants
PI: Yoshihiro Gohda
Collaborators: Zahra Torbatian, Yasutomi Tatetsu, Daisuke Hirai, Shinji Tsuneyuki, andTaisuke Ozaki
5. Publications
Original Papers (total 6 papers):
1) Z. Torbatian, T. Ozaki, S. Tsuneyuki, and Y. Gohda, “Strain effects on the magnetic
anisotropy of Y2Fe14B examined by first-principles calculations”, Appl. Phys. Lett. 105,
242403 (2014).
2) T. Nakamura, A. Yasui, Y. Kotani, T. Fukagawa, T. Nishiuchi, H. Iwai, T. Akiya, T.
Ohkubo, Y. Gohda, K. Hono, and S. Hirosawa, “Direct observation of ferromagnetism in
grain boundary phase of Nd-Fe-B sintered magnet using soft x-ray magnetic circulardichroism”, Appl. Phys. Lett. 105, 202404 (2014).
3) M. Kawamura, Y. Gohda, and S. Tsuneyuki, “Improved tetrahedron method for the
Brillouin-zone integration applicable to response functions”, Phys. Rev. B 89, 094515(2014).
Invited Talks (total 10 talks: 4 in English; 6 in Japanese):
1) Y. Gohda, “First-principles theory on magnetism at interfaces”, Collaborative Conference
on Materials Research 2014 (CCMR2014), (Incheon, Korea, Jun. 24, 2014).
2) Y. Gohda, “First-principles study of microstructure interfaces of Nd-Fe-B magnets”,
ESICMM-G8 Symposium on Next Generation Permanent Magnets, (Tsukuba, Japan, Jun.
2015, confirmed).
3) Y. Gohda, “Microstructure of Fe-based magnetic materials”, The 3rd OpenMX/QMAS
workshop, (Kashiwa, Japan, May 2015, confirmed).

120

A02: Dynamics of electron, spin and heat transport through single atomic
and molecular junctions
Manabu Kiguchi
Tokyo Institute of Technology
1. Introduction
A single molecular junction, where a single molecule bridges between metal electrodes,
is a low-dimensional material having two metal-molecule interfaces (see Fig. 1). Due toits geometrical and electronic characters, single molecular junctions can show novel
properties that are not observed in isolated molecular or bulk systems. The single
molecular junction has also attracted attention due to its potential application in
ultrasmall electronic devices. Since it is hard to see a single molecule bridging between
metal electrodes, we discuss the single
molecular junction by comparing theexperimentally obtained signal from thesingle molecular junction and calculation
results using some models. Here, we study
the dynamic properties of the carriers
Figure 1. Schematic representation of an
through single atomic and molecular atomic scale junction where a benzenedithioljunction by experiments and theoretical molecule is sandwiched between two Au
electrodes.
calculation.
2. Results[1] Single Molecular Resistive Switch
We investigated the possibility to exert dynamic control on molecular anchoring sites
and proposed a single-molecule resistive switch based on “switch of anchor in single
molecule” without altering the binding modes or chemical structures of themetal-molecule contact. For this purpose, we investigated a quaterthiophene-based
molecular wire (QT), bearing two thiophene rings at each end of the molecular
backbone as molecular anchors (T1 and T2) (Figure 2a). These two termini weresymmetrically separated by the“insulating” molecular layer, to prevent theintermolecular interactions among neighboring QT molecules and to promote themeasurement of single QT conductance. Because of this unique QT arrangement, thejunction could be switched between T1-T1, T1-T2, and T2-T2 configurations. This gave121

rise to three different conductance values: 0.05 G 0 (high-conductance state), 0.005 G 0(medium-conductance state), and 0.0005 G 0 (low-conductance state). By estimating the
distance between the contact points as well as tip-Au distance, we calculated that theresulting conductance was in line with the standard length dependence of off-resonant
tunneling, i.e. exponential decay. The most stable molecular structure of each contact
configuration was determined employing first-principles calculations (Figure 2b). The
theoretical studies demonstrated that the resistive switching was caused by varying theeffective length of the conducting π-conjugated backbone via a change in the molecularanchoring points T1 ↔ T2 (Figure 2c). The dynamic control of the three conductancestates was experimentally shown by mechanically controllable modulation. Supported
by DFT calculations we demonstrated that molecular materials bearing multiple
anchoring groups are a viable alternative for the development molecular devices withadvanced electronic functionality eg. resistive switching.

Figure 2. (a) Chemical structure of the quaterthiophene-based molecular wire. (b) One of DFT
models for electronic transport calculation (T2-T2 connection). (c) Schematic illustration ofswitch in the effective length of the conducting π-conjugate backbone via a change in theanchors.

[2] Electron transport through metal-ion wires
In recent years, electron transport though metal-atom wires has attracted the interest of
the scientific community. Their minute size and quantized electric properties makes
them especially interesting for the development of the next generation ultra-dense
electronic devices. Despite their appealing properties, electrical performance of
metal-ion wires has never been studied mainly hindered by the synthetic difficulties.
The preparation of metal-ion arrays with a fixed cross section and length between
nanogap electrodes is nontrivial. More recently, a number of synthetic protocols havebeen developed for the preparation of multinuclear complexes with ordered metal-ion

122arrays. Among them, [3 × n] gold-ion clusters, enclosed in self-assembled cages (Figure3a), are of particular interest because the stacking number n is uniquely determined byheight of the molecular cages. The molecular cages, that keep the metal ions highlyordered in a [3 × n] manner, are excellent frameworks for accommodating π-stacked
molecules and evaluating their electron-transport properties (Figure 3b).
Single-molecule conductance studies performed on the [3 × n] stacks of gold-ion
clusters n (n = 1, 2, or 3) showed a high conductance profile, in the range of 3~9 ×10-3
G 0 regardless of the stacking number. DFT studies further revealed that the metal-ion
arrays are highly conductive with a negligibly small distance dependence and anelectron transport efficiency comparable to that of metal-atom wires. The combinationof atomic scale transport experiments and DFT studies shed light on the potential
benefits of metal-ion molecular wires and on the fundamental aspects of
one-dimensional physics.

Figure 3. (a) [3 × n] gold-ion clusters, enclosed in self-assembled cages. (b) Schematic
illustration of [3 × 3] gold-ion clusters sandwiched by gold electrodes.

[3] Temperature dependence of the thermopower and its variation of the Au atomic
contact
Thermoelectric devices are of great interest because of their ability to efficiently harvest
waste heat and transform it into an electrical signal. While thermopower measurements
have been performed for various metal atomic contacts and molecular junctions, thetemperature dependence of the thermopower effect has not been satisfactorily studied.
Here, we simultaneously measured the thermopower and conductance of Au atomic
contacts during the self-breaking process in air at room temperature. These studies
showed a random fluctuation of the junction thermopower. This fluctuation was foundto be especially large and both positive and negative in sign for the junction in the lower
conductance regime. While the average thermopower was negligibly small, the standard
deviation of thermoelectric voltage linearly increased with the temperature difference123across the Au contacts. The standard deviation of the thermopower was inversely
related to the junction conductance with a dramatic increase below 1 G 0 . The increasein the standard deviation of thermopower was explained by the coherent backscatteringof the electrons near the contact. Above 320 K, the standard deviation of thermopower
was found to decrease, indicating the decrease in the density of defects for the Au
atomic contacts. The Au atomic contacts are low-dimensional nanomaterials, hence theatomic motion could be significantly enhanced above 320 K, leading to the decrease in
the density of defects near the contacts. The present study indicates that the standard
deviation of thermopower provides information about the thermopower of the Au
atomic contact and the atomic structure of the Au atomic contacts (e.g., crystallinity anddistribution of defects).
These works are done with the collaboration with Prof. Fujii and students in Tokyo
Institute of Technology, Prof. Fujita and Prof. Watanabe (Tokyo University), Dr.
Nakamura (AIST), and Prof. Tada (Tokyo Institute of Technology), and Dr. Sugiyasu(NIMS).
3. Selected publications (Total 35 publications)[1] M. Kiguchi, T. Ohto, S. Fujii, K. Sugiyasu, S. Nakajima, M. Takeuchi, and H.
Nakamura “Single molecular resistive switch obtained via sliding multiple
anchoring points and varying effective wire length” J. Am. Chem. Soc., 136,
7327-7332 (2014). Selected as the cover picture and picked up by THE NIKKANKOGYO SHIMBUN.
[2] M. Kiguchi, J. Inatomi, Y. Takahashi, R. Tanaka, T. Osuga, T. Murase, M. Fujita, T.
Tada, S. Watanabe “Highly Conductive [3xn] Gold Ionic Clusters Enclosed within
Self-Assembled Cages” Angew. Chem. Int. Ed., 52, 6202-6205, (2013).
[3] R. Matsushita, S. Kaneko, S. Fujii, H. Nakamura, M. Kiguchi “Temperature
Dependence of the Thermopower and its Variation of the Au Atomic Contact”
Nanotechnology, 26, 045709, (2015).
4. Selected invited talks (Total 12 talks)[1] M. Kiguchi, “Experimental and theoretical study of electron transport throughsingle π-stacked systems and ionic wires” 2nd International Symposium on
Computics: Quantum Simulation and Design, Tokyo, 2014/12/1.

124

Electron-Phonon Coupling Dynamics and High-Sensitive Detection of
Nano-Scale Phonon Wavepackets in Carbon-Based Materials
Jun Takeda
Department of Physics, Graduate School of Engineering, Yokohama National University
1. Introduction
To unveil the non-equilibrium dynamics far from equilibrium states in new functional materials,
ultrafast laser spectroscopy is one of most powerful tools. In this research, we have mainly utilized
three types of new spectroscopic methods we have recently developed: i.e., frequency-resolved
coherent phonon spectroscopy, THz time-domain spectroscopy and single-shot imaging spectroscopy.
Using these methods, we have investigated dynamics of nano-scale optical phonon wavepackets ingraphene, resonant enhancement of first- and second-order coherent phonons due to van-Hove
singularities in metallic single-walled carbon nanotubes (M-SWCNTs) and photoinduced
phase-changes of chalcogenide glasses. Some other related photoinduced phenomena were also
revealed.

2. Results
2.1 Nano-scale optical phonon wavepacket in graphene
Coherent nano-scale optical phonon wavepacket was excited using ultrashort laser pulses via
impulsive inter-Dirac-cone scattering at graphene edges and defects. In contrast to the normal coherent
phonons generated with first-order Raman process such as G-mode, the D-mode generated with thesecond-order Raman process showed the phonon chirp, and the amplitude and frequency strongly
changed depending on the probe wavelength. The results can be discussed on the basis of doubly
resonant Raman scattering process; the D-mode couples strongly with the electronic states near the
Fermi level, and thereby spatiotemporally modulates the electronic states with nanometer-spatial andfemtosecond-temporal resolutions. As a result, the generated D-mode coherent phonons can propagate
as nanoscale optical-phonon wavepackets. (See, Fig. 1)

2-2. Resonant properties of metallic single-walled carbon nanotubes
High-frequency coherent phonons resonantly excited in M-SWCNTs were investigated via
spectrally resolved pump-probe spectroscopy. In addition to first-order coherent phonons such as
RBM, M and G modes, we clearly observed second-order high-frequency coherent phonons of 2D and
2G modes, which can be regarded as squeezed phonons. We found that the amplitudes of the RBM, G
and 2D modes were resonantly enhanced at specific wavelengths: the maximum resides at a125

wavelength whose energy is smaller than that of the van-Hove singularities in M-SWCNTs by anamount corresponding to the phonon energy. Further, the 2D mode has stronger enhancement than theother first-order Raman modes. These results indicate that the enhancement originates from a
Stokes-stimulated Raman-scattering process at van-Hove singularities, and that efficient resonanceenhancement occurs for the 2D mode, possibly through double resonance due to the trigonal warping
effect and strong electron-phonon coupling due to the Kohn anomaly. (See, Fig. 2)

2-3. Ultrafast amorphization in chalcogenide glasses

Using broadband single-shot imaging spectroscopy, we observed an irreversible ultrafast
crystalline-to-amorphous phase transition in chalcogenide alloy thin films. The absorbancechange that accompanied the ultrafast amorphization was measured via single-shot detection
even for laser fluences above the critical value, where a permanent amorphized mark wasformed. The observed rise time to reach the amorphization was found to be 130~200 fs,
which was in good agreement with the half period of the A1 phonon frequency in theoctahedral GeTe6 structure. This result strongly suggests that the ultrafast amorphization canbe attributed to the rearrangement of Ge atoms from an octahedral structure to a tetrahedral
structure. Based on the dependence of the absorbance change on the laser fluence, we could
also discuss the stability of the photoinduced amorphous phase. (See, Fig. 3)

126

3. Members of This Research Project
Jun Takeda (Project Reader; Professor, Yokohama National University)
Maki Suemitsu (Member; Professor, Tohoku University)
Ikufumi Katayama (Member; Associate Professor, Yokohama National University)
Masahiro Kitajima (Member; Senior Researcher, LxRay Co. Ltd.)
Yauo Minami (Member; Assistant Professor, Yokohama National University)

4. Research Activities
Original Papers (total: 11 papers)
1) K. Yoshioka, Y. Minami, K. Shudo, T. D. Dao, T. Nagao, M. Kitajima, J. Takeda, and I.
Katayama, "Terahertz-Field-Induced Nonlinear Electron Delocalization in Au Nanostructures",
Nano Lett. 15, pp. 1036-1040 (2015).
2) K. Sato, K. Tahara, Y. Minami, I. Katayama, M. Kitajima, H. Kawai, K. Yanagi, and J. Takeda,"Resonant Enhancement of First- and Second-Order Coherent Phonons in Metallic Single-Walled
Carbon Nanotubes", Phys. Rev. B 90, 235435: pp.1-5 (2014).
3) Y. Minami, J. Takeda, T. D. Dao, T. Nagao, M. Kitajima, and I. Katayama, "Nonlinear Electron
Dynamics of Gold Ultrathin Films Induced by Intense Terahertz Waves", Appl. Phys. Lett. 105,
241107: pp.1-4 (2014).
4) J. Takeda, W. Oba, Y. Minami, T. Saiki, and I. Katayama, "Ultrafast Crystalline-to-Amorphous
Phase Transition in Ge2Sb2Te5 Chalcogenide Alloy Thin Film Using Single-Shot Imaging
Spectroscopy", Appl. Phys. Lett. 104, 261903: pp. 1-4 (2014).
5) Y. Minami, M. Yamaki, I. Katayama, and J. Takeda, "Broadband Pump-Probe Imaging
Spectroscopy Applicable to Ultrafast Single-Shot Events", Appl. Phys. Exp. 7, 022402: pp. 1-3(2014).

127

6) J. Takeda, I. Katayama, K. Shudo, and M. Kitajima, (Invited paper) "Coherent Phonon and
Surface-Enhanced Raman Scattering Dynamics in Solids", J. Lumin. 152, pp. 23-27 (2014).
7) T. Kawano, I. Katayama, J. Ohara, M. Ashida, and J. Takeda, "Intermolecular THz Vibrations
Relevant to Optically and Thermally Induced Magnetic Phase Transitions in the Strongly
Correlated Organic Radical TTTA", J. Phys. Soc. Jpn. 83, 014713: pp. 1-6 (2014).
8) I. Katayama, K. Sato, S. Koga, J. Takeda, S. Hishita, H. Fukidome, M. Suemitsu, and M. Kitajima,"Coherent Nanoscale Optical-Phonon Wavepacket in Graphene Layers", Phys. Rev. B 88, 245406:
pp. 1-5 (2013).
9) K. Shudo, T. Satake, T. Shimatsu, M. Uchiyama, and J. Takeda, "Electronic Excited States of aStrongly Correlated Organic Radical 1,3,5-trithia-2,4,6-triazapentalenyl (TTTA) Adsorbed on aSi(001) Surface", Phys. Chem. Chem. Phys. 15, pp. 19213-19220 (2013).
10) Y. Minami, Y. Hayashi, J. Takeda, and I. Katayama, "Single-Shot Measurement of a Terahertz
Electric-Field Waveform Using a Reflective Echelon Mirror", Appl. Phys. Lett. 103, 051103: pp.
1-4 (2013).
11) I. Katayama, S. Koga, J. Takeda, S. Hishita, D. Fujita, and M. Kitajima, "Electron-Phonon
Coupling and Defect Scatterings in Ar+-ion Implanted Graphite", J. Cer. Soc. Jpn. 121, pp.
291-294 (2013).
Selected Proceedings & Reviews (total: 13 papers)
1) J. Takeda, Y. Minami, and I. Katayama, (Review article in Japanese) "Time-Frequency
Two-Dimensional Imaging Spectroscopy Using a Reflective Echelon Mirror", Rev. Laser Eng. 43,(2015) to be published.
2) K. Maekawa, K. Sato, Y. Minami, I. Katayama, J. Takeda, K. Yanagi, and M. Kitajima,"Electrochemical Control of Coherent Phonon Generations in Single-walled Metallic Carbon
Nanotubes", Ultrafast Phenomena XIX: Proceedings of the 19th International Conference on
Ultrafast Phenomena, (Springer Proceedings in Physics), pp. 356-359 (2015).
3) W. Oba, I. Katayama, Y. Minami, T. Saiki, and J. Takeda, "Single-shot Real-time Observation of
Ultrafast Amorphization in Ge2Sb2Te5 Thin Film", Ultrafast Phenomena XIX: Proceedings of the
19th International Conference on Ultrafast Phenomena, (Springer Proceedings in Physics), pp.
352-355 (2015).
4) Y. Minami, T. D. Dao, T. Nagao, J. Takeda, M. Kitajima, and I. Katayama, "Nonlinear Carrier
Responses in Gold Thin Films Induced by Intense Terahertz Waves", Ultrafast Phenomena XIX:
Proceedings of the 19th International Conference on Ultrafast Phenomena, (Springer Proceedingsin Physics), pp. 626-629 (2015).
5) K. Araki, Y. Minami, T. D. Dao, T. Nagao, J. Takeda, M. Kitajima, and I. Katayama, "Nonlinear
Carrier Dynamics in Semi-Metal Bismuth Induced by Intense Terahertz Field", Ultrafast
Phenomena XIX: Proceedings of the 19th International Conference on Ultrafast Phenomena,(Springer Proceedings in Physics), pp. 633-636 (2015).
6) I. Katayama, Y. Minami, and J. Takeda, (Review article in Japanese) "Terahertz Time-Domain
Spectroscopy of Thin Films and Surfaces", J. Surf. Sci. Jpn. 35, pp. 680-685 (2014).
Selected Invited Talks (total: 7 talks)
1) I. Katayama, W. Oba, Y. Minami, T. Saiki, and J. Takeda (Invited Talk), "Ultrafast Dynamics ofPhase-Change Materials Observed Using Echelon-Based Single-Shot Setups", The 26th
Symposium on Phase Change Oriented Science (PCOS2014), (December 4-5, Hamamatsu, Japan).
2) I. Katayama, K. Shudo, J. Takeda, and M. Kitajima (Invited Talk), "Nanoscale and Femtosecond
Phonon Dynamics Observed with Surface Enhanced Raman Scattering", Ultrafast Surface
Dynamics 2013, (May 28-31, 2013, Colorado, USA).
3) J. Takeda, I. Katayama, K. Shudo, and M. Kitajima (Invited Talk), "Coherent Phonon and
Surface-Enhanced Raman Scattering Dynamics in Solids", 18th International Conference on
Dynamical Processes in Excited States of Solids, II-1 (August 4-9, 2013, Fuzhou, China).

128

“First-principle molecular dynamics analyses on nonlinear optical properties
change of Si nano-dot”
Yasuteru Shigeta1,3,4, Kenji Shiraishi 2, Masayoshi Nakano3
1

Graduate School of Pure and Applied Sciences, University of Tsukuba
2Graduate School of Engineering, Nagoya University
3
Graduate School of Engineering Sciences, Osaka University
4
CREST, JST

Reviewed Paper (28 including 2 review articles)
1. H. Harada, Y. Takano, T. Baba, Y. Shigeta, “Simple, Yet Powerful Methodologies forConformational Sampling of Proteins”, Physical Chemistry Chemical Physics (invited featurearticle) 17, 6155-6173 (2015).
2. H. Ando, Y. Shigeta, T. Baba, C. Watanabe, Y. Okiyama, Y. Mochizuki, M. Nakano, “Hydration
Effects on Enzyme-Substrate Complex of Nylon Oligomer Hydrolase: Inter-Fragment Interaction
Energy Study by the Fragment Molecular Orbital Method”, Molecular Physics 113, 319-326(2015).
3. R. Harada, Y. Takano, Y. Shigeta, "Enhanced Conformational Sampling Method for Proteins
Based on the TaBoo SeArch (TBSA) Algorithm: Application to the Folding of a Mini-protein,
Chignolin", Journal of Computational Chemistry 36, 763-772 (2015).
4. K. Okuno, Y. Shigeta, R. Kishi, M. Nakano, “Theoretical design of solvatochromism switching by
photochromic reactions using donor-acceptor disubstituted diarylethene derivatives with oxidized
thiophene rings”, Physical Chemistry Chemical Physics 17, 6484-6494 (2015).
5. T. Baba, M. Boero, K. Kamiya, H. Ando, S. Negoro, M. Nakano, Y. Shigeta, “Unraveling thedegradation of artificial amide bonds in nylon oligomer hydrolase: from induced-fit to acylation
processes”, Physical Chemistry Chemical Physics, 17, 4492-4504(2015).
6. T. Matsui, Y. Kitagawa, M. Okumura, Y. Shigeta, “Accurate Standard Hydrogen Electrode
Potential and Applications to the Redox Potentials of Vitamin C and NAD/NADH”, Journal of
Physical Chemistry A, 119, 369-376 (2015).
7. Y. Shigeta, H. Harada, M. Kayanuma, M. Shoji, “Quantal cumulant dynamics for real-time
simulations of quantum many-body systems”, International Journal of Quantum Chemistry(invited review), 115, 300-308 (2015). DOI: 10.1002/qua.24820
8. R. Harada, T. Nakamura, Y. Takano, Y. Shigeta, “Protein folding pathways extracted by Outlier
FLOODing method (OFLOOD)”, Journal of Computational Chemistry, 36, 97-102 (2015).
9. R. Nakamura, Y. Shigeta, K. Okuno, M. Fukushima, M. Hasegawa, S. Suzuki, M. Kozaki, K.
Okada, M. Nakano, “Substitution Effects on Optical Properties of Iminonitroxide- substituted
Iminonitroxide Diradical”, Molecular Physics 113, 267-273 (2015).
10. T. Baba, T. Matsui, K. Kamiya, M. Nakano, Y. Shigeta, "A density functional study on pKa of
small polyprotic molecules”, International Journal of Quantum Chemistry, 114, 1128–1134(2014).
11. T. Baba, R. Harada, M. Nakano, Y. Shigeta, "On the induced-fit mechanism of substrate-enzyme
binding structures of Nylon-oligomerhydrolase”, Journal of Computational Chemistry 35,
1240-47 (2014).
12. K. Kamiya, T. Baba, M. Boero, T. Matsui, S. Negoro, Y. Shigeta, “A Nylon-oligomer Hydrolase
Promoting Cleavage Reactions in Unnatural Amide Compounds”, Journal of Physical Chemistry
Letters 5, 1210-1216 (2014).

12913. R. Harada, Y. Takano, Y. Shigeta, “Fluctuation Flooding Method (FFM) for accelerating
conformational transitions of proteins”, Journal of Chemical Physics 140, 125103 (2014).
14. Y. Kitagawa, T. Matsui, Y. Nakanishi, Y. Shigeta, T. Kawakami, M. Okumura, K. Yamaguchi,“Theoretical studies of electronic structures, magnetic properties and electron conductivities of
one-dimensional Nin (n = 3, 5, 7) complexes”, Dalton Transaction, 42, 16200-16208 (2013).
15. S. Muhammad, H. Xu, Z. Su, K. Fukuda, R. Kishi, Y. Shigeta, M. Nakano, “A new type oforganic-inorganic hybrid NLO-phore with large off-diagonal first hyperpolarizability tensors: A
two-dimensional approach”, Dalton Transaction, 42, 15053-15062 (2013).
16. K. Koyanagi, Y. Kita, Y. Shigeta, M. Tachikawa, “Binding of a positron to nucleic-acid basemolecules and their pairs”, ChemPhysChem. 14, 3458-3462 (2013).
17. T. Inui, Y. Shigeta, K. Okuno, T. Baba, R. Kishi, M. Nakano, “Finite-Field Method with Unbiased
Polarizable Continuum Model for Evaluation of the Second Hyperpolarizability of an Open-Shell
Singlet Molecule in Solvents”, Journal of Computational Chemistry, 34, 2345-2352 (2013).
18. K. Okuno, Y. Shigeta, R. Kishi, M. Nakano, “Non-empirical tuning of CAM-B3LYP functional intime-dependent density functional theory for excitation energies of diarylethene derivatives”,
Chemical Physics Letters, 585, 201-206 (2013).
19. K. Okuno, Y. Shigeta, R. Kishi, M. Nakano, “Photochromic Switching of Diradical Character:
Design of Efficient Nonlinear Optical Switches”, Journal of Physical Chemistry Letters, 4,2418-2422 (2013).
20. T. Matsui, Y. Kitagawa, Y. Shigeta, M. Okumura, ”A Density Functional Theory Based Protocol
to Compute the Redox Potential of Transition Metal Complex with the Correction of
Pseudo-Counterion: General Theory and Applications”, Journal of Chemical Theory andComputation, 9, 2974–2980 (2013).
21. K. Koizumi, M. Boero, Y. Shigeta, A. Oshiyama, “Atom-Scale Reaction Pathways andFree-Energy Landscapes in Oxygen Plasma Etching of Graphene”, Journal of Physical Chemistry
Letters, 4, 1592–1596 (2013).
22. Y. Shigeta, T. Inui, T. Baba, K. Okuno, H. Kuwabara, R. Kishi, M. Nakano, “Qantal Cumulant
Mechanics and Dynamics for Multi-Dimensional Quantum Many-Body Clusters”, International
Journal of Quantum Chemistry, 113, 348-355 (2013).
23. T. Matsui, Y. Kitagawa, M. Okumura, Y. Shigeta, S. Sakaki, “Consistent Scheme for ComputingStandard Hydrogen Electrode and Redox Potentials”, Journal of Computational Chemistry, 34,
21-26 (2013).
24. M. Nakano, R. Kishi, H. Fukui, T. Minami, K. Yoneda, S. Minamide, Y. Inoue, T. Yamada, S. Ito,
S. Muhammad, Y. Shigeta, T. Kubo, B. Champagne, “Diradicalology in Third-Order Nonlinear
Optical Systems: Second Hyperpolarizabilities of Acetylene-Linked Phenalenyl Based
Superpolyenes”, International Journal of Quantum Chemistry, 113, 585-591 (2013).
25. S. Muhammad, T. Minami, H. Fukui, K. Yoneda, S. Minamide, R. Kishi, Y. Shigeta, M. Nakano,“Comparative Study of Diradical Characters and Third-Order Nonlinear Optical Properties ofLinear/Cyclic Acenes versus Phenylenes”, International Journal of Quantum Chemistry, 113,592-598 (2013).
26. K. Yoneda, S. Minamide, T. Yamada, S. Ito, T. Minami, R. Kishi, Y. Shigeta, M. Nakano,“Antidot Effects on the Open-Shell Characters and Second Hyperpolarizabilities of Rectangular
Graphene Nanoflakes”, International Journal of Quantum Chemistry, 113, 605-611 (2013).
27. A. Konishi, Y. Hirao, K. Matsumoto, H. Kurata, R. Kishi, Y. Shigeta, M. Nakano, K. Tokunaga,
K. Kamada, T. Kubo, “Synthesis and Characterization of Quarteranthene: Elucidating theCharacteristics of the Edge State of Graphene Nanoribbons at the Molecular Level”, Journal of
the American Chemistry Society, 135, 1430-1437 (2013).
28. S. Muhammad, K. Fukuda, T. Minami, R. Kishi, Y. Shigeta, M. Nakano, “Interplay between the130Diradical Character and Third-Order Nonlinear Optical Properties in Fullerene Systems,
Chemistry-A European Journal, 19, 1677-1685 (2013).
Invited talks (8 International conferences)
1. Y. Shigeta, “A Molecular Design of Nonlinear Optical Properties and Conductivity Switches on
the Basis of Open-shell Nature”, 2014 Workshop on Innovative Nanoscale Devices and Systems(WINDS), Nov. 30th-Dec. 5th 2014, Hawaii, USA.
2. Y. Shigeta, “Protein Folding Processes Detected by Enhanced Sampling Techniques”, 19thInternational Workshop on Quantum Systems in Chemistry and Physics, Nov. 11th-17th 2014,
Tamsui, Taiwan.
3. Y. Shigeta, “Towards Theoretical Design of Catalytic Activities of Enzymes”, The Asia Hub for
e-Drug Discovery Symposium 2014, Nov. 11th-12th 2014, in Sichuan University (SCU), Chengdu,
China.
4. Y. Shigeta, “First-principles analysis on enzymatic degradation of nylon”, 2nd World Congress on
Petrochemistry and Chemical Engineering, Oct. 27th-29th 2014, Las Vegas, USA.
5. Y. Shigeta, “Computational Studies on Redox Potential of Metal Complexes and Model
Cofactors”, International Conference on Synthetic Metals, June 30th-July 5th 2014, Turku, Finland.
6. Y. Shigeta, “Theoretical Studies on Reaction Mechanism of An Enzyme Cleaving Unnatural
Amide Compound”, Symposium on Molecular Science and Synthesis of Functional Molecules for
Next Generation, Mar.10th-11th 2014, Hiroshima, Japan.
7. Y. Shigeta, “Quantum chemical studies for estimation of pKa and redox potential”, World
Congress on Petrochemistry and Chemical Engineering, Nov. 17th-19th 2013, San Antonio USA.
8. Y. Shigeta, “Free Energy Analyses on Cluster Deformations by Cumulant Mechanics”, 8thCongress of the International Society for Theoretical Chemical Physics, Aug. 25th-31st 2013,
Budapest, Hungary.
Award and others (1 award and 1 press release)
1. [Award] Yasuteru Shigeta, May 2013, 2nd Young researcher’ award of Osaka University,“Dynamics Theory Including Fluctuation” (Osaka University, Japan).
2. [Press release] K. Kamiya, T. Baba, M. Boero, S. Negoro, Y. Shigeta, Nikkei Biotechnology
Onlinehttps://bio.nikkeibp.co.jp/article/news/20140419/175556/

131

Letter(1) Photochromic Switching of Diradical Character: Design of Efficient Nonlinear Opticalpubs.acs.org/JPCL
Switches [19]In conventional designs for enhancing second hyperpolarizabilities, γ, one adopted an elongationof pi-conjugation, introducing donor-acceptor substituents on the both side of molecules, and so on. On
Photochromic
Switchingof Diradical
Character:
Design ofEﬃcient
theother hand, Open-shell
singlet diradical
molecules havebeen widely investigatedbecausethey arekey
tounderstanding
thenatureofchemicalbonds.
We
have
proposedadesignprinciple
forγ, byNonlinear Optical Switches
referring a diradical character, y (0<y<1), which an index of the instability of chemical bonds. In our
Katsukiprinciple,
Okuno, Yasuteru
Shigeta,
Ryohei
Kishi, andMasayoshi
designγ becomes
large atintermediatey region.
ThisNakano*y-γ correlation has been theoreticallyconfirmed
variousmolecules
as School
H2 moleculeandScience,p-quinodimethane.
This theoretical
Department offor
Materials
Engineering
Science,such
Graduateof Engineering
Osaka University, Toyonaka,
Osaka
560-8531, Japanprediction
has subsequently been verified by the experimental results of gigantic two-photon absorption
S Supporting
cross
section, which
is proportional to the imaginary part of γ, of the polycyclic aromatic hydrocarbons*
Information
with intermediate y.
ABSTRACT:
diradical
molecules
been widely
investigatedWe haveOpen-shell
proposedsinglet
a newconcept
for havereversible
switchingofbecause they
are key to understanding
the nature ofbonds. We proposea newdiradical
characterof a moleculebychemical
photochromic
reaction.
concept for reversible switching of diradical characteran index of the instability of
Photochromic
diarylethenewithvariousopen-shell
singletchemical bondsofa molecule derivatives
by photochromic
reaction.
Photochromic
diarylethene
diradical
charactersaretheoreticallydesigned,andtheirphotochromic
derivatives with various open-shell singlet diradical characters are theoretically designed,
and their photochromic
diradical characterswitchingare
behaviors
are clariﬁed.
Theseresults
results
diradical
character switching
behaviorsclarified.
These
contributetodesigninghighlyeﬃcient
third-ordernonlinear
opticalswitchingsubstances
contribute to designing highly efficient third-order nonlinear optical
based on the correlation between the diradical character and second hyperpolarizability.
switching substances based on the correlation between the diradical
character, y, and second hyperpolarizability, γ. With this idea, we
achieve switching ratio γc/γo ~ 25 for model compounds withintermediate y for a closed-ring isomer and y=1 for open-ring isomer.
SECTION: Molecular Structure, Quantum Chemistry, and General Theory(2) Theoretical design of solvatochromism switching by photochromic reactions using
donor-acceptor disubstituted diarylethene derivatives with oxidized thiopehene rings [4]Singlet diradical compounds have been widely investigated“photochromism”. Photochromism is deﬁned as reversible

isomerization
betweenisomers, where
where atatleastone
from both theoreticaland experimental
becausePhotochromism
is defined
as a researchers
reversible
isomerization
betweentwotwoisomers,
least
21−23
directionofreversible
reactionisaphotoinduced
one.
they
have
peculiar
electronic
structures
and
are
key
to
one direction of the reversible reaction
is induced by light irradiation. Among the various
Among various photochromic compounds, diarylethene deunderstanding the nature of chemical bonds.1−10 The diradical
photochromic
compounds,
diarylethene
derivatives
undergocycloreversion
rivatives
undergo cyclization
cyclization and andcycloreversion
in 6π-electronature of a molecule is quantiﬁed by the diradical character
24−26
cyclic openreactions.and closed-ringThrough these reactions,reactionsreactions.
Since bothisomersπ-conjugationof the(y)an indexinof 6π-electrocyclic
the instability of chemicalbondsranging
11,12
lengthinamoleculecanbereversibly
changed,where
typicallyfrom0(closed-shell)
to
1
(pure
diradical),
so
that,
in
other
diarylethene derivatives are thermally stable, the photochromic reactions occur only via lightthe open-ring isomer with less extended π-conjugation iswords, 1 − y indicates an “eﬀective bond order”. Interesting
irradiation.
The various properties of the diarylethene
derivatives
such asisomer
absorption
colorless,while the closed-ringwith more spectra,
extended πphysicochemical properties of singlet diradical compounds are
refractive
indices,dielectric
constants,
redox
potentials,
and
so
on
are
quickly
andprinciple,
reversibly
conjugation is colored. Following thisvariousfound to originate in the diradical character. Among them, wediarylethene
derivatives
with modiﬁcationof aryl rings haveswitched
throughthesereactions.
We havedesigned
severaldiarylethene
have theoreticallyclariﬁedthe photochromic
correlations betweentheIn regard to thewhich
open-shelldiarylethene
been synthesized.24−26
diradical character
excitation thiophene
energies/properties,
derivatives
withandoxidizedrings andand donor-acceptor
substituents,show
thederivatives, Matsuda and Irie reported that the intramolecularthus have revealed a remarkable enhancement of the second
solvatochromism
switching by photochromic magnetic
reactions,usingin time-dependent
densityinteraction
high spin states of diarylethene
hyperpolarizability (γ), which is the microscopic origin of thirdfunctional
theory(TD-DFT)
methodwiththepolarizable
continuum
model(PCM).
It isderivatives with nitronyl nitroxide radical sources can reversiblyorder nonlinear optical (NLO) property in the intermediate
27
13−15
In
changedalong
thecyclizationandcycloreversion
reactions.
foundonly the This
open-ringisomers
exhibitthe solvatochromism, while the closed-ring
diradicalthat
character.
theoreticalprediction
has
sharp
contrasttothistraditional
π-conjugation
based
principle,
subsequentlydo
beennot
veriﬁedexperimental resultsofisomers
in bythethe examined
UV-vis.
spectral region. The mechanism of the
we here propose a reversible change scheme of diradical
gigantic two-photon absorption cross section, which issolvatochromism
behavior and its switching process
are y,clarified
from bondthe viewpoint
theof
characteri.e., the eﬀective
order, within aofrange
proportional to the imaginary part of γ, of the polycyclic
16−18charge-transfer(CT)
excitation
from
the
donor
to
the
acceptor
substituents.
We
demonstrate
closed-shell
(y
=
0)
and
pure
diradical(y=1)
statesbythearomatic hydrocarbons with intermediate y.
Although thischangeofmolecular
structurethroughthephotochromic
thatthis CT between
excitation
be controlled
by choosing
new correlationy and γcan
is expected
to open a new
F F
F F
1,14
F
F
F to γ Fswitching
reactions, and thus
apply
this y tuning method
paradigm for highly
eﬃcient
still acceptor
appropriate
pairs
ofNLO
thematerials,
donor there
andarethe
F
F
F
F
Vis.
by light irradiation.
OO
no reports on the methods for reversible switching or tuning ofOO
substituents
on
the
basis
of
the
orbitalcorrelation
O S the photochromic
S O
O Sand the SdiradicalO
By
combining
reaction
the y and γ values except for switching external strong staticUV
19
diagrambetweenthediarylethene
derivatives
andthecharacter
basedprinciple,
we
proposeanewconcept
forlightelectric ﬁeld applied to symmetric open-shell singlet systems,
D
D
A
A
20
inducedreversible
switchingofγ
valuesthroughtuningydonor-acceptor
substituents,whichisconstructedfromwhich leads to asymmetric open-shell singlet systems. The(Figure 1). There ON
can be nine classes (A−I) of OFF
systems, whichdevelopmentofsuchreversible
tuningmethod
as
well
as
of
the topologies and the orbital energies of the molecularmolecular design of practical compounds would contribute toD: donorA: acceptor
orbitalsprimarily contributing to the excitations.
opening up a new class of highly eﬃcient NLO switching.
As a typical example of such tuning methods of properties
through the change of molecular structures, we can consider© XXXX American Chemical Society

132
2418

Received: June 13, 2013
Accepted: July 11, 2013 : π conjugation

dx.doi.org/10.1021/jz401228c | J. Phys. Chem. Lett. 2013, 4, 2418−2422

Articlepubs.acs.org/JCTC(3) Density
A Density Functional
TheoryBased Protocol
to Compute
the Redox
Potential of theTransition
A
Functional
TheoryBasedProtocolto Compute
Redox
Metal Complex
the Correctionof Pseudo-Counterion:General
Theory and Applications[20]
Potential
ofwithTransition
Metal
Complex with the Correctionof PseudoOxidation andGeneral
reduction processes
one Applicationsof the most fundamental reactions in the fields ofCounterion:
Theoryareandbiology, physics, and chemistry. In this regard, it is necessary to theoretically estimate the redox

,†
†
‡,§
†
Toru
Matsui,*Yasutaka
Kitagawa,toward
Yasuteru
Shigeta,and Mitsutaka
Okumurapotentialof a given
compounddesignof compoundswith desired
redox potential. Many†computational
works employed
this ofstandardhydrogenelectrode(SHE) potential
to evaluate
the redox
Department of Chemistry,
Graduate School
Science, Osaka
University,
1-1 Machikaneyama,
Toyonaka,
Osaka 560-0043,
Japan
potential.
are many
to employ
this SHE
potential,
i.e. (1) the
SHE potential
Department However,of Materials there
Engineering
Science,problems
Graduate Schoolof Engineering
Science,
Osaka University,
1-3 Machikaneyama,
Toyonaka, Osaka
560-8531,
Japan
depends
on themodelsand the computational methods employed, (2) many experimental papers§
Core Research
for Evolutional
Science andranging
Technology(CREST),
andremoveTechnology(JST)
Agency, Kawaguchi,
reported
different
SHE potentialfrom4.05 toJapan
4.73Science
V. To
these
difficulties,
we have
332-0012
Saitama,Japan
proposed a new method to evaluate the standard hydrogen electrode potential in aqueous solution. We
S Supporting*have
added aInformation
charge correction term for transition metal complexes called “Pseudo Counter Ion
Solvation” (PCIS) model. We here evaluate the redox potential of a wide variety of transition metal
ABSTRACT:The
We meanproposeabsolutean accurateerror
schemeto evaluate
the redox
potential of redox
a widecomplexes.
(MAE)
towardexperimental
varietyoftransition
metal
complexesbyaddingacharge-dependent
correction
term
potentials of charged complexes is considerably reduced from 0.81forVa
counterion around the charged complexes, which is based on Generalized Born theory, to(maximum1.22TheV) mean
to 0.22
V (maximum
error
0.50 experimental
V). We foundthe solvationerror
energy.
absolute
error (MAE)
towardredoxaremarkableexchange-correlation
functionalon theerror
results
potentials of charged
complexes is considerably
reduced dependence
from 0.81 V (maximum
1.22V) to than
0.22 Vthe(maximumerror ones.
0.50 V).
We combination
found a remarkableexchange-correlation
ratherbasis setTheof Wachters+f(forfunctionalon the results
rather atoms)
than the basis
ones.
The combinationofmetal)
anddependence
6-31++G(d,p)(for otherwithsetthe
B3LYP
functional
Wachters+f (for metal) and 6-31++G(d,p) (for other atoms) with the B3LYP functional
gives
the least
least
MAE
0.15
V testfor complexes.
the test This
complexes.
schemeisgives theMAE
0.15 Vfor thescheme is Thisapplicableto otherapplicable
other transition
solvents,metal
andcomplexes
heavier suchtransition
metal
complexes(M1 = Cr,suchMo,solvents, andtoheavier
as M1(CO)5(pycn)Pd, Pt),=and(bpy)W),= Fe, Ru, Os)2 withthe samequality.
M2(mnt)(pycn)
Cr,M3Mo,(M(II)
= Ni,
Pd,
asW),
M(I)
(CO)2 5(M2 = Ni, (M(I)
3 (M3 M(II)(mnt)
Pt), and M(III)(bpy)3 (M(III) = Fe, Ru, Os) with the same quality.
‡(4) Unraveling the degradation of artificial amide bonds in nylon oligomer hydrolase: fromobserved by Trasatti.9 Many theoretical10−12 and experimental
1. INTRODUCTION
induced-fit to acylation processes [5]papers13−17 reported various values ranging from 4.1 to 4.8 V
Toward the design of new materials on a molecular level, it isdepending on their experimental conditions, such as thehelpfulto understandthe propertiesof a single
molecule.
Nylon-6
is a synthetic
artificial
fiber
and Themade in40’s. It was unprecedented discovery thattemperature, electrode, and so on. From the view of
redox
potential is onehydrolasesof the most that
fundamental
quantities
ofNylon-oligomer
degradethe artificial
fiber
were found in 70’.
This itis isa also
mystery
thattotheoretical/computational
chemistry,diﬃculttransition metal complexes (TMCs), which clariﬁes whether
bacteria
aquires canan easily
ability
the nylon
i.e. just
30 years.
Tobecause
elucidate
how theneeds
nylon
the SHEpotential,
this property
thethetarget molecule
taketoaneat
electronor not. so
Theearly,obtain
Gibbs
free
energy(GFE)ofaprotonintheaqueousphase.
In
oligomerhydrolase(NylB)
acquiresits
peculiar
degradation
activity
towards
non-biologicalamide
importance of measuring the redox potential of TMCs is not
order
to
avoid
such
problems,
wealso
reportedanewschemelimitedinorganic
chemistry.
redox potential
is one ofprocesses going from the induced-fit upon substratebonds,toweinspected
the Theunderlying
enzymatic
computing the SHE potential by combining experimental
thekey pieces
of informationneeded for awe
metal
cluster to theformutationalbindingto acylation.
Specificallyinvestigatedeffects of two mutants, Y170F andpK
a values of alcohol molecules and theoretical GFE with theunderstand its catalytic property. It is also meaningful to
D181G, theindicatedin former
experiments
as crucial polarizable
systems continuum
because oftheir
specific
model
(PCM)
on theamino
basis ofacid
ourinvestigate
redox potentialof signiﬁcantmetallo-proteins,
18−21previous
work
for
pK
values.
We
havecalculatedanaresidues.
Therefore,byadopting
first-principles
molecular
dynamics
complementedwithwhich play a crucial role in the electron transfer process inabsolute SHE potential of 4.52 V when we employed a goldbiomolecules.
For
example,acytochrome
c,
which
is
a
heme
metadynamics we provide a detailed insight into the underlying
acylation mechanism. Our results show
starndard calculation level as CCSD(T)/aug-cc-pVDZ. Our
compoundan ironatom,
is a mediaterof electronthat whileandincontains
the wild
type(WT)
the Tyr170
residue points
theNH
themolecules
proton-acceptor
previous study
for group
organic towards
and inorganic
consiting oftransport processes in several biologically important systems. In
site
of
an
artificial
amide
bond,
hence
ready
to
react,
in
the
Y170F
this
does
not
occur.
The
reasonistypicalelements
indicated
that
the
value
of
the
SHE
potential
order to understand the correlation between the structure andon thebycomputational
methods,
and isweunable
proposedtoa
ascribed
to the absenceof Tyr170 in
the mutant,is replaced
phenylalanine,
which
the
redox potentialof the metallo-enzymes,
it requires
accuratewhichdepends
necessity for tuning of the SHE potential. The mean absolutecomputation
schemes
fromwithﬁrst principles.
form hydrogenbond
the amide bond; thus,1−8resulting
in anerror of -1the redox potentials became 0.1 V when we used a
Although inmanyreportedoftheir
methods. potential by using a popular density functional
increase
thetheoreticians
activationhave
barriermorethan 10 kcal
mol
tuned SHE
to compute the redox potential, the error from the experimental
Nonetheless,despitethelackofhydrogenbonding
betweenthetheory (DFT) such as B3LYP.
value is quite large, such as 0.20−0.30 V in mean absolute error
Y170FThe
anddeviation
the substrate,
highestvalues
free energy
forspitethe of the success with the typical molecules, our(MAE).
from the theexperimental
originatesbarrier In
previous
scheme cannot describe TMCs with an excess charge.
mainlyfrom theistreatmenthydrogenelectrodeto suggest thatinduced-fit
similar of
to thethatstandard
of WT.
This seems
Thereremaintwo high barriers to achieve the accurate(SHE)
which process,has been treatedas a isconstantin thepotential,
induced-fit
kinetics
little andaffected by the
independent of the computational method in many theoreticalmutation.
On the basis
structural homology
papers.
The international
union ofof
pureadditional
and applied chemistry
Received: April 2, 2013
analysesrecommendson the enzymes
theV samefamily,we suggestthat May 23, 2013
(IUPAC)the valueof4.44
as a SHE
potential
Published:natural selection is responsible for the development of the© 2013 American Chemical Society
2974
dx.doi.org/10.1021/ct4002653 | J. Chem. Theory Comput. 2013, 9, 2974−2980
peculiar hydrolytic activityof Arthrobacter sp. KI72.

133

De novo design of the functional transition metal proteinsusing a coiled-coil scaffold
Yu Takano
Institute for Protein Research, Osaka University
Introduction
De novo design of proteins is becoming increasingly attractive, because it can create novel
materials that adopt a designed function at will. Additionally, it enables us to understand what isrequired for the protein function. In particular, the creation of a metal binding site bearing
unusual characteristic features within a scaffold protein would improve our understanding of thecorrelation between the coordination structures and the spectral features and facilitate design
strategies for metal binding sites. We performed quantum mechanical and molecular mechanical(QM/MM) computations for metalloprotein and found that protein environment enhances theintrinsic character of metal active sites and provides function of metalloproteins. It indicates that
it is required to investigate the electronic structures of the metal active site itself and the effectof protein environments on the active site. Our ultimate goal is to propose a guideline ofcomputational design of de novo metalloproteins as a protocol of a priori design of novel
molecules by using computics. In this project, we performed QM and QM/MM calculations on
the CuA center of cytochrome c oxidase (CcO) and nitrous oxide reductase (N2OR) to elucidate
the quantitative structure-function relationship in copper protein. Based on the computational
results, we conducted de novo designs of copper proteins and successfully created the purple
copper CuA site within a de novo scaffold coiled-coil protein. We also have investigatedmolecular interactions in the α-helices and β-sheets composed of alanine residues, using
quantum chemical (QC) methods (B97D/6-31+G(d)) and molecular mechanics (MM)(AMBER99-SB) to improve the force field for accurate description of protein behavior, whichis required for the reliable computation of de novo metalloproteins.
Computational studies of the electronic structure of the CuA site in cytochrome c oxidase
and nitrous oxide reductase
CcO and N2OR belong to purple copper proteins and possess electron transfer functions. Bothproteins receive the electron from cytochrome c at an unusual copper electron transfer center
called CuA. The CuA site contains a binuclear copper center, which is composed of two bridging
Cys residues forming a Cu2S2 core. The CuA site has two additional His residues interacting
equatorially and several weakly interacting axial groups, such as Met, His, and the main chain
carbonyl of Glu. These sites have been characterized as the binuclear copper center with mixed
valence [Cu(1.5)−Cu(1.5)] species in various spectroscopic studies. The ground state of theoxidized CuA site is strongly related to the singly occupied orbital, since the oxidized state of
the CuA site is doublet. Spectroscopic studies have suggested the σu* ground state of the CuA

134

site due to the direct CuCu interaction, facilitating a rapid electron transfer over long distanceswith low driving forces.
First of all, in order to elucidate what
factors are responsible for the stabilization ofthe σu* ground state of the CuA site, theelectronic structures of the Cu2S2 core of the
CuA site in CcO have been studied, using thedensity functional theory (M06) to elucidate
the origin of the σu* ground state of the CuA
site. Our computational results show that the
ground state of the Cu2S2 core is the πu stateeven with a shorter Cu–Cu distance, and thatan addition of electrostatic and orbital
interactions by ligand coordination stabilizes
Figure 1. Orbital diagram of the Cu2S2 core

the σu* state rather than the πu state, because
electrostatic repulsion and antibonding

orbital interaction between the dσ* orbital of the Cu–Cu dimer and the lone pairs of thecoordinating His residues rise the orbital energy of the dσu* orbital of the CuA site, as shown in
Figure 1.
Next, we have studied the coordinating effects of His, Met, and the peptide carbonyl group of
Glu of the CuA site, using the density functional theory. We constructed models to examine thecoordination effect. The models of the CuA site are shown in Figure 2. In the models, weconsidered the one-ligand removed models of the CuA site, in which each coordinating residue
was removed (the non-H161, non-E198, non-H204, and non-M207 models). The one-ligand
removed model nullifies the orbital and electrostatic effects of the removed ligand. To
investigate the electrostatic effects of the removed ligand, we also built point charge models (thepc-H161, pc-E198, pc-H204, and pc-M207 models), in which we replaced each coordinating
ligand with point charges, as illustrated in Figure 2.

Figure 2. Models of the CuA site.
Comparison between the one-ligand removed and the point charge models suggests that theorbital interactions between the Cu2S2 core and the coordinating ligands dominate   of135

the CuA site. Our computation has shown that the electrostatic effect, due to the ligation of
His161, His204, and Glu198, is more influential on the IP than the orbital effect (Figure 3).
With regard to the coordination of Met207, however, the orbital and electrostatic effects
equivalently affect the IP. These results indicate that the strong orbital and electrostatic effectsof the His198 and His204 residues both contribute to the stabilization of the σu* ground state
and the decrease in the ionization potential, as compared to the Cu2S2 core itself. The peptidegroup of Glu198 weakly coordinates to the Cu2S2 core of the CuA site and only electrostatically
influences the ionization potential of the CuA site. The coordination of Met207 provides weakorbital and electrostatic effects, which are useful for the fine-tuning of the electronic structure of
the CuA site.

Figure 3. Energy gaps between the σu* and πu states (  ) and ionization potentials (IPs)of the models of the CuA site. Black bold lines represent the core and full models. Bluedashed-dotted and red dotted lines represent the one-ligand (non-) and the point charge (pc-)
models, respectively. Black thin line is drawn on those of the full model as a reference.
We also found that three dimensional X-ray crystallographic structures of the CuA site showthat the structural variety of the Cu2S2 core of the CuA site with a short Cu–Cu distance of 2.34–
2.59 Å and a S–S
distance of 3.76–4.24
Å.

We addressed theorigin of the structuralvariety of the Cu2S2
core of the CuA site,
using

the

densityfunctional theory. Our
computation
demonstrated that we
found

the

dependence

structuraloftheenergy gap between the

Figure 4. Scatter plots of the energy gap between the σu* oxidized
and reduced states () as functions of the Cu–Cu bond length(A), the S–S bond length (B). Black lines represent the regression
lines, and R denotes the Pearson's correlations.

136

σu* and πu oxidized states (Figure 4), the ionization potentials, and total energies, though thestructural variety is insensitive to the shape and symmetry of the σu* and πu redox active
molecular orbitals (RAMOs) and the πu ground state of the Cu2S2 core of the CuA site. In
addition, the σu* oxidized state is more robust to the structural changes in ionization potentials
and total energies than the πu oxidized state (Figure 5). This robustness indicates the small
reorganization energy in the redox reaction between σu* oxidized and reduced states of the CuA
site. The optimized Cu2S2 cores have almost same structures in all the states, indicating thatsurrounding proteins are responsible for the structural variety of the Cu2S2 core of the CuA site.

Figure 5. Schematic pictures of the Cu–Cu indirect interaction via the Cu–S bonds for IP u (A)and the Cu–Cu direct interaction for IP * (B).
uCreation of a binuclear purple copper site within a de novo four-helical coiled-coil protein
The creation of a metal binding site bearing unusual characteristic features within a scaffold
protein would improve our understanding of the correlation between the coordination structures
and the spectral features and facilitate design strategies for metal binding sites. On the basis of
the blue copper site described previously, we successfully created the purple copper CuA site
within a de novo scaffold tetrameric α-helical coiled-coil protein, amino acid sequence of whichis shown in Figure 6. This protein is referred to as bi-AM2C.

Figure 6. Amino acid sequences of bi-AM2C (X = H, Y = C) for the purple copper protein.
The protein bound two copper ions and exhibited a purple color, with relatively intense
absorption bands at 488 and 530 nm in the UV−vis spectrum (Figures 7A and B). The EPR
spectrum displayed unresolved hyperfine splittings in the g region, which was similar to thenative or engineered Cus site with an A∼480/A∼530 > 1. The extended X-ray absorption structure
analyses of the protein revealed the presence of the Cu2S2 core structure, with two typicalN(His)−Cu bonds per core at 1.90 Å, two S(Cys)−Cu bonds at 2.21 Å, and the Cu−Cu bond at137

2.51 Å, which are also characteristic structures of a purple copper site. We also investigated thecoordination geometry of the active site in the model of the bi-AM2C-(Cu)2 complex. Since nocrystal structure of bi-AM2C has been reported, the model of the bi-AM2C-(Cu)2 complex wasconstructed and minimized by the conjugate gradient method (Figure 7C). The active site of thebi-AM2C-(Cu)2 complex is located in the vicinity of the second α-helix of the coiled coil. The
active site contains copper ions bridged with the deprotonated cysteinyl groups of Cys47 andCys50. The copper ions and cysteinyl sulfur atoms form an almost planar Cu2S2 core with adihedral angle of 163.2°. The dihedral angle of the Cu2S2 core of the bi-AM2C-(Cu)2 complex isclose to those of the nitrous oxide reductase from Paracoccus denitrificans (164.0°− 165.9° ).
The copper ions are strongly and equatorially coordinated with two His residues (His13 andHis81). The peptide carbonyl groups of Glu46 and Cys47 reside within 2.7 Å from the copper
ion and might axially interact with it. His16 and His84 are located over the Cu2S2 core, atdistances of 4.07 and 2.70 Å from the copper ions, respectively, implying that these His
residues might also form axial interactions with the copper ions.

Figure 7. Purple color of the solution of bi-AM2C and Cu2+ ion (A), UV-vis spectrum of bi-AM2C and Cu2+ ion (B), and computational model of bi-AM2C-(Cu)2 complex (C).

Theoretical study of the molecular interactions in secondary structures of proteins
The three-dimensional structure of a protein
determines

its

functionsandchemicalproperties. The secondary structures such as
α-helicesandβ-sheets

areimportantcomponents for the protein architecture, as
shown in Figure 8. Molecular interactions, inparticular hydrogen bonding, play significant
roles in the formation of the secondary
structures. Quantitative estimate of theseinteractions is required to understand theprincipleoftheformationofthethree-dimensional protein structure. Here, toFigure 8. Three-dimensional structure andsecondary structure of proteins

improve the force field for accurate description of protein behavior, we have investigated138

molecular interactions in the α-helices and β-sheets composed of alanine residues, using
quantum chemical (QC) methods (B97D/6-31+G(d)) and molecular mechanics (MM)(AMBER99-SB) . The characteristic interactions essential for forming the secondary structures
are discussed quantitatively.
The MM force field stabilizes the α-helix much more than the QC method. In addition, theMM force field overestimates the interaction energies of the α-helix, compared to those of theβ-sheet, implying that the improvement of the electrostatic interaction is required (Figure 9).
We also have evaluated each hydrogen bond (H-bond) energy in the α-helices and theβ-sheets by using the molecular tailoring approach (MTA), which is one of the fragment-based
methods. Within MTA, a large
system is divided into somesubsystems.

Theelectronic

properties of the subsystems arecomputationally evaluated. The
properties of the large system are
obtainedbypatching

thecorresponding properties of all thesubsystems. MTA shows that the
H-bonds
strongerin the β-sheets
than those in

are
the

α-helices due to the structural
constraint in α-helix and that the
antiparallel β-sheet can form
linear H-bonds rather than theparallel β-sheet, resulting in thestronger

H-bonds

intheantiparallel β-sheet. In addition,

Figure 9. Relationship between hydrogen-bond energies(interaction energies or relative energies) and the number ofhydrogen bonds

the H-bonds in the α-helices canbeclassified

threetypes:

N-terminal, C-terminal, and inner H-bonds, and inner H-bonds are weaker than other H-bonds
at the same H-bond distances. The differences of the H-bond energies between the parallel andantiparallel β-sheets and between the N-terminal, C-terminal, and inner H-bonds in the
α-helices indicate secondary-structure dependence on the H-bond energies.
Members of this research project
Co-Investigator (renkei-ken kyusha): Prof. Toshiki Tanaka (Nagoya Institute of Technology)
Research Collaborator: Prof. Yasuteru Shigeta (Tsukuba University)
Research Collaborator: Ms. Makoto Kita (Osaka University)
Research Collaborator: Ms. Orio Okuyama (Osaka University)
Research Collaborator: Mr. Tomoyuki Iwamoto (Osaka University)

139Publications and Invited lectures
Papers1)

K. Koizumi, Y. Shigeta, O. Okuyama, H. Nakamura, Y. Takano, "Coordination effects on theelectronic structure of the CuA site of cytochrome c oxidase", Chem. Phys. Lett. 531,
197-201 (2012).

2)

Y. Takano, Y. Shigeta, K. Koizumi, H. Nakamura, "Electronic structures of the Cu2S2 core of
the CuA site in cytochrome c oxidase and nitrous oxide reductase", Int. J. Quantum Chem.
112, 208-218 (2012).

3)

D. Shiga, Y. Funahashi, H. Masuda, A. Kikuchi, M. Noda, S. Uchiyama, K. Fukui, K.
Kanaori, K. Tajima, Y. Takano, H. Nakamura, M. Kamei, T. Tanaka, " Creation of a binuclear
purple copper site within a de novo coiled-coil protein", Biochemistry 51, 7901-7907 (2012).

4)

Y. Takano, O. Okuyama, H. Nakamura, Y. Takano, "Density functional studies of thestructural variety of the Cu2S2 core of the CuA site", Int. J. Quantum Chem. 112, 3583-3591(2012).

5)

Y. Takano, K. Yamaguchi, H. Nakamura, "Theoretical studies of effect of protein
environment on electronic structures and chemical indices of the active site of oxygenatedand deoxygenated hemerythrin", Int. J. Quantum Chem. 113, 497-503 (2013).

6)

Y. Takano, O. Okuyama, Y. Shigeta, H. Nakamura, "Density functional study of the origin of
the strongly delocalized electronic structure of the CuA site in cytochrome c oxidase",
Quantum Systems in Chemistry and Physics -Progress in Methods and Applications-,
Progress in Theoretical Chemistry and Physics Vol. 26, edited by Kiyoshi Nishikawa, Jean
Maruani, Erkki J. Brändas, Gerardo Delgado-Barrio, and Piotr Piecuch, Chapter 30, Springer
2013, pp. 513-524.
Y. Takano, M. Kita, H. Nakamura, "Quantum chemical study of conformational preferences
in intermediates and transition states of the hydrolysis of dimethyl phosphate.", J. Phys.
Conf. Ser. 454, 208-218 (2013).
M. Kita, H. Nakamura, Y. Takano, "Density functional study of the phosphodiesterhydrolysis of RNA in RNA/DNA hybrid by RNase H.", Mol. Phys. 112, 355-364 (2014).

7)

8)Invited lecture
1)

Yu Takano, "How to use PCM (Polarizable Continuum Model) for description the molecularenvironment, and application of the PCM to the protein environment of metalloprotein.", The
13rd Dai Tsukuba Physical Chemistry Seminar (December 7, 2012, Josai University)

2)

Yu Takano, "Towards de novo design of functional metalloproteins by using computer
simulation", Mie University CUTE seminar (January 15, 2013, Mie University)

3)

Yu Takano, "Relationship between molecular structure, electronic structure, and function of
metal cofactors in proteins", The 6th symposium (November, 20, 2013, Kobe University)

4)

Yu Takano, "Computational study of relationship between molecular structure, electronic
structure, and function of metal cofactors in proteins", IPR seminar (February 8, 2014, Osaka
University)

5)

Yu Takano, "Theoretical study of electronic structures of metal cofactors in proteins", The
3rd computational life science WG lecture (March 17, 2014, Tsukuba University)

6)

Yu Takano, "Electronic structures of metal cofactors in proteins -QM and QM/MM
analysis-" Kinki Chemical Society Japan lecture (May 9, 2014, RIKEN Advanced Institutefor Computational Science)

7)

Yu Takano, "Theoretical study of the molecular interactions to form secondary structure inproteins", The 14th Protein Science Society Japan annual meeting (June 26, 2014, Work Pier
Yokohama)

140

Construction of Material Design System on
Next-Generation Density Functional Theory
Yutaka Imamura
RIKEN Advanced Institute for Computational Science
1.

Introduction
Density functional theory (DFT) has been one of widely used methodologies because of

relatively accurate description and low computational cost comparable to that of Hartree-Fock(HF) method. However, the exact DFT total energy expression has not been known but has been
rigorously proposed. Hybrid functionals have become mainstream among a wide range of
exchange-correlation (XC) functionals; for instance, B3LYP, one of representative hybrid
functionals, has been applied to many systems such as molecules and nanomaterials not only by
quantum chemists but also by experimentalists. Recently, more sophisticated XC functionals suchas range-separated and local as well as orbital-specific (OS) hybrid functionals proposed by us
have been developed. The OS hybrid functionals determine portions of HF exchange interaction
by the linearity condition for orbital energies (LCOE) without resorting to numerical fitting toexperimental standard formation of enthalpies. The OS hybrid functionals have succeeded in
reproducing core and valence ionization potentials (IPs). However, the conventional scheme for
the OS hybrid functionals is not applicable to multiple LCOEs that become important for excitedstates and multiple orbital-involved chemical reactions. This study explored a new scheme todevelop new OS hybrid functionals that satisfy multiple LCOEs based on a multiply-range
separated functional [1].
2.

TheoryThis section explains how to build OS hybrid functionals that satisfy multiple LCOEs. Janak’s

theorem in Kohn-Sham DFT (KS-DFT) states that the first derivative of the total energy withrespect to the occupation number fi of the ith KS orbital is equivalent to its orbital energy:

E i .
f i(1)

Almbladh and von Barth proved the relation for the highest occupied molecular orbital (HOMO):

E
f HOMO  HOMO  IP ,(2)which indicates that the HOMO energy is equivalent to the first IP with an opposite sign. Since
the IP must be constant, the following LCOE is derived:

141


2E
 i
0.
2
f i 0 f 1 f i 0 f 1
i

(3)

i

Although Eq. (3) is exact for the HOMO, it may not necessarily hold for other orbitals; however,Eq. (3) provides a better description through the removal of the self-interaction errors.
The multiply range-separated hybrid XC functional is assumed to have the following form in
the case of N LCOEs:
N 1





SR
SR
MRi
MRi
LR
Exc  (1  SR ) EDFTx SR EHFx  (1   MRi ) EDFTx
  MRi EHFx
EHFx
 EDFTc . (4)i 1

HF/DFT
HF/DFT
HF/DFT
EX,SR, EX,MRi , and EX,LR
are the short-, middle-, and long-range HF/DFT exchange

energies, respectively, with corresponding fractions of SR ,  MRi , and  LR .

When the number of parameters of SR ,  MRi , and  LR are equal to or less than that of

LCOEs, optimization of the parameters is basically executable. In the present study,  wasdetermined using the Brent’s method, which combines the bisection and secant methods, andinverse quadratic interpolation so as to exactly satisfy LCOEs. A one-to-one correspondence
between LCOE and  was also assumed.
3.

Results and Discussion
We numerically assessed the new OS hybrid functionals by computing four types of excitedstates for CO and N2 molecules: core-valence, core-Rydberg, valence-valence, and valenceRydberg transitions using time-dependent density functional theory (TDDFT). The OS hybrid
functional with two LCOEs of HOMO and the deeper 1s orbital was adopted. XC functionals
such as BLYP, B3LYP, and LC-BLYP were examined for comparison. Table lists the excitationenergies for conventional and OS hybrid XC functionals. Deviations from experimental values
and percent deviations are listed in parentheses and square brackets. Mean absolute deviation(MADs) and mean absolute percent deviations (MAPDs) of four excitation energies fromexperimental values and overall MADs and MAPDs are also given.
The valence-valence excitation energies of CO and N2 were reproduced within 1 eV using
BLYP and B3LYP. The MADs were (0.31 and 0.37) eV and (3.51% and 4.17%) for (BLYP,
B3LYP), respectively. As for valence-Rydberg transitions, excitation energies were
underestimated especially for BLYP. The MADs and MAPDs were (1.88 and 1.01) eV and(15.72% and 8.43%) for (BLYP, B3LYP), respectively. The long-range corrected hybrid
functional LC-BLYP, which was designed to improve behavior in the asymptotic region,
accurately estimated excitation energies for valence-valence and valence-Rydberg excitations
with MADs of less than 1.0 eV.

142

Table Excitation energies (eV) calculated using conventional and OS hybrid functionals.
BLYP
B3LYP
LC-BLYP
Valence→Valence
CO molecule* 35.87 (-0.45) [-7.15] 5.90 (-0.42) [-6.69]
6.00 (-0.32) [-5.11]6.30
σ→π ( Π)* 1
8.29 (-0.22) [-2.62] 8.45 (-0.06) [-0.67]
8.63 (0.12) [1.37]
8.98
σ→π ( Π)* 3 +8.14 (-0.37) [-4.38] 7.97 (-0.54) [-6.31]
8.11 (-0.40) [-4.75]
8.24
π→π ( Σ )* 3
8.73 (-0.63) [-6.70] 8.67 (-0.69) [-7.32]
8.84 (-0.52) [-5.54]8.90
π→π ( Δ)
* 1 9.84 (-0.04) [-0.41] 9.78 (-0.10) [-1.03]
9.79 (-0.09) [-0.95]9.61
π→π ( Σ )* 1
10.08 (-0.15) [-1.47] 10.10 (-0.13) [-1.31] 10.27 (0.04) [0.41] 10.48
π→π ( Δ)
N2 molecule3+7.63 (-0.12) [-1.55] 7.23 (-0.52) [-6.72]
7.09 (-0.66) [-8.55]
6.81
πu→πg ( Σu )3
7.56 (-0.48) [-6.03] 7.68 (-0.36) [-4.47]
7.97 (-0.07) [-0.81]
8.11
σg→πg ( Πg)3
8.37 (-0.51) [-5.70] 8.11 (-0.77) [-8.63]
8.04 (-0.84) [-9.44]
7.78
πu→πg ( Δu)1
9.21 (-0.10) [-1.11] 9.38 (0.07) [0.79]
9.65 (0.34) [3.63]
9.93
σg→πg ( Πg)
3 9.74 (0.07) [0.69] 9.49 (-0.18) [-1.91]
9.24 (-0.43) [-4.43]
8.81
πu→πg ( Σu )1 9.74 (-0.18) [-1.85] 9.49 (-0.43) [-4.38]
9.24 (-0.68) [-6.84]
8.81
πu→πg ( Σu )1
10.02 (-0.25) [-2.40] 9.88 (-0.39) [-3.81]
9.84 (-0.43) [-4.16]9.87
πu→πg ( Δu)3
10.40 (-0.79) [-7.09] 10.70 (-0.49) [-4.34] 10.99 (-0.20) [-1.78] 11.60
σu→πg ( Πu)Valence→Rydberg
CO molecule
1 +9.06 (-1.72) [-15.97] 9.90 (-0.88) [-8.20] 11.10 (0.32) [2.93] 11.17σ→3sσ ( Σ )1 +9.53 (-1.87) [-16.42] 10.31 (-1.09) [-9.59] 11.64 (0.24) [2.08] 11.70σ→3pσ ( Σ )1
9.67 (-1.86) [-16.10] 10.44 (-1.09) [-9.47] 11.74 (0.21) [1.79] 11.79
σ→3pπ ( Π)N2 molecule
3 +σ→3sσ ( Σg ) 10.24 (-1.76) [-14.68] 11.08 (-0.92) [-7.66] 12.06 (0.06) [0.48] 12.41
1 +σ→3sσ ( Σg ) 10.59 (-1.61) [-13.17] 11.49 (-0.71) [-5.82] 12.64 (0.44) [3.57] 13.08
1σ→3pπ ( Πu) 10.82 (-2.08) [-16.10] 11.80 (-1.10) [-8.52] 13.04 (0.14) [1.06] 13.60
1 +σ→3pσ ( Σu ) 10.69 (-2.29) [-17.61] 11.71 (-1.27) [-9.78] 13.11 (0.13) [1.03] 13.57
Core→
CO moleculeC1s→π*
271.36 (-16.04) [-5.58] 276.23 (-11.17) [-3.89] 271.30 (-16.10) [-5.60] 285.76
O1s→π*
512.60 (-21.60) [-4.04] 520.14 (-14.06) [-2.63] 513.16 (-21.04) [-3.94] 535.46
N2 moleculeN1s→π*
382.63 (-18.37) [-4.58] 388.72 (-12.28) [-3.06] 382.73 (-18.27) [-4.56] 403.63
Core→
CO moleculeC1s→3sσ272.81 (-19.59) [-6.70] 279.41 (-12.99) [-4.44] 275.86 (-16.54) [-5.66] 292.15
C1s→3pπ273.55 (-19.75) [-6.73] 280.12 (-13.18) [-4.49] 276.83 (-16.47) [-5.62] 293.16
C1s→3pσ273.40 (-20.10) [-6.85] 280.06 (-13.44) [-4.58] 276.76 (-16.74) [-5.70] 293.24
O1s→3sσ514.08 (-24.82) [-4.61] 522.97 (-15.93) [-2.96] 516.45 (-22.45) [-4.17] 540.49
O1s→3pπ514.76 (-25.14) [-4.66] 523.79 (-16.11) [-2.98] 517.64 (-22.26) [-4.12] 541.83
N2 moleculeN1s→3sσ384.17 (-22.03) [-5.42] 391.89 (-14.31) [-3.52] 386.76 (-19.44) [-4.79] 410.24
N1s→3pπ384.63 (-22.47) [-5.52] 392.49 (-14.61) [-3.59] 387.72 (-19.38) [-4.76] 411.40
N1s→3pσ384.53 (-22.77) [-5.59] 392.50 (-14.80) [-3.63] 387.99 (-19.31) [-4.74] 411.74
MAD/MAPD (V→V)(0.31) [3.51](0.37) [4.17](0.37) [4.13]
MAD/MAPD (V→R)(1.88) [15.72](1.01) [8.43](0.22) [1.85]
MAD/MAPD (C→V)(18.67) [4.74](12.50) [3.19](18.47) [4.70]
MAD/MAPD (C→R)(22.08) [5.76](14.42) [3.77](19.07) [4.94]
MAD/MAPD (Total)(7.82) [6.86](5.16) [4.91](6.71) [3.89]

OS(-0.02)(0.47)(-0.27)(-0.46)(-0.27)(0.25)

Exptl.
[-0.39] 6.32
[5.54] 8.51[-3.17] 8.51[-4.91] 9.36[-2.75] 9.88
[2.42] 10.23(-0.94) [-12.13] 7.75
(0.07) [0.92] 8.04(-1.10) [-12.34] 8.88
(0.62) [6.64] 9.31(-0.86) [-8.87] 9.67(-1.11) [-11.16] 9.92(-0.40) [-3.90] 10.27(0.41) [3.67] 11.19(0.39)(0.30)(0.26)[3.66] 10.78[2.64] 11.4
[2.24] 11.53(0.41)(0.88)(0.70)(0.59)[3.38]
12[7.25] 12.2
[5.44] 12.9
[4.54] 12.98(-1.64)(1.26)[-0.57] 287.4[0.24] 534.2(2.63)[0.65]

401(-0.25)(-0.14)(-0.26)(1.59)(1.93)[-0.09][-0.05][-0.09][0.29][0.36]

292.4
293.3
293.5
538.9
539.9

(4.04)(4.30)(4.44)(0.52)(0.51)(1.84)(2.12)(1.04)[0.99] 406.2[1.06] 407.1[1.09] 407.3[5.63][4.17][0.49][0.50][3.55]

We turn to the core-valence and core-Rydberg transitions in CO and N2 molecules. The three
conventional functionals severely underestimated excitation energies. The MADs and MAPDs
increased in the order of B3LYP < LC-BLYP < BLYP. The values were 12.50, 18.47, and 18.67
eV for core-valence transitions, and 14.42, 19.07, and 22.08 eV for core-Rydberg transitions. For
core-valence transitions, short-range HFx is the most important; thus, the 20% HFx in B3LYP in
the core region reduced deviations. For core-Rydberg transitions, long-range HFx is also
important for the description of Rydberg orbitals. Thus, a performance difference was observed
between BLYP and LC-BLYP for core-Rydberg transitions.
Finally, we investigated the performance of the OS hybrid functional. Excitation energies are

143well reproduced for the four types of transitions, i.e., valence-valence, valence-Rydberg, corevalence, and core-Rydberg. As for core excitations, MADs and MAPDs were drastically reduced
from those of the three conventional functionals; in particular, core excitation energies of CO
were remarkably improved. As for valence excitations, reasonable MADs were obtained. The
overall MAD, 1.04 eV, and MAPD, 3.55%, were smaller than those of BLYP, B3LYP, and LCBLYP. Only the OS hybrid functionals offered reliable descriptions of all types of excitations.
4. Concluding Remarks
We proposed and numerically assessed a new scheme for building an OS hybrid functional thatsatisfies multiple LCOEs. In terms of excitation energies, the OS hybrid functional has ableddescription of three different types of orbitals, i.e., core, valence, and Rydberg, for single
calculations. This opens applications toward chemical phenomena involving not only core, butalso valence electrons.
5.

Publication list

1) Yutaka Imamura, Kensei Suzuki, Takeshi Iizuka, and Hiromi Nakai, Chem. Phys. Lett. 618,
30 (2015).
6.

Seminar1) Yutaka Imamura, "Development of Orbital-Specific Density Functional Theory", Peralta
Group Seminar at Central Michigan University, Mount Pleasant, USA, Jan. 16, 2014
2) Yutaka Imamura, "Development of Orbital-Specific Density Functional Theory", Perdew
Group Seminar at Temple University, Philadelphia, USA, Jan. 13, 2014.

144

First Principles Simulation of Nonequilibrim Transport at Nano-contact
and Design of Functional Materials
Hisao Nakamura
Nano System Research Institute, “RICS”, AIST

１ ． 	
 

Introduction

In the present research report, we summarize our recent studies of nonequilibrium transport
theory and first principles simulation of “nan-contact”. Here, we use the terminology“nan-contact” as an extended definition, which includes single/SAM molecular junctions, hetero
junction of the different semiconductor materials connected to electrodes etc. In particular, we
focus on developing a computational theory of electric/thermal current as well as interaction
between conducting electrons and ions, and we evaluate IV characteristics, local (vibrational)
heating by inelastic scatterings, current-induced conformation change, and thermoelectric
properties of nano-contact systems quantitatively. Then we analyze detail mechanism of thesenonequiribrium processes in the functional devises and make guiding principles to design
functional nano-contact devises.
We adopt mainly the framework of nonequilibrium Green’s function method combined withdensity functional theory (NEGF-DFT) to develop theoretical scheme. Our concrete target ismolecular electronics to bridge molecular synthesis/ fabrication technologies and devise physics,
design of high ZT (figure of merit) thermoelectric material, organic-inorganic hybrid
thermoelectric materials, and functional materials of nanoelecronics like nonvolatile memories.
Our aim is to establish computational electronics applicable to sub 10 nm scale device based on
quantum transport via first principles calculation.
２ ． 	
 Outline of our theoretical and computational approaches

On the transport calculation at nano-contact, electronic states in the semi-infinite electrodes
and in the contact (bridge molecule/resistance layer) have to be calculated on the equal footing.
Since the Fermi level in each electrode is shifted due to applied bias, calculation of Hartree
potential is obtained by Poisson equation with different boundary condition form that of an
equilribrium system. As stated in Sec. 1, we adopt NEGF-DFT framework. The mosttime-consuming step in the standard algorithm is self-consistent field (SCF) calculation ofelectron density (and other self-energy terms like electron-phonon interactions) since it requires
caluclation of many Green’s function matrices and accurate numerical intergration. In order toreduce heavy computational facility, we introduce “restricted SCF space” defined by thereference equilibrium state (i.e., zero-bias electronic state) and correct difference of electron
density between equilibrium and nonequilibrium conditions by updating only electron
occupation in the molecular orbitals (MOs) in the restricted SCF space. Note that this algorithm
includes the entire orbital relaxation effect by SCF. When the SCF result converges to an
unphysical state, the restricted SCF space violates smoothness (continuous change) by
iterating the SCF cycle. Hence, we can check unphysical result and/or validity of theapproximation and model directly.

145

The other important point in the present study is developing methods to simulate inelastictransport process by electron-phonon interaction and phonon thermal transport as well as
current-induced vibrational excitation. For calculation of inelastic electric current, we use therigorous formulation of the lowest order Born expansion (LOE) then introduce the following
approximations: (1) Energy dependence of the retarded (advanced) Green’s function is
calculated on the coarse-grained energy grids, and averaged over the bias window (2) subpart ofthe lesser (greater) Green’s function is integrated analytically as a function of energy and derive
analytical expression. By combining (1) and (2), computational task to evaluate electron-phonon
self-energy is much less than standard approach, which includes numerical convolution integral.
Furthermore, this conventional LOE does not require the wide-band-limit approximation. We
will apply our algorithm to study inelastic electron tunneling spectroscopy (IETS) of molecularjunctions, temperature dependence of resistance in resistance random access memory (ReRAM)
etc. In the present algorithm of LOE calculations, each phonon Green’s function is defined bynormal mode coordinate. Thermal dissipation of each intrinsic phonon mode to electrodes is
simply represented by related damping parameter. However, thermal conductance obtained bythis simplification is not as accurate as electric conductance by full NEGF-DFT. This is notsuitable to calculate ZT. Hence, separately, we implemented the program to calculate phonon
Green’s function in real space to apply thermoelectric materials. Real space phonon Green’s
function includes phonons of semi-infinite electrodes as the same with electric Green’s function,i.e., electric and thermal transport properties cab be calculated by the same level theory.

３ ． 	
 

Research Activities

3-1. Single molecular electronics by first principles transport calculation
In the ballistic electric transport at nano-contact, the transport properties such as IV
characteristics are determined by both of the intrinsic properties of the conductor (resistancelayer) and interaction between the conductor and electrodes. The former relates to “bandengineering” or “MO engineering”, and the latter is called “contact engineering” or “contact
chemistry”. A single molecular junction is the most fundamental devise and suitable to comparea theory and precise measurement. In addition, analysis of transport properties in the views of
MO engineering and contact chemistry is quite useful to design molecular device. We develop
effective Hamiltonian theory constructed by the molecular projected state Hamiltonian (MPSH)
and Green’s function by NEGF-DFT results. The complex energy of the effective Hamiltonian
spanned by each MO tells us the shifted MO energy and electronic coupling strength of
conducting MO and electrodes.
To realize single molecular electronics, mechanical stability of the junction, i.e., robustness ofanchoring of molecule to the electrode surface is required. It is usually considered thatrobustness of chemical bonds by anchor moieties relates
electric conductance positively.

Figure 1

Table 1146

However, recently, Tao’s group measured electric conductance of several molecular junctions
consisting of Au-C anchoring systematically. They found that the above relation between
mechanical robustness and conductance is sometimes violated, and not a reliable guiding
principles to synthesize molecular devise. We performed calculation of effective Hamiltonian
followed by NEGF-DFT. By using the results, we identified the conducting MO and electric
coupling strength with states in electrodes. Then we compared the coupling strength, adsorptionenergy and conductance. We found that, only when the conducting MO is the reactive MO andwhen Bloch states in the electrodes interact with reactive MO, the positive correlation statedabove is satisfied. Furthermore, the above matching is not so general. (Table 1)
3-2. Novel thermoelectric property of organometallic complex molecular materials
Efficiency of thermoelectric energy conversion is estimated by figure of merit, ZT. In order toincrease ZT, we can consider the following strategies (1) increasing Seebeck coefficient Se, (2)increasing electric conductance, (3) suppression of thermal conductance. However,simultaneous enhancement of (1) and (2) is generally difficult. Additionally, Wiedemann-Franz(WF) law tells us that enhancing electric conductance also increase electric thermal conductance.
Therefore, we have to consider validity of WF law, and possibly, design material (or devise
structure) violating WF law. When quantum dot is adopted as a device, one can expect giant
Seebeck coefficient by quantum confinement effect, interference effect of conducting MOs, andso on. In this sense, use of molecular junctions seems to be promising. However, electrictransport is sometimes tunneling in molecular conductors, in particularly, when quantum
interference effect is sufficiently large. It is well known that conductance reduces exponentiallyas a function of conductor length for tunneling transport. Since sufficient length of devise isrequired to make well-defined tangent of temperature in the devise, ZT (more strictly, power
factor) rapidly reduced in a practical devise even when thermal conductance is suppressed.
Recently, we found that long-range coherent tunneling is possible in organometallic
molecular wire theoretically. Thus, we focus on this long-range transport property and examine
design of thermoelectric devise via organometallic wire. For this purpose, we adopted
terpyridine complex of Ru connected to the electrode by S tom (Ru-complex) and analyzed
length dependence of ZT etc by calculations. Ru-complex has four nearly degenerated MOs
close to Fermi level. When we lengthen the wire (or stacking the layer) as the unit of
Ru-complex, the electric transmission coefficient spectra forms narrow overlapping resonancestructure. Since Se is determined by slope of transmission spectra as a function of energy,
narrow overlapping resonance enhances Se. Note that overlapping resonance also suggests
possibility of long-range coherent transport in thelow~middle temperature regime (~ 300 K). The above
enhancement mechanism of Se can be understand as theresult of quantum interference between four
quasi-degenerate tunneling orbitals. However, this isdifferent from giant Se by Fano resonance, which was
obtained by introducing side-chain group to a main
body organic conductor. In the view of a practicalmaterial design, the present approach has advantage
since chemical modification of the conductor is not
Figure 2
necessary and quantum interference effect is morerobust against molecular conformation change thanFano resonance.

147

Analysis of temperature dependence of Se is interesting because sign of Se relates to carrier.
As shown in Fig. 2, the sign changes by varying temperature, i.e., carrier (electron or hole)changes by varying temperature. This temperature dependence of sign of Se is characteristic
phenomenon of overlapping resonance of electric transport. On the contrary, the calculated ZT
is still much lower than the level expected as efficient thermoelectric devise in the low~middle
temperature regime (Figure 3) because electric conductance is relatively smaller than phonon
thermal conductance. To develop more practical organometallic molecular device, increasing
electric conductance (carrier concentration) by designing anchor moieties (contact chemistry)
will be promising with the strategy (3).

Figure 33-3. Inelastic transport and current-induced force at nano-contact
In sections 3-1 and 3-2, we showed studies of ballistic transport processes, where
electron-phonon interactions are omitted. In this section, we introduce our recent study of
inelastic electric current and current-induced ion migration (vibrational excitation) byelectron-phonon interactions. For this purpose, we focus on a current-induced reaction of
Melamine adsorbed on Cu(001) surface because it consists of fundamental processes to trigger
ion (atom) migration of energy-exchange by current-induced force. The force is often
considered as sum of field-driven and current-driven mean-field forces, and vibrational heating.
Although the system is simple, it can be a template to understand mechanism ofnano-fabrication, SET/REST nonvolatile memories like ReRAM, and so on.
Three conformations of melamine on Cu(001) can be switched by STM, which are labeled as
GS, C1, and C2, respectively, as shown in Fig. 4. In order to trigger reaction from theenergetically lowest state (GS) to C1, the required threshold voltage is 2.6V. This thresholdvoltage is much higher than height of reaction barrier and it is known that the reverse reaction(C1 to GS) is forbidden. On the contrary, transition between C1 and C2 states is reversible
reaction and the required voltage is in the broad range as [-1.0,1.-1.0] volt. Furthermore, thethreshold voltage to trigger the reaction is verylow.
In the present study, we focus on the three
mysteries: (1) why C1 → GS reaction isforbidden ? (2) what kind of force is driving force
for C1/C2 switching ? (3) what is physical originof the observed power law of reaction (C1 →C2) rate as a function of current ? We examine
our theoretical model and first principles
Figure 4
calculations to attack these mysteries.
Usually, the nonequilrbrium force acting on ions is classified to the field-driven force, wind148

force (current-driven force), and vibrational heating. The former two forces are summarized as
current-induced force and defined as Ehrenfest force by mean-field theory. Note thatcurrent-induced force can be defined rigorously; however, it cannot be divided into filed-driven
and wind force terms generally. Vibrational heating is non mean-field type force, i.e., kinetic
energy transfer by inelastic scatterings between conducting electron and ions.
First, we calculated current-induced force by NEGF-DFT with applying bias voltage andsearch the reaction path by nudged-elastic-band (NEB) method. By integrating current-induced
force along the reaction path, we construct nonequilibrium potential energy surface (PES).
The resulting bias-dependence of barrier height of C1/C2 transition is lower than 0.1 V evenwhen the applied bias is larger than 1.0 V. We applied the same procedure to GS/C1 andobtained similar result, i.e., change of barrier height by applying bias is very small. Hence, themain mechanism is vibrational heating.
As the next step, we calculated inelastic electron tunneling (IET) spectra for GS, C1, and C2
states, respectively. We found that one of strong peak positions of IETS of GS relates to thepromoting mode of GS→C1 transition. On the contrary, the promoting mode of C1→GS
transition is IETS inactive. This is direct evidence that C1→GS is forbidden as STM-induced
reaction. For C1/C2, we found that both of promoting modes of forward/backward reactions are
IETS active. Hence, we conclude that single mode vibrational heating model is sufficient toanalyze the current-induced mechanism and introduce rate equation approach, which is firstly
introduced as theoretical model of DIET/DIMET processes.

Figure 6

The vibrational excitation rate is evaluated by
NEGF-DFT results and parameterized to fit themodel. In Fig. 5 the reaction rate is plotted as a function of electric current for several fixed bias,
and clear power low behavior i.e., R(I) ∝ IN is obtained. Then we evaluated the power-order,N, and plotted as a function of bias voltage in Figure 6. The resulting bias dependence agrees
very well with those of experiment. The tendency of N (V) is explained by simple picture as
follows. When bias is increased, kinetic energy per one electron also increase: thus vibrational
excitation required for the reaction can be achieved by the scattering events with small numberof electrons. However, it is not straightforward that N relates to the average number of electrons
per reaction event because N is non-integer in most range and continuous function of V. Thisimplies that vibrational heating is quantum phenomenon not classical.
Figure 5149

3-4. Design of nonvolatile memory cell by first principles calculations
Finally, we show application of first principles transport calculations to design ReRAM cellconsisting of metal oxides. In the ReRAM device, the most simple cell structure consists ofmetal oxide (resistance layer) and metallic electrodes. The switching of SET/REST is controlledcontrolling by concentration of oxygen vacancy (Vo) in the resistance layer. To realize low
power cost for resistive switch, the device scale is now in sub 10 nm, which is close tomolecular electronics; thus quantum transport approach is applicable. Recently, someexperimental groups reported that that insertion of metal buffer layer improves SET/RESET(ON/OFF) ratio. Since distinct SET/RESET states are also key to power-saving read/write
processes as the memory, material design of buffer layer is highly desired. We adopted TiNelectrodes and HfO2 resistance layer. As buffer layer, we examined Ta, Hf, and Ti as well as nobuffer cell. We present schematic figure in Figure 7. When Vo is created in HfO2 layer, oxygenions (Ox) moved to the interface, i.e., the interface is oxidized. We include this oxidized effectto make ON state structure. The calculated IV curve of ON state is very close to that of OFF
state, i.e., not good performance as ReRAM. According to the calculated density of states of
HfO2 system (without electrodes), ON structure of the present model is sufficiently metallic.
Contact resistance at TiN surface and moved Ox ions frustrate electric transport through HfO2
(with Vo) conductor. IV characteristics of ReRAM cell consisting buffer layer is shown in
Figure. 8, where the metal buffer consists of only two atomic mono layers. Drastic improvementof distinct SET/REST states is found for all cases while ON/OFF ratio depends on materials.
Although the expected role of metal buffer layer was oxygen scavenger (i.e., reducing resistanceof SET state by rapid and high concentration of Vo), the present theoretical results suggest thatmetal buffer layer reduce contact resistance caused by Ox migration to the interface and thisimproves ReRAM performance. Recall that concentration of Vo in HfO2 is fixed for the all cellsof SET state, and that inserted buffer is only two atomic mono layers. The theoretical prediction
will change guiding principles to design ReRAM cell, in particular, sub 10nm scale device.

Figure 7
Figure 8

４ ． 	
 Selected publication list

1)

H. Nakamura, Y. Asai, J. Hihath, C. Bruot, and NJ. Tao, “Switch of conducting orbital bybias-induced electronic contact asymmetry in bipyrimidinyl-biphenyl diblock molecule:Mechanism to achieve pn directional molecular diode”, J. Phys. Chem. C, 115, 19931-19938(2011)

2)

K. Terada, H. Nakamura, K. Kanaizuka, M. Haga, Y. Asai, and T. Ishida, “Long Range
Electron Transport of Ruthenium-Centered Multilayer Films via Stepping-Stone Mechanism”，
ACS Nano 6，pp.1988-1999 (2012)

150

3)

R. Matsushita, M. Horikawa, Y. Naitoh, H. Nakamura, and M. Kiguchi, “Conductance and
SERS Measurement of Benzenedithiol Molecules Bridging Between Au Electrodes”, J. Phys.
Chem. C. 117, pp. 1791-1795 (2013)

4)

T. Hines, I. Diez-Perez, H. Nakamura, T. Shimazaki, Y. Asai, and N. Tao, “Controllingformation of single-molecule junctions by electrochemical reduction of diazonium terminalgroups”, J. Am. Chem. Soc.135 pp3319-3322 (2013)

5)

M. Kiguchi, Y. Takahashi, S. Fujii, M. Takase, T. Narita, M. Iyoda, M. Horikawa, Y. Naitoh,and H. Nakamura, “Additive Electron Pathway and Nonadditive Molecular Conductance byUsing a Multipodal Bridging Compound”，J. Phys. Chem. C 118，pp.5275-5283 (2014)

6)

H. Nakamura, T. Ohto, T. Ishida, and Y. Asai, “Thermoelectric Efficiency of Organometallic
Complex Wires via Quantum Resonance Effect and Long-Range Electric Transport
Property”, J. Am. Chem. Soc. 135, pp16545-16552 (2013)

7)

T. Ohto, I. Runger, K. Yamashita, H. Nakamura, and S. Sanvito, “Ab initio theory forcurrent-induced molecular switching: Melamine on Cu(001)”, Phys. Rev. B. 87, pp.
205439-1-7 (2013)

8)

K. Nishio, T. Miyazaki, and H. Nakamura, “Universal medium-range order of amorphous
metal oxide”, Phys. Rev. Lett. 111 pp155502 (2013)9)

Single Molecular Resistive Switch Obtained via Sliding Multiple Anchoring Points andVarying Effective Wire Length”M. Kiguchi, T. Ohto, K. Sugiyasu, S. Nakajima, M.
Takeuchi, and H. Nakamura, J. Am. Chem. Soc. 136, 7327-7332 (2014)10) “Additive Electron Pathway and Nonadditive Molecular Conductance by Using Multipodal
Bridging Compound”, M. Kiguchi, Y. Takahashi, S. Fujii, M. Takase, T. Narita, M. Iyoda, M.
Horikawa, Y. Naitoh, H. Nakamura, J. Phys. Chem. C., 118, 5275-5283 (2014)11) Design of ReRAM Cell Structure by Metal Buffer and Contact Engineering via
First-Principles Transport Calculation”Transport Calculation”, H. Nakamura, T.
Miyazaki, K. Nishio, H. Shima, H. Akinaga, and Y. Asai，Proceedings of IWCE 2014
(IEEE) (Paris 2014)151

Reactions analysis of plasma with self-assembled monolayer
Masanori Shinohara
Graduate School of Engineering, Nagasaki University
Abstract
We investigated an interaction between octadecyltrichlorosilane self-assembled monolayer and hydrogenplasma, with in-situ infrared absorption spectroscopy in multiple internal reflection geometry. The
molecular frame of –CH2–CH2– is etched by hydrogen plasma with the formation of CH3. The etching was
automatically stopped with the polymerization.
1．
．Introduction
Plasma processing is an inevitable tool for material processing such as semiconductor etching, surface
coatings and so on. Plasma is very complicated state, because it is a quasi-stable state and generates a lot ofkinds of radicals, ions and non-ionized neutral molecules. Then, the plasma–induced reaction is also very
complicated. In these days, much attention has been paid to the organic materials, including organic
semiconductor devices and bio-applications. These organic materials have complicated molecular
structures, compared with inorganic materials such as silicon, gallium arsenide. In the near future, plasma
processing will be required to these organic–materials processing. The reaction should be cleared forrealizing high-precision processing technology. Hence, the reaction analysis is difficult, because organic
materials and plasma are both complicated. Then, octadecyltrichlorosilane self-assembled monolayer(OTS–SAM) is adopted as the model materials for organic materials. Hydrogen plasma is adopted in thisstudy because hydrogen plasma is interesting for basic science and applied technology. In hydrogen plasmahydrogen radicals and ions are generated and they play an important role for various process. We
investigate and interaction between OTS–SAM and hydrogen plasma with in-situ infrared absorption
spectroscopy in multiple internal reflection geometry (MIR–IRAS). MIR-IRAS is powerful tool because it
can be used in-situ and real-time monitoring during plasma processing.

2. Experiment

RF power
Matching supply
network (13.56 MHz)A stainless steel chamber, maintained atless than 6.7 × 10–4 Pa was equippedwith an RF (13.56 MHz) plasma source,

Sample holder
H2

Plasmagas–feeding system and in-situ
MIR-IRAS measurement system. An
interferometer (JASCO), a HgCdTe

MIR prism as
the substrate

RF plasma
Glass tube

Thermocouple

Vacuum Pump system
IR(MCT) detector were used for
MIR-IRAS. The IR path was purged
with CO2–free dry air to reduce noise in

Fig. 1

152

IR
Concave
MirrorInterferometer

Experimental setup

MCT
detectorthe IR spectra because of H2O. The MIR prisms used in this study were made from n-type P-doped Si(100)wafers polished on both sides and had resistivity of approximately 10 Ω cm. The size of the prism was 10mm × 40 mm × 0.5 mm. After OTS–SAM was deposited on the chemical oxide covered prism, by
dropping 2.5 mM of OTS (Tokyo Chemical Industry) in toluene, the MIR prism was installed into thechamber. The prim was exposed to hydrogen plasma, generated by applying 30W of RF power at thehydrogen pressure of 6.7 Pa. The substrate was not applied to additional negative bias during plasma
exposure. Then, hydrogen ions should be annihilated near the substrate while the hydrogen radicals would
react with the SAM on the prism. The changes in chemical states were in-situ observed with MIR–IRAS.

3. Results and discussions
Fig. 2 shows the IR absorption spectra of the OTS–SAM exposed to hydrogen plasma in the C-H stretchingvibration region, compared with the as–deposited OTS–SAM. Then, the vertical axis indicates thedifference absorbance. Every spectrum shows three valleys (peaks in the negative direction) at
approximately 2850, 2915,and 2960 cm–1, and small shoulder in the negative direction at approximately
2900 cm–1. The peaks at 2850 and 2915 cm–1 correspond with CH2 components, and the peaks at 2960 cm–1
and the shoulder at 2900 cm–1 correspond with CH3. These indicates that the OTS–SAM were etched byhydrogen plasma exposure. Each peak was decreased with plasma exposure time.
In order to investigate decreasing feature of the peaks, the evolution of each maximum peak intensity wasplotted as a function of the plasma exposure time in Fig. 3, in a comparison with the as–deposited
OTS–SAM. The vertical axis indicates the residual peak intensity: Iresidual (%), which is denoted as follows;
Iresidual (%) = (Io – I decrease) / Io

(1)

where Io expresses the peak intensity acquired
1.5 min.

peak intensity in spectra acquired after plasma

5.5 min.

exposure. Iresidual indicates the residual portion

59.5 min.

Differencer absrobance

before plasma exposure and Idecrease expresses the(%) of each components on the prism. Each plot
indicates the decreasing features stop at a long–exposure time region and that the peak intensity
due to CH2 decreased more progressively than that
due to CH3. The peak intensity due to CH2
decreased into approximately 65 %. It is suggested
that each monomer become stable after

3200

approximately 6 species of CH2 per monomer isetched, considering OTS–SAM has 17 CH2
components per a monomer. The CH3 componentof the monomer should be preferentially etchedbecause it is located at the top of the monomer.
However, the CH3 component remains during

Fig. 2

3000
2800
-1
Frequency (cm )2600

IR absorption spectra of OTS-SAM

exposed to hydrogen plasma at different exposure
time in the C-H stretching vibration region,
compared with the spectrum of the as–deposited
OTS–SAM.

153

100

plasma exposure, indicating that the CH3-1component is generated by the hydrogen plasma: I(2960cm ), CH390

CH3–CH2–CH2– … –CH2–Si– + H →CH4 + –CH2–CH2– … –CH2–Si– (CH3
etching)(2)

Iresuidal (%)

exposure as follows:
80
70
60

–CH2–CH2– … –CH2–Si– + H →-1: I(2850cm ), CH2-1: I(2915cm ), CH 2

50

CH3–CH2– … –CH2–Si– (CH3 generation) (3)

40
0
10
20
30
40
50
60
Hydrogen plasama exposure time (min.)

These indicate that the after CH3 is hydrogenatedinto CH4 by hydrogen radical generated inhydrogen plasma, indicating that CH3 is etched

Fig.3

Evolution of each peak intensity as a(the model (2)). At the same moment remaining

function of hydrogen plasma exposure time

CH2 is hydrogenated into CH3 (the model (3)).
The stop of the etching cannot be explained withthe above model, which indicate the etching proceeds continuously.
Fig. 4 shows the IR absorption spectra in the low vibration region, in a comparison with the spectrum of the
as–deposited OTS–SAM. The peaks are observed at approximately at 1700 cm–1. The peak could bepossibility that the peaks are attributed to C=Obonds. However, we observed much broad peaks
due to C=O in the spectra acquired at the
OTS-SAM exposed to oxygen plasma. The peakintensity increased with the exposure time,
indicating that the polymerization would beoccurred with the formation of C=C bonds
during plasma exposure. the =CH, alkene
components would be generated by hydrogen

Difference absorbance (arb. units)

attributed to C=C bonds, although it has

59.5 min.

5.5 min.
1.5 min.

abstraction from the CH2 alkyl components of1900

the monomer, as follows;

1800–CH2–CH2– … –CH2–Si– + H →=CH–CH2– … –CH2–Si– + H2

(4)

Two =CH would coalesce to form C=C bonds inthe monomer. The polymerization would prevent
the OTS–SAM from being etched any more.

Fig.4

IR

1700
1600
-1
Frequency (cm )1500

absorption spectra of OTS-SAM

exposed to hydrogen plasma at different exposure
time in low vibration region (C=C region),compared with the spectrum of the as–deposited
OTS–SAM.

Conclusion154

We investigate an interaction between OTS–SAM and hydrogen plasma in order to understand a reaction
between organic molecules and plasmas, with MIR–IRAS. The IR absorption data showed OTS–SAM wasetched by hydrogen plasma with hydrogenation. The etching would automatically stopped with thepolymerization of the monomer during plasma exposure.

Reviewed Journal Paper
1) M. Shinohara, N. Maruno, Y. Yoshida, Y. Taniguchi, K. Ito, Y. Nakano, Y. Matsuda, H. Fujiyama,"Infrared spectroscopic study on reactions of octadecyltrichlorosilane self-assembled monolayer with
Plasmas," Jpn. J. Appl. Phys. Vol. 54 No. 1 (2015) 01AB04-1 (6 pages).
2) M. Shinohara, K. Amano, N. Maruno, Y. Yoshida, Y. Matsuda and H. Fujiyama, "Infrared Absorption
Spectroscopic Study on Reaction between Self-Assembled Monolayers and Atmospheric-Pressure Plasma,"J. Spectroscopy Vol. 2015 (2015) 417024 (7 pages).
3) M. Shinohara, Y. Taniguchi, Y. Takaki, Y. Matsuda, H. Fujiyama, "Adsorption of hydrocarbon
components generated in deuterated benzene plasma studied by in situ real-time infrared absorption
spectroscopy," Jpn. J. Appl. Phys. Vol.53 No.9 (2014) 096202 (6 pages).
4) M. Shinohara, T. Kawakami, K. Hara, S. Yagi, Y. Matsuda, H. Fujiyama, "Reactions of Surface
Hydrogen on Amorphous Carbon Films with Hydrogen Plasma," Jpn. J. Appl. Phys. Vol. 53 (2014) 010204(6 pages).
5) M. Shinohara, A. Fuakae, K. Amano, Y. Yoshida, Y. Matsuda, H. Fujiyama, "Development of an in-situ
infrared spectroscopic measurement of plasma-induced reactions in ethanol," Jpn. J. Appl. Phys. Vol. 52(2013) 11NC06 (5 pages).

Invited Tutorial Paper
1) M. Shinohara, Y. Matsuda, H. Fujiyama, "Deposition Process and Its Surface Reactions of an
Amorphous Carbon Film," The Journal of The Institute of Electrical Engineers of Japan
Vol. 135 No. 3 (2015) pp.139-141.
2) M. Shinohara, J. Kim, "Recent Trends of Plasma Engineering: –Infrared Spectroscopy of Plasma Process
and Plasma Processing in Graphene Research–," IEEJ Transactions on Fundamentals and Materials Vol.
134 No. 1 (2014) pp. 41-46.

Invited Talk
1) M. Shinohara, S. Yagi, Y. Taniguchi, N. Maruno, Y. Yoshida, Y. Matsuda, H. Fujiyama, "Infrared
spectroscopy of reactions induced by plasma," Atomic and molecular data forum (2013.10), National
Institute for Fusion Science.

155