Running tight-binding calculations
How to run tight-binding calculations using the pre-fit tight-binding coefficients. Note, only elemental and binary systems are currently supported.
Running a julia function for the first time will compile the function. Future runs will be much faster.
Create a crystal object
Use makecrys to create a crystal from lattice vectors, atomic positions, and atom types:
using ThreeBodyTB
A = [2.1 2.1 0.0;2.1 0.0 2.1;0.0 2.1 2.1];
pos = [0.0 0.0 0.0];
types = ["Al"];
fcc_al = makecrys(A, pos, types)Units: Å A1= 2.10000 2.10000 0.00000 A2= 2.10000 0.00000 2.10000 A3= 0.00000 2.10000 2.10000 Al 0.00000 0.00000 0.00000
Current default units are Angstrom and eV. You can change the global units to atomic units with set_units(both="atomic") if you prefer.
Alternatively, you can read the positions from a simple POSCAR or Quantum Espresso input file.
rbcl = makecrys("../src/POSCAR_rbcl")Units: Å A1= 3.90963 -0.00000 0.00000 A2= -0.00000 3.90963 0.00000 A3= 0.00000 -0.00000 3.90963 Rb 0.00000 0.00000 0.00000 Cl 0.50000 0.50000 0.50000
Do a self-consistent calculation.
Gets the energy and charge density with scf_energy:
alp = makecrys("../src/POSCAR_alp")
energy, tbc_alp = scf_energy(alp);
println("The energy is $energy eV")Begin scf_energy------------- Load TB parameters from file prepare atoms ["Al", "P"] add_to_database Set([:P]) found /juliateam/.julia/packages/ThreeBodyTB/djQTk/src/../dats/pbesol/v1.2/els/coef.el.2bdy.P.xml.gz version 2 added to cache (:P, :P) found /juliateam/.julia/packages/ThreeBodyTB/djQTk/src/../dats/pbesol/v1.2/els/coef.el.3bdy.P.xml.gz version 2 added to cache (:P, :P, :P) add_to_database Set([:P, :Al]) found /juliateam/.julia/packages/ThreeBodyTB/djQTk/src/../dats/pbesol/v1.2/binary/coef.el.2bdy.Al.P.xml.gz version 2 added to cache (:P, :Al) twobody found /juliateam/.julia/packages/ThreeBodyTB/djQTk/src/../dats/pbesol/v1.2/binary/coef.el.3bdy.Al.P.xml.gz version 2 added to cache (:P, :Al) threebody add_to_database Set([:Al]) found /juliateam/.julia/packages/ThreeBodyTB/djQTk/src/../dats/pbesol/v1.2/els/coef.el.2bdy.Al.xml.gz version 2 added to cache (:Al, :Al) found /juliateam/.julia/packages/ThreeBodyTB/djQTk/src/../dats/pbesol/v1.2/els/coef.el.3bdy.Al.xml.gz version 2 added to cache (:Al, :Al, :Al) construct ----- Construct tight-binding model from crystal structure distances CHECK FRONTIER - everything fine 2body 3body onsite ----- 9.168010 seconds (23.29 M allocations: 1.201 GiB, 8.54% gc time, 94.08% compilation time) ------ Do SCF Mixing mode: pulay Get initial guess from tbc DQ: [0.0, -0.0] Parameters: smearing = 0.01 conv_thr = 1.0e-5, iters = 75, mix = 0.8, grid = missing START SCF ---------------- SCF CALC 0001 energy -10.42350021 SCF CALC 0002 energy -10.39235446 en_diff: 3.114575E-02 dq_diff: 9.077615E-01 SCF CALC 0003 energy -10.37239696 en_diff: 1.995750E-02 dq_diff: 7.231044E-01 SCF CALC 0004 energy -10.28933097 en_diff: 8.306599E-02 dq_diff: 1.452555E-01 SCF CALC 0005 energy -10.26727897 en_diff: 2.205201E-02 dq_diff: 2.901010E-02 SCF CALC 0006 energy -10.26272388 en_diff: 4.555090E-03 dq_diff: 5.804629E-03 SCF CALC 0007 energy -10.26180215 en_diff: 9.217278E-04 dq_diff: 1.160782E-03 SCF CALC 0008 energy -10.26161767 en_diff: 1.844756E-04 dq_diff: 2.321653E-04 YES convergence in 8 iters, energy -10.261617674454563 eV dq = [-1.163, 1.163] END SCF ------------------ scf_energy success, done The energy is -10.261617674454563 eV
This returns the (non-magnetic) atomization energy, and a tight-binding object with the TB matrix elements and SCF electron density calculated for post-processing.
Plot the band structure.
Using the tb_crys tight-binding object tbc_alp from above. Note: SCF must be done first.
plot_bandstr(tbc_alp, do_display=false);plot_bandstr align vbm color = blue markersize = 0 no display
Use do_display=true (the default) to produce an interactive plot. Here do_display is set to false because we are saving a static figure with savefig for the docs.
The default plot_bandstr just picks some random kpoints, but you can add your own kpath. We can also project onto the s orbital of Al.
kpath=[0.0 0.0 0.0; 0.5 0.5 0.5; 0.0 0.5 0.5];
knames=["Γ", "X", "V"];
plot_bandstr(tbc_alp, kpath=kpath, names=knames, npts=100, proj_orbs=[:s], proj_types=["Al"], do_display=false);proj_inds [1] plot_bandstr align vbm no display
You can also plot the DOS using dos
dos(tbc_alp, do_display=false);grid [14, 14, 14]
Projection type: atomic
PROJ
("P", [5, 6, 7, 8], 4)
("Al", [1, 2, 3, 4], 4)
Project onto orbitals instead with proj_type=:orbs
Calculate force / stress
energy, force, stress, tbc = scf_energy_force_stress(tbc_alp);
println("energy $energy")
println()
println("Forces")
show(stdout, "text/plain", force)
println()
println("Stress")
show(stdout, "text/plain", stress)prepare atoms ["Al", "P"]
Calculate Force, Stress (no scf)
get_energy_force_stress_fft
21.477943 seconds (57.72 M allocations: 3.094 GiB, 9.72% gc time, 90.21% compilation time)
done
----
energy -10.261580756479562
Forces
2×3 Matrix{Float64}:
-1.37139e-11 -1.37135e-11 -1.37137e-11
0.0 0.0 0.0
Stress
3×3 Matrix{Float64}:
0.138528 0.0 0.0
0.0 0.138528 0.0
0.0 0.0 0.138528Can also be called directly on a new crystal structure instead of a tb_crys object.
Relax structure
Using relax_structure
crys_new, tbc_updated, energy, force, stress = relax_structure(alp);
println("Energy new $energy")
println()
println("Force")
show(stdout, "text/plain", force)
println()
println("Stress")
show(stdout, "text/plain", stress)prepare atoms ["Al", "P"]
relax_structure conv_thr 0.002 energy_conv_thr (Ryd) 0.0002
-----
Construct tight-binding model from crystal structure
distances
CHECK FRONTIER - everything fine
2body
3body
onsite
-----
------
Do SCF
Mixing mode: pulay
Get initial guess from tbc
DQ: [0.0, -0.0]
Parameters:
smearing = 0.01 conv_thr = 1.0e-7, iters = 75, mix = 0.8, grid = missing
START SCF ----------------
SCF CALC 0001 energy -10.42350021
SCF CALC 0002 energy -10.39235446 en_diff: 3.114575E-02 dq_diff: 9.077615E-01
SCF CALC 0003 energy -10.37239696 en_diff: 1.995750E-02 dq_diff: 7.231044E-01
SCF CALC 0004 energy -10.28933097 en_diff: 8.306599E-02 dq_diff: 1.452555E-01
SCF CALC 0005 energy -10.26727897 en_diff: 2.205201E-02 dq_diff: 2.901010E-02
SCF CALC 0006 energy -10.26272388 en_diff: 4.555090E-03 dq_diff: 5.804629E-03
SCF CALC 0007 energy -10.26180215 en_diff: 9.217278E-04 dq_diff: 1.160782E-03
SCF CALC 0008 energy -10.26161767 en_diff: 1.844756E-04 dq_diff: 2.321653E-04
SCF CALC 0009 energy -10.26158076 en_diff: 3.691798E-05 dq_diff: 4.643256E-05
SCF CALC 0010 energy -10.26157337 en_diff: 7.383469E-06 dq_diff: 9.286542E-06
SCF CALC 0011 energy -10.26157190 en_diff: 1.476750E-06 dq_diff: 1.857307E-06
YES convergence in 11 iters, energy -10.261571896261032 eV dq = [-1.163, 1.163]
END SCF ------------------
starting vec
[0.0, 0.0, 0.0, 0.25, 0.25, 0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Conj Grad START
START SCF ----------------
SCF CALC 0001 energy -10.26157160
SCF CALC 0002 energy -10.26157157 en_diff: 3.244721E-08 dq_diff: 2.203338E-07
YES convergence in 2 iters, energy -10.261571568465202 eV dq = [-1.163, 1.163]
END SCF ------------------
get_energy_force_stress_fft
2.088345 seconds (8.98 M allocations: 704.666 MiB, 12.76% gc time)
FCALL 1 en: -0.754211325193509 (Ryd) fsum: 9.238894857355681e-13 ssum: 0.0026132419734498577 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
A1= +2.60000 +2.60000 +0.00000 | +0.13853 +0.00000 +0.00000
A2= +2.60000 +0.00000 +2.60000 | +0.00000 +0.13853 +0.00000
A3= +0.00000 +2.60000 +2.60000 | +0.00000 +0.00000 +0.13853
Al +0.00000 +0.00000 +0.00000 | -0.00000 -0.00000 -0.00000
P +0.25000 +0.25000 +0.25000 | +0.00000 +0.00000 +0.00000
START SCF ----------------
SCF CALC 0001 energy -10.77406918
SCF CALC 0002 energy -10.73232516 en_diff: 4.174402E-02 dq_diff: 1.948531E-01
SCF CALC 0003 energy -10.72003551 en_diff: 1.228965E-02 dq_diff: 1.505436E-01
SCF CALC 0004 energy -10.68614796 en_diff: 3.388756E-02 dq_diff: 3.139261E-02
SCF CALC 0005 energy -10.67841771 en_diff: 7.730242E-03 dq_diff: 6.172964E-03
SCF CALC 0006 energy -10.67689772 en_diff: 1.519994E-03 dq_diff: 1.242025E-03
SCF CALC 0007 energy -10.67658984 en_diff: 3.078795E-04 dq_diff: 2.478466E-04
SCF CALC 0008 energy -10.67652847 en_diff: 6.136700E-05 dq_diff: 4.961009E-05
SCF CALC 0009 energy -10.67651618 en_diff: 1.229194E-05 dq_diff: 9.918999E-06
YES convergence in 9 iters, energy -10.676516181896682 eV dq = [-0.949, 0.949]
END SCF ------------------
get_energy_force_stress_fft
1.562539 seconds (6.98 M allocations: 571.404 MiB, 12.73% gc time, 0.95% compilation time)
FCALL 2 en: -0.7847091806818252 (Ryd) fsum: 9.227596667164204e-9 ssum: 0.00019373959579528408 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
A1= +2.80318 +2.80318 +0.00000 | -0.01027 -0.00000 -0.00000
A2= +2.80318 +0.00000 +2.80318 | -0.00000 -0.01027 -0.00000
A3= +0.00000 +2.80318 +2.80318 | -0.00000 -0.00000 -0.01027
Al +0.00000 +0.00000 +0.00000 | -0.00000 -0.00000 -0.00000
P +0.25000 +0.25000 +0.25000 | +0.00000 +0.00000 +0.00000
MY CG quadratic linesearch iter 1 fn val: -0.7847091806818252 | rms grad: 0.057596919388920365 within_stepsize: false , old step: 0.3275188899474979 , new step_size: 0.10565125482177351
MY CG quadratic linesearch iter 2 fn val: -0.7847091806818252 | rms grad: 0.057596919388920365 within_stepsize: true , old step: 0.10565125482177351 , new step_size: 0.15847688223266027
START SCF ----------------
SCF CALC 0001 energy -10.68005595
SCF CALC 0002 energy -10.68005582 en_diff: 1.305953E-07 dq_diff: 4.971398E-07
YES convergence in 2 iters, energy -10.680055819817516 eV dq = [-0.959, 0.959]
END SCF ------------------
get_energy_force_stress_fft
1.512667 seconds (6.96 M allocations: 570.267 MiB, 12.19% gc time)
FCALL 3 en: -0.7849693391759861 (Ryd) fsum: 1.2197893514326893e-8 ssum: 0.00011613992097199466 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
A1= +2.79405 +2.79405 -0.00000 | -0.00616 -0.00000 -0.00000
A2= +2.79405 -0.00000 +2.79405 | -0.00000 -0.00616 -0.00000
A3= -0.00000 +2.79405 +2.79405 | -0.00000 -0.00000 -0.00616
Al -0.00000 -0.00000 -0.00000 | +0.00000 +0.00000 +0.00000
P +0.25000 +0.25000 +0.25000 | -0.00000 -0.00000 -0.00000
MY CG quadratic linesearch iter 3 fn val: -0.7849693391759861 | rms grad: 0.0341908455769504 within_stepsize: true , old step: 0.15847688223266027 , new step_size: 0.2377153233489904
START SCF ----------------
SCF CALC 0001 energy -10.68164733
SCF CALC 0002 energy -10.68164722 en_diff: 1.148513E-07 dq_diff: 4.443982E-07
YES convergence in 2 iters, energy -10.681647215759517 eV dq = [-0.968, 0.968]
END SCF ------------------
get_energy_force_stress_fft
1.550751 seconds (6.96 M allocations: 570.267 MiB, 13.28% gc time)
FCALL 4 en: -0.7850863045778562 (Ryd) fsum: 5.897271768735629e-8 ssum: 4.108858643914421e-5 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
A1= +2.78591 +2.78591 -0.00000 | -0.00218 -0.00000 -0.00000
A2= +2.78591 -0.00000 +2.78591 | -0.00000 -0.00218 -0.00000
A3= -0.00000 +2.78591 +2.78591 | -0.00000 -0.00000 -0.00218
Al +0.00000 +0.00000 +0.00000 | -0.00000 -0.00000 -0.00000
P +0.25000 +0.25000 +0.25000 | +0.00000 +0.00000 +0.00000
MY CG quadratic linesearch iter 4 fn val: -0.7850863045778562 | rms grad: 0.011990883456665921 within_stepsize: true , old step: 0.2377153233489904 , new step_size: 0.354933648697078
START SCF ----------------
SCF CALC 0001 energy -10.68185680
SCF CALC 0002 energy -10.68186058 en_diff: 3.772649E-06 dq_diff: 1.426638E-05
YES convergence in 2 iters, energy -10.681860575934548 eV dq = [-0.973, 0.973]
END SCF ------------------
get_energy_force_stress_fft
1.520687 seconds (6.96 M allocations: 570.273 MiB, 12.46% gc time)
FCALL 5 en: -0.7851019862557825 (Ryd) fsum: 3.663487348203222e-7 ssum: 2.4639528839960134e-7 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
A1= +2.78166 +2.78166 -0.00000 | -0.00001 -0.00000 -0.00000
A2= +2.78166 -0.00000 +2.78166 | -0.00000 -0.00001 -0.00000
A3= -0.00000 +2.78166 +2.78166 | -0.00000 -0.00000 -0.00001
Al -0.00000 -0.00000 -0.00000 | +0.00000 +0.00000 +0.00000
P +0.25000 +0.25000 +0.25000 | -0.00000 -0.00000 -0.00000
yes conv sum_grad 7.168101918451657e-5 fn_diff_1 0.00011696540187011717 fn_diff_2 1.5681677926315807e-5
Relax done
---------------------------------
Final Energy -10.681860575934548
A1= +2.78166 +2.78166 -0.00000 | -0.00120 -0.00000 -0.00000
A2= +2.78166 -0.00000 +2.78166 | -0.00000 -0.00120 -0.00000
A3= -0.00000 +2.78166 +2.78166 | -0.00000 -0.00000 -0.00120
Al -0.00000 -0.00000 -0.00000 | +0.00010 +0.00010 +0.00010
P +0.25000 +0.25000 +0.25000 | -0.00010 -0.00010 -0.00010
Energy new -10.681860575934548
Force
2×3 Matrix{Float64}:
3.84538e-6 3.84538e-6 3.84538e-6
-3.84538e-6 -3.84538e-6 -3.84538e-6
Stress
3×3 Matrix{Float64}:
-1.30612e-5 -5.37793e-8 -5.37793e-8
-5.37793e-8 -1.30612e-5 -5.37793e-8
-5.37793e-8 -5.37793e-8 -1.30612e-5Energy is lower, stress is near zero, forces are zero by symmetry in Zinc Blende structure.
Force/Stress defaults are eV/Ang and eV/Ang^3.