{
"cells": [
{
"cell_type": "code",
"execution_count": 164,
"metadata": {},
"outputs": [],
"source": [
"import plotly.express as px\n",
"import numpy as np\n",
"import pandas as pd\n",
"import random\n",
"from tqdm import tqdm"
]
},
{
"cell_type": "code",
"execution_count": 165,
"metadata": {},
"outputs": [],
"source": [
"import plotly.io as pio\n",
"pio.renderers.default = \"notebook_connected\""
]
},
{
"cell_type": "code",
"execution_count": 150,
"metadata": {},
"outputs": [],
"source": [
"def rho(r, a0=1, Z=1):\n",
" return Z*r/a0\n",
"\n",
"def wf1s(sphericalPos, a0=1, Z=1):\n",
" r, theta, phi = sphericalPos\n",
" normConst = (Z/a0)**(3/2)*(1/np.sqrt(np.pi))\n",
" return normConst*np.exp(-rho(r, a0, Z))\n",
"\n",
"def wf2s(sphericalPos, a0=1, Z=1):\n",
" r, theta, phi = sphericalPos\n",
" normConst = (Z/a0)**(3/2)*(1/4/np.sqrt(2*np.pi))\n",
" return normConst*(2-rho(r, a0, Z))*np.exp(-rho(r, a0, Z)/2)\n",
"\n",
"def wf2pz(sphericalPos, a0=1, Z=1):\n",
" r, theta, phi = sphericalPos\n",
" normConst = (Z/a0)**(3/2)*(1/4/np.sqrt(2*np.pi))\n",
" result = normConst*rho(r, a0, Z)*np.exp(-rho(r, a0, Z)/2)*np.cos(theta)\n",
" return result\n",
"\n",
"def wf2px(sphericalPos, a0=1, Z=1):\n",
" r, theta, phi = sphericalPos\n",
" normConst = (Z/a0)**(3/2)*(1/4/np.sqrt(2*np.pi))\n",
" result = normConst*rho(r, a0, Z)*np.exp(-rho(r, a0, Z)/2)*np.sin(theta)*np.cos(phi)\n",
" return result\n",
"\n",
"def wf2py(sphericalPos, a0=1, Z=1):\n",
" r, theta, phi = sphericalPos\n",
" normConst = (Z/a0)**(3/2)*(1/4/np.sqrt(2*np.pi))\n",
" result = normConst*rho(r, a0, Z)*np.exp(-rho(r, a0, Z)/2)*np.sin(theta)*np.sin(phi)\n",
" return result\n",
"\n",
"def wf3d(sphericalPos, a0=1, Z=1):\n",
" r, theta, phi = sphericalPos\n",
" normConst = (Z/a0)**(3/2)*(np.sqrt(2)/81/np.sqrt(np.pi))\n",
" result = normConst*rho(r, a0, Z)*(6-rho(r, a0, Z))*np.exp(-rho(r, a0, Z)/3)*np.cos(theta)\n",
" return result\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 166,
"metadata": {},
"outputs": [],
"source": [
"def spheric2cartesian(sphericalCoordinatesPoint):\n",
" r, theta, phi = sphericalCoordinatesPoint\n",
" x = r*np.sin(theta)*np.cos(phi)\n",
" y = r*np.sin(theta)*np.sin(phi)\n",
" z = r*np.cos(theta)\n",
" return np.array([x, y, z])\n",
"\n",
"from scipy.special import gammainc\n",
"def sample(center,radius,n_per_sphere):\n",
" r = radius\n",
" ndim = center.size\n",
" x = np.random.normal(size=(n_per_sphere, ndim))\n",
" ssq = np.sum(x**2,axis=1)\n",
" fr = r*gammainc(ndim/2,ssq/2)**(1/ndim)/np.sqrt(ssq)\n",
" frtiled = np.tile(fr.reshape(n_per_sphere,1),(1,ndim))\n",
" p = center + np.multiply(x,frtiled)\n",
" return p\n",
"\n",
"def cartesian2spheric(cartesianCoordinate):\n",
" x, y, z = cartesianCoordinate\n",
" x_2 = x**2\n",
" y_2 = y**2\n",
" z_2 = z**2\n",
" r = np.sqrt(x_2 + y_2 + z_2)\n",
" phi = np.arctan2(y, x)\n",
" theta = np.arccos(z/r)\n",
" return np.array([r, theta, phi])\n",
"\n",
"def filter_point(point, distributionFunction, power=1):\n",
" if random.uniform(0, 1) <= distributionFunction(point)**power:\n",
" return True\n",
" else:\n",
" return False\n",
"\n",
"def simulate_particle_position(wavefunction, nPoints = 200, radius = 10):\n",
" resultingDataset = []\n",
" while len(resultingDataset) < nPoints:\n",
" #generar puntos\n",
" seedPoints = sample(np.zeros(3), radius=radius, n_per_sphere=50)\n",
" #transformarlos a polar\n",
" sphericalPoints = []\n",
" for point in seedPoints:\n",
" sphericalPoints.append(cartesian2spheric(point))\n",
" #filtrar con función de distribución\n",
" for point in sphericalPoints:\n",
" if filter_point(point, wavefunction, power=2):\n",
" resultingDataset.append(point)\n",
" return resultingDataset"
]
},
{
"cell_type": "code",
"execution_count": 169,
"metadata": {},
"outputs": [],
"source": [
"def simulateWavefunction(wavefunction, orbitalName, radius = 12, nPoints = 1000):\n",
" cartesianSimulatedPoints = []\n",
" for point in simulate_particle_position(wavefunction, nPoints=nPoints, radius=radius):\n",
" cartesianSimulatedPoints.append(spheric2cartesian(point))\n",
" simulatedPointsDF = pd.DataFrame(cartesianSimulatedPoints, columns=[\"x\", \"y\", \"z\"])\n",
" simulatedPointsDF[\"Orbital\"] = [orbitalName for i in range(len(simulatedPointsDF))]\n",
" return simulatedPointsDF\n",
"\n",
"DFList = []\n",
"\n",
"# Simular 1s\n",
"DFList.append(simulateWavefunction(wf1s, \"1s\"))\n",
"\n",
"# Simular 2s\n",
"#DFList.append(simulateWavefunction(wf2s, \"2s\"))\n",
"\n",
"# Simular 2pz\n",
"#DFList.append(simulateWavefunction(wf2pz, \"2pz\"))\n",
"\n",
"# Simular 2px\n",
"#DFList.append(simulateWavefunction(wf2px, \"2px\"))\n",
"\n",
"# Simular 2py\n",
"#DFList.append(simulateWavefunction(wf2py, \"2py\"))\n",
"\n",
"# Simular 3d\n",
"DFList.append(simulateWavefunction(wf3d, \"3d\"))"
]
},
{
"cell_type": "code",
"execution_count": 170,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = px.scatter_3d(pd.concat(DFList), x='x', y='y', z='z', color=\"Orbital\",\n",
" opacity=0.7)\n",
"fig.update_traces(marker={'size':2})\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": 154,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import plotly.graph_objects as go \n",
"\n",
"resolution = 240\n",
"phi = np.linspace(0, 2*np.pi, resolution)\n",
"theta = np.linspace(0, np.pi, resolution)\n",
"\n",
"phiGrid, thetaGrid = np.meshgrid(phi, theta)\n",
"\n",
"r = np.cos(thetaGrid)**2 # 2p_z\n",
"\n",
"x = r*np.sin(thetaGrid)*np.cos(phiGrid)\n",
"y = r*np.sin(thetaGrid)*np.sin(phiGrid)\n",
"z = r*np.cos(thetaGrid)\n",
"\n",
"surface = go.Surface(x=x, y=y, z=z)\n",
"figure = go.Figure(data=[surface])\n",
"\n",
"figure.show()"
]
},
{
"cell_type": "code",
"execution_count": 163,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import plotly.graph_objects as go \n",
"\n",
"resolution = 240\n",
"phi = np.linspace(0, 2*np.pi, resolution)\n",
"theta = np.linspace(0, np.pi, resolution)\n",
"\n",
"phiGrid, thetaGrid = np.meshgrid(phi, theta)\n",
"\n",
"r = (np.sin(thetaGrid)*np.sin(phiGrid))**2 # 2p_\n",
"\n",
"x = r*np.sin(thetaGrid)*np.cos(phiGrid)\n",
"y = r*np.sin(thetaGrid)*np.sin(phiGrid)\n",
"z = r*np.cos(thetaGrid)\n",
"\n",
"surface = go.Surface(x=x, y=y, z=z)\n",
"figure = go.Figure(data=[surface])\n",
"\n",
"figure.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.10.7 ('plotting')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.8"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "9ffe10734b8d1e9fe694779c050d4b860d6cc9a014de59011a2246c64fa2ec74"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}