diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..cc2489f --- /dev/null +++ b/.flake8 @@ -0,0 +1,2 @@ +[flake8] +ignore=E741,E743 diff --git a/README.md b/README.md index 2d5f7c1..438279e 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,62 @@ -# quantum-entanglement-simulation -Simulate an entangled state of a quantum computer +# Quantum mechanics spin simulation + +Welcome to the Quantum Mechanics Spin Simulation repository. This Python-based project is designed to simulate quantum spin dynamics, providing an interactive platform to explore the fundamental principles of quantum mechanics through computational methods. Whether you're a student, educator, or researcher, this repository offers valuable tools to simulate and analyze the behavior of quantum spins under various conditions. + +## Single spin simulation + +The single spin simulation script allows users to simulate the quantum mechanics of a single spin system. Users can specify the orientation of a measurement apparatus in three-dimensional space and observe the behavior of a quantum spin subjected to this configuration. The simulation outputs include the probabilities of the spin aligning in either the up or down state and can be visualized in real-time to aid in understanding quantum superposition and measurement collapse. + +### Features: + +- **Dynamic Apparatus Orientation:** Users can input coordinates to align the measurement apparatus in any arbitrary direction. +- **Randomized Spin Measurement:** Utilizes a random number generator to simulate the probabilistic nature of quantum measurements. +- **Measurement Repeatability:** Simulate consecutive measurements to observe the quantum mechanical property of state collapse. + +## Two spin simulation + +The two spin simulation script extends the capabilities of the single spin simulator by introducing interactions between two spins. This script includes various predefined states such as entangled states and product states, allowing users to simulate complex quantum phenomena like entanglement and spin correlation. + +### Simulation types: + +- **Product States:** Simulate two independent spins in various configurations. +- **Entangled States:** Explore the intriguing properties of entangled spins with several types of entangled states available, including the singlet and triplet states. +- **Measurement Options:** Perform measurements on both spins individually or simultaneously to observe correlations directly resulting from quantum entanglement. + +### Output: + +- **Statistical Analysis:** The script provides detailed statistics on spin states and their correlations after multiple measurements, essential for understanding entanglement. + +## Getting Started + +To get started with these simulations: +1. Clone the repository: + ``` + git clone https://github.com/azimonti/quantum-entanglement-simulation.git + ``` +2. Navigate to the repository directory: + ``` + cd quantum-entanglement-simulation + ``` +3. Install required dependencies: + ``` + pip install -r requirements.txt + ``` +4. Run the simulation scripts: + ``` + python single_spin_sim.py + python two_spin_sim.py + ``` + +## Contributing + +Contributions to the Quantum Mechanics Spin Simulation project are welcome. Whether it's through submitting bug reports, proposing new features, or contributing to the code, your help is appreciated. For major changes, please open an issue first to discuss what you would like to change. + +## License + +This project is licensed under the MIT License - see the [LICENSE](LICENSE.md) file for details. + +## Contact + +If you have any questions or want to get in touch regarding the project, please open an issue or contact the repository maintainers directly through GitHub. + +Thank you for exploring the quantum mechanics of spins with us! diff --git a/mod_spin_operators.py b/mod_spin_operators.py new file mode 100644 index 0000000..c99e203 --- /dev/null +++ b/mod_spin_operators.py @@ -0,0 +1,215 @@ +#!/usr/bin/env python3 +''' +/************************/ +/* mod_spin_operators */ +/* Version 1.0 */ +/* 2024/05/11 */ +/************************/ +''' +import cmath +import math +import numpy as np +import sys + + +class SingleSpin: + + def __init__(self, basis: str = 'ud'): + self.__basis = basis + self.__state = None + if (basis == 'ud'): + self.__u = np.array([[1], [0]], dtype=complex) + self.__d = np.array([[0], [1]], dtype=complex) + self.__r = np.array( + [[1 / np.sqrt(2)], [1 / np.sqrt(2)]], dtype=complex) + self.__l = np.array( + [[1 / np.sqrt(2)], [-1 / np.sqrt(2)]], dtype=complex) + self.__i = np.array( + [[1 / np.sqrt(2)], [1j / np.sqrt(2)]], dtype=complex) + self.__o = np.array( + [[1 / np.sqrt(2)], [-1j / np.sqrt(2)]], dtype=complex) + self.__sz = np.array([[1, 0], [0, -1]], dtype=complex) + self.__sx = np.array([[0, 1], [1, 0]], dtype=complex) + self.__sy = np.array([[0, -1j], [1j, 0]], dtype=complex) + elif (basis == 'rl'): + self.__u = np.array( + [[1 / np.sqrt(2)], [1 / np.sqrt(2)]], dtype=complex) + self.__d = np.array( + [[1 / np.sqrt(2)], [-1 / np.sqrt(2)]], dtype=complex) + self.__r = np.array([[1], [0]], dtype=complex) + self.__l = np.array([[0], [1]], dtype=complex) + self.__i = np.array( + [[(1 + 1j) / 2], [(1 - 1j) / 2]], dtype=complex) + self.__o = np.array( + [[(1 - 1j) / 2], [(1 + 1j) / 2]], dtype=complex) + elif (basis == 'io'): + self.__u = np.array( + [[1 / np.sqrt(2)], [1 / np.sqrt(2)]], dtype=complex) + self.__d = np.array( + [[1j / np.sqrt(2)], [-1j / np.sqrt(2)]], dtype=complex) + self.__r = np.array( + [[(1 + 1j) / 2], [(1 - 1j) / 2]], dtype=complex) + self.__l = np.array( + [[(1 - 1j) / 2], [(1 + 1j) / 2]], dtype=complex) + self.__o = np.array([[1], [0]], dtype=complex) + self.__i = np.array([[0], [1]], dtype=complex) + else: + raise NotImplementedError( + "Basis " + basis + "not Implemented") + + @property + def u(self): + return self.__u + + @property + def d(self): + return self.__d + + @property + def r(self): + return self.__r + + @property + def l(self): + return self.__l + + @property + def i(self): + return self.__i + + @property + def o(self): + return self.__o + + @property + def s_z(self): + if (self.__basis != 'ud'): + raise NotImplementedError( + "S_z for basis " + self.__basis + "not Implemented") + return self.__sz + + @property + def s_x(self): + if (self.__basis != 'ud'): + raise NotImplementedError( + "S_x for basis " + self.__basis + "not Implemented") + return self.__sx + + @property + def s_y(self): + if (self.__basis != 'ud'): + raise NotImplementedError( + "S_y for basis " + self.__basis + "not Implemented") + return self.__sy + + @property + def psi(self): + return self.__state + + @psi.setter + def psi(self, value): + # check that length is unitary + assert math.isclose(np.linalg.norm(value), 1) + self.__state = value + + def theta(self): + assert self.__state + return 2 * np.arccos(np.linalg.norm(self.__state[0][0])) + + def phi(self): + assert self.__state + if self.__state[1][0] != 0: + return cmath.phase(self.__state[1][0]) + else: + return 0 + + def angles(self): + return [self.theta(), self.phi()] + + +class TwolSpins: + + def __init__(self, basis: str = 'ud'): + self.__basis = basis + self.__state = None + if (basis == 'ud'): + u = np.array([[1], [0]], dtype=complex) + d = np.array([[0], [1]], dtype=complex) + self.__b = [ + np.kron(u, u), np.kron(u, d), + np.kron(d, u), np.kron(d, d) + ] + self.__bmap = {'uu': 0, 'ud': 1, 'du': 2, 'dd': 3} + self.__s = [ + np.array([[1, 0], [0, -1]], dtype=complex), + np.array([[0, 1], [1, 0]], dtype=complex), + np.array([[0, -1j], [1j, 0]], dtype=complex), + np.array([[1, 0], [0, 1]], dtype=complex) + ] + self.__smap = {'z': 0, 'x': 1, 'y': 2, 'I': 3} + else: + raise NotImplementedError( + "Basis " + basis + "not Implemented") + + @property + def psi(self): + return self.__state + + @ psi.setter + def psi(self, value): + # check that length is unitary + assert math.isclose(np.linalg.norm(value), 1) + self.__state = value + + def BasisVector(self, s): + return self.__b[self.__bmap[s]] + + def Sigma(self, sA: str, sB: str): + return np.kron( + self.__s[self.__smap[sA]], self.__s[self.__smap[sB]]) + + def Sigma_A(self, s: str): + return np.kron(self.__s[self.__smap[s]], self.__s[3]) + + def Sigma_B(self, s: str): + return np.kron(self.__s[3], self.__s[self.__smap[s]]) + + def Expectation(self, sA: str, sB: str): + if (sA == 'I'): + # expectation of System B + return np.linalg.multi_dot([ + self.__state.conj().T, self.Sigma_B(sB), self.__state])[0, 0] + elif (sB == 'I'): + # expectation of System A + return np.linalg.multi_dot([ + self.__state.conj().T, self.Sigma_A(sA), self.__state])[0, 0] + else: + # expectation of the composite system + return np.linalg.multi_dot([ + self.__state.conj().T, self.Sigma(sA, sB), self.__state])[0, 0] + + def ProductState(self, A: np.array, B: np.array): + # check that length is unitary + assert math.isclose(np.linalg.norm(A), 1) + assert math.isclose(np.linalg.norm(B), 1) + self.psi = np.kron(A, B) + + def Singlet(self): + self.psi = 1 / np.sqrt(2) * (self.__b[1] - self.__b[2]) + + def Triplet(self, i: int): + match i: + case 1: + self.psi = 1 / np.sqrt(2) * (self.__b[1] + self.__b[2]) + case 2: + self.psi = 1 / np.sqrt(2) * (self.__b[0] + self.__b[3]) + case 3: + self.psi = 1 / np.sqrt(2) * (self.__b[0] - self.__b[3]) + case _: + raise ValueError("Incorrect index " + str(i)) + + +if __name__ == '__main__': + if sys.version_info[0] < 3: + raise 'Must be using Python 3' + pass diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..0f9fddf --- /dev/null +++ b/requirements.txt @@ -0,0 +1,2 @@ +pyqt6==6.5.2 +numpy==1.25.2 diff --git a/single_spin_sim.py b/single_spin_sim.py new file mode 100644 index 0000000..39a7ad1 --- /dev/null +++ b/single_spin_sim.py @@ -0,0 +1,379 @@ +#!/usr/bin/env python3 +''' +/************************/ +/* single_spin_sim.py */ +/* Version 1.0 */ +/* 2024/05/11 */ +/************************/ +''' +import argparse +import cmath +import math +import numpy as np +import random +import sys +import time +from PyQt6 import QtWidgets +from PyQt6.QtWidgets import QPushButton, QSlider, QLabel +from PyQt6.QtWidgets import QVBoxLayout, QGridLayout +from PyQt6.QtWidgets import QWidget +from PyQt6.QtOpenGLWidgets import QOpenGLWidget +from PyQt6.QtGui import QPainter, QFont, QColor +from PyQt6.QtCore import Qt, QThread, pyqtSignal, QRect +from OpenGL.GL import ( + glClear, glClearColor, glEnable, glPushMatrix, glPopMatrix, glRotatef, + glTranslatef, glBegin, glEnd, glVertex3f, glViewport, glMatrixMode, + glLoadIdentity, glColor3f) +from OpenGL.GLU import gluPerspective, gluLookAt +from OpenGL.GL import ( + GL_DEPTH_TEST, GL_COLOR_BUFFER_BIT, GL_DEPTH_BUFFER_BIT, + GL_QUADS, GL_LINES, glFlush, GL_PROJECTION, GL_MODELVIEW) +from mod_spin_operators import SingleSpin + +description = ''' +This script simulate a single spin following quantum mechanics principles. + +Simulation types available (-t SIMUL_TYPE, --simul_type SIMUL_TYPE): + +NO Time evolution +1 - The spin is prepared to be always in the up direction (reset at +each time step) [DEFAULT] +2 - The spin is prepared to be always in the left direction (reset at +each time step) +2 - The spin is prepared to be always in the inner direction (reset at +each time step) +4 - The spin is collapsed in the direction of the measurement +''' + + +class SimulationThread(QThread): + result = pyqtSignal(np.ndarray) + + def __init__(self, simul_type: int): + super().__init__() + self.simul_type = simul_type + self.current_state = None + + def run(self): + simul_spin = SingleSpin() + + # initial condition for the case needed + match self.simul_type: + case 1: + self.current_state = simul_spin.u + case 2: + self.current_state = simul_spin.l + case 3: + self.current_state = simul_spin.i + case 4: + self.current_state = simul_spin.u + case _: + raise ValueError( + f"Incorrect simulation type {self.simul_type}") + while True: + self.result.emit(self.current_state) + time.sleep(1) # slows down the loop for demonstration purposes + + def collapse_wave_function(self, state: np.ndarray): + self.current_state = state + + +class OpenGLWidget(QOpenGLWidget): + + def __init__(self, parent, simul_type: int): + super(OpenGLWidget, self).__init__(parent) + self.a_theta = 0 + self.a_phi = 0 + self.measurement = None + self.count_p1 = 0 + self.count_m1 = 0 + self.current_state = None + self.num_measurements = 0 + self.spin = SingleSpin() + self.simul_type = simul_type + + @property + def apparatus_direction(self): + return np.array([ + [np.cos(self.a_theta / 2)], + [np.exp(1j * self.a_phi) * np.sin(self.a_theta / 2)]]) + + @property + def apparatus_opposite_direction(self): + return np.array([ + [np.cos(np.pi / 2 + self.a_theta / 2)], + [np.exp(1j * self.a_phi) * np.sin( + np.pi / 2 + self.a_theta / 2)]]) + + def initializeGL(self): + glClearColor(0.0, 0.0, 0.0, 1.0) + glEnable(GL_DEPTH_TEST) + + def paintGL(self): + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + glLoadIdentity() + glColor3f(1.0, 1.0, 1.0) + # Adjust the camera view + gluLookAt(0.0, -5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0) + glPushMatrix() + glMatrixMode(GL_MODELVIEW) + # Rotate the rectangle along Y-axis + a_theta_deg = math.degrees(self.a_theta) + a_phi_deg = math.degrees(self.a_phi) + # print(f"apparatus theta = {a_theta_deg}\n phi = {a_phi_deg}") + + # Apply the rotation + glRotatef(-a_phi_deg, 0.0, 0.0, 1.0) + glRotatef(a_theta_deg, 0.0, 1.0, 0.0) + # Draw a rectangle on ZX plane + glBegin(GL_QUADS) + glVertex3f(-0.5, 0, -0.5) + glVertex3f(0.5, 0, -0.5) + glVertex3f(0.5, 0, 0.5) + glVertex3f(-0.5, 0, 0.5) + glEnd() + + # Draw an arrow + glBegin(GL_LINES) + glVertex3f(0., 0, -1) + glVertex3f(0., 0, 1) + glVertex3f(0., 0, 1) + glVertex3f(-.2, 0, 0.75) + glVertex3f(0., 0, 1) + glVertex3f(.2, 0, 0.75) + glVertex3f(-1, 0, 0) + glVertex3f(1, 0, 0) + glEnd() + glPopMatrix() + + if self.measurement is not None: + painter = QPainter(self) + painter.setFont(QFont('Arial', 14)) + painter.setPen(QColor(255, 255, 255)) + y = int(0.25 * self.height() + 20) + rect = QRect(0, y, self.width(), self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Measurement: {self.measurement}") + y = int(0.25 * self.height() + 55) + rect = QRect(0, y, self.width(), self.height() - y) + prob_p1 = self.count_p1 / self.num_measurements + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Total Measurements: {self.num_measurements}") + y = int(0.25 * self.height() + 90) + rect = QRect(0, y, self.width(), self.height() - y) + prob_p1 = self.count_p1 / self.num_measurements + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< +1 > = {prob_p1:.2f}") + y = int(0.25 * self.height() + 125) + rect = QRect(0, y, self.width(), self.height() - y) + prob_m1 = self.count_m1 / self.num_measurements + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< -1 > = {prob_m1:.2f}") + + y = int(0.75 * self.height() - 40) + x = int(0.25 * self.width()) + rect = QRect(x, 0, self.width() - x, y) + painter.setPen(QColor(255, 255, 0)) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + "current state") + # display the current spin + glPushMatrix() + glMatrixMode(GL_MODELVIEW) + glTranslatef(1.0, 0., 0.7) + glColor3f(1.0, 1.0, 0.0) + + # compute the spin angle + s_theta = 2 * np.arctan2( + abs(self.current_state[1][0]), abs(self.current_state[0][0])) + if self.current_state[1][0] != 0: + s_phi = cmath.phase(self.current_state[1][0]) + else: + s_phi = 0 + # If the real part is negative, set it > pi + if self.current_state[0][0].real < 0: + s_theta = 2 * np.pi - s_theta + + if s_phi < 0: + s_phi += 2 * np.pi + # Convert theta and phi to degrees + s_theta_deg = math.degrees(s_theta) + s_phi_deg = math.degrees(s_phi) + # print(f"spin theta = {s_theta_deg}\nspin phi = {s_phi_deg}") + + y = int(0.25 * self.height() + 160) + rect = QRect(0, y, self.width(), self.height() - y) + prob_m1 = self.count_m1 / self.num_measurements + + # calculate the dot product + dot_product = np.sin(self.a_theta) * np.sin(s_theta) * \ + np.cos(self.a_phi - s_phi) + \ + np.cos(self.a_theta) * np.cos(s_theta) + + # calculate the angle difference + angle_difference = np.arccos(dot_product) + # convert angle difference from radians to degrees + angle_difference_degrees = np.degrees(angle_difference) + cos_half_alpha_2 = np.cos( + np.deg2rad(angle_difference_degrees / 2))**2 + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"cos(θ_mn/2)^2: {cos_half_alpha_2:.2f}") + painter.end() + + # Apply the rotation - Order is important + glRotatef(s_phi_deg, 0., 0., 1.) + glRotatef(s_theta_deg, 0., 1., 0.) + + # Draw an arrow + glBegin(GL_LINES) + glVertex3f(0., 0, -0.1) + glVertex3f(0., 0, 0.1) + glVertex3f(0., 0, 0.1) + glVertex3f(-0.1, 0, 0.) + glVertex3f(0., 0, 0.1) + glVertex3f(0.1, 0, 0.) + glEnd() + glPopMatrix() + + glFlush() + + def resizeGL(self, w: int, h: int): + glViewport(0, 0, w, h) + glMatrixMode(GL_PROJECTION) + glLoadIdentity() + gluPerspective(60.0, w / h, 0.1, 100.0) + glMatrixMode(GL_MODELVIEW) + + def update_rotation_theta(self, value: float): + self.a_theta = np.deg2rad(value) + self.update() + + def update_rotation_phi(self, value: float): + self.a_phi = np.deg2rad(value) + self.update() + + def measure(self, current_state: np.ndarray): + self.current_state = current_state + # get the measurement direction + direction = np.array([ + np.cos(self.a_theta / 2), + np.exp(1j * self.a_phi) * np.sin(self.a_theta / 2)]) + prob_p1 = np.abs(np.vdot(direction, self.current_state)) ** 2 + # Generate a random number between 0 and 1 + random_number = random.uniform(0, 1) + + # Perform the measurement in apparatus direction + self.num_measurements += 1 + if random_number < prob_p1: + self.count_p1 += 1 + self.measurement = 1 + else: + self.count_m1 += 1 + self.measurement = -1 + self.update() + + +class MainWindow(QWidget): + + def __init__(self, simul_type: int): + super(MainWindow, self).__init__() + self.simul_type = simul_type + + self.initUI() + + self.current_state = None + self.simulation_thread = SimulationThread(simul_type) + self.simulation_thread.result.connect(self.store_simul_spin) + self.simulation_thread.start() + + def initUI(self): + self.setGeometry(300, 300, 800, 600) + match self.simul_type: + case 1: + desc = 'spin always |up>' + case 2: + desc = 'spin always |left>' + case 3: + desc = 'spin always |inner>' + case 4: + desc = 'spin takes measurement direction' + self.setWindowTitle(f"Single quantum spin simulation: {desc}") + + self.opengl_widget = OpenGLWidget(self, self.simul_type) + + self.slider1 = QSlider(Qt.Orientation.Horizontal, self) + self.slider1.setRange(0, 360) + self.slider1.valueChanged.connect(self.update_rotation_theta) + self.label1 = QLabel("Apparatus Rotation θ: 0", self) + + self.slider2 = QSlider(Qt.Orientation.Horizontal, self) + self.slider2.setRange(0, 360) + self.slider2.valueChanged.connect(self.update_rotation_phi) + self.label2 = QLabel("Apparatus Rotation 𝜙: 0", self) + + self.button = QPushButton('Make measurement', self) + self.button.clicked.connect(self.on_button_clicked) + + self.layout = QVBoxLayout(self) + self.layout.addWidget(self.opengl_widget) + self.gridlayout = QGridLayout() + self.gridlayout.addWidget(self.slider1, 0, 0) + self.gridlayout.addWidget(self.label1, 0, 1) + self.gridlayout.addWidget(self.slider2, 1, 0) + self.gridlayout.addWidget(self.label2, 1, 1) + self.gridlayout.addWidget(self.button, 2, 0, 1, 2) + self.layout.addLayout(self.gridlayout) + + def store_simul_spin(self, current_simul_state: np.ndarray): + self.current_state = current_simul_state + + def update_rotation_theta(self, value: float): + self.opengl_widget.update_rotation_theta(value) + self.label1.setText(f"Apparatus Rotation θ: {value}") + + def update_rotation_phi(self, value: float): + self.opengl_widget.update_rotation_phi(value) + self.label2.setText(f"Apparatus Rotation 𝜙: {value}") + + def on_button_clicked(self): + self.opengl_widget.measure(self.current_state) + match self.simul_type: + case 1 | 2 | 3: + pass + case 4: + if self.opengl_widget.measurement == 1: + state = self.opengl_widget.apparatus_direction + else: + state = self.opengl_widget.apparatus_opposite_direction + self.simulation_thread.collapse_wave_function(state) + + +class CustomHelpFormatter(argparse.HelpFormatter): + def _fill_text(self, text, width, indent): + # Preserve line breaks by not wrapping text + return "\n".join([indent + line for line in text.splitlines()]) + + +def main(): + # Set a fixed seed value + seed_value = 5692 + random.seed(seed_value) + parser = argparse.ArgumentParser(description=description, + formatter_class=CustomHelpFormatter) + parser.add_argument('-t', '--simul_type', help='simulation type', + required=False) + + args = parser.parse_args() + if (args.simul_type): + simul_type = int(args.simul_type) + else: + simul_type = 1 + app = QtWidgets.QApplication(sys.argv) + window = MainWindow(simul_type) + window.show() + sys.exit(app.exec()) + + +if __name__ == '__main__': + if sys.version_info[0] < 3: + raise 'Must be using Python 3' + main() diff --git a/two_spin_sim.py b/two_spin_sim.py new file mode 100644 index 0000000..fe409b8 --- /dev/null +++ b/two_spin_sim.py @@ -0,0 +1,631 @@ +#!/usr/bin/env python3 +''' +/************************/ +/* two_spin_sim.py */ +/* Version 1.0 */ +/* 2024/05/11 */ +/************************/ +''' +import argparse +import math +import numpy as np +import random +import sys +import time +from PyQt6 import QtWidgets +from PyQt6.QtWidgets import QPushButton, QSlider, QLabel +from PyQt6.QtWidgets import QVBoxLayout, QGridLayout +from PyQt6.QtWidgets import QWidget +from PyQt6.QtOpenGLWidgets import QOpenGLWidget +from PyQt6.QtGui import QPainter, QFont, QColor +from PyQt6.QtCore import Qt, QThread, pyqtSignal, QRect +from OpenGL.GL import ( + glClear, glClearColor, glEnable, glPushMatrix, glPopMatrix, glRotatef, + glTranslatef, glBegin, glEnd, glVertex3f, glViewport, glMatrixMode, + glLoadIdentity, glColor3f) +from OpenGL.GLU import gluPerspective, gluLookAt +from OpenGL.GL import ( + GL_DEPTH_TEST, GL_COLOR_BUFFER_BIT, GL_DEPTH_BUFFER_BIT, + GL_QUADS, GL_LINES, glFlush, GL_PROJECTION, GL_MODELVIEW) +from mod_spin_operators import SingleSpin, TwolSpins + +description = ''' +This script simulate a two spins following quantum mechanics principles. + +Simulation types available (-t SIMUL_TYPE, --simul_type SIMUL_TYPE): + +NO Time evolution +1 - Product state + A = | u > + B = | d > +2 - Product state + A = 1 / sqrt(2) * (| u > + | d >) + B = | u > / 2 + sqrt(3) / 2 * | d > +3 - Singlet state + | Psi > = 1 / sqrt(2) ( | ud > - | du > ) [DEFAULT] +4 - Triplet state I + | Psi > = 1 / sqrt(2) ( | ud > + | du > ) +5 - Triplet state II + | Psi > = 1 / sqrt(2) ( | uu > + | dd > ) +6 - Triplet state III + | Psi > = 1 / sqrt(2) ( | uu > - | dd > ) +7 - Partially entangled state + | Psi > = sqrt(0.6) | ud > - sqrt(0.4) | du > ) + It requires a few hundrends measurements to show the correct + correlation σ^z = 1.00, σ^x = σ^y = 0.98 + + +The program is measuring each component separately when any of the two +measure button is pressed and resetting the status. it is possible to +measure both system as the same time when one of the measure button is +pressed and the correlation between these measurement is shown setting +the command line option "-m, --measure_both" (default = False) + +''' + + +class SimulationThread(QThread): + result = pyqtSignal(np.ndarray) + + def __init__(self, simul_type: int): + super().__init__() + self.simul_type = simul_type + self.current_state = None + + def run(self): + spin = TwolSpins() + s = SingleSpin() + + # initial condition for the case needed + match self.simul_type: + case 1: + A = s.u + B = s.d + spin.ProductState(A, B) + case 2: + A = 1 / np.sqrt(2) * (s.u + s.d) + B = s.u / 2 + np.sqrt(3) / 2 * s.d + spin.ProductState(A, B) + case 3: + spin.Singlet() + case 4: + spin.Triplet(1) + case 5: + spin.Triplet(2) + case 6: + spin.Triplet(3) + case 7: + spin.psi = np.sqrt(0.6) * spin.BasisVector('ud') - \ + np.sqrt(0.4) * spin.BasisVector('du') + case _: + raise ValueError( + f"Incorrect simulation type {self.simul_type}") + while True: + self.current_state = spin.psi + self.result.emit(self.current_state) + time.sleep(1) # slows down the loop for demonstration purposes + + +class OpenGLWidget(QOpenGLWidget): + + def __init__(self, parent, simul_type: int, measure_both: bool): + super(OpenGLWidget, self).__init__(parent) + self.simul_type = simul_type + self.measure_both = measure_both + self.a_thetaA = 0 + self.a_phiA = 0 + self.measurementA = None + self.count_p1A = 0 + self.count_m1A = 0 + self.a_thetaB = 0 + self.a_phiB = 0 + self.measurementB = None + self.count_p1B = 0 + self.count_m1B = 0 + self.current_state = None + self.sigma = {'A': {'x': [], 'y': [], 'z': [], 'th_ph': []}, + 'B': {'x': [], 'y': [], 'z': [], 'th_ph': []} + } + self.directions = { + 'z': np.array([[1, 0], + [0, 1]]), + 'x': np.array([[1 / np.sqrt(2), 1 / np.sqrt(2)], + [1 / np.sqrt(2), -1 / np.sqrt(2)]]), + 'y': np.array([[1 / np.sqrt(2), 1j / np.sqrt(2)], + [1 / np.sqrt(2), -1j / np.sqrt(2)]]) + } + + @property + def apparatusA_direction(self): + return np.array([ + [np.cos(self.a_thetaA / 2)], + [np.exp(1j * self.a_phiA) * np.sin(self.a_thetaA / 2)]]) + + def initializeGL(self): + glClearColor(0.0, 0.0, 0.0, 1.0) + glEnable(GL_DEPTH_TEST) + + def paintGL(self): + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) + glLoadIdentity() + glColor3f(0.0, 1.0, 0.0) # Set color to green + # Adjust the camera view + gluLookAt(0.0, -5.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0) + glTranslatef(-3.0, 0.0, 0.0) # Move to the left + glPushMatrix() + glMatrixMode(GL_MODELVIEW) + # Rotate the rectangle along Y-axis + a_theta_deg = math.degrees(self.a_thetaA) + a_phi_deg = math.degrees(self.a_phiA) + # Apply the rotation + glRotatef(-a_phi_deg, 0.0, 0.0, 1.0) + glRotatef(a_theta_deg, 0.0, 1.0, 0.0) + self.drawRectangleAndArrow() + glPopMatrix() + glTranslatef(6.0, 0.0, 0.0) # Move to the right + glColor3f(1.0, 0.0, 0.0) # Set color to red + glPushMatrix() + glMatrixMode(GL_MODELVIEW) + # Rotate the rectangle along Y-axis + a_theta_deg = math.degrees(self.a_thetaB) + a_phi_deg = math.degrees(self.a_phiB) + # Apply the rotation + glRotatef(-a_phi_deg, 0.0, 0.0, 1.0) + glRotatef(a_theta_deg, 0.0, 1.0, 0.0) + self.drawRectangleAndArrow() + glPopMatrix() + + if len(self.sigma['A']['th_ph']) > 0: + saz = np.mean(self.sigma['A']['z']) + sax = np.mean(self.sigma['A']['x']) + say = np.mean(self.sigma['A']['y']) + sai = saz**2 + sax**2 + say**2 + painter = QPainter(self) + painter.setFont(QFont('Arial', 14)) + painter.setPen(QColor(255, 255, 255)) + half_width = int(self.width() / 2) + y = int(0.25 * self.height() + 55) + rect = QRect(0, 0, half_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Σ< σ_Ai > = {sai:.2f}") + y = int(0.25 * self.height() + 90) + rect = QRect(0, 0, half_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< σ^Ay > = {say:.2f}") + y = int(0.25 * self.height() + 125) + rect = QRect(0, 0, half_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< σ^Ax > = {sax:.2f}") + y = int(0.25 * self.height() + 160) + rect = QRect(0, 0, half_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< σ^Az > = {saz:.2f}") + y = int(0.25 * self.height() + 55) + rect = QRect(0, y, half_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Measurement: {self.measurementA}") + y = int(0.25 * self.height() + 90) + rect = QRect(0, y, half_width, self.height() - y) + painter.drawText( + rect, Qt.AlignmentFlag.AlignCenter, + f"Total Measurements: {len(self.sigma['A']['th_ph'])}") + y = int(0.25 * self.height() + 125) + rect = QRect(0, y, half_width, self.height() - y) + prob_p1 = self.count_p1A / len(self.sigma['A']['th_ph']) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f" < +1 > = {prob_p1:.2f}") + y = int(0.25 * self.height() + 160) + rect = QRect(0, y, half_width, self.height() - y) + prob_m1 = self.count_m1A / len(self.sigma['A']['th_ph']) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< -1 > = {prob_m1:.2f}") + painter.end() + + if len(self.sigma['B']['th_ph']) > 0: + sbz = np.mean(self.sigma['B']['z']) + sbx = np.mean(self.sigma['B']['x']) + sby = np.mean(self.sigma['B']['y']) + sbi = sbz**2 + sbx**2 + sby**2 + painter = QPainter(self) + painter.setFont(QFont('Arial', 14)) + painter.setPen(QColor(255, 255, 255)) + tq_width = int(self.width() * 3 / 2) + y = int(0.25 * self.height() + 55) + rect = QRect(0, 0, tq_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Σ< σ_Bi > = {sbi:.2f}") + y = int(0.25 * self.height() + 90) + rect = QRect(0, 0, tq_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< σ^By > = {sby:.2f}") + y = int(0.25 * self.height() + 125) + rect = QRect(0, 0, tq_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< σ^Bx > = {sbx:.2f}") + y = int(0.25 * self.height() + 160) + rect = QRect(0, 0, tq_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< σ^Bz > = {sbz:.2f}") + + y = int(0.25 * self.height() + 55) + rect = QRect(0, y, tq_width, self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Measurement: {self.measurementB}") + y = int(0.25 * self.height() + 90) + rect = QRect(0, y, tq_width, self.height() - y) + painter.drawText( + rect, Qt.AlignmentFlag.AlignCenter, + f"Total Measurements: {len(self.sigma['B']['th_ph'])}") + y = int(0.25 * self.height() + 125) + rect = QRect(0, y, tq_width, self.height() - y) + prob_p1 = self.count_p1B / len(self.sigma['B']['th_ph']) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< +1 > = {prob_p1:.2f}") + y = int(0.25 * self.height() + 160) + rect = QRect(0, y, tq_width, self.height() - y) + prob_m1 = self.count_m1B / len(self.sigma['B']['th_ph']) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"< -1 > = {prob_m1:.2f}") + painter.end() + + if (len(self.sigma['A']['th_ph']) > 0) and \ + (len(self.sigma['B']['th_ph']) > 0) and \ + (len(self.sigma['A']['z']) > 2): + corrz = np.corrcoef(self.sigma['A']['z'], + self.sigma['B']['z'])[0, 1] + corrx = np.corrcoef(self.sigma['A']['x'], + self.sigma['B']['x'])[0, 1] + corry = np.corrcoef(self.sigma['A']['y'], + self.sigma['B']['y'])[0, 1] + painter = QPainter(self) + painter.setFont(QFont('Arial', 14)) + painter.setPen(QColor(255, 255, 255)) + y = int(0.25 * self.height() + 90) + rect = QRect(0, 0, self.width(), self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Correlation <σ^Ay> <σ^By> = {corry:.2f}") + y = int(0.25 * self.height() + 125) + rect = QRect(0, 0, self.width(), self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Correlation <σ^Ax> <σ^Bx> = {corrx:.2f}") + y = int(0.25 * self.height() + 160) + rect = QRect(0, 0, self.width(), self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Correlation <σ^Az> <σ^Bz> = {corrz:.2f}") + if (self.measure_both): + y = int(0.25 * self.height() + 55) + corrm = np.corrcoef(self.sigma['A']['th_ph'], + self.sigma['B']['th_ph'])[0, 1] + rect = QRect(0, y, self.width(), self.height() - y) + painter.drawText(rect, Qt.AlignmentFlag.AlignCenter, + f"Correlation <σ^Am> <σ^Bm> = {corrm:.2f}") + painter.end() + + glFlush() + + def drawRectangleAndArrow(self): + glBegin(GL_QUADS) + glVertex3f(-0.5, 0, -0.5) + glVertex3f(0.5, 0, -0.5) + glVertex3f(0.5, 0, 0.5) + glVertex3f(-0.5, 0, 0.5) + glEnd() + + # Draw an arrow + glBegin(GL_LINES) + glVertex3f(0., 0, -1) + glVertex3f(0., 0, 1) + glVertex3f(0., 0, 1) + glVertex3f(-.2, 0, 0.75) + glVertex3f(0., 0, 1) + glVertex3f(.2, 0, 0.75) + glVertex3f(-1, 0, 0) + glVertex3f(1, 0, 0) + glEnd() + + def resizeGL(self, w: int, h: int): + glViewport(0, 0, w, h) + glMatrixMode(GL_PROJECTION) + glLoadIdentity() + gluPerspective(60.0, w / h, 0.1, 100.0) + glMatrixMode(GL_MODELVIEW) + + def update_rotation_thetaA(self, value: float): + self.a_thetaA = np.deg2rad(value) + self.update() + + def update_rotation_phiA(self, value: float): + self.a_phiA = np.deg2rad(value) + self.update() + + def measureA(self, current_state: np.ndarray): + self.measureAB(current_state, True) + + def updateCountA(self, value: int): + if value == 1: + self.count_p1A += 1 + self.measurementA = 1 + else: + self.count_m1A += 1 + self.measurementA = -1 + + def measureAB(self, current_state: np.ndarray, measureA: bool): + self.current_state = current_state + for axis in ['z', 'x', 'y']: + sp = self.measureSpin(self.directions[axis], + self.directions[axis], True) + self.sigma['A'][axis].append(sp[0]) + self.sigma['B'][axis].append(sp[1]) + + # get the measurement direction for A + directionAp = np.array([ + np.cos(self.a_thetaA / 2), + np.exp(1j * self.a_phiA) * np.sin(self.a_thetaA / 2)]) + # get the opposite measurement direction + directionAm = np.array([ + np.cos(np.pi / 2 + self.a_thetaA / 2), + np.exp(1j * self.a_phiA) * np.sin(np.pi / 2 + self.a_thetaA / 2)]) + + # get the measurement direction for B + directionBp = np.array([ + np.cos(self.a_thetaB / 2), + np.exp(1j * self.a_phiB) * np.sin(self.a_thetaB / 2)]) + # get the opposite measurement direction + directionBm = np.array([ + np.cos(np.pi / 2 + self.a_thetaB / 2), + np.exp(1j * self.a_phiB) * np.sin(np.pi / 2 + self.a_thetaB / 2)]) + + if measureA: + sp = self.measureSpin(np.array([directionAp, directionAm]), + np.array([directionBp, directionBm]), True) + self.sigma['A']['th_ph'].append(sp[0]) + self.updateCountA(sp[0]) + if self.measure_both: + self.sigma['B']['th_ph'].append(sp[1]) + self.updateCountB(sp[1]) + else: + sp = self.measureSpin(np.array([directionAp, directionAm]), + np.array([directionBp, directionBm]), False) + self.sigma['B']['th_ph'].append(sp[1]) + self.updateCountB(sp[1]) + if self.measure_both: + self.sigma['A']['th_ph'].append(sp[0]) + self.updateCountA(sp[0]) + + self.update() + + def update_rotation_thetaB(self, value: float): + self.a_thetaB = np.deg2rad(value) + self.update() + + def update_rotation_phiB(self, value: float): + self.a_phiB = np.deg2rad(value) + self.update() + + def measureB(self, current_state: np.ndarray): + self.measureAB(current_state, False) + + def updateCountB(self, value: int): + if value == 1: + self.count_p1B += 1 + self.measurementB = 1 + else: + self.count_m1B += 1 + self.measurementB = -1 + + def measureSpin(self, directionsA: np.ndarray, + directionsB: np.ndarray, simulate_A: bool): + ''' + Perform the measurement of the spin of two systems A and B, + which one to simulate is decided by the caller. + ''' + def RhoA(psi: np.ndarray): + return np.outer(psi, psi.conj()).reshape( + (2, 2, 2, 2)).trace(axis1=1, axis2=3) + + def RhoB(psi: np.ndarray): + return np.outer(psi, psi.conj()).reshape(( + 2, 2, 2, 2)).trace(axis1=0, axis2=2) + + psi = self.current_state + # Define the projector operator for the "+1" state + directionA_p1 = np.array(directionsA[0]) + directionA_m1 = np.array(directionsA[1]) + directionB_p1 = np.array(directionsB[0]) + directionB_m1 = np.array(directionsB[1]) + + projectorA_p1 = np.outer(directionA_p1, directionA_p1.conj()) + projectorA_m1 = np.outer(directionA_m1, directionA_m1.conj()) + projectorB_p1 = np.outer(directionB_p1, directionB_p1.conj()) + projectorB_m1 = np.outer(directionB_m1, directionB_m1.conj()) + + if simulate_A: + # Create 4x4 projectors for the two-spin system + projector_p1_s = np.kron(projectorA_p1, np.eye(2)) + projector_m1_s = np.kron(projectorA_m1, np.eye(2)) + projector_p11 = projectorA_p1 + projector_p21 = projectorB_p1 + + # Calculate the reduced density matrix for the first spin + # which is measured + rho1 = RhoA(psi) + else: + projector_p1_s = np.kron(np.eye(2), projectorB_p1) + projector_m1_s = np.kron(np.eye(2), projectorB_m1) + projector_p11 = projectorB_p1 + projector_p21 = projectorA_p1 + rho1 = RhoB(psi) + + # Calculate the probability of this first spin being "+1" + prob_p11 = np.linalg.norm(np.trace(np.dot(projector_p11, rho1))) + # Generate a random number between 0 and 1 + random_number1 = random.uniform(0, 1) + + # Perform the measurement in apparatus direction + sp1 = 1 if random_number1 < prob_p11 else -1 + + # reduce psi projecting on the direction + psi_r = np.dot(projector_p1_s, psi) if sp1 == 1 else \ + np.dot(projector_m1_s, psi) + # Normalize psi_r + psi_r = psi_r / np.linalg.norm(psi_r) + # Calculate the reduced density matrix for the second spin + # with the collapsed wave function for the first system + rho2 = RhoB(psi_r) if simulate_A else RhoA(psi_r) + + # Calculate the probability of "system 2" being "+1" + prob_p12 = np.linalg.norm(np.trace(np.dot(projector_p21, rho2))) + + # Generate a random number between 0 and 1 + random_number2 = random.uniform(0, 1) + + sp2 = 1 if random_number2 < prob_p12 else -1 + return (sp1, sp2) if simulate_A else (sp2, sp1) + + +class MainWindow(QWidget): + + def __init__(self, simul_type: int, measure_both: bool): + super(MainWindow, self).__init__() + self.simul_type = simul_type + self.measure_both = measure_both + + self.initUI() + + self.current_state = None + self.simulation_thread = SimulationThread(simul_type) + self.simulation_thread.result.connect(self.store_simul_spin) + self.simulation_thread.start() + + def initUI(self): + self.setGeometry(300, 300, 800, 600) + match self.simul_type: + case 1: + desc = 'Product state: A = | u > B = | d >' + case 2: + desc = 'Product state: A = 1 / sqrt(2) * (| u > + | d >) ' \ + 'B = | u > / 2 + sqrt(3) / 2 * | d >' + case 3: + desc = 'Singlet state: '\ + '| Psi > = 1 / sqrt(2) ( | ud > - | du > )' + case 4: + desc = 'Triplet state I: '\ + '| Psi > = 1 / sqrt(2) ( | ud > + | du > )' + case 5: + desc = 'Triplet state II: '\ + '| Psi > = 1 / sqrt(2) ( | uu > + | dd > )' + case 6: + desc = 'Triplet state III: '\ + '| Psi > = 1 / sqrt(2) ( | uu > - | dd > )' + case 7: + desc = 'Partially entangled state: '\ + '| Psi > = sqrt(0.6) | ud > - sqrt(0.4) | du >' + + self.setWindowTitle(f"Two quantum spins simulation - {desc}") + + self.opengl_widget = OpenGLWidget( + self, self.simul_type, self.measure_both) + + self.slider1A = QSlider(Qt.Orientation.Horizontal, self) + self.slider1A.setRange(0, 360) + self.slider1A.valueChanged.connect(self.update_rotation_thetaA) + self.label1A = QLabel("Apparatus A Rotation θ: 0", self) + + self.slider2A = QSlider(Qt.Orientation.Horizontal, self) + self.slider2A.setRange(0, 360) + self.slider2A.valueChanged.connect(self.update_rotation_phiA) + self.label2A = QLabel("Apparatus A Rotation 𝜙: 0", self) + + self.buttonA = QPushButton('Make measurement system A', self) + self.buttonA.clicked.connect(self.on_buttonA_clicked) + + self.slider1B = QSlider(Qt.Orientation.Horizontal, self) + self.slider1B.setRange(0, 360) + self.slider1B.valueChanged.connect(self.update_rotation_thetaB) + self.label1B = QLabel("Apparatus B Rotation θ: 0", self) + + self.slider2B = QSlider(Qt.Orientation.Horizontal, self) + self.slider2B.setRange(0, 360) + self.slider2B.valueChanged.connect(self.update_rotation_phiB) + self.label2B = QLabel("Apparatus B Rotation 𝜙: 0", self) + + self.buttonB = QPushButton('Make measurement system B', self) + self.buttonB.clicked.connect(self.on_buttonB_clicked) + + self.layout = QVBoxLayout(self) + self.layout.addWidget(self.opengl_widget) + self.gridlayout = QGridLayout() + self.gridlayout.addWidget(self.slider1A, 0, 0) + self.gridlayout.addWidget(self.label1A, 0, 1) + self.gridlayout.addWidget(self.slider2A, 1, 0) + self.gridlayout.addWidget(self.label2A, 1, 1) + self.gridlayout.addWidget(self.buttonA, 2, 0, 1, 2) + self.gridlayout.addWidget(self.slider1B, 3, 0) + self.gridlayout.addWidget(self.label1B, 3, 1) + self.gridlayout.addWidget(self.slider2B, 4, 0) + self.gridlayout.addWidget(self.label2B, 4, 1) + self.gridlayout.addWidget(self.buttonB, 5, 0, 1, 2) + self.layout.addLayout(self.gridlayout) + + def store_simul_spin(self, current_simul_state: np.ndarray): + self.current_state = current_simul_state + + def update_rotation_thetaA(self, value: float): + self.opengl_widget.update_rotation_thetaA(value) + self.label1A.setText(f"Apparatus A Rotation θ: {value}") + + def update_rotation_phiA(self, value: float): + self.opengl_widget.update_rotation_phiA(value) + self.label2A.setText(f"Apparatus A Rotation 𝜙: {value}") + + def on_buttonA_clicked(self): + self.opengl_widget.measureA(self.current_state) + + def update_rotation_thetaB(self, value: float): + self.opengl_widget.update_rotation_thetaB(value) + self.label1B.setText(f"Apparatus B Rotation θ: {value}") + + def update_rotation_phiB(self, value: float): + self.opengl_widget.update_rotation_phiB(value) + self.label2B.setText(f"Apparatus B Rotation 𝜙: {value}") + + def on_buttonB_clicked(self): + self.opengl_widget.measureB(self.current_state) + + +class CustomHelpFormatter(argparse.HelpFormatter): + def _fill_text(self, text, width, indent): + # Preserve line breaks by not wrapping text + return "\n".join([indent + line for line in text.splitlines()]) + + +def main(): + # Set a fixed seed value + seed_value = 5948 + random.seed(seed_value) + parser = argparse.ArgumentParser(description=description, + formatter_class=CustomHelpFormatter) + parser.add_argument('-t', '--simul_type', help='simulation type', + required=False) + parser.add_argument('-m', '--measure_both', action='store_true', + help='Measure both systems', required=False) + + args = parser.parse_args() + if (args.simul_type): + simul_type = int(args.simul_type) + else: + simul_type = 3 + if (args.measure_both): + measure_both = True + else: + measure_both = False + app = QtWidgets.QApplication(sys.argv) + window = MainWindow(simul_type, measure_both) + window.show() + sys.exit(app.exec()) + + +if __name__ == '__main__': + if sys.version_info[0] < 3: + raise 'Must be using Python 3' + main()