New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Inaccurate results #17
Comments
Here is a multipole expansion that does give accurate results, within the theoretical error bounds: import math
import numpy as np
from scipy.special import lpmv
# Prepare the charge distribution dict for the MultipoleExpansion object:
charge_dist = {
'discrete': True, # point charges are discrete charge distributions
'charges': [
{'q': 1, 'xyz': (0, 0, 1)},
{'q': -1.53, 'xyz': (0, 0, -1)},
{'q': -1, 'xyz': (0, 3.2, -1)}
]
}
charges = np.array([c['q'] for c in charge_dist['charges']])
positions = np.array([c['xyz'] for c in charge_dist['charges']])
def true_potential(x, y, z):
sum_ = 0.0
for c in charge_dist['charges']:
p = np.array(c['xyz'])
r = np.linalg.norm(np.array([x,y,z]) - p)
sum_ += c['q']/r
return sum_
l_max = 8
lm_indices = np.array([(l, m) for l in range(0, l_max+1) for m in range(-l, l+1)])
def harmonic(l, m, phi, theta):
f = math.factorial
return math.sqrt( f(l - abs(m)) / f(l + abs(m)) ) * lpmv(abs(m), l, np.cos(theta)) * np.exp(1j * m * phi)
def cartesian_to_spherical(x, y, z):
r = np.linalg.norm([x,y,z])
if(np.isclose(r, 0.)):
return 0., 0., np.pi/2
theta = math.acos(z/r)
phi = math.atan2(y, x)
return r, phi, theta
def multipole_coeff(charges, positions, l, m):
N = len(charges)
assert charges.shape == (N,)
assert positions.shape == (N, 3)
coeff = 0.
for q, p in zip(charges, positions):
r, phi, theta = cartesian_to_spherical(*p)
coeff += q * r**l * harmonic(l, -m, phi, theta)
return coeff
def multipole_expansion(charges, positions):
exp = np.zeros( len(lm_indices), dtype=np.cdouble )
for i, (l, m) in enumerate(lm_indices):
exp[i] = multipole_coeff(charges, positions, l, m)
return exp
def expansion_to_voltage(expansion, target):
voltage = 0.
for i, (l, m) in enumerate(lm_indices):
r, phi, theta = cartesian_to_spherical(*target)
voltage += harmonic(l, m, phi, theta) * expansion[i] / (r**(l+1))
return np.real(voltage)
expansion = multipole_expansion(charges, positions)
x, y, z = 30.5, 30.6, 30.7
correct = true_potential(x, y, z)
print('Accuracy: ', expansion_to_voltage(expansion, (x, y, z))/true_potential(x, y, z) - 1)
max_x = np.max(np.linalg.norm(positions, axis=1))
y = np.linalg.norm( (x, y, z) )
print('Error bound: ', np.sum(np.abs(charges))/(y - max_x) * (max_x/y)**l_max)
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Describe the bug
The packages give inaccurate answers for some charge distributions.
To Reproduce
Consider the following code.
The accuracy is only 5e-3, and does not improve when increasing
l_max
.Expected behavior
I expect an accuracy down to floating point limit (1e-15) for large value of
l_max
.Desktop (please complete the following information):
Linux, Debian, Python 3.9
Additional context
I believe there is a sign error related to the spherical harmonics functions used by this package. To be specific, I believe for some values of l, m the sign of the imaginary part of the spherical harmonics functions used is incorrect. I'm still investigating.. I will let you know when I figure it out.
The text was updated successfully, but these errors were encountered: