This repository has been archived by the owner on Mar 5, 2018. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathTightBinding.jl
1743 lines (1448 loc) · 58.5 KB
/
TightBinding.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
`module TightBinding`
### Summary
Implements some functionality for tight binding models
### Notes
* at the moment we assume availability of ASE, MatSciPy
"""
module TightBinding
using AtomsInterface
importall AtomsInterface
using Potentials, ASE, MatSciPy, Prototypes, SparseTools
import MatSciPy.potential_energy
import Potentials.evaluate, Potentials.evaluate_d
export AbstractTBModel, SimpleFunction
export TBModel, FermiDiracSmearing, potential_energy, forces, evaluate,
potential_energy_d, site_energy, band_structure_all, band_structure_near_eF
abstract AbstractTBModel <: AbstractCalculator
abstract SmearingFunction <: SimpleFunction
"""`TBModel`: basic non-self consistent tight binding calculator. This
type admits TB models where the Hamiltonian is of the form
H_ij = f_hop(r_ij) if i ≠ j
H_ii = f_os( {r_ij}_j )
i.e. the hopping terms is a pair potential while the on-site terms are
more general; this is consistent in particular with the NRL TB model.
This model can only descibe a single species of atoms.
"""
type TBModel{P_os, P_hop, P_ol, P_p} <: AbstractTBModel
# Hamiltonian parameters
onsite::P_os
hop::P_hop
overlap::P_ol
# repulsive potential
pair::P_p
# HJ: add a parameter Rcut
# since the functions "cutoff" in Potentials.jl and NRLTB.jl may conflict
Rcut::Float64
# remaining model parameters
smearing::SmearingFunction
norbitals::Int
# WHERE DOES THIS GO?
fixed_eF::Bool
eF::Float64
# beta::Float64
# eF::Float64
# k-point sampling information:
# 0 = open boundary
# 1 = Gamma point
nkpoints::Tuple{Int, Int, Int}
# internals
hfd::Float64 # step used for finite-difference approximations
needupdate::Bool # tells whether hamiltonian and spectrum are up-to-date
arrays::Dict{Any, Any} # storage for various
end
typealias TightBindingModel TBModel
TBModel(;onsite = ZeroSitePotential(),
hop = ZeroPairPotential(),
overlap = ZeroPairPotential(),
pair = ZeroPairPotential(),
Rcut = 0.0,
smearing = ZeroTemperature(),
norbitals = 0,
fixed_eF = true,
eF = 0.0,
nkpoints = (0,0,0),
hfd = 0.0,
needupdate = true,
arrays = Dict()) =
TBModel(onsite, hop, overlap, pair, Rcut, smearing, norbitals,
fixed_eF, eF, nkpoints, hfd, needupdate, arrays)
############################################################
#### UTILITY FUNCTIONS
# import Potentials.cutoff
cutoff(tbm::TBModel) = tbm.Rcut
# HJ: not sure this returns right Rcut for NRL ----------------------------------
# max(cutoff(tbm.hop), cutoff(tbm.onsite), cutoff(tbm.pair))
# -------------------------------------------- ----------------------------------
"""`indexblock`:
a little auxiliary function to compute indices for several orbitals
"""
# indexblock(n::Vector, tbm::TBModel) =
# (n-1) * tbm.norbitals .+ [1:tbm.norbitals;]'
indexblock(n::Integer, tbm::TBModel) =
Int[(n-1) * tbm.norbitals + j for j = 1:tbm.norbitals]
"""`FermiDiracSmearing`:
f(e) = ( 1 + exp( beta (e - eF) ) )^{-1}
"""
type FermiDiracSmearing <: SmearingFunction
beta
eF
end
FermiDiracSmearing(beta;eF=0.0) = FermiDiracSmearing(beta, eF)
# FD distribution and its derivative. both are vectorised implementations
fermi_dirac(eF, beta, epsn) =
2.0 ./ ( 1.0 + exp((epsn - eF) * beta) )
fermi_dirac_d(eF, beta, epsn) =
- fermi_dirac(eF, beta, epsn).^2 .* exp((epsn - eF) * beta) * beta / 2.0
fermi_dirac_d2(eF, beta, epsn) =
fermi_dirac(eF, beta, epsn).^3 .* exp((epsn - eF) * 2.0 * beta) * beta^2 / 2.0 -
fermi_dirac(eF, beta, epsn).^2 .* exp((epsn - eF) * beta) * beta^2 / 2.0
fermi_dirac_d3(eF, beta, epsn) =
- 12.0 * fermi_dirac(eF, beta, epsn).^4 .* exp((epsn - eF) * 3.0 * beta) * beta^3 / 16.0 +
12.0 * fermi_dirac(eF, beta, epsn).^3 .* exp((epsn - eF) * 2.0 * beta) * beta^3 / 8.0 -
fermi_dirac(eF, beta, epsn).^2 .* exp((epsn - eF) * beta) * beta^3 / 2.0
# Boilerplate to work with the FermiDiracSmearing type
@inline evaluate(fd::FermiDiracSmearing, epsn) =
fermi_dirac(fd.eF, fd.beta, epsn)
@inline evaluate_d(fd::FermiDiracSmearing, epsn) =
fermi_dirac_d(fd.eF, fd.beta, epsn)
# Boilerplate to work with the FermiDiracSmearing type
@inline evaluate(fd::FermiDiracSmearing, epsn, eF) =
fermi_dirac(eF, fd.beta, epsn)
@inline evaluate_d(fd::FermiDiracSmearing, epsn, eF) =
fermi_dirac_d(eF, fd.beta, epsn)
function set_eF!(fd::FermiDiracSmearing, eF)
fd.eF = eF
end
"""`ZeroTemperature`:
TODO
"""
type ZeroTemperature <: SmearingFunction
eF
end
function full_hermitian(A)
A = 0.5 * (A + A')
A[diagind(A)] = real(A[diagind(A)])
if vecnorm(conj(A) - A, Inf) < 1e-10
A = Symmetric(full(real(A)))
else
A = Hermitian(full(A))
end
return A
end
"""`sorted_eig`: helper function to compute eigenvalues, then sort them
in ascending order and sort the eig-vectors as well."""
function sorted_eig(H, M)
epsn, C = eig(full_hermitian(H), full_hermitian(M))
Isort = sortperm(epsn)
return epsn[Isort], C[:, Isort]
end
import Base.setindex!
setindex!(tbm::TBModel, val, symbol) = set_array!(tbm, symbol, val)
import Base.getindex
getindex(tbm::TBModel, symbol) = get_array(tbm, symbol)
"store an array"
function set_array!(tbm::TBModel, key, val)
tbm.arrays[key] = val
end
"""retrieve an array; instead of raising an exception if `key` does not exist,
this function returns `nothing`"""
function get_array(tbm::TBModel, key)
if haskey(tbm.arrays, key)
return tbm.arrays[key]
else
return nothing
end
end
""" store k-point dependent arrays"""
set_k_array!(tbm, q, symbol, k) = set_array!(tbm, (symbol, k), q)
"""retrieve k-point dependent arrays"""
get_k_array(tbm, symbol, k) = get_array(tbm, (symbol, k))
"""`monkhorstpackgrid(cell, nkpoints)` : constructs an MP grid for the
computational cell defined by `cell` and `nkpoints`.
MonkhorstPack: K = {b/kz * j + shift}_{j=-kz/2+1,...,kz/2} with shift = 0.0.
Returns
### Parameters
* 'cell' : 3 × 1 array of lattice vector for (super)cell
* 'nkpoints' : 3 × 1 array of number of k-points in each direction. Now
it can only be (0, 0, kz::Int).
### Output
* `K`: 3 × Nk array of k-points
* `weights`: integration weights; scalar (uniform grid) or Nk array.
"""
function monkhorstpackgrid(cell::Matrix{Float64},
nkpoints::Tuple{Int64, Int64, Int64})
kx, ky, kz = nkpoints
## We need to check somewhere that 'nkpoints' and 'pbc' are compatable,
## e.g., if pbc[1]==false, then kx!=0 should return an error.
# if kx != 0 || ky != 0
# error("This boundary condition has not been implemented yet!")
# end
## We want to sample the Γ-point (which is not really necessary?)
if mod(kx,2) == 1 || mod(ky,2) == 1 || mod(kz,2) == 1
throw(ArgumentError("k should be an even number in Monkhorst-Pack
grid so that the Γ-point can be sampled!"))
end
# compute the lattice vector of reciprocal space
v1 = cell[1,:][:]
v2 = cell[2,:][:]
v3 = cell[3,:][:]
c12 = cross(v1,v2)
c23 = cross(v2,v3)
c31 = cross(v3,v1)
b1 = 2 * π * c23 / dot(v1,c23)
b2 = 2 * π * c31 / dot(v2,c31)
b3 = 2 * π * c12 / dot(v3,c12)
# We can exploit the symmetry of the BZ.
# TODO: NOTE THAT this is not necessarily first BZ
# and THE SYMMETRY HAS NOT BEEN FULLY EXPLOITED YET!!
nx = (kx==0? 1:kx)
ny = (ky==0? 1:ky)
nz = (kz==0? 1:kz)
N = nx * ny * nz
K = zeros(3, N)
weight = zeros(N)
kx_step = b1 / nx
ky_step = b2 / ny
kz_step = b3 / nz
w_step = 1.0 / ( nx * ny * nz )
# evaluate K and weight
for k1 = 1:nx, k2 = 1:ny, k3 = 1:nz
k = k1 + (k2-1) * nx + (k3-1) * nx * ny
# check when kx==0 or ky==0 or kz==0
# K[:,k] = (k1-(kx/2)) * kx_step + (k2-(ky/2)) * ky_step + (k3-(kz/2)) * kz_step
K[:,k] = (k1-(kx==0? nx:(kx/2))) * kx_step +
(k2-(ky==0? ny:(ky/2))) * ky_step + (k3-(kz==0? nz:(kz/2))) * kz_step
weight[k] = w_step
end
#print(K); println("\n"); print(weight); println("\n")
#println("sum_weight = "); print(sum(weight))
return K, weight
end
"""`monkhorstpackgrid(atm::ASEAtoms, tbm::TBModel)` : extracts cell and grid
information and returns an MP grid.
"""
monkhorstpackgrid(atm::ASEAtoms, tbm::TBModel) =
monkhorstpackgrid(cell(atm), tbm.nkpoints)
############################################################
##### update functions
"""`update_eig!(atm::ASEAtoms, tbm::TBModel)` : updates the hamiltonians
and spectral decompositions on the MP grid.
"""
function update_eig!(atm::ASEAtoms, tbm::TBModel)
K, weight = monkhorstpackgrid(atm, tbm)
nlist = NeighbourList(cutoff(tbm), atm)
nnz_est = length(nlist) * tbm.norbitals^2 + length(atm) * tbm.norbitals^2
It = zeros(Int32, nnz_est)
Jt = zeros(Int32, nnz_est)
Ht = zeros(Complex{Float64}, nnz_est)
Mt = zeros(Complex{Float64}, nnz_est)
X = positions(atm)
for n = 1:size(K, 2)
k = K[:,n]
H, M = hamiltonian!(tbm, k, It, Jt, Ht, Mt, nlist, X)
epsn, C = sorted_eig(H, M)
set_k_array!(tbm, epsn, :epsn, k)
set_k_array!(tbm, C, :C, k)
end
end
"""`update!(atm::ASEAtoms, tbm:TBModel)`: checks whether the precomputed
data stored in `tbm` needs to be updated (by comparing atom positions) and
if so, does all necessary updates. At the moment, the following are updated:
* spectral decompositions (`update_eig!`)
* the fermi-level (`update_eF!`)
"""
function update!(atm::ASEAtoms, tbm::TBModel)
Xnew = positions(atm)
Xold = tbm[:X] # (returns nothing if X has not been stored previously)
if Xnew != Xold
tbm[:X] = Xnew
# do all the updates
update_eig!(atm, tbm)
update_eF!(atm, tbm)
end
end
"""`update_eF!(tbm::TBModel)`: recompute the correct
fermi-level; using the precomputed data in `tbm.arrays`
"""
function update_eF!(atm::ASEAtoms, tbm::TBModel)
if tbm.fixed_eF
set_eF!(tbm.smearing, tbm.eF)
return
end
# the following algorithm works for Fermi-Dirac, not general Smearing
K, weight = monkhorstpackgrid(atm, tbm)
Ne = tbm.norbitals * length(atm)
nf = round(Int, ceil(Ne/2))
# update_eig!(atm, tbm)
# set an initial eF
μ = 0.0
for n = 1:size(K, 2)
k = K[:, n]
epsn_k = get_k_array(tbm, :epsn, k)
μ += weight[n] * (epsn_k[nf] + epsn_k[nf+1]) /2
end
# iteration by Newton algorithm
err = 1.0
while abs(err) > 1.0e-8
Ni = 0.0
gi = 0.0
for n = 1:size(K,2)
k = K[:, n]
epsn_k = get_k_array(tbm, :epsn, k)
Ni += weight[n] * r_sum( tbm.smearing(epsn_k, μ) )
gi += weight[n] * r_sum( @D tbm.smearing(epsn_k, μ) )
end
err = Ne - Ni
#println("\n err="); print(err)
μ = μ - err / gi
end
tbm.eF = μ
set_eF!(tbm.smearing, tbm.eF)
end
############################################################
### Hamiltonian Evaluation
"""`hamiltonian`: computes the hamiltonitan and overlap matrix for a tight
binding model.
#### Parameters:
* `atm::ASEAtoms`
* `tbm::TBModel`
* `k = k=[0.;0.;0.]` : k-point at which the hamiltonian is evaluated
### Output: H, M
* `H` : hamiltoian in CSC format
* `M` : overlap matrix in CSC format
"""
function hamiltonian(atm::ASEAtoms, tbm::TBModel, k)
nlist = NeighbourList(cutoff(tbm), atm)
nnz_est = length(nlist) * tbm.norbitals^2 + length(atm) * tbm.norbitals^2
It = zeros(Int32, nnz_est)
Jt = zeros(Int32, nnz_est)
Ht = zeros(Complex{Float64}, nnz_est)
Mt = zeros(Complex{Float64}, nnz_est)
X = positions(atm)
return hamiltonian!( tbm, k, It, Jt, Ht, Mt, nlist, X)
end
hamiltonian(atm::ASEAtoms, tbm::TBModel) =
hamiltonian(atm::ASEAtoms, tbm::TBModel, [0.;0.;0.])
function hamiltonian!(tbm::TBModel, k,
It, Jt, Ht, Mt, nlist, X)
idx = 0 # index ito triplet format
H_nm = zeros(tbm.norbitals, tbm.norbitals) # temporary arrays
M_nm = zeros(tbm.norbitals, tbm.norbitals)
temp = zeros(10)
# loop through sites
for (n, neigs, r, R, _) in Sites(nlist)
In = indexblock(n, tbm) # index-block for atom index n
exp_i_kR = exp(im * (k' * (R - (X[:,neigs] .- X[:,n]))))
# loop through the neighbours of the current atom
for m = 1:length(neigs)
Im = TightBinding.indexblock(neigs[m], tbm)
# compute hamiltonian block
H_nm = evaluate!(tbm.hop, r[m], R[:, m], H_nm)
# compute overlap block
M_nm = evaluate!(tbm.overlap, r[m], R[:, m], M_nm)
# add new indices into the sparse matrix
@inbounds for i = 1:tbm.norbitals, j = 1:tbm.norbitals
idx += 1
It[idx] = In[i]
Jt[idx] = Im[j]
Ht[idx] = H_nm[i,j] * exp_i_kR[m]
Mt[idx] = M_nm[i,j] * exp_i_kR[m]
end
end
# now compute the on-site terms (we could move these to be done in-place)
H_nn = tbm.onsite(r, R)
M_nn = tbm.overlap(0.0)
# add into sparse matrix
for i = 1:tbm.norbitals, j = 1:tbm.norbitals
idx += 1
It[idx] = In[i]
Jt[idx] = In[j]
Ht[idx] = H_nn[i,j]
Mt[idx] = M_nn[i,j]
end
end
# convert M, H into Sparse CCS and return
# TODO: The conversion to sparse format accounts for about 1/2 of the
# total cost. Since It, Jt are in an ordered format, it should be
# possible to write a specialised code that converts it to
# CCS much faster, possibly with less additional allocation?
# another option would be to store a single It, Jt somewhere
# for ALL the hamiltonians, and store multiple Ht, Mt and convert
# these "on-the-fly", depending on whether full or sparse is needed.
# but at the moment, eigfact cost MUCH more than the assembly,
# so we could choose to stop here.
return sparse(It, Jt, Ht), sparse(It, Jt, Mt)
end
"""`densitymatrix(at::ASEAtoms, tbm::TBModel) -> rho`:
### Input
* `at::ASEAtoms` : configuration
* `tbm::TBModel` : calculator
### Output
* `rho::Matrix{Float64}`: density matrix,
ρ = ∑_s f(ϵ_s) ψ_s ⊗ ψ_s
where `f` is given by `tbm.SmearingFunction`. With BZ integration, it becomes
ρ = ∑_k w^k ∑_s f(ϵ_s^k) ψ_s^k ⊗ ψ_s^k
"""
function densitymatrix(at::ASEAtoms, tbm::TBModel)
update!(at, tbm)
K, weight = monkhorstpackgrid(atm, tbm)
rho = 0.0
for n = 1:size(K, 2)
k = K[:, n]
epsn_k = get_k_array(tbm, :epsn, k)
C_k = get_k_array(tbm, :C, k)
f = tbm.smearing(epsn_k, tbm.eF)
# TODO: should eF be passed or should it be sotred in SmearingFunction?
for m = 1:length(epsn_k)
rho += weight[n] * f[m] * C_k[:,m] * C_k[:,m]'
end
end
return rho
end
############################################################
### Standard Calculator Functions
function potential_energy(at::ASEAtoms, tbm::TBModel)
update!(at, tbm)
K, weight = monkhorstpackgrid(at, tbm)
E = 0.0
for n = 1:size(K, 2)
k = K[:, n]
epsn_k = get_k_array(tbm, :epsn, k)
E += weight[n] * r_sum(tbm.smearing(epsn_k, tbm.eF) .* epsn_k)
# TODO: pass eF?
end
return E
end
function band_structure_all(at::ASEAtoms, tbm::TBModel)
update!(at, tbm)
# tbm.fixed_eF = false
# TightBinding.update_eF!(at, tbm)
na = length(at) * tbm.norbitals
K, weight = monkhorstpackgrid(at, tbm)
E = zeros(na, size(K,2))
Ne = tbm.norbitals * length(at)
nf = round(Int, ceil(Ne/2))
for n = 1:size(K, 2)
k = K[:, n]
epsn_k = get_k_array(tbm, :epsn, k)
for j = 1:na
E[j,n] = epsn_k[j]
end
end
return K, E
end
# get 2*Nb+1 bands around the fermi level
function band_structure_near_eF(Nb, at::ASEAtoms, tbm::TBModel)
update!(at, tbm)
# tbm.fixed_eF = false
# TightBinding.update_eF!(at, tbm)
K, weight = monkhorstpackgrid(at, tbm)
E = zeros(2*Nb+1, size(K,2))
Ne = tbm.norbitals * length(at)
nf = round(Int, ceil(Ne/2))
for n = 1:size(K, 2)
k = K[:, n]
epsn_k = get_k_array(tbm, :epsn, k)
E[Nb+1,n] = epsn_k[nf]
for j = 1:Nb
E[Nb+1-j,n] = epsn_k[nf-j]
E[Nb+1+j,n] = epsn_k[nf+j]
end
end
return K, E
end
function forces_k(X::Matrix{Float64}, tbm::TBModel, nlist, k::Vector{Float64})
# obtain the precomputed arrays
epsn = get_k_array(tbm, :epsn, k)
C = get_k_array(tbm, :C, k)
df = tbm.smearing(epsn, tbm.eF) + epsn .* (@D tbm.smearing(epsn, tbm.eF))
# precompute some products
const C_df_Ct = (C * (df' .* C)')::Matrix{Complex{Float64}}
const C_dfepsn_Ct = (C * ((df.*epsn)' .* C)')::Matrix{Complex{Float64}}
# allocate forces
const frc = zeros(Complex{Float64}, 3, size(X,2))
# pre-allocate dH, with a (dumb) initial guess for the size
const dH_nn = zeros(3, tbm.norbitals, tbm.norbitals, 6)
const dH_nm = zeros(3, tbm.norbitals, tbm.norbitals)
const dM_nm = zeros(3, tbm.norbitals, tbm.norbitals)
# loop through all atoms, to compute the force on atm[n]
for (n, neigs, r, R) in Sites(nlist)
neigs::Vector{Int}
R::Matrix{Float64}
# compute the block of indices for the orbitals belonging to n
In = indexblock(n, tbm)
# compute ∂H_mm/∂y_n (onsite terms) M_nn = const ⇒ dM_nn = 0
# dH_nn should be 3 x norbitals x norbitals x nneigs
# dH_nn = (@D tbm.onsite(r, R))::Array{Float64,4}
# in-place version
if length(neigs) > size(dH_nn, 4)
dH_nn = zeros(3, tbm.norbitals, tbm.norbitals,
ceil(Int, 1.5*length(neigs)))
end
evaluate_d!(tbm.onsite, r, R, dH_nn)
for i_n = 1:length(neigs)
m = neigs[i_n]
Im = indexblock(m, tbm)
kR = dot(R[:,i_n] - (X[:,neigs[i_n]] - X[:,n]), k)
eikr = exp(im * kR)::Complex{Float64}
# compute ∂H_nm/∂y_n (hopping terms) and ∂M_nm/∂y_n
grad!(tbm.hop, r[i_n], -R[:,i_n], dH_nm)
grad!(tbm.overlap, r[i_n], -R[:,i_n], dM_nm)
# the following is a hack to put the on-site assembly into the
# innermost loop
# F_n = - ∑_s f'(ϵ_s) < ψ_s | H,n - ϵ_s * M,n | ψ_s >
for a = 1:tbm.norbitals, b = 1:tbm.norbitals
t1 = 2.0 * real(C_df_Ct[Im[a], In[b]] * eikr)
t2 = 2.0 * real(C_dfepsn_Ct[Im[a], In[b]] * eikr)
t3 = C_df_Ct[In[a],In[b]]
# add contributions to the force
for j = 1:3
frc[j,n] = frc[j,n] - dH_nm[j,a,b] * t1 +
dM_nm[j,a,b] * t2 + dH_nn[j,a,b,i_n] * t3
frc[j,m] = frc[j,m] - t3 * dH_nn[j,a,b,i_n]
end
end
end # m in neigs-loop
end # sites-loop
return frc
end
function forces(atm::ASEAtoms, tbm::TBModel)
# tell tbm to update the spectral decompositions
update!(atm, tbm)
# precompute neighbourlist
nlist = NeighbourList(cutoff(tbm), atm)
X = positions(atm)
# BZ integration loop
K, weight = monkhorstpackgrid(atm, tbm)
# allocate output
frc = zeros(3, length(atm))
for iK = 1:size(K,2)
frc += weight[iK] * real(forces_k(X, tbm, nlist, K[:,iK]))
end
return frc
end
# compute all forces on all the atoms
function forces_debug(atm::ASEAtoms, tbm)
# tell tbm to update the spectral decompositions
update!(atm, tbm)
# allocate output
frc = zeros(3, length(atm))
# precompute neighbourlist
nlist = NeighbourList(cutoff(tbm), atm)
X = positions(atm)
@code_warntype forces_k(X, tbm, nlist, zeros(3))
end
potential_energy_d(atm::ASEAtoms, tbm::TBModel) = -forces(atm, tbm)
############################################################
### Site Energy Stuff
function site_energy(l::Integer, atm::ASEAtoms, tbm::TBModel)
# tell tbm to update the spectral decompositions
update!(atm, tbm)
# BZ integration loop
K, weight = monkhorstpackgrid(atm, tbm)
# use the following parameters as those in update_eig!
nlist = NeighbourList(cutoff(tbm), atm)
nnz_est = length(nlist) * tbm.norbitals^2 + length(atm) * tbm.norbitals^2
It = zeros(Int32, nnz_est)
Jt = zeros(Int32, nnz_est)
Ht = zeros(Complex{Float64}, nnz_est)
Mt = zeros(Complex{Float64}, nnz_est)
X = positions(atm)
Es = 0.0
for n = 1:size(K, 2)
k = K[:, n]
epsn = get_k_array(tbm, :epsn, k)
C = get_k_array(tbm, :C, k)::Matrix{Complex{Float64}}
# precompute electron distribution function
f = tbm.smearing(epsn, tbm.eF) .* epsn
# overlap matrix is needed in this calculation
# ([M^{1/2}*ψ]_i)^2 → [M*ψ]_i*[ψ]_i
H, M = hamiltonian!(tbm, k, It, Jt, Ht, Mt, nlist, X)
MC = M * C::Matrix{Complex{Float64}}
I = indexblock(l, tbm)
for j = 1:tbm.norbitals
# the first component of the following line should be conjugate
Es += weight[n] * r_sum(real(f .* slice(C, I[j], :) .* slice(MC, I[j], :)))
# Es += weight[n] * r_sum( f .* (slice(C, I[j], :) .* slice(MC, I[j], :)) )
end
end
return Es
end
site_energy(nn::Array{Int}, atm::ASEAtoms, tbm::TBModel) =
reshape(Float64[ site_energy(n, atm, tbm) for n in nn ], size(nn))
# site_forces always returns a complete gradient, i.e. dEs = d x Natm
# When idx is an array, then the return-value is the gradient of \sum_{i ∈ idx} E_i
function site_forces(idx::Array{Int,1}, atm::ASEAtoms, tbm::TBModel)
# tell tbm to update the spectral decompositions
update!(atm, tbm)
# BZ integration loop
K, weight = monkhorstpackgrid(atm, tbm)
# allocate output
sfrc = zeros(Float64, 3, length(atm))
# precompute neighbourlist
nlist = NeighbourList(cutoff(tbm), atm)
Nneig = 1
for (n, neigs, r, R) in Sites(nlist)
if length(neigs) > Nneig
Nneig = length(neigs)
end
end
X = positions(atm)
# assemble the site forces for k-points
for iK = 1:size(K,2)
sfrc += weight[iK] *
real(site_forces_k(idx, X, tbm, nlist, Nneig, K[:,iK]))
end
return sfrc
end
# scalar index: just wraps the vector version
site_forces(n::Int, atm::ASEAtoms, tbm::TBModel) = site_forces([n;], atm, tbm)
# With a given k-point, compute the site force by loop through eigenpairs (index s)
# note that in the old version, we loop through through atoms
# E_l = ∑_s f(ɛ_s)⋅[ψ_s]_l^2
# E_l = ∑_s f(ɛ_s)⋅[ψ_s]_l⋅[M⋅ψ_s]_l
# E_{l,n} = ∑_s ( f'(ɛ_s)⋅ɛ_{s,n}⋅[ψ_s]_l⋅[M⋅ψ_s]_l + 2.0⋅f(ɛ_s)⋅[ψ_s]_{l,n}⋅[M⋅ψ_s]_l
# + f(ɛ_s)⋅[ψ_s]⋅[M_{,n}⋅ψ_s]_l )
# We loop through eigenpair s to compute the first two parts and through atom n for the third part
#
function site_forces_k(idx::Array{Int,1}, X::Matrix{Float64},
tbm::TBModel, nlist, Nneig, k::Vector{Float64};
beta = ones(size(X,2)))
# obtain the precomputed arrays
epsn = get_k_array(tbm, :epsn, k)
C = get_k_array(tbm, :C, k)::Matrix{Complex{Float64}}
# some constant parameters
Nelc = length(epsn)
Natm = size(X,2)
Norb = tbm.norbitals
# overlap matrix is needed in this calculation
# use the following parameters as those in update_eig!
nnz_est = length(nlist) * Norb^2 + Natm * Norb^2
It = zeros(Int32, nnz_est)
Jt = zeros(Int32, nnz_est)
Ht = zeros(Complex{Float64}, nnz_est)
Mt = zeros(Complex{Float64}, nnz_est)
~, M = hamiltonian!(tbm, k, It, Jt, Ht, Mt, nlist, X)
MC = M * C::Matrix{Complex{Float64}}
# allocate output
const dEs = zeros(Complex{Float64}, 3, Natm)
# pre-allocate dM
const dM_nm = zeros(3, Norb, Norb)
const Msn = zeros(Complex{Float64}, Nelc)
# const eps_s_n = zeros(Float64, 3, Natm)
# const psi_s_n = zeros(Float64, 3, Natm, Nelc)
# precompute electron distribution function
f = tbm.smearing(epsn, tbm.eF) .* epsn
df = tbm.smearing(epsn, tbm.eF) + epsn .* (@D tbm.smearing(epsn, tbm.eF))
# loop through all eigenstates to compute the hessian
for s = 1 : Nelc
# compute ϵ_{s,n} and ψ_{s,n}
eps_s_n, psi_s_n = d_eigenstate_k(s, tbm, X, nlist, Nneig, k)
# loop for the first part
for d = 1:3
for n = 1:Natm
Msn = M * psi_s_n[d, n, :][:]
for id in idx
# in this iteration of the loop we compute the contributions
# that come from the site i. hence multiply everything with beta[i]
Ii = indexblock(id, tbm)
dEs[d, n] -= beta[id] * df[s] * eps_s_n[d, n] * sum( C[Ii, s] .* MC[Ii, s] )
dEs[d, n] -= beta[id] * f[s] * sum( MC[Ii, s] .* psi_s_n[d, n, Ii][:] )
dEs[d, n] -= beta[id] * f[s] * sum( C[Ii, s] .* Msn[Ii] )
end # loop of id
end # loop of d
end # loop of n
end # loop of s
# loop through all atoms, to compute the last part
for (n, neigs, r, R) in Sites(nlist)
# consider only the rows related to site idx
if n in idx
# compute the block of indices for the orbitals belonging to n
In = indexblock(n, tbm)
exp_i_kR = exp(im * (k' * (R - (X[:, neigs] .- X[:, n]))))
for i_n = 1:length(neigs)
m = neigs[i_n]
Im = indexblock(m, tbm)
eikr = exp_i_kR[i_n]
# compute and store ∂M_mn/∂y_n
grad!(tbm.overlap, r[i_n], R[:,i_n], dM_nm)
# sum over all eigenpairs
for s = 1:Nelc
for d = 1:3
dEs[d, n] -= beta[n] * f[s] * sum( C[In, s] .* ( - slice(dM_nm, d, :, :) * C[Im, s] ) ) * eikr
dEs[d, m] -= beta[n] * f[s] * sum( C[In, s] .* ( slice(dM_nm, d, :, :) * C[Im, s] ) ) * eikr
end # loop for d
end # loop for s
end # loop for neighbour i_n
end # end of if
end # loop for atom n
return dEs
# note that this is in fact the site force, -dEs
end
###################### Hessian and Higher-oerder derivatives ##########################
# For a given s and a given k-point, returns ψ_{s,n} and ϵ_{s,n} for all n∈{1,⋯,d×Natm}
# Input
# s : which eigenstate
# k : k-point
# Output
# psi_s_n : ψ_{s,n} for all n, a 3 × Natm × Nelc matrix
# eps_s_n : ϵ_{s,n} for all n, a 3 × Natm matrix
#
# Algorithm
# ϵ_{s,n} = < ψ_s | H_{,n} - ϵ_s * M,n | ψ_s >
#
# ψ_{s,n} = ∑_{t,ϵ_t≠ϵ_s} ψ_t < ψ_t | ϵ_s⋅M_{,n} - H_{,n} | ψ_s > / (ϵ_t-ϵ_s)
# - 1/2 ∑_{t,ϵ_t=ϵ_s} ψ_t < ψ_t | M_{,n} | ψ_s >
#
# Step 1. compute g_s_n = (ϵ_s⋅M_{,n} - H_{,n}) ⋅ ψ_s
# and f_s_n = M_{,n} ⋅ ψ_s
# Step 2. (C' * g) ./ (epsilon - epsilon[s])
# with the second part added in the loop for ϵ_t =≠ ϵ_s
function d_eigenstate_k(s::Int, tbm::TBModel, X::Matrix{Float64}, nlist, Nneig::Int,
k::Vector{Float64})
# obtain the precomputed arrays
epsn = get_k_array(tbm, :epsn, k)
C = get_k_array(tbm, :C, k)::Matrix{Complex{Float64}}
# some constant parameters
Nelc = length(epsn)
Natm = size(X,2)
Norb = tbm.norbitals
# allocate memory
psi_s_n = zeros(Complex{Float64}, 3*Natm, Nelc)
eps_s_n = zeros(Complex{Float64}, 3*Natm)
g_s_n = zeros(Complex{Float64}, 3*Natm, Nelc)
f_s_n = zeros(Complex{Float64}, 3*Natm, Nelc)
const dH_nn = zeros(3, Norb, Norb, Nneig)
const dH_nm = zeros(3, Norb, Norb)
const dM_nm = zeros(3, Norb, Norb)
# Step 1. loop through all atoms to compute g_s_n and f_s_n for all n
for (n, neigs, r, R) in Sites(nlist)
In = indexblock(n, tbm)
exp_i_kR = ( exp(im * (k' * (R - (X[:, neigs] .- X[:, n])))) )' # conjugate here for later use
# compute and store ∂H_nn/∂y_n (onsite terms)
evaluate_d!(tbm.onsite, r, R, dH_nn)
for i_n = 1:length(neigs)
m = neigs[i_n]
Im = indexblock(m, tbm)
eikr = exp_i_kR[i_n]
# compute and store ∂H_nm/∂y_m (hopping terms) and ∂M_nm/∂y_m
grad!(tbm.hop, r[i_n], R[:,i_n], dH_nm)
grad!(tbm.overlap, r[i_n], R[:,i_n], dM_nm)
for d = 1:3
md = d + 3*(m-1)
nd = d + 3*(n-1)
g_s_n[md, In] += ( slice(dH_nn, d, :, :, i_n) * C[In, s] )'
g_s_n[nd, In] -= ( slice(dH_nn, d, :, :, i_n) * C[In, s] )'
g_s_n[md, In] += ( slice(dH_nm, d, :, :) * C[Im, s] )' * eikr
g_s_n[nd, In] -= ( slice(dH_nm, d, :, :) * C[Im, s] )' * eikr
f_s_n[md, In] += ( slice(dM_nm, d, :, :) * C[Im, s] )' * eikr
f_s_n[nd, In] -= ( slice(dM_nm, d, :, :) * C[Im, s] )' * eikr
end # loop for dimension
end # loop for neighbours
end # loop for atomic sites
g_s_n = epsn[s] * f_s_n - g_s_n
# Step 2. compute eps_s_n and psi_s_n for all n
# TODO: use BLAS for matrix-matrix/vector multiplication?
# compute ϵ_{s,n}
eps_s_n = real( - g_s_n * C[:,s] )
diff_eps_inv = zeros(Float64, Nelc)
# loop through all orbitals to compute 1/(ϵ_t-ϵ_s) and add the second part of ψ_{s,n}
for t = 1:Nelc
if abs(epsn[t]-epsn[s]) > 1e-10
diff_eps_inv[t] = 1.0/(epsn[t]-epsn[s])
else
diff_eps_inv[t] = 0.0
# psi_s_n -= 0.5 * ( C[:,t] * (f_s_n * C[:,t])' )'
psi_s_n -= 0.5 * transpose( C[:,t] * (f_s_n * C[:,t])' )
end
end # loop for orbitals
# g = - (C' * gsn) ./ (epsilon - epsilon[s])
# use BLAS here!! gemm!
g_s_n = g_s_n * C
for jj = 1 : Nelc
@simd for ii = 1 : 3*Natm
@inbounds g_s_n[ii,jj] *= diff_eps_inv[jj]
end
end
# add the first part of ψ_{s,n}
# psi_s_n += ( C * g_s_n' )'
psi_s_n += transpose( C * g_s_n' )
# return eps_s_n, psi_s_n
return reshape(eps_s_n, 3, Natm), reshape(psi_s_n, 3, Natm, Nelc)
end
# hessian always returns a complete hessian, i.e. hessian = ( d × Natm )^2
function hessian(atm::ASEAtoms, tbm::TBModel)
# tell tbm to update the spectral decompositions
update!(atm, tbm)
# BZ integration loop
K, weight = monkhorstpackgrid(atm, tbm)
# allocate output
hessian = zeros(Float64, 3, length(atm), 3, length(atm))
# precompute neighbourlist
nlist = NeighbourList(cutoff(tbm), atm)
Nneig = 1
for (n, neigs, r, R) in Sites(nlist)
if length(neigs) > Nneig
Nneig = length(neigs)
end
end
X = positions(atm)
# loop for all k-points
for iK = 1:size(K,2)
Hess_k, ~ = hessian_k(X, tbm, nlist, Nneig, K[:,iK])