#!/usr/bin/env python3 """ LEECH LATTICE 24D MINIMAL-SHELL PHONON DYNAMICS ================================================ Exact finite-shell harmonic dynamics for the Leech lattice. Normalization ------------- Minimal Leech vectors satisfy ||R||^2 = 4 and there are exactly 196,560 such vectors. Potential --------- Central pair potential F(U) = U^(-2), U = ||R||^2. At U = 4, F'(4) = -1/32 F''(4) = 3/128 giving kappa_L = 5/16 kappa_T = -1/16. For a minimal vector R, H(R) = kappa_T I + ((kappa_L-kappa_T)/4) R R^T. Bloch dynamical matrix ---------------------- Using the Bloch character exp(2 pi i ), the real symmetric dynamical matrix is D(q) = sum_R [1-cos(2*pi*q.R)] [ ((kappa_L-kappa_T)/4) R R^T + kappa_T I ]. Because kappa_T < 0, this finite-shell model is mechanically unstable. Negative eigenvalues are physical results of this chosen model and correspond to imaginary harmonic frequencies: omega = i sqrt(-lambda). Deep-hole representatives -------------------------- 24A1: q = (1,1,...,1)/sqrt(8) D24: q = (4,0,...,0)/sqrt(8) The program verifies: * Extended binary Golay code size = 4096 * 759 octads * Family A = 1104 * Family B = 97152 * Family C = 98304 * Total shell = 196560 * Every minimal vector has norm squared 4 * No duplicate vectors * Inversion symmetry * Full-shell second moment * Family-C second moment * Gamma-point D(0)=0 * 24A1 phase structure * 24A1 scalar dynamical matrix * D24 axial symmetry * Acoustic-limit behavior * Deep-hole spectra * Gamma -> 24A1 -> D24 -> Gamma dispersion path IMPORTANT --------- This is an exact minimal-shell model, not a full all-shell Leech-lattice interaction model. """ import math import time import numpy as np class LeechMinimalShellDynamics: # ================================================================== # INITIALIZATION # ================================================================== def __init__(self, path_points=25, verbose=True): self.path_points = int(path_points) self.verbose = bool(verbose) self.n = 24 # Minimal Leech vectors have squared norm 4. self.shell_norm2 = 4.0 # Integer Leech construction vectors have squared norm 32. self.scale = 1.0 / math.sqrt(8.0) # -------------------------------------------------------------- # Pair potential F(U) = U^(-2) # -------------------------------------------------------------- U = self.shell_norm2 self.F_prime = -2.0 / U**3 self.F_double_prime = 6.0 / U**4 # Central-force Hessian coefficients: # # kappa_L = 4 U F''(U) + 2 F'(U) # kappa_T = 2 F'(U) self.kappa_L = ( 4.0 * U * self.F_double_prime + 2.0 * self.F_prime ) self.kappa_T = ( 2.0 * self.F_prime ) # For ||R||^2 = 4: # # H(R) # = kappa_T I # + ((kappa_L-kappa_T)/4) R R^T self.RR_coefficient = ( self.kappa_L - self.kappa_T ) / 4.0 self.I_coefficient = self.kappa_T # -------------------------------------------------------------- # Header # -------------------------------------------------------------- self._print("=" * 80) self._print( "LEECH LATTICE 24D MINIMAL-SHELL PHONON DYNAMICS" ) self._print("=" * 80) # -------------------------------------------------------------- # 1. Golay code # -------------------------------------------------------------- self._print( "\n[1/5] Constructing extended binary Golay code..." ) t0 = time.time() self.codewords = self._build_golay_code() self._print( f" Codewords: {len(self.codewords):,} " f"({time.time() - t0:.3f} s)" ) # -------------------------------------------------------------- # 2. Minimal shell # -------------------------------------------------------------- self._print( "\n[2/5] Enumerating all 196,560 minimal Leech vectors..." ) t0 = time.time() ( self.shell1, self.family_A, self.family_B, self.family_C, ) = self._build_shell1() self._print( f" Family A: {len(self.family_A):,} vectors" ) self._print( f" Family B: {len(self.family_B):,} vectors" ) self._print( f" Family C: {len(self.family_C):,} vectors" ) self._print( f" Shell size: {len(self.shell1):,} vectors " f"({time.time() - t0:.3f} s)" ) # -------------------------------------------------------------- # 3. Shell geometry # -------------------------------------------------------------- self._print( "\n[3/5] Verifying shell geometry..." ) self._verify_shell() # -------------------------------------------------------------- # 4. Deep-hole representatives # -------------------------------------------------------------- self._print( "\n[4/5] Defining Conway-Sloane deep-hole points..." ) self.k_gamma = np.zeros( 24, dtype=np.float64 ) # 24A1 representative. self.k_24A1 = ( np.ones( 24, dtype=np.float64 ) / math.sqrt(8.0) ) # D24 representative. self.k_D24 = np.zeros( 24, dtype=np.float64 ) self.k_D24[0] = ( 4.0 / math.sqrt(8.0) ) self._print( f" ||q_24A1||^2 = " f"{np.dot(self.k_24A1, self.k_24A1):.12f}" ) self._print( f" ||q_D24||^2 = " f"{np.dot(self.k_D24, self.k_D24):.12f}" ) # -------------------------------------------------------------- # 5. Dynamical model # -------------------------------------------------------------- self._print( "\n[5/5] Verifying dynamical model..." ) self._verify_dynamical_model() self._print( "\nInitialization complete." ) # ================================================================== # PRINT # ================================================================== def _print(self, message): if self.verbose: print(message) # ================================================================== # GOLAY CODE # ================================================================== def _build_golay_code(self): """ Construct the extended binary Golay [24,12,8] code. Generator matrix is written as G = [I | B]. Returns ------- codewords : ndarray, shape (4096,24) Binary codewords. """ b_row = [ 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0 ] B = np.zeros( (12, 12), dtype=np.uint8 ) # 11 x 11 circulant block. for i in range(11): for j in range(11): B[i, j] = ( b_row[(j - i) % 11] ) B[i, 11] = 1 # Final row. B[11, :11] = 1 B[11, 11] = 0 G = np.hstack( [ np.eye( 12, dtype=np.uint8 ), B ] ) # All 4096 messages. messages = np.array( [ [ (m >> i) & 1 for i in range(12) ] for m in range(4096) ], dtype=np.uint8 ) codewords = ( messages @ G ) % 2 if codewords.shape != (4096, 24): raise RuntimeError( "Golay code construction produced " f"shape {codewords.shape}, expected (4096,24)." ) # Verify minimum-distance structure. weights = np.sum( codewords, axis=1 ) nonzero_weights = sorted( set( int(w) for w in weights if w != 0 ) ) expected_weights = [ 8, 12, 16, 24 ] if nonzero_weights != expected_weights: raise RuntimeError( "Unexpected extended Golay weight spectrum: " f"{nonzero_weights}" ) return codewords # ================================================================== # SHELL CONSTRUCTION # ================================================================== def _build_shell1(self): """ Construct the 196,560 minimal Leech vectors. Family A: (+/-4,+/-4,0^22)/sqrt(8) 1104 vectors. Family B: (+/-2^8,0^16)/sqrt(8) supported on Golay octads with product of signs +1. 97152 vectors. Family C: one coordinate +/-3 and 23 coordinates +/-1, generated from Golay codewords. 98304 vectors. The critical Family-C construction is v[i] = base[i] - 4*base[i], which maps +1 -> -3 -1 -> +3. Simply doing v[i] -= 4 is incorrect because it can produce -5 when base[i] = -1. """ family_A = [] family_B = [] family_C = [] scale = self.scale # ============================================================== # FAMILY A # # (+/-4,+/-4,0^22)/sqrt(8) # # 4 * C(24,2) = 1104. # ============================================================== for i in range(24): for j in range(i + 1, 24): for s1 in (-4, 4): for s2 in (-4, 4): v = np.zeros( 24, dtype=np.float64 ) v[i] = s1 v[j] = s2 family_A.append( v * scale ) # ============================================================== # FAMILY B # # (+/-2^8,0^16)/sqrt(8) # # Support is a Golay octad. # Signs have product +1. # # 759 * 128 = 97152. # ============================================================== octads = [ cw for cw in self.codewords if int(np.sum(cw)) == 8 ] if len(octads) != 759: raise RuntimeError( f"Expected 759 octads, found {len(octads)}." ) for octad in octads: idx = np.where( octad == 1 )[0] if len(idx) != 8: raise RuntimeError( "Invalid Golay octad." ) # 2^7 sign patterns with even number of minus signs. for p in range(128): signs = np.ones( 8, dtype=np.float64 ) minus_parity = 0 for b in range(7): if (p >> b) & 1: signs[b] = -1.0 minus_parity ^= 1 # Enforce product +1. if minus_parity: signs[7] = -1.0 v = np.zeros( 24, dtype=np.float64 ) v[idx] = 2.0 * signs family_B.append( v * scale ) # ============================================================== # FAMILY C # # (+/-3,+/-1^23)/sqrt(8) # # 4096 * 24 = 98304. # # IMPORTANT: # # base[i] is either +1 or -1. # # Distinguished coordinate must transform as: # # +1 -> -3 # -1 -> +3 # # Therefore: # # v[i] = base[i] - 4*base[i]. # # The previous erroneous expression # # v[i] = base[i] - 4 # # produces -5 whenever base[i] = -1. # ============================================================== for c in self.codewords: base = ( 1.0 - 2.0 * c.astype( np.float64 ) ) for i in range(24): v = base.copy() # Correct Family-C transformation. v[i] -= 4.0 * base[i] family_C.append( v * scale ) # Convert to arrays. family_A = np.asarray( family_A, dtype=np.float64 ) family_B = np.asarray( family_B, dtype=np.float64 ) family_C = np.asarray( family_C, dtype=np.float64 ) shell1 = np.vstack( [ family_A, family_B, family_C ] ) return ( shell1, family_A, family_B, family_C ) # ================================================================== # SHELL VERIFICATION # ================================================================== def _verify_shell(self): shell = self.shell1 expected = 196560 # -------------------------------------------------------------- # Count # -------------------------------------------------------------- if len(shell) != expected: raise RuntimeError( f"Expected {expected} shell vectors, " f"found {len(shell)}." ) # -------------------------------------------------------------- # Norms # -------------------------------------------------------------- norms2 = np.einsum( "ij,ij->i", shell, shell ) max_norm_error = float( np.max( np.abs( norms2 - 4.0 ) ) ) self._print( f" max | ||R||^2 - 4 | = " f"{max_norm_error:.3e}" ) if max_norm_error > 1e-10: # Helpful diagnostic. bad = np.where( np.abs(norms2 - 4.0) > 1e-10 )[0] if len(bad) > 0: sample = shell[ bad[0] ] self._print( " First bad vector:" ) self._print( f" {sample}" ) self._print( f" norm^2 = {norms2[bad[0]]}" ) raise RuntimeError( "Minimal-shell enumeration contains " "vectors with incorrect norm." ) # -------------------------------------------------------------- # Uniqueness # -------------------------------------------------------------- integer_shell = np.rint( shell / self.scale ).astype(np.int16) unique_vectors = np.unique( integer_shell, axis=0 ) self._print( f" Unique shell vectors: " f"{len(unique_vectors):,}" ) if len(unique_vectors) != expected: raise RuntimeError( "Duplicate minimal vectors detected." ) # -------------------------------------------------------------- # Inversion symmetry # -------------------------------------------------------------- integer_set = { tuple(row) for row in integer_shell } inversion_ok = all( tuple(-row) in integer_set for row in integer_shell ) self._print( " Inversion symmetry: " + ( "verified." if inversion_ok else "FAILED." ) ) if not inversion_ok: raise RuntimeError( "Shell is not inversion symmetric." ) # -------------------------------------------------------------- # Full-shell second moment # # Sum_R R R^T = # # (N * 4 / 24) I # # = 32760 I. # -------------------------------------------------------------- full_second = ( shell.T @ shell ) expected_full = ( len(shell) * 4.0 / 24.0 ) full_error = float( np.max( np.abs( full_second - expected_full * np.eye(24) ) ) ) self._print( f" full-shell second-moment error: " f"{full_error:.3e}" ) if full_error > 1e-8: raise RuntimeError( "Full-shell second moment is not isotropic." ) # -------------------------------------------------------------- # Family-C second moment # # 98304 * 4 / 24 = 16384. # -------------------------------------------------------------- C_second = ( self.family_C.T @ self.family_C ) expected_C = ( len(self.family_C) * 4.0 / 24.0 ) C_error = float( np.max( np.abs( C_second - expected_C * np.eye(24) ) ) ) self._print( f" Family-C second-moment error: " f"{C_error:.3e}" ) if C_error > 1e-8: raise RuntimeError( "Family-C second moment is not isotropic." ) self._print( " Shell geometry: VERIFIED." ) # ================================================================== # DYNAMICAL MATRIX # ================================================================== def compute_dynamical_matrix(self, q): """ Compute D(q) = sum_R [1-cos(2*pi*q.R)] [ ((kappa_L-kappa_T)/4) R R^T + kappa_T I ]. The result is a real symmetric 24x24 matrix. """ q = np.asarray( q, dtype=np.float64 ) if q.shape != (24,): raise ValueError( "q must have shape (24,)." ) # Bloch phase. phase = ( 2.0 * math.pi * (self.shell1 @ q) ) modulation = ( 1.0 - np.cos(phase) ) # Weighted second moment: # # sum_R modulation_R R R^T R_outer_sum = ( self.shell1.T @ ( self.shell1 * modulation[:, None] ) ) # Scalar modulation. scalar_sum = float( np.sum(modulation) ) # IMPORTANT: # # No /24 factor appears in the identity contribution. # # D = # c_RR * R_outer_sum # + kappa_T * scalar_sum * I D = ( self.RR_coefficient * R_outer_sum + self.I_coefficient * scalar_sum * np.eye(24) ) # Remove tiny floating-point antisymmetry. D = 0.5 * ( D + D.T ) return D # ================================================================== # PHASE MODULATION # ================================================================== def phase_modulation(self, q): q = np.asarray( q, dtype=np.float64 ) phase = ( 2.0 * math.pi * (self.shell1 @ q) ) return ( 1.0 - np.cos(phase) ) # ================================================================== # DYNAMICAL MODEL VERIFICATION # ================================================================== def _verify_dynamical_model(self): self._print( f" kappa_L = " f"{self.kappa_L:.12f}" ) self._print( f" kappa_T = " f"{self.kappa_T:.12f}" ) self._print( f" RR coefficient = " f"{self.RR_coefficient:.12f}" ) self._print( f" I coefficient = " f"{self.I_coefficient:.12f}" ) # -------------------------------------------------------------- # Gamma # -------------------------------------------------------------- D0 = self.compute_dynamical_matrix( self.k_gamma ) gamma_error = float( np.max( np.abs(D0) ) ) self._print( f" Gamma-point ||D(0)||_max = " f"{gamma_error:.3e}" ) if gamma_error > 1e-12: raise RuntimeError( "Gamma point is not zero." ) # -------------------------------------------------------------- # 24A1 phases # -------------------------------------------------------------- modulation = ( self.phase_modulation( self.k_24A1 ) ) nA = len( self.family_A ) nB = len( self.family_B ) mod_A = modulation[ :nA ] mod_B = modulation[ nA:nA + nB ] mod_C = modulation[ nA + nB: ] A_error = float( np.max( np.abs(mod_A) ) ) B_error = float( np.max( np.abs(mod_B) ) ) C_error = float( np.max( np.abs( mod_C - 2.0 ) ) ) self._print( f" 24A1 Family-A cancellation error: " f"{A_error:.3e}" ) self._print( f" 24A1 Family-B cancellation error: " f"{B_error:.3e}" ) self._print( f" 24A1 Family-C phase error: " f"{C_error:.3e}" ) if A_error > 1e-12: raise RuntimeError( "Family A does not cancel at 24A1." ) if B_error > 1e-12: raise RuntimeError( "Family B does not cancel at 24A1." ) if C_error > 1e-12: raise RuntimeError( "Family C does not have phase factor 2." ) # -------------------------------------------------------------- # 24A1 scalar matrix # -------------------------------------------------------------- D24A1 = ( self.compute_dynamical_matrix( self.k_24A1 ) ) offdiag = ( D24A1 - np.diag( np.diag(D24A1) ) ) offdiag_error = float( np.max( np.abs(offdiag) ) ) diagonal = np.diag( D24A1 ) scalar_block_error = float( np.max( np.abs( diagonal - np.mean(diagonal) ) ) ) self._print( f" 24A1 off-diagonal error: " f"{offdiag_error:.3e}" ) self._print( f" 24A1 scalar-block error: " f"{scalar_block_error:.3e}" ) if offdiag_error > 1e-8: raise RuntimeError( "24A1 matrix is not diagonal." ) if scalar_block_error > 1e-8: raise RuntimeError( "24A1 matrix is not scalar." ) self._print( " Dynamical model: VERIFIED." ) # ================================================================== # EIGENVALUES # ================================================================== def eigenvalues(self, q): D = self.compute_dynamical_matrix( q ) return np.linalg.eigvalsh(D) # ================================================================== # FREQUENCY CLASSIFICATION # ================================================================== @staticmethod def classify_frequencies(eigenvalues): eigenvalues = np.asarray( eigenvalues, dtype=np.float64 ) tolerance = 1e-10 real_frequencies = np.sqrt( np.maximum( eigenvalues, 0.0 ) ) imaginary_frequencies = np.sqrt( np.maximum( -eigenvalues, 0.0 ) ) # Exact numerical zeros. zero_mask = ( np.abs(eigenvalues) < tolerance ) real_frequencies[ zero_mask ] = 0.0 imaginary_frequencies[ zero_mask ] = 0.0 return ( real_frequencies, imaginary_frequencies ) # ================================================================== # DEEP-HOLE REPORT # ================================================================== def report_deep_hole( self, name, q ): eigs = self.eigenvalues( q ) ( real_freqs, imag_freqs ) = self.classify_frequencies( eigs ) self._print( f"\n {name}" ) self._print( f" q norm^2 = " f"{np.dot(q, q):.12f}" ) self._print( f" eigenvalue minimum = " f"{np.min(eigs):.12f}" ) self._print( f" eigenvalue maximum = " f"{np.max(eigs):.12f}" ) if np.min(eigs) < -1e-10: self._print( " imaginary frequency magnitude = " f"[{np.min(imag_freqs):.12f}, " f"{np.max(imag_freqs):.12f}]" ) elif np.max(eigs) > 1e-10: self._print( " real frequency range = " f"[{np.min(real_freqs):.12f}, " f"{np.max(real_freqs):.12f}]" ) else: self._print( " all eigenvalues are numerically zero." ) # -------------------------------------------------------------- # Eigenvalue multiplicities # -------------------------------------------------------------- sorted_eigs = np.sort( eigs ) groups = [] current = [ sorted_eigs[0] ] for value in sorted_eigs[1:]: if abs( value - current[-1] ) < 1e-7: current.append( value ) else: groups.append( current ) current = [ value ] groups.append( current ) multiplicities = [ len(group) for group in groups ] self._print( f" approximate eigenvalue multiplicities: " f"{multiplicities}" ) return eigs # ================================================================== # EXACT 24A1 CHECK # ================================================================== def exact_24A1_check(self): D = self.compute_dynamical_matrix( self.k_24A1 ) lambda_scalar = ( np.trace(D) / 24.0 ) expected = -9216.0 error = abs( lambda_scalar - expected ) self._print( "\n[3] Exact 24A1 scalar check" ) self._print( f" lambda_scalar = " f"{lambda_scalar:.12f}" ) self._print( f" expected exact value = " f"{expected:.12f}" ) self._print( f" absolute error = " f"{error:.3e}" ) if error > 1e-8: raise RuntimeError( "Exact 24A1 scalar check failed." ) # ================================================================== # D24 AXIAL SYMMETRY # ================================================================== def axial_symmetry_check(self): D = self.compute_dynamical_matrix( self.k_D24 ) offdiag = ( D - np.diag( np.diag(D) ) ) offdiag_error = float( np.max( np.abs(offdiag) ) ) transverse = np.diag( D )[1:] transverse_spread = ( np.max(transverse) - np.min(transverse) ) longitudinal = D[0, 0] transverse_value = np.mean( transverse ) self._print( "\n[4] D24 axial symmetry check" ) self._print( f" max |D_ij|, i != j = " f"{offdiag_error:.3e}" ) self._print( f" transverse diagonal spread = " f"{transverse_spread:.3e}" ) self._print( f" longitudinal eigenvalue = " f"{longitudinal:.12f}" ) self._print( f" transverse eigenvalue = " f"{transverse_value:.12f}" ) # ================================================================== # ACOUSTIC LIMIT # ================================================================== def acoustic_limit( self, dk=1e-5 ): """ Examine q = dk * e_1. For lambda ~ c * k^2: c > 0: omega/k = sqrt(c) c < 0: omega/k = i sqrt(-c) The current U^(-2) minimal-shell model has c < 0. """ q = np.zeros( 24, dtype=np.float64 ) q[0] = dk eigs = self.eigenvalues( q ) # Along a coordinate axis: # # 23 transverse eigenvalues # 1 longitudinal eigenvalue transverse_lambda = float( np.mean( eigs[:23] ) ) longitudinal_lambda = float( eigs[23] ) transverse_coeff = ( transverse_lambda / dk**2 ) longitudinal_coeff = ( longitudinal_lambda / dk**2 ) self._print( "\n[1] Acoustic limit near Gamma" ) self._print( f" transverse lambda/k^2 = " f"{transverse_coeff:.12f}" ) self._print( f" longitudinal lambda/k^2 = " f"{longitudinal_coeff:.12f}" ) # -------------------------------------------------------------- # Transverse # -------------------------------------------------------------- if transverse_coeff < 0: vT_imag = math.sqrt( -transverse_coeff ) self._print( f" transverse mode is unstable: " f"omega/k = i {vT_imag:.12f}" ) else: vT = math.sqrt( transverse_coeff ) self._print( f" transverse sound velocity = " f"{vT:.12f}" ) # -------------------------------------------------------------- # Longitudinal # -------------------------------------------------------------- if longitudinal_coeff < 0: vL_imag = math.sqrt( -longitudinal_coeff ) self._print( f" longitudinal mode is unstable: " f"omega/k = i {vL_imag:.12f}" ) else: vL = math.sqrt( longitudinal_coeff ) self._print( f" longitudinal sound velocity = " f"{vL:.12f}" ) # -------------------------------------------------------------- # Ratio # -------------------------------------------------------------- if ( transverse_coeff < 0 and longitudinal_coeff < 0 ): ratio = math.sqrt( longitudinal_coeff / transverse_coeff ) self._print( f" |v_L/v_T| = " f"{ratio:.12f}" ) self._print( f" sqrt(3) = " f"{math.sqrt(3.0):.12f}" ) return eigs # ================================================================== # DISPERSION PATH # ================================================================== def compute_dispersion_path( self, n_points=None ): """ Compute the straight-line path Gamma -> 24A1 -> D24 -> Gamma. Note: This is a selected path through the specified q representatives. It is not, by itself, a construction of the full first Brillouin zone. """ if n_points is None: n_points = self.path_points n_points = int( n_points ) if n_points < 1: raise ValueError( "n_points must be >= 1." ) segments = [ ( "Gamma -> 24A1", self.k_gamma, self.k_24A1 ), ( "24A1 -> D24", self.k_24A1, self.k_D24 ), ( "D24 -> Gamma", self.k_D24, self.k_gamma ) ] all_eigenvalues = [] all_real_frequencies = [] all_imaginary_frequencies = [] positions = [] current_distance = 0.0 self._print( "\n[5] Complete dispersion path" ) self._print( f"\nComputing dispersion path " f"({3 * n_points} segment samples)..." ) for ( name, q0, q1 ) in segments: length = float( np.linalg.norm( q1 - q0 ) ) self._print( f" {name}: length = " f"{length:.8f}" ) for step in range( n_points ): t = ( step / float(n_points) ) q = ( (1.0 - t) * q0 + t * q1 ) D = ( self.compute_dynamical_matrix( q ) ) eigs = np.linalg.eigvalsh( D ) ( real_freqs, imag_freqs ) = self.classify_frequencies( eigs ) all_eigenvalues.append( eigs ) all_real_frequencies.append( real_freqs ) all_imaginary_frequencies.append( imag_freqs ) positions.append( current_distance + t * length ) current_distance += length # -------------------------------------------------------------- # Final Gamma endpoint # -------------------------------------------------------------- D = ( self.compute_dynamical_matrix( self.k_gamma ) ) eigs = np.linalg.eigvalsh( D ) ( real_freqs, imag_freqs ) = self.classify_frequencies( eigs ) all_eigenvalues.append( eigs ) all_real_frequencies.append( real_freqs ) all_imaginary_frequencies.append( imag_freqs ) positions.append( current_distance ) eigenvalues = np.asarray( all_eigenvalues, dtype=np.float64 ) real_frequencies = np.asarray( all_real_frequencies, dtype=np.float64 ) imaginary_frequencies = np.asarray( all_imaginary_frequencies, dtype=np.float64 ) positions = np.asarray( positions, dtype=np.float64 ) self._print( f" Number of path points: " f"{len(positions)}" ) return ( positions, eigenvalues, real_frequencies, imaginary_frequencies ) # ================================================================== # STABILITY DIAGNOSTIC # ================================================================== def stability_diagnostic( self, eigenvalues ): minimum = float( np.min(eigenvalues) ) self._print( "\n[6] Stability diagnostic" ) if minimum < -1e-10: self._print( " WARNING: the minimal-shell model " "has negative dynamical-matrix eigenvalues." ) self._print( " These are mechanical instabilities, " "not zero-frequency acoustic modes." ) self._print( " For lambda < 0:" ) self._print( " omega = i sqrt(-lambda)." ) elif minimum > 1e-10: self._print( " Dynamical matrix is positive definite " "on the sampled path." ) else: self._print( " Zero eigenvalues are present." ) # ================================================================== # FULL ANALYSIS # ================================================================== def run_analysis(self): self._print( "\n" + "=" * 80 ) self._print( "DETAILED PHONON ANALYSIS" ) self._print( "=" * 80 ) # -------------------------------------------------------------- # 1. Acoustic limit # -------------------------------------------------------------- self.acoustic_limit( dk=1e-5 ) # -------------------------------------------------------------- # 2. Deep-hole spectra # -------------------------------------------------------------- self._print( "\n[2] Deep-hole spectra" ) eigs_24A1 = ( self.report_deep_hole( "24A1", self.k_24A1 ) ) eigs_D24 = ( self.report_deep_hole( "D24", self.k_D24 ) ) # -------------------------------------------------------------- # 3. Exact 24A1 check # -------------------------------------------------------------- self.exact_24A1_check() # -------------------------------------------------------------- # 4. D24 symmetry # -------------------------------------------------------------- self.axial_symmetry_check() # -------------------------------------------------------------- # 5. Complete dispersion path # -------------------------------------------------------------- ( positions, eigenvalues, real_frequencies, imaginary_frequencies ) = self.compute_dispersion_path() # -------------------------------------------------------------- # Global statistics # -------------------------------------------------------------- global_min_eigenvalue = float( np.min(eigenvalues) ) global_max_eigenvalue = float( np.max(eigenvalues) ) maximum_real_frequency = float( np.max(real_frequencies) ) maximum_imaginary_frequency = float( np.max(imaginary_frequencies) ) self._print( f" Global minimum eigenvalue: " f"{global_min_eigenvalue:.12f}" ) self._print( f" Global maximum eigenvalue: " f"{global_max_eigenvalue:.12f}" ) self._print( f" Maximum real frequency: " f"{maximum_real_frequency:.12f}" ) self._print( f" Maximum imaginary frequency: " f"{maximum_imaginary_frequency:.12f}" ) # -------------------------------------------------------------- # Deep-hole imaginary frequencies # -------------------------------------------------------------- if np.min( eigs_24A1 ) < 0: imag_24A1 = math.sqrt( -np.min( eigs_24A1 ) ) self._print( f" 24A1 imaginary frequency magnitude: " f"{imag_24A1:.12f}" ) if np.min( eigs_D24 ) < 0: imag_D24 = math.sqrt( -np.min( eigs_D24 ) ) self._print( f" D24 imaginary frequency magnitude: " f"{imag_D24:.12f}" ) # -------------------------------------------------------------- # 6. Stability # -------------------------------------------------------------- self.stability_diagnostic( eigenvalues ) self._print( "\n" + "=" * 80 ) self._print( "ANALYSIS COMPLETE" ) self._print( "=" * 80 ) return { "positions": positions, "eigenvalues": eigenvalues, "real_frequencies": real_frequencies, "imaginary_frequencies": imaginary_frequencies, } # ====================================================================== # MAIN # ====================================================================== if __name__ == "__main__": engine = LeechMinimalShellDynamics( path_points=25, verbose=True ) results = engine.run_analysis()