Hi, is anyone there? Antony Padilla Morales El jue, 13 mar 2025, 02:41, antony padilla <antonypamo@gmail.com> escribió:
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