The 35th CMD workshop”Practice of Esopt”
Introduction
Koichi KusakabeGrad. Sch. Eng. Sci, Osaka University

Let’s start the X window system◼ We use machines operated by a Linux system .
◼ Each user is identified by “user id” & “password”.
◼ Your terminal might start by pushing the monitor button.

Use the id of “cmd” for ubuntsu with pw, “c.m.dworkshop”.
cmdパスワード

Access to the PC cluster in phys.sci.
◼ Login the X terminal machine:◼ Account : stud** or cmd◼ Password : none or “cmdworkshop”◼ Open a terminal window, and access `rl’.
[stud** ~]$ ssh –Y –l stud** rl.phys.sci.osaka-u.ac.jp

Enter passphrase for key :…
[stud**@rl ~]$Passphrase for ‘rl’ is ‘cmdstud**’ .

Preparation for the program package◼ Copy the program package.
Type the red characters in a single line without “return” in the middle.

[stud**@rl ~] cp /home/CMD/teac01/ESopt-CMD35-fw.tgz .
“a space and a period” are requested.

◼ Unpack the package by tar command.
[stud**@rl ~] tar zxvf ESopt-CMD35-fw.tgz

How to compile the program?
◼ Change the directry to the main directry of ESopt, checkthe header of makefile, then `make’.
[stud**@rl ~]$ cd ESopt-CMD35-fw[stud**@rl ~/ESopt-CMD35-fw]$ ls
a0ene.pl
getll.f90
mkvel.f90 ….
….
[stud**@rl ~/ESopt-CMD35-fw]$ head makefile# The makefile for ES_opt#DESTDIR=Si-band#

Please rewrite this line!!!

<= this is the destination…
[stud**@rl ~/ESopt-CMD35-fw]$ cp INPUT_DATA/CORD-Si-b
INPUT_DATA/CORD[stud**@rl ~/ESopt-CMD35-fw]$ make < input_makeRun the program◼ Change the directory to the destination written

in the header of the makefile
DESTDIR=Si-band <= you may rewrite this name.
◼ If you have a working directory named `Si-band’,

you may submit a job using the queing system,[stud**@rl ~/ESopt-CMD35-fw]$ cd Si-band[stud**@rl Si-band]$ qsub qsub
[stud**@rl Si-band]$ qstat

When you are doing “bandstructure calculation”,
after finishing the job, you can type,$ ../utility/vbpef2gp –s –e fort.50$ gnuplotWhat happens, when you compile it?
◼ At beginning, make creates executable file,“opt”.
◼ Once “opt” is made, make will create a working
directry for each project given in
INPUT_DATA/CORD.
◼ The name of the working directry is given in theheader of makefile.
◼ Some temporary data files are created in a sub
directry “mknon”.

Bulk modulus: Second derivative of E(V)Vdp Vd 2 E=,
Bulk modulus at T=0: B0 = −2dVdV
B ' B0V0
B0V  (V0 / V ) 0Murnaghan eq. of states: E (V ) = E0 ++ 1 −.
B0 '  B0 '−1B'−10
calc_mode: bas
band_calc: metal
Ecut = 40 Ry.
“Total energy” is seen
by looking atMONIT/energy.

In this calculation,“murnfit” tells that
B0=89GPa.

A non-linear fitting problem◼ Murnaghan’s equation of states
B0 '
B0V (V0 / V )+CE (V ) =+1
B0 '  (B0 '−1)= V 1− B0 ' + V + B 'V B '== V  +  V + 000B0 ' ( B0 '−1), =B0,  = C ,  = 1 − B0 'B0 'We let four sets of (V, E(V)) be (xi,yi) (i=1,4).

yx
 1
   1 y 2  =  x2
 y   x
 3  3x1 1   x2 1     ( ),  ( ),  ( )x3 1  ()g ( ) = y4 −  ( )x4 +  ( )x4 +  ( )
Finding  which gives
g()=0 is important.

Trial 1: Bulk modulus of c-BN
◼ Change the directory to the destination in “makefile” as,DESTDIR=cBN <= you may rewrite this name.
◼ If you see utility/cBN_scripts, CORD files for the test as well assome utility scripts, chev.sh, sub.sh, and summ.sh are stored.
[stud**@rl ~/ESopt-CMD35-fw]$ cd utility/cBN_scripts[stud**@rl cBN_scripts]$ cp CORD_5 ../../INPUT_DATA/CORD[stud**@rl cBN_scripts]$ cd ../../[stud**@rl ~/ESopt-CMD35-fw]$ head makefile[stud**@rl ~/ESopt-CMD35-fw]$ make < input_make[stud**@rl ~/ESopt-CMD35-fw]$ cd cBN[stud**@rl cBN]$ qsub qsub

◼ Create a0ene.dat and a0ene.dati. (An example is in utility/cBN_scripts)◼ You can do the fitting by asking murnfit in fit_esopt.
[stud**@rl cBN]$ pushd ../fit_esopt
Equilibrium volume: V0[stud**@rl fit_esopt]$ make
Bulk modulus: B0[stud**@rl fit_esopt]$ popd
Pressure deriv. of bulk modulus: B0’[stud**@rl cBN]$ ../fit_esopt/murnfit

Flow diagram to check accuracy
Select the target system
GoodP.P.
NoAre P.P.s available？

Create P.P.

Convergence?
Logarithmic derivative?

Not good

Yes
Test of P.P. OK？

No

Test calc.
Estimation of params.

Accuracy？
Band str.? Ghost?

Yes
Estimation of Ecutetc. OK？

GoodNoNo

YesCalc. of the target!

Accuracy OK?
P.P. good？

YesInterpretationof the result

Ecut dependence of the total energy
calc_mode: bas
band_calc: metal

For an ordinary use, 40 Ry may be allowed.
To have a converged stress tensor, Ecut > 60 Ry would be required.

Trial 2: Cutoff check of PP◼ Select an element among “Al, Ca, Cu, Mg, Mn, Si”.
◼ You can see an example of “CORD” for them in INPUT_DATA.
◼ Change the directory to the destination in “makefile” like,

DESTDIR=Al <= you may rewrite this name.
◼ Using an example, try the cutoff test.
[stud**@rl ~/ESopt-CMD35-fw]$ cp INPUT_DATA/CORD-Al INPUT_DATA/CORD[stud**@rl ~/ESopt-CMD35-fw]$ emacs INPUT_DATA/CORD… please “edit” the input file “CORD”. …
[stud**@rl ~/ESopt-CMD35-fw]$ head makefile[stud**@rl ~/ESopt-CMD35-fw]$ make < input_make[stud**@rl ~/ESopt-CMD35-fw]$ cd Al[stud**@rl Al]$ qsub qsub

◼ Create a data file of “the lattice constant a” v.s. “the total energy”.
◼ Total energy may be seen in MONIT/energy, RESULTS/Htot001 (for basicmode), or STRUCT/result0001 (for extended mode).

Some post-processing commands◼ For writing a band structure calculation,◼ Select band_calc = ‘bands’ in INPUT_DATA/CORD,◼ Do ‘make’ after checking
INPUT_DATA/kpoint◼ Top of the ‘makefile’◼ Do run your job by ‘qsub qsub’ in the workingdirectry.
◼ After finishing the job, you can type,◼

$ ../utility/vbpef2gp –s –e fort.50◼ $ gnuplot◼ We now see the functions of ‘xcrysden’ by seeing ‘axsf’ file inyour ‘Si-band’.
◼

Kohn-Sham band structure of Si
calc_mode: bas
band_calc: bands

For this band calculation, Ecut = 60 Ry is used.
ncband=6 in CORD is used to obtain unoccupied levels.

Band structures of SWNT
Tight binding modelwith p orbitals only

DFT-LDA using optimizedstructure1/3-rule in a simple TBM band structure

Saito, Dresselhaus, Dresselhaus, ‘Physical Properties of Carbon Nanotubes’,Imperial College Press.

LDA band structures for thin SWNTMachón, et al. Phys. Rev. B 66, 155410 (2002).

Trial 3: band structure of CNT◼ We can use prepared inputs as follows.
[stud**@rl ~/ESopt-CMD35-fw]$ cp INPUT_DATA/CORD-50Tube-b
INPUT_DATA/CORD# Check E_cut in CORD
[stud**@rl ~/ESopt-CMD35-fw]$ cp INPUT_DATA/kpoint-b
INPUT_DATA/kpoint[stud**@rl ~/ESopt-CMD35-fw]$ emacs –nw makefile<< rewrite header of makefile so that “DEST-DIR = Tube-band” >>[stud**@rl ~/ESopt-CMD35-fw]$ make < input_make[stud**@rl ~/ESopt-CMD35-fw]$ cd Tube-band[stud**@rl Tube-b]$ qsub qsub
# We need to wait for a while. If we need, please see the result tomorrow.
[stud**@rl Tube-b]$ ../utility/vbpef2gp –s –e fort.50[stud**@rl Tube-b]$ gnuplot
gnuplot> load ‘bandfig.gp’The 35th CMD workshop”Practice of Esopt”
Overview of optimization problems
Koichi KusakabeGrad. Sch. Eng. Sci, Osaka University

Multi-dimensional optimization problem
We regard next quantities as
the variables determining of
the system.
1. Electronic state

2. Atomic position
3. Unit cell (lattice parameters)

Restricted by the periodic boundary
condition, various electronic states
appear depending on atomic positionin the variable unit cell.

Multi-dimensional optimization problem
In the actual simulation, we have

3. Electronic state

2. Atomic positions
1. Unit cellto be optimized in an order ofthe number.

Multi-dimensional optimization problem
Structural change given by modification of lattice vectors
induced by difference between internal stress and theexternal pressure.

a3a3

a2

a2
a1

a1

9 - 3 = 6 degrees of freedom after removing total rotation of the system.
1. Lattice constants
2. Elastic constants
3. Structural transformation← optimization of all parameters← deformation of stable structure← molecular dynamics

Multi-dimensional optimization problem
Structural change and electronic structure of the structure
obtained by shift in atomic position due to the interatomic force.

Optimization ofElectronic state

Interatomic force
ESopt
Shift in the 3×Natom structural parameters
Change in the crystal symmetry
We can determine,
1. Stable atomic configuration
2. Intermediate state (transition state), reaction energy

Multi-dimensional optimization problem◼ Perform optimization of essential parameters tomake the energy be the lowest. Theoptimization is done with respect to◼ Lattice constants◼ Atomic configuration◼ Electronic structure◼ The optimization process is very similar,

because all of the procedure is regarded as
motion of a particle on an energy surface.

Multi-dimensional optimization problem◼ Technique of the optimization◼ Global search ： most powerful. Global structure of theenergy surface should be known before the optimization.
◼ Searching methods utilizing a gradient vector ： very fast.
Optimization follows a line toward the minimal point.

P2

L2

L1

P1

Multi-dimensional optimization problem◼ Available methods◼ Metadynamics ： search of potential minimal points
beyond a potential barrier◼ Search of saddle points ： force inversion technique
Making the change in thepotential surface, a motion go
beyond the barrier is realized.

P2

L2

L1

Force inversion make the motion
going up the hill.

P1

Methods to find the transition states
Nudged elastic band method：• Replicas （yellow points）feel repulsion.
• They give a path through TS.

Force inversion：• If the landscape is simple,• The transition state is fixed.

BO surface

L2

P2

P2

L2

L1

L1

P1

P128

Minimum search on the Born-Oppenheimer(BO) surface

If one look for the minimum locally around this point, one cannot
find existence of another true minimum.
• Electronic structure at this point shows only local stability.
• Even if one make small vibration around this local minimum, onlystable oscillation will appear.
In this situation, we need global search on the BO surface.
But, search in the whole 3×Natom-dimensional space ispractically impossible.
: At present, for the first-order transition, the first-principles
molecular dynamics can do little without wisdom on materials.

Search for the saddle point on the BO surface
If one has some locally stable points, the saddle point of the BO surface may
be found. But, we could have difficulty even in this case.

One dimensional search:
Optimization methods which cansearch a minimum point in asection perpendicular to thereaction coordinate may give thesaddle point.

In this case, this saddle point
given by the one dimensional
search does not tell us the true
potential barrier, since we haveanother way to go without thebarrier.

Unstable fixed points

A saddle point
If we optimize all the coordinate
other than the reaction coordinate,
the saddle point is given.

But, we may have restriction due toundesired symmetry in thesimulation due to1) Finite simulation cell2) Assumed symmetry in the run.

A unstable fixedpoint which may
be obtained as asaddle point in a
A section of BOS
restricted search.
At the unstable fixed point, interatomic
forces are balanced. If one look at thesecond-order derivative (the Hessianmatrix), we may have several paths tosome stable points.

An example in the structural transformationof carbon systems
by Y. Tateyama et al.

The force inversion techniquedetermines two saddle points fromthe graphite to two diamondstructures.

Saddle point

Unstable fixed point

The meta-dynamics simulation
Following the motion of a system inthe molecular dynamics, we add
artificial potential to the true potential.

Vg ( x) = W exp(− x(t ) − x(ti ) 2h 2 )2ti tUtilizing this technique, we can reach 1)
the lowest saddle point, and 2) a stablestructure beyond the saddle point. If weproceed to perform this simulation without
removing Vg, we should reach all of thelocal minima on the BO surface.

A. Laio, M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (2002).
R. Martonak, et al., Phys. Rev. Lett. 90, 075503 (2003).

Recommended procedure of the structural
optimization◼ The case of bulk crystal structure◼ Rough estimation of the stable structure with theconstant-pressure MD.
◼ Optimization of the lattice constant by accurate
electronic structure calculation.
◼ Estimation of the bulk modulus and the elastic
constants. Comparison with the experimentallyobserved values.
◼ Confirmation of positivity in the phonon modes.
◼ If an imaginary frequency appears, this may
indicate that the simulation cell is too small. Try
the larger simulation cell.

Recommended procedure of the structuraloptimization◼ The case of bulk crystal structure◼ Rough estimation of the stable structure with theconstant-pressure MD.
◼ Optimization of the lattice constant by accurate
electronic structure calculation.
◼ Estimation of the bulk modulus and the elasticconstants. Comparison with the experimentallyobserved values.
◼ Confirmation of positivity in the phonon modes.
◼ If an imaginary frequency appears, this may
indicate that the simulation cell is too small. Try
the larger simulation cell.

MD simulation of graphite in high pressureStructural optimization method◼ Optimization of bulk crystal structure◼ Estimate roughly by the constant pressure mode.
◼ Improve the accuracy, optimize lattice constantsfor the case with the unit e.
◼ Calculate the bulk modulus and elastic constants
and compare the real observation.
◼ Check that phonon modes have positive realfrequency. (We need the super cell technique.)◼ If an imaginary frequency appears, find aconstraint in the simulation, remove therestriction and re-optimize the parameters.