Coverage for gpaw/lcao/tightbinding.py: 73%
109 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1"""Simple tight-binding add-on for GPAWs LCAO module."""
3from math import pi
5import numpy as np
6import scipy.linalg as sla
8import ase.units as units
10from gpaw.utilities.tools import tri2full
13class TightBinding:
14 """Simple class for tight-binding calculations."""
16 def __init__(self, atoms, calc):
17 """Init with ``Atoms`` and a converged LCAO calculation."""
19 # Store and extract useful attributes from the calculator
20 self.atoms = atoms
21 self.calc = calc
22 self.kd = calc.wfs.kd
23 wfs = calc.wfs
24 kd = wfs.kd
26 # Matrix size
27 self.nao = wfs.setups.nao
28 # K-point info
29 self.gamma = kd.gamma
30 if self.gamma:
31 self.Nk_c = (1, 1, 1)
32 else:
33 self.Nk_c = tuple(kd.N_c)
35 self.ibzk_kc = kd.ibzk_kc
36 self.ibzk_qc = kd.ibzk_qc
37 self.bzk_kc = kd.bzk_kc
39 # Symmetry
40 self.symmetry = kd.symmetry
41 if self.symmetry.point_group:
42 raise NotImplementedError("Only time-reversal symmetry supported.")
44 # Lattice vectors and number of repetitions
45 self.R_cN = None
46 self.N_c = None
48 # Init with default number of real-space cells
49 self.set_num_cells()
51 def set_num_cells(self, N_c=None):
52 """Set number of real-space cells to use.
54 Parameters
55 ----------
56 N_c: tuple or ndarray
57 Number of unit cells in each direction of the basis vectors.
59 """
61 if N_c is None:
62 self.N_c = tuple(self.Nk_c)
63 else:
64 self.N_c = tuple(N_c)
66 if np.any(np.asarray(self.Nk_c) < np.asarray(self.N_c)):
67 print("WARNING: insufficient k-point sampling.")
69 # Lattice vectors
70 R_cN = np.indices(self.N_c).reshape(3, -1)
71 N_c = np.array(self.N_c)[:, np.newaxis]
72 R_cN += N_c // 2
73 R_cN %= N_c
74 R_cN -= N_c // 2
75 self.R_cN = R_cN
77 def lattice_vectors(self):
78 """Return real-space lattice vectors."""
80 return self.R_cN
82 def bloch_to_real_space(self, A_qxMM, R_c=None):
83 """Transform quantity from Bloch to real-space representation.
85 Parameters
86 ----------
87 A_qxMM: ndarray
88 Bloch representation of matrix. May be parallelized over k-points.
89 R_cN: ndarray
90 Cell vectors for which the real-space matrices will be calculated
91 and returned.
93 """
95 # Include all cells per default
96 if R_c is None:
97 R_Nc = self.R_cN.transpose()
98 else:
99 R_Nc = [R_c]
101 # Real-space quantities
102 A_NxMM = []
103 # Reshape input array
104 shape = A_qxMM.shape
105 A_qx = A_qxMM.reshape(shape[0], -1)
107 # Fourier transform to real-space
108 for R_c in R_Nc:
109 # Evaluate fourier sum
110 phase_q = np.exp(2.j * pi * np.dot(self.ibzk_qc, R_c))
111 A_x = np.sum(phase_q[:, np.newaxis] * A_qx, axis=0)
112 self.kd.comm.sum(A_x)
113 A_xMM = A_x.reshape(shape[1:])
115 # Time-reversal symmetry
116 if not len(self.ibzk_kc) == len(self.bzk_kc):
117 # Broadcast Gamma component
118 gamma = np.where(abs(self.ibzk_kc).sum(axis=1) < 1e-13)[0]
119 if gamma.size > 0:
120 rank, myu = self.kd.get_rank_and_index(gamma[0])
122 if self.kd.comm.rank == rank:
123 A0_xMM = A_qxMM[myu]
124 else:
125 A0_xMM = np.zeros_like(A_xMM)
126 self.kd.comm.broadcast(A0_xMM, rank)
127 else:
128 A0_xMM = 0.
129 # Add conjugate and subtract double counted Gamma component
130 A_xMM += A_xMM.conj() - A0_xMM
132 A_xMM /= np.prod(self.Nk_c)
134 try:
135 assert np.all(np.abs(A_xMM.imag) < 1e-10)
136 except AssertionError:
137 raise ValueError("MAX Im(A_MM): % .2e" %
138 np.amax(np.abs(A_xMM.imag)))
140 A_NxMM.append(A_xMM.real)
142 return np.array(A_NxMM)
144 def h_and_s(self):
145 """Return LCAO Hamiltonian and overlap matrix in real-space."""
147 # Extract Bloch Hamiltonian and overlap matrix
148 H_kMM = []
149 S_kMM = []
151 h = self.calc.hamiltonian
152 wfs = self.calc.wfs
153 kpt_u = wfs.kpt_u
155 for kpt in kpt_u:
156 H_MM = wfs.eigensolver.calculate_hamiltonian_matrix(h, wfs, kpt)
157 S_MM = wfs.S_qMM[kpt.q]
158 # XXX Converting to full matrices here
159 tri2full(H_MM)
160 tri2full(S_MM)
161 H_kMM.append(H_MM)
162 S_kMM.append(S_MM)
164 # Convert to arrays
165 H_kMM = np.array(H_kMM)
166 S_kMM = np.array(S_kMM)
168 H_NMM = self.bloch_to_real_space(H_kMM)
169 S_NMM = self.bloch_to_real_space(S_kMM)
171 return H_NMM, S_NMM
173 def band_structure(self, path_kc, blochstates=False):
174 """Calculate dispersion along a path in the Brillouin zone.
176 Parameters
177 ----------
178 path_kc: ndarray
179 List of k-point coordinates (in units of the reciprocal lattice
180 vectors) specifying the path in the Brillouin zone for which the
181 dynamical matrix will be calculated.
182 blochstates: bool
183 Return LCAO expansion coefficients when True (default: False).
185 """
187 # Real-space matrices
188 self.H_NMM, self.S_NMM = self.h_and_s()
190 assert self.H_NMM is not None
191 assert self.S_NMM is not None
193 # Lattice vectors
194 R_cN = self.R_cN
196 # Lists for eigenvalues and eigenvectors along path
197 eps_kn = []
198 psi_kn = []
200 for k_c in path_kc:
201 # Evaluate fourier sum
202 phase_N = np.exp(-2.j * pi * np.dot(k_c, R_cN))
203 H_MM = np.sum(phase_N[:, np.newaxis, np.newaxis] * self.H_NMM,
204 axis=0)
205 S_MM = np.sum(phase_N[:, np.newaxis, np.newaxis] * self.S_NMM,
206 axis=0)
208 if blochstates:
209 eps_n, c_Mn = sla.eigh(H_MM, S_MM)
210 # Sort eigenmodes according to increasing eigenvalues
211 c_nM = c_Mn[:, eps_n.argsort()].transpose()
212 psi_kn.append(c_nM)
213 else:
214 eps_n = sla.eigvalsh(H_MM, S_MM)
216 # Sort eigenvalues in increasing order
217 eps_n.sort()
218 eps_kn.append(eps_n)
220 # Convert to eV
221 eps_kn = np.array(eps_kn) * units.Hartree
223 if blochstates:
224 return eps_kn, np.array(psi_kn)
226 return eps_kn