Simulation quantum part. 1
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
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
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
Yes, we are here. But it wasn’t clear from your email what kind of response you wanted. all the best cjc On 14 Mar 2025, at 07:49, antony padilla <antonypamo@gmail.com> wrote: Hi, is anyone there? Antony Padilla Morales El jue, 13 mar 2025, 02:41, antony padilla <antonypamo@gmail.com<mailto: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 _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
antony padilla
-
Cotter, Colin J