⇞ᩱ㟁ụ㟁ゎ㉁⭷ࡢ୙ᆒ୍ᵓ㐀୰࡟࠾ࡅࡿỈ⣲ࠊ㓟⣲ࡢ㍺㏦ᮾ኱㝔᪂㡿ᇦᒸᓮ 㐍Ọ஭ဴ㑻ᐩᓅྜྠ◊✲఍ 2021.12.9

Collaborators

Dr. Y. Andohat NIMS

Dr. N. Yoshii
at Shin-Etsu Chemical

Dr. K. Fujimotoat Nagoya Univ.

Dr. T. Nagai

Dr. Z. Tangࡑࡢ௚ࡢ୺࡞MODYLASඹྠ㛤Ⓨ⪅
ᆏୗ㐩ဢ
ᒣ⏣⠜ᚿ
Yang Chunwei
Jiachao Zhangᾆ㔝 ㄹᕷᕝ┾୍୕ୖ┈ᘯ
Zhu Zheඹྠ◊✲⪅Current group membersࢺࣚࢱ⮬ື㌴
ྜྷ᳃ ᫂㸦᪂₲኱㸧❧ᕝṇ඾㸦ᶓ὾ᕷ኱㸧
ᚿ㈡ᇶஅ㸦ཎᏊຊᶵᵓ㸧

Substance transportation in the inhomogeneous systems
Fuel cell electric currentcross leakToyota‫ ە‬transportation of hydrogen, oxygen,and proton࣭bulk polymer electrolyte࣭electrode interface‫ ە‬prediction of transportation coefficient‫ ە‬mechanism
࣭free energy topographical map
࣭local and global paths
1,000,000 atoms

100,000,000 atoms

100,000 atoms

polymer electrolytebulk

4-phases interface
platinum, carbon, ionomer, and water

real interface㻶㻲㻱㻌䝔䜽䝜䝸䝃䞊䝏ᰴᘧ఍♫ ⇞ᩱ㟁ụ⏝㻼㼠ゐ፹ศᕸ䛾㼁㻸㼂㻙㻿㻱㻹ほᐹ㻔㻌㼔㼠㼠㼜㼟㻦㻛㻛㼣㼣㼣㻚㼖㼒㼑㻙㼠㼑㼏㻚㼏㼛㻚㼖㼜㻛㼚㼍㼚㼛㻙㼏㻛㼚㼍㼚㼛㼋㼏㼍㼟㼑㻜㻝㼋㼁㻸㼂㻙㻿㻱㻹㻚㼔㼠㼙㼘㻌㻕

MODYLAS余୺僐ィ⟬ᡭἲ兟ᶵ⬟I. ỗ⏝ศᏊືຊᏛィ⟬兟儆兗儙兗儺兏
NVE傎NVT傎NPT(full flexible)兟᫬㛫ⓎᒎRESPA傎multiple time step兟ᣊ᮰ཎᏊ傎ศᏊ㞟ᅋ傎཯ᛂᗙᶆ SHAKE,RATTLE,ROLL兟ຊሙྛ✀ຊሙ傎優兠儺兏㛵ᩘII. 儓兠兑兗┦஫స⏝ FMM傎Ewald傎PMEIII. ⮬⏤儌儱兏儒兠 ⇕ຊᏛⓗ✚ศἲ傎⢏Ꮚᤄධἲ傎overlap distribution method
IV. ୙ᆒ୍⣔఩⨨౫Ꮡ僔ᣑᩓಀᩘ傎኱ᇦⓗ僐ᾐ㏱ಀᩘ
V. 儻兑儬兗⛣ືvehicleᶵᵓ傎Grottussᶵᵓ傎EVB兆儫兏
䞉FugakuVI. ㄏ㟁ⓗᛶ㉁ㄏ㟁⦆࿴ࠊㄏ㟁ᦆኻࠊ㉥እ儝儾儓儬兏䞉PC clusters
VII. ᮦᩱ⣔僞僔ᒎ㛤 㧗ศᏊ傎儏兠兀兗傎◚ቯ傎⇞↝Open source
http://www.modylas.org/

余බ㛤凚♫఍ᐇ⿦凛僑ྥ傷僌兟tutorial, user’s manual, developer’s manual兟免儈償兗儝僔᳨ウ࣭ᐩᓅㄢ㢟࡬ࡢᥦ౪
兟GROMACS஫᥮傎VMD㐃ᦠ࣭௻ᴗ࡬ࡢᥦ౪
兟償兏ศ๭傎୪ิᩘ僔⮬ືタᐃ࣭᝟ሗᇶ┙ࢭࣥࢱ࣮➼࡬ࡢ࢖ࣥࢫࢺ࣮ࣝ
兟100୓ཎᏊ僸㉸傮僱኱つᶍ⣔ィ⟬僑᭱㐺㸦RIST㸧

MODYLAS‫ە‬㧗㏿໬
࣭᪂つ㏻ಙ࢔ࣝࢦࣜࢬ࣒MTDἲࡢ㛤Ⓨ࡜ᐇ⿦Ᏻ⸨࣭᭷ຠSIMD㛗7.2㸦┤᥋ィ⟬㸧ࠊ4.5㸦FMM㸧
࣭FMMἲ࡟࠾ࡅࡿ₇⟬ᩘࡢ๐ῶ ྜྷ஭䚸ᆏୗsolid harmonicsᇶᗏࡢ᥇⏝බ㛤∧࡜ẚ㍑ࡋ࡚
ᴟᗙᶆ㍈ࡢᅇ㌿
2.4ಸࡢ㧗㏿໬ᑐ⛠ᛶࡢ฼⏝§₞໬ᘧࡢ฼⏝࣭┦஫స⏝┤᥋ィ⟬㒊ࡢ༢⢭ᗘ/ಸ⢭ᗘΰྜ₇⟬໬ᙇ࣭FMMィ⟬㒊ࡢ᏶඲༢⢭ᗘ₇⟬໬࣭୺せMPI㏻ಙࡢRDMA㏻ಙ࡟ࡼࡿ⨨᥮

PMMA 6,187,520 atoms(◊✲༠ຊ⪅㸸୰ᑿ,బ⸨ (⌮◊RCCS))࣭EVBィ⟬㒊ࡢ୪ิ໬⸨ᮏ

Dr. Y. Andoh䠄⌧ NIMS䠅‫ە‬⢭ᗘࡢྥୖ
࣭FMMィ⟬ࡢṇ๎໬࡟ࡼࡿ඲࢚ࢿࣝࢠ࣮ಖᏑ
ᮒ

Dr. N. Yoshii䠄⌧ಙ㉺໬Ꮫ䠅

Andoh et al., JCTC 9, 3201(2013)
Yoshii et al. JCC 39, 1192(2018)
Yoshii et al., JCC 41, 940(2020)
Andoh et al., JCC 42, 1073(2021)

Substance transportation in the inhomogeneous systems
Fuel cell electric current

cross leakToyota‫ ە‬transportation of hydrogen, oxygen,and proton࣭bulk polymer electrolyte࣭electrode interface‫ ە‬prediction of transportation coefficient‫ ە‬mechanism
࣭free energy topographical map
࣭local and global paths
1,000,000 atoms

100,000,000 atoms

100,000 atoms

polymer electrolytebulk

4-phases interface
platinum, carbon, ionomer, and water

real interface㻶㻲㻱㻌䝔䜽䝜䝸䝃䞊䝏ᰴᘧ఍♫ ⇞ᩱ㟁ụ⏝㻼㼠ゐ፹ศᕸ䛾㼁㻸㼂㻙㻿㻱㻹ほᐹ㻔㻌㼔㼠㼠㼜㼟㻦㻛㻛㼣㼣㼣㻚㼖㼒㼑㻙㼠㼑㼏㻚㼏㼛㻚㼖㼜㻛㼚㼍㼚㼛㻙㼏㻛㼚㼍㼚㼛㼋㼏㼍㼟㼑㻜㻝㼋㼁㻸㼂㻙㻿㻱㻹㻚㼔㼠㼙㼘㻌㻕

Diffusions in the inhomogeneous systems
QV06'DWRPV䠖Ỉ┦䠖㧗ศᏊ┦W
QPᒁᡤⓗ䛺ᣑᩓQP= 2JXPL HWDO-(&6 132 ,QRUGHUWRGLIIXVHE\QPLWWDNHVQP 2 㽢-FP2V QP 2 㽢-QP2QV§ PV06'

R[\JHQ䚸K\GURJHQ '§㽢-FP2VPV

06'

Scale of inhomogeneity § QP኱ᇦⓗ䛺ᣑᩓ኱ᇦⓗ䛺ᣑᩓ

ᣑᩓ㡿ᇦᒁᡤⓗ䛺ᣑᩓWᒁᡤⓗ䛺▷᫬㛫ᣑᩓ䛾✵㛫ᖹᆒᐩᓅࡣࠊ኱つᶍ⣔ࡣᚓព࡛࠶ࡿࡀࠊ㛗᫬㛫ࡢ㏣㊧ࡣⱞᡭW

Coarse graining to a lattice model
Position-dependent free energy and diffusion coefficient

G(i, j, k)

32㽢32㽢32=32,768 latticecoarse grainingof particle dynamics

Dr. T. Nagai

D(i, j, k)

One-particle dynamics
on the lattice

H2 in Nafion O=6dynamic Monte Carlocoarse grainingof space
G(i) : ⢏ᏊᤄධἲࠊOD
D(i) : flat bottom

statistically correct
diffusive trajectories
on G(r) and D(r) surfaces⮬⏤࢚ࢿࣝࢠ࣮໙㓄ࢆᕪࡋᘬ࠸ࡓ᫬ࡢᣑᩓಀᩘ㸦㧗ศᏊ┦࡜Ỉ┦㸧Free energy calculations'Vinsertion

Particle insertion method'G­° 'V (r N , rN 1 ) ½°
kT ln exp ®
¾kT¯°¿° N , N 1insertion

Overlapping distribution method
finsertion 'V

f removal 'V'V
2kT'V
ln Premoval 'V 2kT
ln Pinsertion 'V 'Vinsertion'VremovalPinsertion SUREDELOLW\WKDW'Vinsertion LV'V
Premoval SUREDELOLW\WKDW'Vremoval LV'V'G

kTf removal 'V  kTfinsertion 'Vinsertion

removalFree energy calculations
Comparison between test particle insertion method
and overlap distribution method
For 1024 positions solvation free energy of hydrogen and oxygen was evaluated via OD and TPI
As to hydrogen, TPI and OD methods agree under 15 kJ/mol for solvation free energy.
As to oxygen, TPI and OD methods agree under 10 kJ/mol for solvation free energy.

For gas diffusion, TPI may be sufficiently accurate

Spatial resolution:4ύ㽢 4ύ㽢 4ύ
Ensemble of solvent for TPI: 5ns

&RPSDULVRQRI+2 DQG22200,000 ཎᏊ䚸32,768 meshes

H2O = 6, 'x=0.4 nmO2

㓟⣲䛸Ỉ⣲ඹ䛻䚸⤒㊰䛻䛺䜚䛭䛖䛺⮬⏤䜶䝛䝹䜼䞊䛾ప䛔ሙᡤ䛿㧗ศᏊ䛾ෆ㒊䛻Ꮡᅾ

1DILRQᵓ㐀࡜㏱㐣ศᏊࡢ⮬⏤࢚ࢿࣝࢠ࣮㧗ศᏊ䛾ศᏊᵓ㐀䛸Ỉ⣲䛾⮬⏤䜶䝛䝹䜼䞊ᆅᙧ

ȴ& [kJ/mol]㧗ศᏊ┦䛾⮬⏤䜶䝛䝹䜼䞊ᆅᙧ䛿㉳అ䛻ᐩ䜐⮬⏤䜶䝛䝹䜼䞊䛾ప䛔ሙᡤ䛿ᇶᮏⓗ䛻㧗ศᏊෆ㒊Ỉ┦䛾⮬⏤䜶䝛䝹䜼䞊ᆅᙧ䛿ᖹᆠ䛷㧗䜑䛾್䜢䜒䛴

12

Percolation cluster analysis
O=6
Identification of diffusion pathG < 5 kJ/molỈ⣲ศᏊࡣ࡝ࡇࢆ㏻ࡿ࠿
O=6

Present workỈ┦Previous studies
through water phase
T. Sakai et al., J. Electrochem. Soc. 133, 88 (1986).
H. Ito et al., Int. J. Hydrogen Energy 36, 10527 (2011).

⏺㠃through interface䠖Z. Ogumi et al., J. Electrochem. Soc. 132, 2601 (1985).
M. Inaba et al., Electrochim. Acta 38, 1727 (1993).

㧗ศᏊ┦through water and polymer phases
K. Borka et al., J. Appl. Electrochem. 27, 117 (1997).
M. Schalenbach et al., J. Phys. Chem. C 119, 25145 (2015).

14

Position dependent diffusion coefficient
Berendsen et al.

Constrained at position r(delta function-like constraint)

0DUULQN DQG%HUHQGVHQJPC 98  DS r

2RT'F r , tf

³ 'F r , t 'F r , 0 dt

Assuming Langevin equation

0Woolf et al.
:RROIDQG5RX[JACS 116 

Lipid bilayers

DS r

Gr 0 Gr 0f0

Constrained at position r by spring

2³ G r t G r 0 dt

F r, t  F r, t

Gr tr t  r t

Assuming generalized Langevin equationwater molecule in bulk water
Error from the assumption
20~30% for the solute
molecule whose size iscomparable to the solvent
molecules

DVFP2V-࣭Small molecules smaller thanthe solvent molecules࣭Small molecules in solids

Berendsen et al.

Woolf et al.

Regular MD

3.5㽢-5㽢-54.2Ң4.4㽢-5࣮࢝࣎ࣥ㧗ศᏊ↓ᶵᮦᩱ

Greater errors

New methodFlat-bottom potential
Generation of fragments of trajectory without perturbations+LJKOLJKWHGSRUWLRQRIWKHWUDMHFWRU\LVDQDO\]HG
㻣㻌䊅x

1RH[WHUQDOILHOGDFWVLQWKHKLJKOLJKWHGSRUWLRQ㻣䊅㻌“size of a box”

㻸㼏“core”
regionWLPH

Flat-bottom potential

Mean squared displacement

㻸㼒flat region

Correction for the calculated MSD)RUPIDFWRUREWDLQHGDVDIXQFWLRQRI'WLV
LQGHSHQGHQWRI'MSD(t)=6Dt (+C)

NOT sampled

MSDFB(t)FF(Dt)06')%LVFRUUHFWHGXVLQJ)) 'W 
VXFKWKDWWKHVORSHRIWKHFRUUHFWHG㻸㼏IXQFWLRQLV'FF(Dt) is optimized such that
MSDFB(t)/FF(Dt)=6Dt +C(HJ %D\HVLDQLQIHUHQFH

DVFP2V-

MSDFB(t)sampledBerendsen et al.

Woolf et al.

Regular MD

Present study

3.5㽢-5㽢-5

4.2Ң4.4㽢-5㽢-5࣭ᣑᩓ㡿ᇦࡲ࡛ࢧࣥࣉࣜࣥࢢ
࣭グ᠈ࡣ↓ࡋࠊStochastic࣭⛣ὶ㡯ࢆᾘཤࡋ࡚DࢆồࡵࡿV total (r ) V molecular (r )  V FB (r )  G (r )wP( x, t )wtw wwU ( x) ½­w® D( x)  D( x) E
¾ P ( x, t )wx wx
wx ¿¯ wxt

Coarse-graining to diffusive dynamics on a latticeMD┤᥋ィ⟬
G(i) : ⢏ᏊᤄධἲࠊOD
D(i) : ㄪ࿴᣺ືᏊࠊflat bottom㧗ศᏊࡣ✵㛫࡟fixࡉࢀ࡚࠸ࡿ
཯ᛂᗙᶆ : (i, j, k)ࠊ(x, y, z)G(i), D(i)orMSM
trajectoryゎᯒTij, wij⢏Ꮚࢲ࢖ࢼ࣑ࢡࢫ⁐ゎᣑᩓࣔࢹࣝ

Kramers modelゎᯒⓗ⾲⌧࣭trajectory࣭ࡍ࡭࡚ࡢ᝟ሗ࣭ከᵝ࡞ゎᯒ࣭┤ឤ࣭ཝᐦ

࣭1ḟඖࣔࢹࣝ
࣭㏆ఝ

࣭㏿ᗘᐃᩘࡢゎᯒⓗ⾲⌧࣭1ḟඖࣔࢹ࣭ࣝ㏆ఝ

࣭MFPTࠊcommitter࣭ཝᐦ⢏Ꮚࢲ࢖ࢼ࣑ࢡࢫࡀ‶㊊ࡍ࡭ࡁ᮲௳࣭ᣑᩓ᪉⛬ᘧ࣭ヲ⣽㔮ࡾྜ࠸wP( x, t )wtw wwU ( x) ½­w® D( x)  D( x) E
¾ P ( x, t )wx wx
wx ¿¯ wx

Diffusive particle dynamics
No MC algorism
for position-dependent D(r)wP( x, t )wtdiffusion equation
Fokker-Plank equation

Metropolis MC
Proof by Kikuchi
for constant Dw w
wU ( x) ½­w® D( x)  D( x) E
¾ P ( x, t )wx wx
wx ¿¯ wxparticle

Monte Carlo

distribution

Doob, Szabo for positiondependent U(r) and D(r)master equation

suitable to both lattice andcontinuous coordinate system

Gillespie algorithm
suitable to lattice but not suitable tocontinuous coordinate system

Langevin equation
Brownian dynamics
suitable to continuous
coordinate system but
not suitable to latticenumerical stability

New MC algorithm whose molecular trajectories satisfy
diffusion equation with position-dependent U(r) and D(r)Ỉ⣲ࡢ㍺㏦
G(i,j,k)Dynamic MC in G and D spaceD(i,j,k)32,768 meshesỈ⣲ࡢ኱ᇦⓗ࡞ᣑᩓ㍺㏦⤒㊰21

Thank you!