pdf Version

32 downloads 469 Views 2MB Size Report
local part and Kleinman-Bylander projectors) in SIESTA and ABINIT. Compare total ... Regarding the ABINIT code, you can download the required version from:.
How to run with the same pseudos in SIESTA and ABINIT -215.78 Siesta (DZP) Abinit (13.005 Ha)

Free Energy (eV)

-215.79

-215.80

-215.81

-215.82

5.30

5.35 5.40 5.45 Lattice constant (Ang)

5.50

Objectives Run examples with the same pseudos (same decomposition in local part and Kleinman-Bylander projectors) in SIESTA and ABINIT. Compare total energies

Download the last versions of both codes, SIESTA and ABINIT Regarding the ABINIT code, you can download the required version from:

http://personales.unican.es/junqueraj/Abinit.tar.gz But the merge of the relevant subroutines into the main trunk will be done soon

Regarding SIESTA, the code is available at the usual web site:

http://www.icmab.es/siesta

A few modifications to be done before running: SIESTA Edit the file atom.F in the Src directory I.

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIELD.

1. Replace “nrval”1by “nrwf” in the call to schro_eq inside the subroutine KBgen (file atom.F) CiteI.. ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC

I. ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIEL if(erefkb(ikb,l).ge.1.0d3) then Cite1 . nnodes=ikb Cite1 . nprin=l+1 call schro_eq(Zval,rofi,vps(1,l),ve,s,drdi, if(erefkb(ikb,l).ge.1.0d3) then if(erefkb(ikb,l).ge.1.0d3) then . nrwf,l,a,b,nnodes,nprin, nnodes=ikb nnodes=ikb erefkb(ikb,l),rphi(1,ikb)) . nprin=l+1 nprin=l+1 call schro_eq(Zval,rofi,vps(1,l),ve,s,drdi, call schro_eq(Zval,rofi,vps(1,l),ve,s,drdi, if(nkbl(l).eq.1) then . nrwf,l,a,b,nnodes,nprin, . nrwf,l,a,b,nnodes,nprin, call in ghost(Zval,rofi,vps(:,l),vlocal, 2. Replace “nrval” by. “nrwf” the call toerefkb(ikb,l),rphi(1,ikb)) ghost inside the subroutine KBgen (file . ve,s,drdi,nrwf,l,a,b,nrwf, . erefkb(ikb,l),rphi(1,ikb)) atom.F) . erefkb(ikb,l),rphi(:,ikb),ighost) if(nkbl(l).eq.1) then else if(nkbl(l).eq.1) then call ghost(Zval,rofi,vps(:,l),vlocal, .call ghost(Zval,rofi,vps(:,l),vlocal, ve,s,drdi,nrwf,l,a,b,nrwf, real(dp), parameter :: Rmax_kb_default = 60.0d0 . erefkb(ikb,l),rphi(:,ikb),ighost) . ve,s,drdi,nrwf,l,a,b,nrwf, else

.

erefkb(ikb,l),rphi(:,ikb),ighost) else parameter :: to Rmax_kb_default = 60.0d0 3. Increase the default ofreal(dp), Rmax_kb_default 60.0 bohrs (file atom.F) 1 2 3

::J. Rmax_kb_default = 60.0d0 G. real(dp), Lee, Y.-H. Shin,parameter and J. Y. Son, Am. Ceram. Soc., 95, 2773 (2012). E. Wiesendanger, Ferroelectrics, 6, 263 (1974). H. P. Sun, W. Tian, X. Q. Pan, J. H. Haeni, and D. G. Schlom, Appl. Phys. Lett., 84, 3298 (2004).

Definition of the Kleinman-Bylander projectors X. Gonze et al., Phys. Rev. B 44, 8503 (1991) The normalized Kleinman-Bylander projectors are given by

where

and

Note that these

are the eigenstates of the semilocal pseudopotential (screened by the pseudovalence charge density).

are the radial part of the wave function multiplied by

,

For the sake of simplicity, we assume here only one projector per angular momentum shell. If more than one is used, they must be orthogonalized

Definition of the Kleinman-Bylander projectors (old choice in SIESTA for the atomic eigenstates) In the standard version of SIESTA, the Schrödinger equation for the isolated atom while generating the KB projectors is solved inside a box whose size is determined by nrval. This is usually a very large radius (of the order of 120 bohrs)

s - shell p - shell d - shell f - shell

ul

PS

(r)

1,0

0,5

0,0

0

10

Rmax_KB

20 30 radius (bohr)

40

50

Then, this wave functions is normalized inside a sphere of a much smaller radius, determined by Rmax_KB (default value = 6.0 bohr)

Definition of the Kleinman-Bylander projectors (new choice in SIESTA for the atomic eigenstates) 0,8

s - shell p - shell d - shell f - shell

0,6

0,4

ul

PS

(r)

In the new version of SIESTA, everything is consistent, and the Schrödinger equation and the normalization are solved with respect the same boundary conditions

0,2

0,0

0

10

20

30 40 radius (bohr)

50

60 Rmax_KB

Almost no change in total energies observed, but the Kleinman-Bylander energies might be very different, specially for unbounded orbitals

A few modifications to be done before running: ABINIT Edit the file src/65_psp/psp5nl.F90 1. Uncomment the last two lines for the sake of comparing I. ENERGY FUNCTIONAL FOR Kleinman-Bylander A DIELECTRIC INSIDE energies and cosines with the ones obtained with SIESTA Cite1 . !DEBUG write(std_out,*) ’EKB=’,(ekb(iq),iq=1,4) write(std_out,*) ’COSKB=’,(ckb(iq),iq=1,4) !ENDDEBUG if(nkbl(l).eq.1) then call ghost(Zval,rofi,vps(:,l),vlocal, 2. Remember to compile the code enabling the FOX library . ve,s,drdi,nrwf,l,a,b,nrwf, . erefkb(ikb,l),rphi(:,ikb),ighost) else

./configure --with-trio-flavor=netcdf+etsf_io+fox

real(dp), parameter :: Rmax_kb_default = 60.0d0

1

AN ELE

Examples to run SIESTA and ABINIT with the same pseudos 1. Visit the web page: http://personales.unican.es/junqueraj and follow the links: Teaching Métodos Computacionales en Estructura de la Materia Hand-on sessions Pseudos 2. Download the Pseudos and input files for both codes 3. Untar the ball file $ tar –xvf Siesta-Abinit.tar This will generate a directory called Comparison-Siesta-Abinit with 4 directories: $ cd Comparison-Siesta-Abinit $ ls -ltr $ Si (example for a covalent semiconductor, LDA) $ Al (example for a sp-metal, LDA) $ Au (example for a noble metal, includes d-orbital, LDA) $ Fe (example for a transition metal, includes NLCC, GGA)

Examples to run SIESTA and ABINIT with the same pseudos

In every subdirectory it can be found: $ cd Si $ ls –ltr $ Pseudo (files to generate and test the pseudopotential) $ Optimized-Basis (files to optimize the basis set) $ Runsiesta (files to run Siesta) $ Runabinit (files to run Abinit)

How to generate and test a norm-conserving pseudopotential Generate the pseudopotential using the ATM code as usual, following the notes in the Tutorial “How to generate a norm conserving pseudopotential” Copy the input file in the corresponding atom/Tutorial/PS_Generation directory I. ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIELD. and run Cite1 . $ ../../Utils/pg.sh Si.tm2.inp ==> Output data in directory Si.tm2 ==> Pseudopotential in Si.tm2.vps and Si.tm2.psf (and maybe in Si.tm2.xml) if(nkbl(l).eq.1) then call ghost(Zval,rofi,vps(:,l),vlocal, The pseudopotentials will be on the same parent directory: . ve,s,drdi,nrwf,l,a,b,nrwf, .vps (required to test the pseudopotential) . (unformatted) erefkb(ikb,l),rphi(:,ikb),ighost) else .psf (formatted)

.xml (inreal(dp), XML format) to run Abinit) parameter(required :: Rmax_kb_default = 60.0d0 Remember to test the pseudopotential using the ATM code as usual, following the notes in the Tutorial “How to test the transferability of a norm conserving pseudopotential” 1 G. Lee, Y.-H. Shin, and J. Y. Son, J. Am. Ceram. Soc., 95, 2773 (2012). 2 3

E. Wiesendanger, Ferroelectrics, 6, 263 (1974). H. P. Sun, W. Tian, X. Q. Pan, J. H. Haeni, and D. G. Schlom, Appl. Phys. Lett., 84, 3298 (2004).

Running the energy versus lattice constant curve in SIESTA Run the energy versus lattice constant curve in SIESTA as usual. You can use both the .psf or the .xml pseudopotential. Follow the rules given in the tutorial “Lattice constant, bulk modulus, and equilibrium energy of solids” The input file has been prepared for you (file Si.fdf). Since we are interested in compare the performance of the basis set, it is important to converge all the rest of approximations (Mesh Cutoff, k-point grid, etc.) as much as possible 2 I. ENERGY FUNCTIONAL FOR to A DIELECTRIC INSIDE AN called ELECTRIC FIELD. At the end, we would be able write a file (here Si.siesta.latcon.dat) that looks like this: Cite1 .

# Lattice Constant (Ang) 5.30 5.32 5.34 5.36 5.38 5.40 5.42 5.44 5.46 5.48

Free Energy (eV) -215.790381 -215.802337 -215.811082 -215.816643 -215.819095 -215.818374 -215.814586 -215.808022 -215.798710 -215.786727

KBgen: Kleinman-Bylander projectors: l= 0 rc= 1.936440 el= -0.796617

Ekb=

These data have been obtained with a double-zeta plus polarization basis set, optimized at the theoretical lattice constant with a pressure of 0.05 GPa

4.661340

kbcos=

0.299549

Running the energy versus cutoff energy in ABINIT To check the equivalent cutoff energy in Abinit: 1.  We run the same system (same lattice vectors and internal coordinates) at the same level of approximations (same exchange and correlation functional, Monkhorst-Pack mesh etc.) at a given lattice constant. Here it has been written for you (file Si.input.convergence)

3 #Number of atoms, chemical species and atom types natom 2 # Number of atoms in the unit cell ntypat 1 # Number of types of atoms typat 1 1 # Type of atoms znucl 14.0 # Gives nuclear charge for each type of #Coordinates and cell variables rprim 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 acell 5.38 5.38 5.38 Angstrom xred -0.125 -0.125 -0.125 0.125 0.125 0.125 #PlaneWave cutoff and k-grid mesh integration ecutsm 0.5 # Energy cutoff smearing (Ha) nband 10 # Number of bands kptopt 1 # Kpoints option # 0 = read directly nkpt, kpt, kptnrm and wtk # 1 = rely on ngkpt or kptrlatt, as well as # on nshiftk and shiftk to set up # the k points. # Full symmetry taken into account. # 2 = 1, but only time reversal symmetry # is taken into account. # 3 = 1, but do not take into account # any symmetry # A negative value = rely on kptbounds, # and ndivk to set up a band structure # calculation along different lines ngkpt 3 3 3 # Number of grid points for k points # generation # This is a 3x3x3 FCC grid, # based on the primitive vectors # of the reciprocal space. # For a FCC real space lattice, # like the present one, # it actually corresponds to the # so-called 6x6x6 Monkhorst-Pack grid, # if the following shifts are used : nshiftk 4 # shiftk 0.5 0.5 0.5 # Shift for k points 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 #Exchange-correlation ixc 2 # Integer for exchange-correlation choice # 0 = No xc # 1 = LDA or LSD, Teter Pade parametrization # 2 = LDA, Perdew-Zunger-Ceperley-Alder # 3 = LDA, old Teter rational polynomial # parametrization, fit to Ceperley-Alder # data (no spin-polarization: no sp) # 4 = LDA, Wigner functional (no sp) # 5 = LDA, Hedin-Lundqvist functional (no sp) # 6 = LDA, "X-alpha" functional (no sp) # 7 = LDA or LSD, Perdew-Wang 92 functional # 8 = LDA or LSD, x-only part of the PW 92 # 9 = LDA or LSD, x- and RPA part of the PW92 # 10= GGA, Perdew-Burke-Ernzerhof

Diamond structure at the lattice constant that minimizes the energy in SIESTA

6 × 6 × 6 Monkhorst-Pack mesh

Ceperley-Alder (LDA) functional

Running the energy versus cutoff energy in ABINIT To check the equivalent cutoff energy in Abinit: 1.  We run the same system (same lattice vectors and internal coordinates) at the same level of approximations (same exchange and correlation functional, Monkhorst-Pack mesh etc.) at a given lattice constant. Here it has been written for you (file Si.input.convergence) ndtset

2. Change the cutoff energy for the plane waves $ more Si.files

$ more Si.files Si.input.convergence Si.input.convergence Si.out 3. Edit the .files file and select Sii Si.out Sio the input file and theSii pseudo t1x file (in XML format) Sio Si.xml t1x Si.xml $ abinit < Si.files > Si.log

4. Run the code

$ abinit < Si.files > Si.log &

1 2 3

ecut1 ecut2 ecut3 ecut4 ecut5 ecut6 ecut7 ecut8 ecut9 ecut10 ecut11 ecut12 ecut13 ecut14 &ecut15 ecut16 ecut17 ecut18 ecut19 ecut20 ecut21

21 4.00 5.00 6.00 7.00 8.00 9.00 10.00 11.00 12.00 13.00 14.00 15.00 16.00 17.00 18.00 19.00 20.00 25.00 30.00 35.00 40.00

G. Lee, Y.-H. Shin, and J. Y. Son, J. Am. Ceram. Soc., 95, 2773 (2012). E. Wiesendanger, Ferroelectrics, 6, 263 (1974). H. P. Sun, W. Tian, X. Q. Pan, J. H. Haeni, and D. G. Schlom, Appl. Phys. Lett., 84, 329

Running the energy versus cutoff energy in ABINIT: bulk Si (covalent semiconductor) Write the total energy as a function of the cutoff energy and edit the corresponding file that should look like this 3 $ grep Total Si.out > Si.abinit.convergence.dat $ more Si.abinit.convergence.dat # Cutoff energy (Ha) Total energy (eV) 4.0 -2.14194449515052E+02 5.0 -2.14976484971340E+02 6.0 -2.15411139789484E+02 7.0 -2.15628684537067E+02 8.0 -2.15720857649398E+02 9.0 -2.15750170913495E+02 10.0 -2.15763328854496E+02 11.0 -2.15779811817589E+02 12.0 -2.15799418159205E+02 13.0 -2.15818945848621E+02 14.0 -2.15839165600763E+02 15.0 -2.15857713873448E+02 16.0 -2.15872186183511E+02 17.0 -2.15882236230850E+02 18.0 -2.15889235327402E+02 19.0 -2.15893738919208E+02 20.0 -2.15896300628071E+02 25.0 -2.15898209214678E+02 30.0 -2.15900079984311E+02 35.0 -2.15901734504665E+02 40.0 -2.15902017447268E+02

$ gnuplot gnuplot> plot "Si.abinit.convergence.dat" u 1:2 w l gnuplot> set terminal postscript Terminal type set to ’postscript’ Options are ’landscape noenhanced defaultplex \ leveldefault monochrome colortext \ dashed dashlength 1.0 linewidth 1.0 butt \ palfuncparam 2000,0.003 \ "Helvetica" 14 ’ gnuplot> set output "Si.abinit.convergence.ps" gnuplot> replot -214

"Si.abinit.convergence.dat" u 1:2

-214.2 -214.4 -214.6

Equivalent PW cutoff a DZP basis set at 5.38 Å

-214.8 -215 -215.2 -215.4 -215.6 -215.8 -216

0

5

10

15

20

25

30

35

40

3 $ $ more Si.abinit.latcon.dat # Lattice constant (Ang) Total Energy (eV) (Cutoff = 13.005 Ha) $ $ more5.30 Si.abinit.latcon.dat-2.15791847116547E+02 # Lattice constant (Ang) Total Energy (eV) (Cutoff = 13.005 Ha) 5.32 -2.15803559248711E+02 5.30 -2.15791847116547E+02 5.34 -2.15811826556959E+02 5.32 -2.15803559248711E+02 5.36 -2.15816906510895E+02 5.38 -2.15819045481337E+02 5.34 -2.15811826556959E+02 5.40 -2.15818274869825E+02 5.36 -2.15816906510895E+02 5.42 -2.15814668146456E+02 5.38 -2.15819045481337E+02 5.44 -2.15808351816212E+02 5.40 -2.15818274869825E+02 5.46 -2.15799487297689E+02 5.42 -2.15814668146456E+02 5.48 -2.15788012963537E+02

Running the energy versus lattice constant curve in ABINIT 1. Same input as before but… 5.44

-2.15808351816212E+02

5.46 -2.15799487297689E+02 #input for bulk Si in the diamond structure. ecut

5.48

-2.15788012963537E+02

13.005

# Energy cutoff (Ha)

#input for bulk Si in the diamond structure. ndtset

10

ecut

13.005

acell1 acell2 ndtset acell3 acell4 acell1 acell5 acell6 acell2 acell7 acell3 acell8 acell4 acell9 acell5 acell10

5.30 5.32 10 5.34 5.36 5.30 5.38 5.40 5.32 5.42 5.34 5.44 5.36 5.46 5.38 5.48

acell6 5.40 acell7 5.42 $ more Si.files acell8 5.44 Si.input.latcon Si.out acell9 5.46 Sii acell10 5.48

5.30 5.32 5.34 5.36 5.30 5.38 5.40 5.32 5.42 5.34 5.44 5.36 5.46 5.38 5.48

5.40 5.42 5.44 5.46 5.48

… setting the plane wave cutoff to the equivalent one to a DZP basis set

# Energy cutoff (Ha)

5.30 Angstrom 5.32 Angstrom 5.34 Angstrom 5.36 Angstrom 5.30 Angstrom Angstrom 5.38 5.40 5.32 Angstrom Angstrom 5.42 5.34 Angstrom Angstrom 5.44 5.36 Angstrom Angstrom 5.46 Angstrom 5.38 Angstrom 5.48 Angstrom

5.40 5.42 5.44 5.46 5.48

… and changing the lattice constant embracing the minimum

Angstrom Angstrom Angstrom Angstrom Angstrom

2. Change in the .files the name of the input file

$ more Si.files Si.input.convergence Si.out Si.out $ gnuplot Sii Sii gnuplot> plot "Si.abinit.convergence.dat" u 1:2 w l Sio gnuplot> set terminal postscriptSio t1x Terminal type set to ’postscript’ t1x Options are ’landscape noenhanced defaultplex \ Si.xml Si.xml leveldefault monochrome colortext \ Sio $t1x more Si.files Si.xml Si.input.latcon

dashed dashlength 1.0 linewidth 1.0 butt \ $ gnuplot palfuncparam 2000,0.003 \ gnuplot> plot "Si.abinit.convergence.dat" 1:2 w l $ abinit < uSi.files "Helvetica" 14 ’ gnuplot> set terminal postscript gnuplot> set output "Si.abinit.convergence.ps" Terminal type set to ’postscript’ gnuplot> replot Options are ’landscape noenhanced defaultplex \ leveldefault monochrome colortext \

3. Run the code

> Si.log &

Running the energy versus lattice constant curve in ABINIT Write the total energy as a function of the lattice constant and edit the corresponding file that should look like this 3 $ grep "Total" Si.out > Si.abinit.latcon.dat $ more Si.abinit.latcon.dat # Lattice constant (Ang) Total Energy (eV) (Cutoff = 13.005 Ha) 5.30 -2.15791847116547E+02 5.32 -2.15803559248711E+02 5.34 -2.15811826556959E+02 5.36 -2.15816906510895E+02 5.38 -2.15819045481337E+02 5.40 -2.15818274869825E+02 5.42 -2.15814668146456E+02 5.44 -2.15808351816212E+02 5.46 -2.15799487297689E+02 5.48 -2.15788012963537E+02 #input for bulk Si in the diamond structure. ecut

13.005

ndtset

10

acell1 acell2 acell3 acell4 acell5 acell6 acell7 acell8 acell9 acell10

5.30 5.32 5.34 5.36 5.38 5.40 5.42 5.44 5.46 5.48

# Energy cutoff (Ha)

5.30 5.32 5.34 5.36 5.38 5.40 5.42 5.44 5.46 5.48

5.30 5.32 5.34 5.36 5.38 5.40 5.42 5.44 5.46 5.48

Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom

Comparing the pseudopotential in SIESTA and ABINIT: bulk Si (covalent semiconductor) To be totally sure that we have run SIESTA and ABINIT with the same peudopotential operator, i.e. with the same decomposition in local part and Kleinman-Bylander projectors: I.

2

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIELD.

1. Edit one of the output files in SIESTA and search for the following lines: 1 Cite .

KBgen: Kleinman-Bylander l= 0 rc= 1.936440 l= 1 rc= 1.936440 l= 2 rc= 1.936440 l= 3 rc= 1.936440

projectors: el= -0.796617 4.661340 kbcos= 0.299549 I. ENERGYEkb= FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIELD. el= -0.307040 Ekb= 1.494238 kbcos= 0.298299 Cite1 . 0.008863 el= Ekb= -2.809034 kbcos= -0.000566 el= 0.013111 Ekb= -0.959387 kbcos= -0.000004 KBgen: Kleinman-Bylander projectors: l= 0

2. Edit the

rc=

1.936440

el= -0.796617

Ekb=

4.661340

kbcos=

0.299549

2

EKB= 2.3306701075631429 l= 1 0.74711916272986667 -1.4045230577882524 rc= 1.936440 el= -0.307040 Ekb= 1.494238 kbcos= 0.298299-0.4796986487264058 l= 2 rc= 1.936440 el= 0.008863 Ekb= -2.809034 kbcos= -0.000566 -3.9422142725424 COSKB= 0.30046489563475459 0.29921145875238120 -5.9527141624513349E-004 log file in ABINIT and l=search for theel=following lines: 3 rc= 1.936440 0.013111 Ekb= -0.959387 kbcos= -0.000004 pspatm: epsatm= 9.82921496 --- l ekb(1:nproj) EKB=--> 2.3306701075631429 The Kleinman-Bylander 0 2.330670 0.74711916272986667 -1.4045230577882524 1 0.747119 energies and cosines should 2 -1.404523COSKB= -0.47969864872640589 0.30046489563475459 be the same upto numerical 3 -0.479699 0.29921145875238120 -5.9527141624513349E-004 -3.9422142725424715E-006 pspatm: epsatm= 9.82921496 --- l ekb(1:nproj) --> 0 2.330670 1 0.747119 2 -1.404523 3 -0.479699

roundoff errors

Note: In SIESTA they are written in Ry and in ABINIT they are in Ha.

Comparing the energy versus lattice constant in SIESTA and ABINIT: bulk Si (covalent semiconductor) 3 $ gnuplot gnuplot> plot "Si.abinit.latcon.dat" u 1:2 w l, "../Runsiesta/Si.siesta.latcon.dat" u 1:2 w l gnuplot> set terminal postscript colorTerminal type set to ’postscript’ Options are ’landscape noenhanced defaultplex \ leveldefault color colortext \ dashed dashlength 1.0 linewidth 1.0 butt \ palfuncparam 2000,0.003 \ "Helvetica" 14 ’ gnuplot> set output "Si.compar.latcon.ps"gnuplot> replot #input for bulk Si in the diamond structure. -215.78

13.005

ndtset

10

acell1 acell2 acell3 acell4 acell5 acell6 acell7 acell8 acell9 acell10

5.30 5.32 5.34 5.36 5.38 5.40 5.42 5.44 5.46 5.48

$ more Si.files Si.input.latcon Si.out Sii Sio t1x Si.xml

# Energy cutoff (Ha)

-215.79

5.30 5.30 5.32 5.32 5.34 5.34 5.36 5.36 5.38 -215.80 5.38 5.40 5.40 5.42 5.42 5.44 5.44 5.46 5.46 5.48 -215.81 5.48

Free Energy (eV)

ecut

-215.82

Siesta (DZP) Abinit (13.005 Ha)

Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom Angstrom

5.30

5.35 5.40 5.45 Lattice constant (Ang)

5.50

Comparing the pseudopotential in SIESTA and ABINIT: bulk Al (sp metal) To be totally sure that we have run SIESTA and ABINIT with the same peudopotential operator, i.e. with the same decomposition in local part and Kleinman-Bylander projectors: I.

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC

Cite1 . # Ekb (in Ha): # SIESTA # l= 0 1.556933 # l= 1 0.439408 # l= 2 -0.891462 # l= 3 -0.308287 # # Coskb # SIESTA # l= 0 0.285004 # l= 1 0.270194 # l= 2 -0.001042 # l= 3 -0.000009

ABINIT 1.556933 0.439408 -0.891467 -0.308292

ABINIT 0.285875 0.271020 -0.001096 -0.000009

ElectronicTemperature 0.02 Ry occopt

3

tsmear

0.01

# Occupation Option # Fermi-Dirac smearing (finite-temperature metal) # Temperature of smearing (in Ha)

Comparing pseudopotential in I. the ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIELD. SIESTA and ABINIT: bulk Al (sp metal) Cite . I.

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRI

Cite1 .

1

# Ekb (in Ha): # SIESTA ABINIT # l= 0 1.556933 1.556933 Ekb (in Ha): # l= 1 0.439408 0.439408 SIESTA ABINIT # l= 2 -0.891462 -0.891467 -0.308292 l= 0 1.556933 1.556933 # l= 3 -0.308287 # l= 1 0.439408 0.439408 # Coskb I. ENERGY FUNCTIONAL SIESTA ABINIT l= 2 -0.891462 -0.891467 # # l= 0 0.285004 0.285875 l= 3 -0.308287 -0.308292 # l= 1 0.270194 IESTA BINIT 0.271020 Cite1 . # l= 2 -0.001042 -0.001096 # l= 3 -0.000009 Coskb I. ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE-0.000009 AN ELECTRIC FIELD.

For #the case of metallic system, besides the k-point sampling we have to pay particular attention to the occupation option # # # # # S # # # SIESTA ABINIT 1 # . l= 0 0.285004 0.285875 Cite Default: Fermi-Dirac # l= 1 0.270194 0.271020 # Ekb (in Ha): # l= 2 -0.001042 -0.001096 # # SIESTA ABINIT l= 3 -0.000009 -0.000009

FOR A DIELE

A

2

# Ekb (in Ha): # SIESTA ABINIT occopt #3 Option l= 0# Occupation 1.432475 1.432472 # Fermi-Dirac smearing (finite-temperature metal) # l= 1# Temperature 0.827346 0.827336 tsmear 0.01 of smearing (in Ha) # l= 2 -3.174442 -3.174442 # l= 3 -1.038935 -1.039329 # l= 0 1.432475 1.432472 # ElectronicTemperature l= 1 0.827346 0.827336 # 0.02 Ry # l= 2 -3.174442 -3.174442 # Coskb # occopt l= 3 -1.038935 # SIESTA ABINIT 3 # -1.039329 Occupation Option # l= 0 0.304116 metal) 0.305046 # Fermi-Dirac smearing # (finite-temperature # Coskb Also, as explained in # the l= Tutorial 1 0.208262 tsmear 0.01 # Temperature of smearing (in Ha) 0.207625 # SIESTA ABINIT l= 2 -0.703574 -0.705726 “Convergence of electronic and structural#properties of a metal with respect to # l= 0 0.304116 0.305046 # l= 3 -0.000012 -0.000012 # l= 1 0.207625 0.208262 # #

l= 2 -0.703574 we should look l= 3 -0.000012

ElectronicTemperature 0.02 Ry

the k-point sampling: bulk Al” Energy and notFreeEng to the *.out Kohn-Sham energy grep > Al.siesta.latcon.dat

at-0.705726 the Free -0.000012

grep FreeEng *.out > Al.siesta.latcon.dat

grep Total Al.out > Al.abinit.latcon.dat

grep Total Al.out > Al.abinit.latcon.dat

occopt

3

# Occupation Option # Fermi-Dirac smearing (finite-te

Running the energy versus cutoff energy in ABINIT: bulk Al (a sp metal) Lattice constant 3.97 Å -56,9

Free Energy (eV)

-57,0

-57,1

DZP

-57,2

0

5

10 20 15 Cutoff energy (Ha)

25

30

Comparing the energy versus lattice constant in SIESTA and ABINIT: bulk Al (sp metal) Basis set of SIESTA: DZP optimized with a pressure of 0.001 GPa at the theoretical lattice constant of 3.97 Å) Plane wave cutoff in ABINIT: 8.97 Ha -57.14 Siesta (DZP) Abinit (8.97 Ha)

Free Energy (eV)

-57.15

-57.16

-57.17

-57.18

3.85

3.90

4.00 3.95 4.05 Lattice constant (Ang)

4.10

Comparing the pseudopotential in SIESTA and ABINIT: bulk Au (a noble metal) To be totally sure that we have run SIESTA and ABINIT with the same peudopotential operator, i.e. with the same decomposition in local part and Kleinman-Bylander projectors: I.

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC F

Cite1 . # Ekb (in Ha): # SIESTA # l= 0 1.432475 # l= 1 0.827346 # l= 2 -3.174442 # l= 3 -1.038935 # # Coskb # SIESTA # l= 0 0.304116 # l= 1 0.207625 # l= 2 -0.703574 # l= 3 -0.000012

ABINIT 1.432472 0.827336 -3.174442 -1.039329

ABINIT 0.305046 0.208262 -0.705726 -0.000012

grep FreeEng *.out > Al.siesta.latcon.dat grep Total Al.out > Al.abinit.latcon.dat occopt

3

tsmear

0.01

# Occupation Option # Fermi-Dirac smearing (finite-temperature metal) # Temperature of smearing (in Ha)

Running the energy versus cutoff energy in ABINIT: bulk Au (a noble metal) Lattice constant 4.08 Å -903.6 -903.8

Free Energy (eV)

-904.0 -904.2 -904.4

DZP

-904.6 -904.8 -905.0

20

30 40 Cutoff energy (Ha)

50

60

Comparing the energy versus lattice constant in SIESTA and ABINIT: bulk Au (a noble metal) Basis set of SIESTA: DZP optimized with a pressure of 0.02 GPa at the theoretical lattice constant of 4.08 Å Plane wave cutoff in ABINIT: 17.432 Ha -904.66 Siesta (DZP) Abinit (17.432 Ha)

Free Energy (eV)

-904.67

-904.68

-904.69

-904.70

4.00

4.10 4.05 Lattice constant (Ang)

4.15

Comparing the pseudopotential in SIESTA and ABINIT: bulk Fe (a magnetic transition metal) To be totally sure that we have run SIESTA and ABINIT with the same peudopotential operator, i.e. with the same decomposition in local part and Kleinman-Bylander projectors: I.

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIE

Cite1 . # Ekb (in Ha): # SIESTA # l= 0 1.693954 # l= 1 0.682661 # l= 2 -7.130128 # l= 3 -0.827854 # # Coskb # SIESTA # l= 0 0.237303 # l= 1 0.207724 # l= 2 -0.711799 # l= 3 -0.000008 XC.functional XC.authors SpinPolarized

ABINIT 1.693938 0.682647 -7.130129 -0.828135

ABINIT 0.238031 0.208363 -0.713975 -0.000008 GGA PBE .true.

nsppol spinat ixc

2 0.0 0.0 1.0 11

occopt

3

tsmear

0.01

# Number of spin polarizations # Spin for atoms # Integer for exchange-correlation choice

# Occupation Option # Fermi-Dirac smearing (finite-temperature metal) # Temperature of smearing (in Ha)

Comparing the pseudopotential in SIESTA and ABINIT: bulk Fe (a magnetic transition metal) I.

2

ENERGY FUNCTIONAL FOR A DIELECTRIC INSIDE AN ELECTRIC FIELD. ENERGY FUNCTIONAL FOR A I.DIELECTRIC INSIDE AN ELECTRIC FIELD.

For the case of metallic system,1 besides the k-point sampling we have to pay Cite . particular attention to the occupation option. Cite1 . # Ekb (in Ha): Now, besides: ABINIT # SIESTA # Ekb (in Ha): 2.129660 - # TheSIESTA system is spin polarized ABINIT # l= 0 2.129660 # l= 1 1.425393 1.425393 l= 0use2.129660 2.129660 # l= 2 -6.283666 - # We a GGA functional -6.283666 l= 1 1.425393 1.425393 # l= 3 -0.825350 -0.825350 - ## We include non-linear partial core corrections in the pseudo l= 2 -6.283666 -6.283666 # # l= 3 -0.825350 # # Coskb # SIESTA # l= 0 S0.260207 IESTA # l= 1 0.174652 # l= 2 -0.683298 # l= 3 -0.000005 XC.functional XC.authors SpinPolarized

-0.825350 # Coskb # # # ABINIT # 0.261003 #

0.175186 -0.685387 XC.functional XC.authors -0.000005 SpinPolarized GGA PBE .true.

nsppol spinat ixc

2 0.0 0.0 1.0 11

occopt

3

tsmear

0.01

l= l= l= l=

SIESTA 0 0.260207 1 0.174652 2 -0.683298 3 -0.000005

ABINIT 0.261003 0.175186 -0.685387 -0.000005ABINIT GGA PBE .true.

nsppol spinat ixc

2 0.0 0.0 1.0 11

occopt

3

# Number of spin polarizations # Spin for atoms # Integer for exchange-correlation choice

# Occupation Option

# Number of spin# polarizations Fermi-Dirac smearing (finite-temperature metal) #tsmear Spin for 0.01 atoms# Temperature of smearing (in Ha) # Integer for exchange-correlation choice

# Occupation Option # Fermi-Dirac smearing (finite-temperature metal) # Temperature of smearing (in Ha)

Running the energy versus cutoff energy in ABINIT: bulk Fe (a magnetic transition metal) Lattice constant 2.87 Å -765.0

Free Energy (eV)

-770.0

-775.0

-780.0

20

DZP 30

40 50 Cutoff energy (Ha)

60

Comparing the energy versus lattice constant in SIESTA and ABINIT: bulk Fe(magnetic transition metal) Basis set of SIESTA: DZP optimized without pressure at the experimental lattice constant of 2.87 Å Plane wave cutoff in ABINIT: 34.82 Ha

-782.48

Siesta (DZP) Abinit (34.82 Ha)

Free Energy (eV)

-782.50

-782.52

-782.54

-782.56

-782.58

-782.60

2.75

2.80

2.90 2.85 2.95 Lattice constant (Ang)

3,00