From c0fd0ac544fcb669b3b4df905d6d97d10c8c3556 Mon Sep 17 00:00:00 2001 From: hjchen1983 Date: Sun, 28 Feb 2016 22:46:48 +0000 Subject: [PATCH] a few updates w.r.t. site force and k-points --- TightBinding.jl | 200 +++++++++++++++++------ notebooks/TightBinding Tests.ipynb | 248 ++++++++++++++--------------- 2 files changed, 272 insertions(+), 176 deletions(-) diff --git a/TightBinding.jl b/TightBinding.jl index 3067ec9..092b0cb 100644 --- a/TightBinding.jl +++ b/TightBinding.jl @@ -453,8 +453,8 @@ function hamiltonian!(tbm::TBModel, k, 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] + 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) @@ -748,7 +748,6 @@ site_energy(nn::Array{Int}, atm::ASEAtoms, tbm::TBModel) = - # 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 @@ -764,10 +763,18 @@ function site_forces(idx::Array{Int,1}, atm::ASEAtoms, tbm::TBModel) # 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, K[:,iK])) + real(site_forces_k(idx, X, tbm, nlist, Nneig, K[:,iK])) end return sfrc @@ -780,7 +787,103 @@ 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 + + + +# an old version for site force +function site_forces_k_old(idx::Array{Int,1}, X::Matrix{Float64}, tbm::TBModel, nlist, k::Vector{Float64}; beta = ones(size(X,2))) # obtain the precomputed arrays @@ -927,10 +1030,8 @@ function site_forces_k(idx::Array{Int,1}, X::Matrix{Float64}, dEs[a,m] += beta[id] * f[s] * sum( C[Ii, s] .* MC_s_m[Ii,a] ) end end - end # loop for s, eigenpairs end # loop for n, atomic sites - return -dEs # , [1:Natm;] end @@ -938,7 +1039,6 @@ end - ###################### Hessian and Higher-oerder derivatives ########################## @@ -976,7 +1076,7 @@ function d_eigenstate_k(s::Int, tbm::TBModel, X::Matrix{Float64}, nlist, Nneig:: # allocate memory psi_s_n = zeros(Complex{Float64}, 3*Natm, Nelc) - eps_s_n = zeros(Float64, 3*Natm) + 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) @@ -985,7 +1085,6 @@ function d_eigenstate_k(s::Int, tbm::TBModel, X::Matrix{Float64}, nlist, Nneig:: # 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])))) @@ -995,6 +1094,7 @@ function d_eigenstate_k(s::Int, tbm::TBModel, X::Matrix{Float64}, nlist, Nneig:: 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) @@ -1006,11 +1106,10 @@ function d_eigenstate_k(s::Int, tbm::TBModel, X::Matrix{Float64}, nlist, Nneig:: 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 + 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 @@ -1053,7 +1152,6 @@ end - # hessian always returns a complete hessian, i.e. hessian = ( d × Natm )^2 function hessian(atm::ASEAtoms, tbm::TBModel) @@ -1103,6 +1201,7 @@ potential_energy_d2(atm::ASEAtoms, tbm::TBModel) = hessian(atm, tbm) # ɛ_{s,mn} ∈ R^{ Nelc × 3 × Natm × 3 × Natm } # note that the output of ɛ_{s,mn} is stored for usage of computing d3E # TODO: have not added e^ikr into the hamiltonian yet +# TODO: we do not need the whole hessian matrix, but only those related to the centred atom function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64}) @@ -1186,6 +1285,7 @@ function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Flo for i_n = 1:length(neigs) m = neigs[i_n] Im = indexblock(m, tbm) + eikr = exp_i_kR[i_n] # compute and store ∂H, ∂^2H and ∂M, ∂^2M # evaluate!(tbm.overlap, r[i_n], R[:, i_n], M_nm) @@ -1206,7 +1306,7 @@ function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Flo + C[In, s]' * ( eps_s_n[d2, l] * slice(dM_nm, d1, :, :) ) * C[Im,s] - )[1] + )[1] * eikr eps_s_mn[s, d1, l, d2, n] += ( C[In, s]' * ( - slice(dH_nm, d2, :, :) + epsn[s] * slice(dM_nm, d2, :, :) @@ -1214,7 +1314,7 @@ function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Flo + C[In, s]' * ( eps_s_n[d1, l] * slice(dM_nm, d2, :, :) ) * C[Im,s] - )[1] + )[1] * eikr eps_s_mn[s, d1, m, d2, l] += ( C[In, s]' * ( slice(dH_nm, d1, :, :) - epsn[s] * slice(dM_nm, d1, :, :) @@ -1222,7 +1322,7 @@ function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Flo - C[In, s]' * ( eps_s_n[d2, l] * slice(dM_nm, d1, :, :) ) * C[Im,s] - )[1] + )[1] * eikr eps_s_mn[s, d1, l, d2, m] += ( C[In, s]' * ( slice(dH_nm, d2, :, :) - epsn[s] * slice(dM_nm, d2, :, :) @@ -1230,7 +1330,7 @@ function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Flo - C[In, s]' * ( eps_s_n[d1, l] * slice(dM_nm, d2, :, :) ) * C[Im,s] - )[1] + )[1] * eikr end # loop for atom l # contributions from hopping terms @@ -1238,19 +1338,19 @@ function hessian_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Flo eps_s_mn[s, d1, n, d2, n] += ( C[In, s]' * ( slice(d2H_nm, d1, d2, :, :) - epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * C[Im,s] - )[1] + )[1] * eikr eps_s_mn[s, d1, m, d2, m] += ( C[In, s]' * ( slice(d2H_nm, d1, d2, :, :) - epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * C[Im,s] - )[1] + )[1] * eikr eps_s_mn[s, d1, m, d2, n] += ( C[In, s]' * ( - slice(d2H_nm, d1, d2, :, :) + epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * C[Im,s] - )[1] + )[1] * eikr eps_s_mn[s, d1, n, d2, m] += ( C[In, s]' * ( - slice(d2H_nm, d1, d2, :, :) + epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * C[Im,s] - )[1] + )[1] * eikr # contributions from onsite terms m1 = 3*(i_n-1) + d1 @@ -1347,7 +1447,6 @@ function d3E(atm::ASEAtoms, tbm::TBModel) end - potential_energy_d3(atm::ASEAtoms, tbm::TBModel) = d3E(atm, tbm) @@ -1473,6 +1572,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 for i_n = 1:length(neigs) m = neigs[i_n] Im = indexblock(m, tbm) + eikr = exp_i_kR[i_n] # compute and store ∂H, ∂^2H and ∂M, ∂^2M evaluate_fd!(tbm.hop, R[:,i_n], dH_nm) @@ -1496,7 +1596,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 + 2.0 * eps_s_n[d3, q] * C[In, s]' * slice(dM_nm, d1, :, :) * psi_s_n[d2, p, Im][:] + 2.0 * psi_s_n[d2, p, In][:]' * ( - slice(dH_nm, d1, :, :) + epsn[s] * slice(dM_nm, d1, :, :) ) * psi_s_n[d3, q, Im][:] - )[1] + )[1] * eikr # 2. mpq D3E[d1, m, d2, p, d3, q] += feps3[s] * ( - eps_s_mn[s, d2, p, d3, q] * C[In, s]' * slice(dM_nm, d1, :, :) * C[Im, s] @@ -1504,7 +1604,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 - 2.0 * eps_s_n[d3, q] * C[In, s]' * slice(dM_nm, d1, :, :) * psi_s_n[d2, p, Im][:] + 2.0 * psi_s_n[d2, p, In][:]' * ( slice(dH_nm, d1, :, :) - epsn[s] * slice(dM_nm, d1, :, :) ) * psi_s_n[d3, q, Im][:] - )[1] + )[1] * eikr # 3. pnq D3E[d1, p, d2, n, d3, q] += feps3[s] * ( eps_s_mn[s, d1, p, d3, q] * C[In, s]' * slice(dM_nm, d2, :, :) * C[Im, s] @@ -1512,7 +1612,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 + 2.0 * eps_s_n[d3, q] * C[In, s]' * slice(dM_nm, d2, :, :) * psi_s_n[d1, p, Im][:] + 2.0 * psi_s_n[d1, p, In][:]' * ( - slice(dH_nm, d2, :, :) + epsn[s] * slice(dM_nm, d2, :, :) ) * psi_s_n[d3, q, Im][:] - )[1] + )[1] * eikr # 4. pmq D3E[d1, p, d2, m, d3, q] += feps3[s] * ( - eps_s_mn[s, d1, p, d3, q] * C[In, s]' * slice(dM_nm, d2, :, :) * C[Im, s] @@ -1520,7 +1620,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 - 2.0 * eps_s_n[d3, q] * C[In, s]' * slice(dM_nm, d2, :, :) * psi_s_n[d1, p, Im][:] + 2.0 * psi_s_n[d1, p, In][:]' * ( slice(dH_nm, d2, :, :) - epsn[s] * slice(dM_nm, d2, :, :) ) * psi_s_n[d3, q, Im][:] - )[1] + )[1] * eikr # 5. pqn D3E[d1, p, d2, q, d3, n] += feps3[s] * ( eps_s_mn[s, d1, p, d2, q] * C[In, s]' * slice(dM_nm, d3, :, :) * C[Im, s] @@ -1528,7 +1628,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 + 2.0 * eps_s_n[d2, q] * C[In, s]' * slice(dM_nm, d3, :, :) * psi_s_n[d1, p, Im][:] + 2.0 * psi_s_n[d1, p, In][:]' * ( - slice(dH_nm, d3, :, :) + epsn[s] * slice(dM_nm, d3, :, :) ) * psi_s_n[d2, q, Im][:] - )[1] + )[1] * eikr # 6. pqm D3E[d1, p, d2, q, d3, m] += feps3[s] * ( - eps_s_mn[s, d1, p, d2, q] * C[In, s]' * slice(dM_nm, d3, :, :) * C[Im, s] @@ -1536,7 +1636,7 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 - 2.0 * eps_s_n[d2, q] * C[In, s]' * slice(dM_nm, d3, :, :) * psi_s_n[d1, p, Im][:] + 2.0 * psi_s_n[d1, p, In][:]' * ( slice(dH_nm, d3, :, :) - epsn[s] * slice(dM_nm, d3, :, :) ) * psi_s_n[d2, q, Im][:] - )[1] + )[1] * eikr end # loop for atom p end # loop for atom q @@ -1549,73 +1649,73 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 - eps_s_n[d3, l] * C[In, s]' * slice(d2M_nm, d1, d2, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( slice(d2H_nm, d1, d2, :, :) - epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * psi_s_n[d3, l, Im][:] - )[1] + )[1] * eikr # 2. mml D3E[d1, m, d2, m, d3, l] += feps3[s] * ( - eps_s_n[d3, l] * C[In, s]' * slice(d2M_nm, d1, d2, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( slice(d2H_nm, d1, d2, :, :) - epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * psi_s_n[d3, l, Im][:] - )[1] + )[1] * eikr # 3. nml D3E[d1, n, d2, m, d3, l] += feps3[s] * ( eps_s_n[d3, l] * C[In, s]' * slice(d2M_nm, d1, d2, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( - slice(d2H_nm, d1, d2, :, :) + epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * psi_s_n[d3, l, Im][:] - )[1] + )[1] * eikr # 4. mnl D3E[d1, m, d2, n, d3, l] += feps3[s] * ( eps_s_n[d3, l] * C[In, s]' * slice(d2M_nm, d1, d2, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( - slice(d2H_nm, d1, d2, :, :) + epsn[s] * slice(d2M_nm, d1, d2, :, :) ) * psi_s_n[d3, l, Im][:] - )[1] + )[1] * eikr # 5. nln D3E[d1, n, d2, l, d3, n] += feps3[s] * ( - eps_s_n[d2, l] * C[In, s]' * slice(d2M_nm, d1, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( slice(d2H_nm, d1, d3, :, :) - epsn[s] * slice(d2M_nm, d1, d3, :, :) ) * psi_s_n[d2, l, Im][:] - )[1] + )[1] * eikr # 6. mlm D3E[d1, m, d2, l, d3, m] += feps3[s] * ( - eps_s_n[d2, l] * C[In, s]' * slice(d2M_nm, d1, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( slice(d2H_nm, d1, d3, :, :) - epsn[s] * slice(d2M_nm, d1, d3, :, :) ) * psi_s_n[d2, l, Im][:] - )[1] + )[1] * eikr # 7. nlm D3E[d1, n, d2, l, d3, m] += feps3[s] * ( eps_s_n[d2, l] * C[In, s]' * slice(d2M_nm, d1, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( - slice(d2H_nm, d1, d3, :, :) + epsn[s] * slice(d2M_nm, d1, d3, :, :) ) * psi_s_n[d2, l, Im][:] - )[1] + )[1] * eikr # 8. mln D3E[d1, m, d2, l, d3, n] += feps3[s] * ( eps_s_n[d2, l] * C[In, s]' * slice(d2M_nm, d1, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( - slice(d2H_nm, d1, d3, :, :) + epsn[s] * slice(d2M_nm, d1, d3, :, :) ) * psi_s_n[d2, l, Im][:] - )[1] + )[1] * eikr # 9. lnn D3E[d1, l, d2, n, d3, n] += feps3[s] * ( - eps_s_n[d1, l] * C[In, s]' * slice(d2M_nm, d2, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( slice(d2H_nm, d2, d3, :, :) - epsn[s] * slice(d2M_nm, d2, d3, :, :) ) * psi_s_n[d1, l, Im][:] - )[1] + )[1] * eikr # 10. lmm D3E[d1, l, d2, m, d3, m] += feps3[s] * ( - eps_s_n[d1, l] * C[In, s]' * slice(d2M_nm, d2, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( slice(d2H_nm, d2, d3, :, :) - epsn[s] * slice(d2M_nm, d2, d3, :, :) ) * psi_s_n[d1, l, Im][:] - )[1] + )[1] * eikr # 11. lnm D3E[d1, l, d2, n, d3, m] += feps3[s] * ( eps_s_n[d1, l] * C[In, s]' * slice(d2M_nm, d2, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( - slice(d2H_nm, d2, d3, :, :) + epsn[s] * slice(d2M_nm, d2, d3, :, :) ) * psi_s_n[d1, l, Im][:] - )[1] + )[1] * eikr # 12. lmn D3E[d1, l, d2, m, d3, n] += feps3[s] * ( eps_s_n[d1, l] * C[In, s]' * slice(d2M_nm, d2, d3, :, :) * C[Im, s] + 2.0 * C[In, s]' * ( - slice(d2H_nm, d2, d3, :, :) + epsn[s] * slice(d2M_nm, d2, d3, :, :) ) * psi_s_n[d1, l, Im][:] - )[1] + )[1] * eikr end # loop for atom l # contributions from hopping terms @@ -1625,42 +1725,42 @@ function d3E_k(X::Matrix{Float64}, tbm::TBModel, nlist, Nneig, k::Vector{Float64 D3E[d1, n, d2, n, d3, n] += feps3[s] * ( C[In, s]' * ( - slice(d3H_nm, d1, d2, d3, :, :) + epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 2. nnm D3E[d1, n, d2, n, d3, m] += feps3[s] * ( C[In, s]' * ( slice(d3H_nm, d1, d2, d3, :, :) - epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 3. nmn D3E[d1, n, d2, m, d3, n] += feps3[s] * ( C[In, s]' * ( slice(d3H_nm, d1, d2, d3, :, :) - epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 4. nmm D3E[d1, n, d2, m, d3, m] += feps3[s] * ( C[In, s]' * ( - slice(d3H_nm, d1, d2, d3, :, :) + epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 5. mmm D3E[d1, m, d2, m, d3, m] += feps3[s] * ( C[In, s]' * ( slice(d3H_nm, d1, d2, d3, :, :) - epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 6. mmn D3E[d1, m, d2, m, d3, n] += feps3[s] * ( C[In, s]' * ( - slice(d3H_nm, d1, d2, d3, :, :) + epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 7. mnm D3E[d1, m, d2, n, d3, m] += feps3[s] * ( C[In, s]' * ( - slice(d3H_nm, d1, d2, d3, :, :) + epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # 8. mnn D3E[d1, m, d2, n, d3, n] += feps3[s] * ( C[In, s]' * ( slice(d3H_nm, d1, d2, d3, :, :) - epsn[s] * slice(d3M_nm, d1, d2, d3, :, :) ) * C[Im, s] - )[1] + )[1] * eikr # contributions from onsite terms diff --git a/notebooks/TightBinding Tests.ipynb b/notebooks/TightBinding Tests.ipynb index 9f8f519..5424f5a 100644 --- a/notebooks/TightBinding Tests.ipynb +++ b/notebooks/TightBinding Tests.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 12, + "execution_count": 3, "metadata": { "collapsed": false }, @@ -178,7 +178,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 40, "metadata": { "collapsed": false }, @@ -200,7 +200,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "metadata": { "collapsed": false, "scrolled": true @@ -223,17 +223,17 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 4.9699125e-04 \n", - " 3 | 4.9673946e-05 \n", - " 4 | 4.9671500e-06 \n", - " 5 | 4.9686804e-07 \n", - " 6 | 5.1582373e-08 \n", - " 7 | 1.8583959e-08 \n", - " 8 | 1.6943414e-08 \n", - " 9 | 1.4380289e-06 \n", - " 10 | 3.3832819e-07 \n", - " 11 | 5.3629033e-05 \n", - " 12 | 9.2336723e-04 \n", + " 2 | 3.9513270e-04 \n", + " 3 | 3.9401826e-05 \n", + " 4 | 3.9390578e-06 \n", + " 5 | 3.9374294e-07 \n", + " 6 | 3.8355864e-08 \n", + " 7 | 2.5033061e-08 \n", + " 8 | 1.7924798e-07 \n", + " 9 | 2.2220583e-06 \n", + " 10 | 2.2646938e-05 \n", + " 11 | 2.6157016e-04 \n", + " 12 | 2.9358722e-03 \n", "-----------------------------\n" ] } @@ -275,7 +275,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -297,17 +297,17 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 4.2598985e-04 \n", - " 3 | 4.2599965e-05 \n", - " 4 | 4.2599976e-06 \n", - " 5 | 4.2599917e-07 \n", - " 6 | 4.2599604e-08 \n", - " 7 | 3.9204751e-09 \n", - " 8 | 4.3368087e-09 \n", - " 9 | 3.8163916e-08 \n", - " 10 | 2.5243684e-07 \n", - " 11 | 2.4286129e-06 \n", - " 12 | 2.5560249e-05 \n", + " 2 | 2.9797294e-03 \n", + " 3 | 2.9797806e-04 \n", + " 4 | 2.9797814e-05 \n", + " 5 | 2.9797817e-06 \n", + " 6 | 2.9793876e-07 \n", + " 7 | 3.4069969e-08 \n", + " 8 | 6.3837824e-08 \n", + " 9 | 3.4694470e-07 \n", + " 10 | 6.1062266e-06 \n", + " 11 | 7.8062556e-05 \n", + " 12 | 5.3776428e-04 \n", "-----------------------------\n" ] } @@ -325,9 +325,8 @@ "tbm = NRLTB.NRLTBModel(elem = NRLTB.Al_spd)\n", "tbm.nkpoints = (0,0,0)\n", "\n", - "s=3\n", + "s=1\n", "\n", - "tbm.nkpoints = (0,0,0)\n", "K, weight = TightBinding.monkhorstpackgrid(at, tbm)\n", "X = positions(at)\n", "TightBinding.update!(at, tbm)\n", @@ -369,7 +368,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 13, "metadata": { "collapsed": false, "scrolled": true @@ -392,17 +391,17 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 1.0444713e-03 \n", - " 3 | 1.0384657e-04 \n", - " 4 | 1.0378725e-05 \n", - " 5 | 1.0378383e-06 \n", - " 6 | 1.0398033e-07 \n", - " 7 | 1.1762433e-08 \n", - " 8 | 5.4592755e-08 \n", - " 9 | 6.0699508e-07 \n", - " 10 | 5.0106858e-06 \n", - " 11 | 5.3865628e-05 \n", - " 12 | 4.5044132e-04 \n", + " 2 | 2.4388635e-04 \n", + " 3 | 2.4180200e-05 \n", + " 4 | 2.4159564e-06 \n", + " 5 | 2.4158272e-07 \n", + " 6 | 2.4305288e-08 \n", + " 7 | 4.2067818e-09 \n", + " 8 | 1.0174231e-08 \n", + " 9 | 1.8135322e-07 \n", + " 10 | 1.0999275e-06 \n", + " 11 | 1.1680686e-05 \n", + " 12 | 6.3722390e-05 \n", "-----------------------------\n" ] } @@ -445,7 +444,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 41, "metadata": { "collapsed": false, "scrolled": true @@ -468,17 +467,17 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 3.5736021e-03 \n", - " 3 | 3.5273328e-04 \n", - " 4 | 3.5227281e-05 \n", - " 5 | 3.5222757e-06 \n", - " 6 | 3.5273431e-07 \n", - " 7 | 3.6121307e-08 \n", - " 8 | 9.3608476e-08 \n", - " 9 | 1.1241399e-06 \n", - " 10 | 5.4816871e-06 \n", - " 11 | 7.0984846e-05 \n", - " 12 | 4.5614270e-04 \n", + " 2 | 8.9340052e-04 \n", + " 3 | 8.8183320e-05 \n", + " 4 | 8.8068202e-06 \n", + " 5 | 8.8056893e-07 \n", + " 6 | 8.8183577e-08 \n", + " 7 | 9.0303267e-09 \n", + " 8 | 2.3402119e-08 \n", + " 9 | 2.8103497e-07 \n", + " 10 | 1.3704218e-06 \n", + " 11 | 1.7746211e-05 \n", + " 12 | 1.1403567e-04 \n", "-----------------------------\n" ] } @@ -497,9 +496,6 @@ "tbm.nkpoints = (0,0,0)\n", "\n", "X = positions(at)\n", - "# X[2,2] += 1.0e-15\n", - "# X[3,2] -= 1.0e-15\n", - "# set_positions!(at, X)\n", "Hess = TightBinding.hessian(at, tbm)\n", "f = reshape(Hess, 3*length(at), 3*length(at))\n", "Tens = TightBinding.d3E(at, tbm)\n", @@ -618,7 +614,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 20, "metadata": { "collapsed": false, "scrolled": true @@ -641,17 +637,17 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 1.5930766e-03 \n", - " 3 | 1.5576992e-04 \n", - " 4 | 1.5541830e-05 \n", - " 5 | 1.5538196e-06 \n", - " 6 | 1.5564084e-07 \n", - " 7 | 1.6554841e-08 \n", - " 8 | 1.6653345e-08 \n", - " 9 | 2.5295189e-07 \n", - " 10 | 1.6344273e-06 \n", - " 11 | 3.1332893e-05 \n", - " 12 | 3.3386867e-04 \n", + " 2 | 9.5899932e-05 \n", + " 3 | 9.5900283e-06 \n", + " 4 | 9.5900611e-07 \n", + " 5 | 9.5903936e-08 \n", + " 6 | 9.7450780e-09 \n", + " 7 | 4.1882209e-09 \n", + " 8 | 4.8154754e-08 \n", + " 9 | 2.9517938e-07 \n", + " 10 | 4.7203659e-06 \n", + " 11 | 3.4545560e-05 \n", + " 12 | 2.3161015e-04 \n", "-----------------------------\n" ] } @@ -710,7 +706,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 21, "metadata": { "collapsed": false, "scrolled": true @@ -733,17 +729,17 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 1.3503005e+02 \n", - " 3 | 1.7916602e-04 \n", - " 4 | 1.7902184e-05 \n", - " 5 | 1.7900218e-06 \n", - " 6 | 1.7952114e-07 \n", - " 7 | 1.5319154e-08 \n", - " 8 | 7.3291674e-08 \n", - " 9 | 6.1528752e-07 \n", - " 10 | 1.3937964e-05 \n", - " 11 | 6.9539720e-05 \n", - " 12 | 8.7991192e-04 \n", + " 2 | 4.9673185e-04 \n", + " 3 | 1.3000268e+03 \n", + " 4 | 1.3001546e+04 \n", + " 5 | 1.3001533e+05 \n", + " 6 | 1.3001548e+06 \n", + " 7 | 1.3001546e+07 \n", + " 8 | 2.4740555e-07 \n", + " 9 | 1.3001546e+09 \n", + " 10 | 2.9972865e-05 \n", + " 11 | 1.3001546e+11 \n", + " 12 | 2.9975990e-03 \n", "-----------------------------\n" ] } @@ -803,7 +799,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 22, "metadata": { "collapsed": false }, @@ -825,15 +821,15 @@ "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 8.0597075e-04 \n", - " 3 | 8.0563563e-05 \n", - " 4 | 8.0563202e-06 \n", - " 5 | 8.0566677e-07 \n", - " 6 | 8.0383845e-08 \n", - " 7 | 8.8819055e-09 \n", - " 8 | 2.2204055e-07 \n", - " 9 | 2.2204443e-06 \n", - " 10 | 1.3322688e-05 \n", + " 2 | 1.9246636e-04 \n", + " 3 | 1.9247712e-05 \n", + " 4 | 1.9247743e-06 \n", + " 5 | 1.9250385e-07 \n", + " 6 | 1.9726388e-08 \n", + " 7 | 3.6830786e-09 \n", + " 8 | 2.3546291e-08 \n", + " 9 | 2.2070277e-07 \n", + " 10 | 3.3293272e-06 \n", "-----------------------------\n" ] } @@ -851,7 +847,7 @@ "\n", "set_pbc!(at, [true, true, true])\n", "tbm = NRLTB.NRLTBModel(elem = NRLTB.Al_spd)\n", - "tbm.nkpoints = (2,8,2)\n", + "tbm.nkpoints = (0,2,4)\n", "\n", "\n", "X = positions(at)\n", @@ -885,7 +881,7 @@ }, { "cell_type": "code", - "execution_count": 63, + "execution_count": 24, "metadata": { "collapsed": false }, @@ -894,22 +890,22 @@ "name": "stdout", "output_type": "stream", "text": [ - "4" + "6" ] } ], "source": [ - "at = bulk(\"Al\")#; cubic=true)\n", - "at = repeat(at, (1, 2, 2))\n", + "at = bulk(\"Al\") #; cubic=true)\n", + "at = repeat(at, (1, 2, 3))\n", "X = positions(at)\n", - "# set_pbc!(at, [false, false, false])\n", + "set_pbc!(at, [false, false, false])\n", "# plot3D(X[1,:][:], X[2,:][:], X[3,:][:], \"b.\")\n", "print(length(at))" ] }, { "cell_type": "code", - "execution_count": 64, + "execution_count": 37, "metadata": { "collapsed": false }, @@ -928,30 +924,30 @@ "name": "stdout", "output_type": "stream", "text": [ - "E - ∑ E_i = 5.329070518200751e-15\n", + "E - ∑ E_i = -1.0658141036401503e-14\n", "------------------------------\n", "Finite-difference test\n", "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 1.8566967e-04 \n", - " 3 | 1.8596655e-05 \n", - " 4 | 1.8600000e-06 \n", - " 5 | 1.8629885e-07 \n", - " 6 | 2.0875614e-08 \n", - " 7 | 3.4353978e-08 \n", - " 8 | 2.5728062e-07 \n", - " 9 | 3.7453059e-06 \n", - " 10 | 2.9216484e-05 \n", - " 11 | 2.1344333e-04 \n", - " 12 | 3.4997035e-03 \n", + " 2 | 1.8723765e-04 \n", + " 3 | 1.8753995e-05 \n", + " 4 | 1.8756845e-06 \n", + " 5 | 1.8734170e-07 \n", + " 6 | 1.6832353e-08 \n", + " 7 | 3.6902442e-08 \n", + " 8 | 3.1859084e-07 \n", + " 9 | 4.5285329e-06 \n", + " 10 | 3.5460745e-05 \n", + " 11 | 4.1085812e-04 \n", + " 12 | 5.0581097e-03 \n", "-----------------------------\n" ] } ], "source": [ "# TEST NRL-TB site energy for Aluminum FCC\n", - "# WITH open BOUNDARY CONDITION ON THIRD DIMENSION\n", + "# WITH open BOUNDARY CONDITION \n", "\n", "using AtomsInterface\n", "reload(\"SparseTools\")\n", @@ -993,7 +989,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 36, "metadata": { "collapsed": false }, @@ -1012,30 +1008,30 @@ "name": "stdout", "output_type": "stream", "text": [ - "E - ∑ E_i = 0.0\n", + "E - ∑ E_i = -1.7763568394002505e-15\n", "------------------------------\n", "Finite-difference test\n", "-----------------------------\n", " p | error \n", "----|------------------------\n", - " 2 | 1.0063240e-03 \n", - " 3 | 1.0064698e-04 \n", - " 4 | 1.0064709e-05 \n", - " 5 | 1.0063949e-06 \n", - " 6 | 9.9476247e-08 \n", - " 7 | 2.2204517e-08 \n", - " 8 | 2.6645355e-07 \n", - " 9 | 1.7763569e-06 \n", - " 10 | 1.3322676e-05 \n", - " 11 | 3.1086245e-04 \n", - " 12 | 2.2204460e-03 \n", + " 2 | 1.0275612e-04 \n", + " 3 | 1.0275015e-05 \n", + " 4 | 1.0274886e-06 \n", + " 5 | 1.0271511e-07 \n", + " 6 | 9.9161716e-09 \n", + " 7 | 8.4521695e-09 \n", + " 8 | 8.2039690e-08 \n", + " 9 | 9.7064306e-07 \n", + " 10 | 9.0142620e-06 \n", + " 11 | 8.5313861e-05 \n", + " 12 | 7.9791681e-04 \n", "-----------------------------\n" ] } ], "source": [ "# TEST NRL-TB site energy for Aluminum FCC\n", - "# WITH periodic BOUNDARY CONDITION ON THIRD DIMENSION\n", + "# WITH periodic BOUNDARY CONDITION ON FIRST DIMENSION\n", "\n", "using AtomsInterface\n", "reload(\"SparseTools\")\n", @@ -1043,9 +1039,9 @@ "reload(\"TightBinding\")\n", "reload(\"NRLTB\") \n", "\n", - "set_pbc!(at, [false, true, true])\n", + "set_pbc!(at, [true, false, false])\n", "tbm = NRLTB.NRLTBModel(elem = NRLTB.Al_spd)\n", - "tbm.nkpoints = (0,0,0)\n", + "tbm.nkpoints = (0,0,4)\n", "\n", "X = positions(at)\n", "\n", @@ -1807,7 +1803,7 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "0.4.0" + "version": "0.4.1" } }, "nbformat": 4,