Antony Padilla Morales ---------- Forwarded message --------- De: antony padilla <antonypamo@gmail.com> Date: jue, 13 mar 2025, 02:41 Subject: Simulation quantum part. 1 To: <firedrake@imperial.ac.uk> import numpy as np import scipy.linalg import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D ############################################################################## 1. Parámetros y configuración ############################################################################## hbar = 1.0 t_hop = 1.0 G = 1.0 M = 1.0 lambda_g = 0.2 dt = 0.02 steps = 50 factor para la evolución cuántica: exp(-iH dt / hbar) factor_evol = -1j * dt / hbar ############################################################################## 2. Definición del icosaedro base (sin deformación) ############################################################################## phi = (1 + np.sqrt(5)) / 2.0 base_vertices = np.array([ [-1, phi, 0], [ 1, phi, 0], [-1, -phi, 0], [ 1, -phi, 0], [ 0, -1, phi], [ 0, 1, phi], [ 0, -1, -phi], [ 0, 1, -phi], [ phi, 0, -1], [ phi, 0, 1], [-phi, 0, -1], [-phi, 0, 1] ], dtype=float) Normalizamos para que un vértice tenga radio 1 r0 = np.linalg.norm(base_vertices[0]) base_vertices /= r0 Lista de aristas del icosaedro (12 vértices, 30 aristas) edges = [ (0,1), (0,5), (0,7), (0,10), (0,11), (1,5), (1,7), (1,8), (1,9), (2,3), (2,6), (2,7), (2,9), (2,11), (3,6), (3,7), (3,8), (3,9), (4,5), (4,6), (4,9), (4,10), (4,11), (5,9), (5,10), (6,7), (6,9), (6,10), (7,8), (8,9), (8,10), (10,11) ] def build_adjacency_matrix(): """ Construye la matriz de adyacencia (12x12) con -t_hop en las posiciones correspondientes a las aristas. """ A = np.zeros((12, 12), dtype=complex) for (i, j) in edges: A[i, j] = -t_hop A[j, i] = -t_hop return A ############################################################################## 3. Potencial logarítmico estático (sin deformar vértices) ############################################################################## def build_potential_log(vertices): """ Calcula el potencial gravitatorio discreto V_G(i) = - sum(G*M / dist(i,j)) y luego aplica la corrección logarítmica V_log(i) = lambda_g * log(1 + V_G(i)). Como la geometría no cambia, este potencial es constante en el tiempo. """ N = len(vertices) V_G = np.zeros(N, dtype=float) for i in range(N): for j in range(N): if i != j: dist_ij = np.linalg.norm(vertices[i] - vertices[j]) V_G[i] += - (G * M / dist_ij) # Evitamos valores no válidos en la log min_val = np.min(V_G) shift = 0.0 if min_val <= -1.0: shift = -(min_val + 0.9999) V_log = lambda_g * np.log(1.0 + V_G + shift) return V_log ############################################################################## 4. Construcción del Hamiltoniano estático ############################################################################## adjacency = build_adjacency_matrix() El icosaedro permanece siempre igual (sin "update_vertices") Por tanto, calculamos el potencial logarítmico solo una vez: V_log_static = build_potential_log(base_vertices) Hamiltoniano total fijo: H = adjacency + diag(V_log_static) H_static = adjacency + np.diag(V_log_static) Operador de evolución para un paso de tiempo dt U_static = scipy.linalg.expm(factor_evol * H_static) ############################################################################## 5. Estado cuántico inicial y bucle de evolución ############################################################################## psi = np.zeros(12, dtype=complex) psi[0] = 1.0 # partícula localizada en el vértice 0 al inicio psi_t_history = [psi.copy()] for _ in range(steps): psi = U_static @ psi # Renormalizamos para evitar errores numéricos norm_psi = np.linalg.norm(psi) if norm_psi > 1e-12: psi /= norm_psi psi_t_history.append(psi.copy()) ############################################################################## 6. Visualización de resultados ############################################################################## prob_t = np.array([np.abs(state)**2 for state in psi_t_history]) Evolución de la probabilidad plt.figure(figsize=(10,6)) for i in range(12): plt.plot(prob_t[:, i], label=f'Vértice {i}') plt.xlabel('Paso de tiempo') plt.ylabel('Probabilidad') plt.title('Evolución de la probabilidad en un icosaedro estático (potencial fijo)') plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left') plt.tight_layout() plt.show() Visualización 3D de la configuración final (los vértices no cambian) final_prob = prob_t[-1] fig = plt.figure(figsize=(7,7)) ax = fig.add_subplot(111, projection='3d') colors = plt.cm.viridis(final_prob / final_prob.max()) for i in range(12): x, y, z = base_vertices[i] ax.scatter(x, y, z, color=colors[i], s=50 + 400*final_prob[i]) ax.set_title('Configuración final en el icosaedro estático') plt.show() Antony Padilla Morales