Atomorbitale darstellen führt zur Fehlermeldung

[BrickBardo]

Hallo zusammen,

ich habe ein Programm zur grafischen Darstellung von Atomorbitalen geschrieben (wobei ich mich auch an Code von anderen bedient habe). Wenn ich das Programm versuche auszuführen, kommt die Fehlermeldung:

TypeError: only size-1 arrays can be converted to Python scalars

Googeln hat ergeben, dass offensichtlich in einer Funktion ein Skalar, anstatt eines Arrays erwartet wird. Ich verstehe allerdings nicht warum. Die Zeilen des Fehlers sind offensichtlich immer die Funktionsaufrufe bzw. die Definitionen selber. Habe das mit #Fehler im Code markiert.

Code: Alles auswählen

import math as ma
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm, genlaguerre

Z   = 1  # Anzahl der Elementarladungen
a_0 = 1 # Bohrscher Atomradius

#n = 1 #Hauptquantenzahl
#l = 0 #Drehimpulsquantenzahl
#m = 0 #Magnetquantenzahl

#Radialfunktion

def R(n, l, r):
  rho = (2*Z*r)/(n*a_0)
  wurzel = ma.sqrt(((rho/r)**3) * (ma.factorial(n-l-1))/(2*n*ma.factorial(n+l)))
  L = np.polyval(genlaguerre(n-l-1, 2*l+1), rho) #Fehler
  return wurzel * ma.e(-rho/2) * rho**l * L

#Wellenfunktion

def psi(n, l, m, r, phi, theta):
  R_nl = R(n,l,r) #Fehler
  Y_lm = sph_harm(m,l,phi,theta)
  return R_nl * Y_lm

#Zeichnen des Orbitals

def Orbital(n, l, m, x, z, r, phi, theta):
  plt.figure()
  s_nlm = np.abs(psi(n,l,m,r,phi,theta)) #Fehler
  plt.pcolormesh(x, z, s_nlm)

n_max = 4

for n in range(1,n_max+1):
  for l in range(n):
    plot_range = (4*n + 4*l) * a_0
    x_1d = np.linspace(-plot_range,plot_range)
    y = 0
    z_1d = np.linspace(-plot_range,plot_range)
    x,z = np.meshgrid(x_1d,z_1d)
    r = np.sqrt(x**2 + y**2 + z**2)
    theta = np.arctan2(np.sqrt(x**2+y**2), z)
    phi = np.arctan2(y, x)
    for m in range(0, l+1):
      Orbital(n,l,m,x,z,r,phi,theta) #Fehler

plt.show()

[simplesimon]

Der folgende Code führt zum Ziel:

Code: Alles auswählen

def R(n, l, r):
  rho = (2*Z*r)/(n*a_0)
  wurzel = np.sqrt(((rho/r)**3) * (ma.factorial(n-l-1))/(2*n*ma.factorial(n+l)))
  L = np.polyval(genlaguerre(n-l-1, 2*l+1), rho) #Fehler
  
  return wurzel * np.exp(-rho/2) * rho**l * L

[BrickBardo]

Ich danke dir vielmals. Das zwei Funktionen, ob wohl sie das gleiche bewirken sollen, eben nicht das gleiche tun, nur weil sie aus unterschiedlichen Bibliotheken stammen, ist schon erstaunlich. Ist das ein bekanntes Problem der math-Bibliothek?

[__blackjack__]

@BrickBardo: Das ist kein Problem der `math`-Bibliothek. Die stellt Funktionen zur Verfügung, die mit `float`-Werten arbeiten bzw. mit Werten die sich in solche umwandeln lassen. Numpy-Arrays sind halt keine `float`-Werte und nur Numpy-Arrays mit nur einem Element lassen sich in einen `float`-Wert umwandeln (wenn sich das eine Element umwandeln lässt). Das ist dann dieses eine Element beziehungsweise dessen `float`-Wert. Genau das sagt auch die Fehlermeldung die Du bekommen hast.

Woher sollte denn die Funktion `math.sqrt()` aus der Standardbibliothek wissen was sie mit einem Numpy-Array anfangen soll, das aus einer ganz anderen Bibliothek von ganz anderen Programmierern stammt?

[BrickBardo]

Ich bin ja noch sehr unerfahren mit Python (und mit Programmierung im allgemeinen), daher bin ich davon ausgegangen, dass die Wurzelfunktion aus Numpy und Math beide ein Skalar zurückgeben und dass diese ununterscheidbar sind. Und wobei ich mich da auch Frage, warum die Fakultät aus der math-Bibliothek nicht das gleiche Problem hervorruft. Wurzel und Fakultät sind beides ja nur Zahlen.

[__blackjack__]

@BrickBardo: Aus dem Beitrag wird jetzt nicht so ganz klar ob die Annahme das beide ein Skalar zurück geben immer noch besteht — `numpy.sqrt()` tut das nicht:

Code: Alles auswählen

In [275]: import numpy as np

In [276]: np.sqrt(np.array([25, 42, 23, 4711]))                                 
Out[276]: array([ 5.        ,  6.4807407 ,  4.79583152, 68.63672486])

Was hättest Du denn angenommen was der skalare Wert bei so einer Wurzel aus einem Array wäre?

Und die Fakultät aus dem `math`-Modul ruft genau das gleiche Problem hervor. Auch die `math.factorial()`-Funktion kann man nicht mit einem Array füttern das mehr als ein Element enthält:

Code: Alles auswählen

In [277]: import math                                                           

In [278]: math.factorial(np.array([25, 42, 23, 4711]))                          
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-278-656a7d87906b> in <module>
----> 1 math.factorial(np.array([25, 42, 23, 4711]))

TypeError: only size-1 arrays can be converted to Python scalars

[BrickBardo]

So, ich bin es wieder. Hatte leider einiges zu tun und konnte nicht sofort antworten. Wenn ich das richtig sehe, dann wird die Wurzelfunktion aus der math-Bibliothek lediglich auf eine Zahl angewendet, während die aus der numpy-Bibliothek auf alle Zahlen eines Arrays angewendet wird. Dementsprechend passte der Datentyp nicht, was zu der Fehlermeldung führte. Was mir zuvor nicht klar war, war die Tatsache, dass da überhaupt eine Array vorliegt. Aber nun gut, das Programm läuft ja nun und ich habe jetzt versucht nach der zweidimensionalen Darstellung auch eine dreidimensionale zu programmieren.

[BrickBardo]

Ich beziehe mich auf diesen Thread, weil ich durch den Code von Magben stark inspiriert wurde (eine schöne Formulierung für: Ich habe dort einiges abgeschrieben). In seiner 2D-Darstellung hat Magben ja auch konkret auf Atomorbitale bezug genommen, bei der der 3D-Darstellung jedoch nur auf Kugelflächenfunktionen, also ohne die Radialfunktion. Daher habe ich versucht beides zu kombinieren. Allerdings ist die 3D-Darstellung meiner Orbitale stark deformiert und woran es genau liegt, kann ich nicht so genau sagen. Ich vermute, dass es bei mir an einer fehlerhaften Deklaration der Variable r liegt. Hat da vielleicht jemand eine Idee?

Code: Alles auswählen

import math as ma
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import sph_harm, genlaguerre

Z = 1  # Anzahl der Elementarladungen
a_0 = 1 # Bohrscher Atomradius

n = 3 #Hauptquantenzahl
l = 2 #Drehimpulsquantenzahl
m = 0 #Magnetquantenzahl

theta_1d = np.linspace(0, np.pi)
phi_1d = np.linspace(0, 2*np.pi)
theta_2d, phi_2d = np.meshgrid(theta_1d, phi_1d)
xyz_2d = np.array([np.sin(theta_2d) * np.sin(phi_2d),
                   np.sin(theta_2d) * np.cos(phi_2d),
                   np.cos(theta_2d)]) 
plot_range = (4*n + 4*l) * a_0
x_1d = np.linspace(-plot_range,plot_range)
y = 0
z_1d = np.linspace(-plot_range,plot_range)
x,z = np.meshgrid(x_1d,z_1d)
r = np.sqrt(x**2 + y**2 + z**2)

#Radialfunktion

def R(n, l, r):
  rho = (2*Z*r)/(n*a_0)
  wurzel = np.sqrt(((rho/r)**3) * (ma.factorial(n-l-1))/(2*n*ma.factorial(n+l)))
  L = np.polyval(genlaguerre(n-l-1, 2*l+1), rho)
  return wurzel * np.exp(-rho/2) * rho**l * L

#Wellenfunktion

def psi(n, l, m, r, phi_2d, theta_2d):
  R_nl = R(n,l,r)
  Y_lm = sph_harm(m,l,phi_2d,theta_2d)
  return R_nl * Y_lm

#Zeichnen des Orbitals

def Orbital(n, l, m, x, z, r, phi_2d, theta_2d):
  plt.figure()
  ax = plt.gca(projection = "3d")
  plt.title(r"$\Psi_{%i%i%i}$" % (n,l,m))
  s_nlm = np.abs(psi(n,l,m,r,phi_2d,theta_2d))*xyz_2d
  ax.plot_surface(s_nlm[0], s_nlm[1], s_nlm[2])

Orbital(n,l,m,x,z,r,phi_2d,theta_2d)
plt.show()