Finite Electric Field Calculations
based on the Berry Phase Method
in OpenMX
Naoya Yamaguchi
Nanomaterials Research Institute(NanoMaRi), Kanazawa University

First-principles LCPAO
approach for insulators
under finite electric fields
with forces
Comput. Phys. Commun. 280, 108487 (2022).
(arXiv:2203.10441)

2

Crystalline potential
Electronic density

EnergyElectronic density𝑛(𝑟)ԦIn a crystal, theelectronic density is
distributed.

𝑉𝐾𝑆
Position

Period (Lattice constant)3

Energy

Crystalline potential
Sawtooth potential

Sawtooth potential
is added to apply
an electric field E.

𝑉 = 𝑒𝑬 ⋅ 𝒓
Position

Elementarycharge

Position4

Potential under auniform electricfield E

EnergyCrystalline potential
Sawtooth potential
Under an electric field

𝑉 = 𝑒𝑬 ⋅ 𝒓
Position

Elementarycharge

Position𝑬
𝑉𝐾𝑆 = 𝑉𝐾𝑆 + 𝑒𝑬 ⋅ 𝒓

5Energy

Crystalline potential
Sawtooth potential
Under an electric field

Periodic boundary
condition gives
unphysicalpotential.
In systems withboundaries and nounphysical regions in thevacuum region, the use ofa sawtooth potential is
not applicable.

𝑉 = 𝑒𝑬 ⋅ 𝒓
𝑬𝑉𝐾𝑆 = 𝑉𝐾𝑆 + 𝑒𝑬 ⋅ 𝒓
Position

6Energy

Crystalline potential

𝑉 = 𝑒𝑬 ⋅ 𝒓

For systems withoutedges such as crystals,
potential that is asingle-valued functionof the position is notappropriate.
→ A multivalued function
seems to be appropriate?

𝑒𝒓Polarization 𝑷 can be→ 𝑷, 𝑉 = Ω𝑬 ⋅ 𝑷 expressed by a
Ω Cell volume
multivalued function.
Position

7

Homogeneous electric fields
in insulating solids based on MTP
R.W. Nunes, D. Vanderbilt, Phys. Rev. Lett. 73, 712 (1994).

Electric fields in first-principles calculations using a plane-wave basis
Periodic part of Bloch orbitals

I. Souza, J. Íñiguez, D. Vanderbilt, Phys. Rev. Lett. 89, 117602 (2002).

P. Umari, A. Pasquarello, Phys. Rev. Lett. 89, 157602 (2002). (Car-Parrinello method)
J.W. Zwanziger, J. Galbraith, Y. Kipouros, M. Torrent, M. Giantomassi, X. Gonze,
Comput. Mater. Sci. 58, 113 (2012). (projector augmented wave framework)
The implementation for finite electric fields with forces in the localized basis set,
particularly in the pseudo-atomic orbital linear combination (LCPAO) method,which is ideal for nanoscale simulations, have not been developed until then. 8

Electric polarization in solids(modern theory of polarization (MTP))(cf. R. D. King-Smith, and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).)
Electronic contributionof polarizationSpindegeneracy

Reciprocallattice vectors

Berry phases
Spin
Overlap matrices between
two successive k-points
Periodic part of Bloch orbitals

e.g. 𝑁𝑏 = 𝑁𝑐 = 2, i.e. , 𝐼𝑏 , 𝐼𝑐 = 0, 0 , 0, 1 , 1, 0 , 1, 1 ;𝑁𝑎 = 5, i.e. , 𝐼𝑎 = 0, 1, 2, 3, 4bcSolving the eigenvalue problem at eachwave number (k-point) to obtaineigenvectors, from which M can becalculated.

9Linear combination of pseudo atomic orbitals（LCPAO）
The Bloch function of band μ can be expressed asfollows:

N: the number of cells considered in the calculation
Rn: lattice vector of cell n. c: LCPAO expansion coefficients.
𝜙: PAO. i: index of the atom. α: index of the orbital.

The previous implementation of Mμν calculations in polB (cf. the technical note).

10Efficient LCPAO approach for MTP
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).

The matrix form is optimal in terms of efficiency and simplicity incalculations, and it is particularly effective for iterative calculations like
SCF calculations.

T(k) need preparing only onceat the beginning of each SCF calculation.

11

LCPAO approachfor finite electric field in insulating solids
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).

Extract the matrix
elements A of theelectric field potential V

12

LCPAO approachfor finite electric field in insulating solids
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).

i−
lattice vectori−Constructed from information at the three kpoints including the points before and after

Adding A to the Hamiltonian results in a Hamiltonian under a uniform electric field.
→We implemented it in OpenMX

13

Force acting on an atom
under an electric field
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).
Ionic contribution to
the electric polarization

Ionic
Position ofatom i

Electronic

Ionic charge

In LCPAO, unlike a complete basis set like plane waves, theelectronic part remains as a Pulay force.
𝝏𝐓→ Only 𝒄† 𝝏𝝉 𝐜 needs to be explicitly computed.
𝒊Hamiltonian𝜕𝑐 †𝜕𝑐（Other terms are included in 𝜕𝝉 Hc + 𝑐 † H 𝜕𝝉 ）𝑖𝑖

Evaluation of Forces → effective charges & static dielectric constant.
14

Electronic dielectric constants
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).

This work

1 𝜕𝑃𝑥
𝜀𝑥𝑥 = 1 +𝜀0 𝜕𝐸𝑥
Dielectric
constant

Electricconstant

LDA: [5] Phys. Rev. Lett. 89, 117602 (2002).
GGA: [9] Comput. Mater. Sci. 58, 113 (2012);[20] Phys. Rev. B 95, 075302 (2017).
expt.: [21] Properties of Group-IV, III-V and II-VI Semiconductors,
15John Wiley & Sons, Ltd, 2005, pp.211–281.

Born effective charges
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).

This work

Born effectivecharge ∗[8] Phys. Rev. Lett. 89, 157602 (2002); [22] J. Appl. Phys. 121, 125701 (2017).
[23] Phys. Rev. B 75, 115116 (2007); [24] Phys. Rev. B 50, 10715 (1994).
[25] AIP Conf. Proc. 1653, 020065 (2015).; [26] Phys. Rev. B 50, 14125 (1994).
[27] J. Vac. Sci. Technol. B 39, 052201 (2021);
Force[28] Solid State Commun. 57, 483 (1986).
𝑥[29] Phys. Rev. B 4, 1367 (1971).
[30] AIP Adv. B 4, 067138 (2014).
𝑥
The values in [9] and [30] were those in GGA,
16
and the other reference values were those in LDA.

1 𝜕𝐹
𝑍𝑥𝑥 =𝑒 𝜕𝐸

Static dielectric constants
N. Yamaguchi, F. Ishii, Comput. Phys. Commun. 280, 108487 (2022).

System LDA
AlP
10.26

GGA
9.94

AlAs

10.90

10.56

GaPMgO

11.96
7.92

11.69
8.79
9.56 This work(for expt. lattice constant )Ref.(theor.)10.2(LDA[5]),10.26(GGA[9])11.5(LDA[5]),11.05(GGA[9])11.2(LDA[5])7.93(LDA[8])Ref.(expt.)
9.6[21]
10.6[21]11.0[21]9.8[21]The static dielectric constants were accurately determined.
→ the electronic states and forces acting on atoms under anelectric field were calculated with sufficient accuracy.
→ First-principles molecular dynamics under an applied electricfield in LCPAO can be achieved at a sufficient level.
17

Summary• We developed a method to construct the effective
potential for applying an electric field using theBerry phase method within the framework of theLCPAO method, including the forces.
• In our implemented code, the calculated values forelectronic and static dielectric constants and Born
effective charges match the reference values.
• This method enables large-scale first-principles
calculations under an electric field, making it
possible to study the behavior of electronic states
under electric fields even in complex systems.
This section is based on
Comput. Phys. Commun. 280, 108487 (2022). (arXiv:2203.10441)
18

New Functionsfor Finite Electric Fields
in OpenMX19

Implementation• Output_CompTime.c• HWC.c• Force.c• HWCfunc.c• Make_InputFile_with_FinalCoord.c• Cluster_DFT_Col_FEF.c• Total_Energy.c• Hamiltonian_Cluster_Hs_FEF.c• Set_XC_Grid.c• Cluster_DFT_NonCol_FEF.c• DFT.c• Hamiltonian_Cluster_NC_Hs2_FEF.c• Input_std.c• mimic_read_scfout.h• FEF.c• Band_DFT_Col_FEF.c• Band_DFT_NonCol_FEF.c

(1)

DFT

• HWCfunc.h• openmx_common.h(Band,Cluster)_DFT_(NonCol,Col)_FEF

In 1 SCF step(1) Eigen value problem(2) Effective potential A(k)

(2)
FEF

HWF (Berry Phase, P(polB) + A)

20

Input keywords• scf.Electric.Field.Type

BerryP # sawtooth|BerryP• scf.Electric.Field
0.1 0.0 0.0 # in GV/m
# x, y, z in BerryP, while a, b, c in sawtooth• FEF.Polarization.Switch
On # Off|On# To evaluate electric polarization P• FEF.Continue.Mode
Update # Off|Zero|Update# Too strong an electric field is(not acceptable|not applied|reduced and applied).
• FEF.Born.Effective.Charges.Switch# To evaluate F_(x,y,z)/EOn # Off|On21

Calculation process
1. Get the SCF convergence w/o electric field
2. Break the convergence flag
3. Calculate the k-independent part of T
4. Calculate A(k)5. Start the next SCF step:
In setting the Hamiltonian, A(k) is added to theglobal form of the Hamiltonian (k-dependent).
6. Check the SCF convergence7. Convergence → finish; otherwise → 422

Example: Diamond23