Coverage for gpaw/fdtd/polarizable_material.py: 86%
279 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1"""Part of the module for electrodynamic simulations
3"""
5# flake8: noqa
7from ase.units import Hartree, Bohr
8from gpaw.fd_operators import Gradient
9import numpy as np
11# in atomic units, 1/(4*pi*e_0) = 1
12_eps0_au = 1.0 / (4.0 * np.pi)
14# Base class for the classical polarizable material:
15# -holds various functions at each point in the calculation grid:
16# 1) the dielectric function (permittivity)
17# 2) electric field
18# 3) classical polarization charge density
19# -contains routines for calculating them from each other and/or external potential
21class PolarizableMaterial():
22 def __init__(self, *, sign=-1.0):
23 self.gd = None
24 self.initialized = False
25 self.sign = sign
26 self.messages = []
27 self.components = []
29 def add_component(self, component):
30 self.components.append(component)
32 def permittivity_value(self, omega=0.0):
33 return self.eps_infty + _eps0_au * np.sum(self.beta / (self.bar_omega**2.0 - 1J * self.alpha * omega - omega**2.0), axis=0)
35 def get_static_permittivity(self):
36 return self.gd.collect(self.permittivity_value(0.0))
38 def initialize(self, gd):
39 if self.initialized: # double initialization leads to problems
40 return
41 self.initialized = True
42 self.messages.append("Polarizable Material:")
44 try:
45 self.Nj = max(component.permittivity.Nj for component in self.components)
46 except:
47 self.Nj = 0
49 self.gd = gd
51 # 3-dimensional scalar array: rho, eps_infty
52 self.charge_density = self.gd.zeros()
53 self.eps_infty = np.ones(self.gd.empty().shape) * _eps0_au
55 # 3-dimensional vector arrays:
56 # electric field, total polarization density
57 dims = [3] + list(self.gd.empty().shape)
58 self.electric_field = np.zeros(dims)
59 self.polarization_total = np.zeros(dims)
61 # 4-dimensional vector arrays:
62 # currents, polarizations
63 dims = [3, self.Nj] + list(self.gd.empty().shape)
64 self.currents = np.zeros(dims)
65 self.polarizations = np.zeros(dims)
67 # 4-dimensional scalar arrays:
68 # oscillator parameters alpha, beta, bar_omega, eps_infty
69 dims = [self.Nj] + list(self.gd.empty().shape)
70 self.alpha = np.zeros(dims)
71 self.beta = np.zeros(dims)
72 self.bar_omega = np.ones(dims)
74 # Set the permittivity for each grid point
75 for component in self.components:
76 self.apply_mask(mask = component.get_mask(self.gd),
77 permittivity = component.permittivity)
79 # Restart by regenerating the structures
80 # TODO: is everything always fully reproduced? Should there be a test for it (e.g. checksum[mask])
81 def read(self, reader):
82 for component in reader.classmat.components:
83 name = component['name']
84 arguments = component['arguments']
85 eps_infty = component['eps_infty']
86 eps = component['eps'] # Data as array
87 eval(f'self.add_component({name}(permittivity = Permittivity(data=eps), {arguments}))')
89 # Save information on the structures for restarting
90 def write(self, writer):
91 writer.write(components=[
92 {'name': component.name,
93 'arguments': component.arguments,
94 'eps_infty': component.permittivity.eps_infty,
95 'eps': component.permittivity.data_eVA()}
96 for component in self.components])
98 # Here the 3D-arrays are filled with material-specific information
99 def apply_mask(self, mask, permittivity):
100 for j in range(permittivity.Nj):
101 self.bar_omega[j] = np.logical_not(mask) * self.bar_omega[j] + mask * permittivity.oscillators[j].bar_omega
102 self.alpha [j] = np.logical_not(mask) * self.alpha[j] + mask * permittivity.oscillators[j].alpha
103 self.beta [j] = np.logical_not(mask) * self.beta[j] + mask * permittivity.oscillators[j].beta
105 # Add dummy oscillators if needed
106 for j in range(permittivity.Nj, self.Nj):
107 self.bar_omega [j] = np.logical_not(mask) * self.bar_omega[j] + mask * 1.0
108 self.alpha [j] = np.logical_not(mask) * self.alpha[j] + mask * 0.0
109 self.beta [j] = np.logical_not(mask) * self.beta[j] + mask * 0.0
111 # Print the permittivity information
112 self.messages.append(" Permittivity data:")
113 self.messages.append(" bar_omega alpha beta")
114 self.messages.append(" ----------------------------------------")
115 for j in range(permittivity.Nj):
116 self.messages.append("{:12.6f} {:12.6f} {:12.6f}".format(permittivity.oscillators[j].bar_omega,
117 permittivity.oscillators[j].alpha,
118 permittivity.oscillators[j].beta))
119 self.messages.append(" ----------------------------------------")
120 self.messages.append("...done initializing Polarizable Material")
121 masksum = self.gd.comm.sum_scalar(int(np.sum(mask)))
122 masksize = self.gd.comm.sum_scalar(int(np.size(mask)))
123 self.messages.append("Fill ratio: %f percent" % (100.0 * float(masksum)/float(masksize)))
126 # E(r) = -Grad V(r)
127 # NB: Here -V(r) is used (if/when sign=-1), because in GPAW the
128 # electron unit charge is +1 so that the calculated V(r) is
129 # positive around negative charge. In order to get the correct
130 # direction of the electric field, the sign must be changed.
131 def solve_electric_field(self, phi):
132 for v in range(3):
133 Gradient(self.gd, v, n=3).apply(-1.0 * self.sign * phi, self.electric_field[v])
135 # n(r) = -Div P(r)
136 def solve_rho(self):
137 self.charge_density *= 0.0
138 dmy = self.gd.empty()
139 for v in range(3):
140 Gradient(self.gd, v, n=3).apply(self.polarization_total[v], dmy)
141 self.charge_density -= dmy
143 # P(r, omega) = [eps(r, omega) - eps0] E(r, omega)
144 # P0(r) = [eps_inf(r) - eps0] E0(r) + sum_j P0_j(r) // Gao2012, Eq. 10
145 def solve_polarizations(self):
146 for v in range(3):
147 self.polarizations[v] = _eps0_au * self.beta / (self.bar_omega**2.0) * self.electric_field[v]
148 self.polarization_total = np.sum(self.polarizations, axis=1) + (self.eps_infty - _eps0_au ) * self.electric_field
150 def propagate_polarizations(self, timestep):
151 for v in range(3):
152 self.polarizations[v] = self.polarizations[v] + timestep * self.currents[v]
153 self.polarization_total = np.sum(self.polarizations, axis=1)
155 def propagate_currents(self, timestep):
156 c1 = (1.0 - 0.5 * self.alpha*timestep)/(1.0 + 0.5 * self.alpha*timestep)
157 c2 = - timestep / (1.0 + 0.5 * self.alpha*timestep) * (self.bar_omega**2.0)
158 c3 = - timestep / (1.0 + 0.5 * self.alpha*timestep) * (-1.0) * _eps0_au * self.beta
159 for v in range(3):
160 self.currents[v] = c1 * self.currents[v] + c2 * self.polarizations[v] + c3 * self.electric_field[v]
162 def kick_electric_field(self, timestep, kick):
163 for v in range(3):
164 self.electric_field[v] = self.electric_field[v] + kick[v] / timestep
166# Box-shaped classical material
167class PolarizableBox():
168 def __init__(self, *, corner1, corner2, permittivity):
169 # sanity check
170 assert(len(corner1)==3)
171 assert(len(corner2)==3)
173 self.corner1 = np.array(corner1)/Bohr # from Angstroms to atomic units
174 self.corner2 = np.array(corner2)/Bohr # from Angstroms to atomic units
175 self.permittivity = permittivity
177 self.name = 'PolarizableBox'
178 self.arguments = 'corner1=[{:f}, {:f}, {:f}], corner2=[{:f}, {:f}, {:f}]'.format(corner1[0], corner1[1], corner1[2],
179 corner2[0], corner2[1], corner2[2])
181 # Setup grid descriptor and the permittivity values inside the box
182 def get_mask(self, gd):
184 # 3D coordinates at each grid point
185 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
187 # inside or outside
188 return np.logical_and(np.logical_and( # z
189 np.logical_and(np.logical_and( # y
190 np.logical_and( #x
191 r_gv[:, :, :, 0] > self.corner1[0],
192 r_gv[:, :, :, 0] < self.corner2[0]),
193 r_gv[:, :, :, 1] > self.corner1[1]),
194 r_gv[:, :, :, 1] < self.corner2[1]),
195 r_gv[:, :, :, 2] > self.corner1[2]),
196 r_gv[:, :, :, 2] < self.corner2[2])
199# Shape from atom positions (surrounding region)
200class PolarizableAtomisticRegion():
201 def __init__(self, *, atoms=None, atom_positions=None, distance=0.0, permittivity=None):
203 if atoms is not None:
204 self.atom_positions = np.array(atoms.get_positions())
205 else:
206 self.atom_positions = np.array(atom_positions)
208 # sanity check
209 assert(len(self.atom_positions)>1)
211 self.permittivity = permittivity
212 self.name = 'PolarizableAtomisticRegion'
214 # use the minimum interatomic distance
215 if distance<1e-10:
216 self.distance = np.sqrt(np.min([np.sort(np.sum((self.atom_positions-ap)**2, axis=1))[1] for ap in self.atom_positions]))/Bohr
217 else:
218 self.distance = distance/Bohr # from Angstroms to atomic units
220 dbohr = self.distance*Bohr
221 self.arguments = 'distance = %20.12e, atom_positions=[' % dbohr
222 for ap in self.atom_positions:
223 self.arguments += f'[{ap[0]:20.12e}, {ap[1]:20.12e}, {ap[2]:20.12e}],'
224 self.arguments = self.arguments[:-1] + ']'
226 # Setup grid descriptor and the permittivity values inside the box
227 def get_mask(self, gd):
229 # 3D coordinates at each grid point
230 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
232 # inside or outside
233 mask_tbl = False * np.ones(r_gv.shape[:-1])
234 for ap in self.atom_positions:
235 mask_tbl = np.logical_or(mask_tbl, np.array( (r_gv[:, :, :, 0] - ap[0]/Bohr)**2.0 +
236 (r_gv[:, :, :, 1] - ap[1]/Bohr)**2.0 +
237 (r_gv[:, :, :, 2] - ap[2]/Bohr)**2.0 < self.distance**2))
238 return mask_tbl
240# Sphere-shaped classical material
241class PolarizableSphere():
242 def __init__(self, *, center, radius, permittivity):
243 self.permittivity = permittivity
244 self.center = np.array(center)/Bohr # from Angstroms to atomic units
245 self.radius = radius/Bohr # from Angstroms to atomic units
247 # sanity check
248 assert(len(self.center)==3)
250 self.name = 'PolarizableSphere'
251 self.arguments = f'center=[{center[0]:20.12e}, {center[1]:20.12e}, {center[2]:20.12e}], radius={radius:20.12e}'
253 def get_mask(self, gd):
255 # 3D coordinates at each grid point
256 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
258 # inside or outside
259 return np.array( (r_gv[:, :, :, 0] - self.center[0])**2.0 +
260 (r_gv[:, :, :, 1] - self.center[1])**2.0 +
261 (r_gv[:, :, :, 2] - self.center[2])**2.0 < self.radius**2)
263# Sphere-shaped classical material
264class PolarizableEllipsoid():
265 def __init__(self, *, center, radii, permittivity):
266 # sanity check
267 assert(len(center)==3)
268 assert(len(radii)==3)
270 self.center = np.array(center)/Bohr # from Angstroms to atomic units
271 self.radii = np.array(radii)/Bohr # from Angstroms to atomic units
272 self.permittivity = permittivity
274 self.name = 'PolarizableEllipsoid'
275 self.arguments = f'center=[{center[0]:20.12e}, {center[1]:20.12e}, {center[2]:20.12e}], radii=[{radii[0]:20.12e}, {radii[1]:20.12e}, {radii[2]:20.12e}]'
278 def get_mask(self, gd):
280 # 3D coordinates at each grid point
281 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
283 # inside or outside
284 return np.array( (r_gv[:, :, :, 0] - self.center[0])**2.0/self.radii[0]**2.0 +
285 (r_gv[:, :, :, 1] - self.center[1])**2.0/self.radii[1]**2.0 +
286 (r_gv[:, :, :, 2] - self.center[2])**2.0/self.radii[2]**2.0 < 1.0)
288 # Rod-shaped classical material
289class PolarizableRod():
290 def __init__(self, *, corners, radius, round_corners, permittivity):
291 # sanity check
292 assert(np.array(corners).shape[0]>1) # at least two points
293 assert(np.array(corners).shape[1]==3) # 3D
295 self.name = 'PolarizableRod'
296 self.arguments = 'radius = %20.12e, corners=[' % radius
297 for c in corners:
298 self.arguments += f'[{c[0]:20.12e}, {c[1]:20.12e}, {c[2]:20.12e}],'
299 self.arguments = self.arguments[:-1] + ']'
300 self.arguments += f', round_corners={round_corners}'
302 self.corners = np.array(corners)/Bohr # from Angstroms to atomic units
303 self.radius = radius/Bohr # from Angstroms to atomic units
304 self.round_corners = round_corners
305 self.permittivity = permittivity
307 def get_mask(self, gd):
309 # 3D coordinates at each grid point
310 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
311 ng = r_gv.shape[0:-1]
312 ngv = r_gv.shape
314 a = self.corners[0]
316 mask = False * np.ones(ng)
318 for p in self.corners[1:]:
319 # https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line:
320 # d = |(a-p)-((a-p).n)n| point p, line a+tn (|n|=1)
321 n = (p-a)/np.sqrt((p-a).dot(p-a))
322 v1 = np.array([a[w]-r_gv[:, :, :, w] for w in range(3)]).transpose((1, 2, 3, 0)) # a-p
324 v2 = np.sum(np.array([v1[:, :, :, w]*n[w] for w in range(3)]), axis=0) # (a-p).n
326 v3 = np.array([v2*n[w] for w in range(3)]).transpose(1, 2, 3, 0) # ((a-p).n)n
328 d = np.zeros(ng)
329 for ind, idx in np.ndenumerate(v3[:, :, :, 0]):
330 d[ind] = np.array(v1[ind]-v3[ind]).dot(v1[ind]-v3[ind])
332 # angle between (p-a) and (r-a):
333 pa = p-a # (3)
334 ra = np.array([r_gv[:, :, :, w]-a[w] for w in range(3)]).transpose((1, 2, 3, 0)) # (ng1, ng2, ng3, 3)
335 para = np.sum([pa[w]*ra[:, :, :, w] for w in range(3)], axis=0)
336 ll2 = pa.dot(pa)*np.sum([ra[:, :, :, w]*ra[:, :, :, w] for w in range(3)], axis=0)
337 angle1 = np.arccos(para/(1.0e-9+np.sqrt(ll2)))
339 # angle between (a-p) and (r-p):
340 ap = a-p # (3)
341 rp = np.array([r_gv[:, :, :, w]-p[w] for w in range(3)]).transpose((1, 2, 3, 0)) # (ng1, ng2, ng3, 3)
342 aprp = np.sum([ap[w]*rp[:, :, :, w] for w in range(3)], axis=0)
343 ll2 = ap.dot(ap)*np.sum([rp[:, :, :, w]*rp[:, :, :, w] for w in range(3)], axis=0)
344 angle2 = np.arccos(aprp/(1.0e-9+np.sqrt(ll2)))
346 # Include in the mask
347 this_mask = np.logical_and(np.logical_and(angle1 < 0.5*np.pi, angle2 < 0.5*np.pi),
348 d < self.radius**2.0 )
350 # Add spheres around current end points
351 if self.round_corners:
352 # |r-a| and |r-p|
353 raDist = np.sum([ra[:, :, :, w]*ra[:, :, :, w] for w in range(3)], axis=0)
354 rpDist = np.sum([rp[:, :, :, w]*rp[:, :, :, w] for w in range(3)], axis=0)
355 this_mask = np.logical_or(this_mask,
356 np.logical_or(raDist < self.radius**2.0, rpDist < self.radius**2.0))
358 mask = np.logical_or(mask, this_mask)
360 # move to next point
361 a = p
363 return mask
365class PolarizableTetrahedron():
366 #https://steve.hollasch.net/cgindex/geometry/ptintet.html
367 #obrecht@imagen.com (Doug Obrecht) writes:
368 #
369 # Can someone point me to an algorithm that determines if a point is within a tetrahedron?
370 #
371 # Let the tetrahedron have vertices
372 # V1 = (x1, y1, z1)
373 # V2 = (x2, y2, z2)
374 # V3 = (x3, y3, z3)
375 # V4 = (x4, y4, z4)
376 #
377 #and your test point be
378 #
379 # P = (x, y, z).
380 #Then the point P is in the tetrahedron if following five determinants all have the same sign.
381 #
382 # |x1 y1 z1 1|
383 # D0 = |x2 y2 z2 1|
384 # |x3 y3 z3 1|
385 # |x4 y4 z4 1|
386 #
387 # |x y z 1|
388 # D1 = |x2 y2 z2 1|
389 # |x3 y3 z3 1|
390 # |x4 y4 z4 1|
391 #
392 # |x1 y1 z1 1|
393 # D2 = |x y z 1|
394 # |x3 y3 z3 1|
395 # |x4 y4 z4 1|
396 #
397 # |x1 y1 z1 1|
398 # D3 = |x2 y2 z2 1|
399 # |x y z 1|
400 # |x4 y4 z4 1|
401 #
402 # |x1 y1 z1 1|
403 # D4 = |x2 y2 z2 1|
404 # |x3 y3 z3 1|
405 # |x y z 1|
406 #
407 # Some additional notes:
408 #
409 # If by chance the D0=0, then your tetrahedron is degenerate (the points are coplanar).
410 # If any other Di=0, then P lies on boundary i (boundary i being that boundary formed by the three points other than Vi).
411 # If the sign of any Di differs from that of D0 then P is outside boundary i.
412 # If the sign of any Di equals that of D0 then P is inside boundary i.
413 # If P is inside all 4 boundaries, then it is inside the tetrahedron.
414 # As a check, it must be that D0 = D1+D2+D3+D4.
415 # The pattern here should be clear; the computations can be extended to simplicies of any dimension. (The 2D and 3D case are the triangle and the tetrahedron).
416 # If it is meaningful to you, the quantities bi = Di/D0 are the usual barycentric coordinates.
417 # Comparing signs of Di and D0 is only a check that P and Vi are on the same side of boundary i.
419 def __init__(self, *, corners, permittivity):
420 # sanity check
421 assert(len(corners)==4) # exactly 4 points
422 assert(len(corners[0])==3) # 3D
424 self.name = 'PolarizableTetrahedron'
425 self.arguments = 'corners=['
426 for c in corners:
427 self.arguments += f'[{c[0]:20.12e}, {c[1]:20.12e}, {c[2]:20.12e}],'
428 self.arguments += ']'
430 self.corners = np.array(corners)/Bohr # from Angstroms to atomic units
431 self.permittivity = permittivity
433 def determinant_value(self, x, y, z,
434 x1, y1, z1,
435 x2, y2, z2,
436 x3, y3, z3,
437 x4, y4, z4, ind):
438 mat = np.array([[x1, y1, z1, 1], [x2, y2, z2, 1], [x3, y3, z3, 1], [x4, y4, z4, 1]])
439 mat[ind][:] = np.array([x, y, z, 1])
440 return np.linalg.det(mat)
442 def get_mask(self, gd):
443 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
444 ng = r_gv.shape[0:-1]
445 ngv = r_gv.shape
446 x1, y1, z1 = self.corners[0]
447 x2, y2, z2 = self.corners[1]
448 x3, y3, z3 = self.corners[2]
449 x4, y4, z4 = self.corners[3]
451 mask = np.ones(ng)==0
453 # TODO: associate a determinant for each point, and use numpy tools to determine
454 # the mask without the *very slow* loop over grid points
455 for ind, pt in np.ndenumerate(r_gv[:, :, :, 0]):
456 x, y, z = r_gv[ind][:]
457 d0 = np.array([[x1, y1, z1, 1],
458 [x2, y2, z2, 1],
459 [x3, y3, z3, 1],
460 [x4, y4, z4, 1]])
461 d1 = np.array([[x, y, z, 1],
462 [x2, y2, z2, 1],
463 [x3, y3, z3, 1],
464 [x4, y4, z4, 1]])
465 d2 = np.array([[x1, y1, z1, 1],
466 [x, y, z, 1],
467 [x3, y3, z3, 1],
468 [x4, y4, z4, 1]])
469 d3 = np.array([[x1, y1, z1, 1],
470 [x2, y2, z2, 1],
471 [x, y, z, 1],
472 [x4, y4, z4, 1]])
473 d4 = np.array([[x1, y1, z1, 1],
474 [x2, y2, z2, 1],
475 [x3, y3, z3, 1],
476 [x, y, z, 1]])
477 s0 = np.linalg.det(d0)
478 s1 = np.linalg.det(d1)
479 s2 = np.linalg.det(d2)
480 s3 = np.linalg.det(d3)
481 s4 = np.linalg.det(d4)
483 if (np.sign(s0) == np.sign(s1) or abs(s1) < 1e-12) and \
484 (np.sign(s0) == np.sign(s2) or abs(s2) < 1e-12) and \
485 (np.sign(s0) == np.sign(s3) or abs(s3) < 1e-12) and \
486 (np.sign(s0) == np.sign(s4) or abs(s4) < 1e-12):
487 mask[ind] = True
488 return mask
490# Lorentzian oscillator function: L(omega) = eps0 * beta / (w**2 - i*alpha*omega - omega**2) // Coomar2011, Eq. 2
491class LorentzOscillator:
492 def __init__(self, bar_omega, alpha, beta):
493 self.bar_omega = bar_omega
494 self.alpha = alpha
495 self.beta = beta
497 def value(self, omega):
498 return _eps0_au * self.beta / (self.bar_omega**2 - 1J * self.alpha * omega - omega**2)
500# Dieletric function: e(omega) = eps_inf + sum_j L_j(omega) // Coomar2011, Eq. 2
501class Permittivity:
502 def __init__(self, fname=None, data=None, eps_infty = _eps0_au ):
504 # Initialize to vacuum permittivity
505 self.eps_infty = eps_infty
506 self.Nj = 0
507 self.oscillators = []
509 # Input as data array
510 if data is not None:
511 assert fname is None
512 self.Nj = len(data)
513 for v in range(self.Nj):
514 self.oscillators.append(LorentzOscillator(data[v][0]/Hartree, data[v][1]/Hartree, data[v][2]/Hartree/Hartree))
515 return
517 # Input as filename
518 if fname != None: # read permittivity from a 3-column file
519 with open(fname) as fp:
520 lines = fp.readlines()
522 self.Nj = len(lines)
524 for line in lines:
525 bar_omega = float(line.split()[0]) / Hartree
526 alpha = float(line.split()[1]) / Hartree
527 beta = float(line.split()[2]) / Hartree / Hartree
528 self.oscillators.append(LorentzOscillator(bar_omega, alpha, beta))
529 return
533 def value(self, omega = 0):
534 return self.eps_infty + sum([osc.value(omega) for osc in self.oscillators])
536 def data(self):
537 return [[osc.bar_omega, osc.alpha, osc.beta] for osc in self.oscillators]
539 def data_eVA(self):
540 return [[osc.bar_omega*Hartree, osc.alpha*Hartree, osc.beta*Hartree*Hartree] for osc in self.oscillators]
542# Dieletric function that renormalizes the static permittivity to the requested value (usually epsZero)
543class PermittivityPlus(Permittivity):
544 def __init__(self, fname=None, data=None, eps_infty = _eps0_au, epsZero = _eps0_au, newbar_omega = 0.01, new_alpha = 0.10, **kwargs):
545 Permittivity.__init__(self, fname=fname, data=data, eps_infty=eps_infty)
547 # Convert given values from eVs to Hartrees
548 _newbar_omega = newbar_omega / Hartree
549 _new_alpha = new_alpha / Hartree
551 # Evaluate the new value
552 _new_beta = ((epsZero - self.value(0.0))*_newbar_omega**2.0/_eps0_au).real
553 self.oscillators.append(LorentzOscillator(_newbar_omega, _new_alpha, _new_beta))
554 self.Nj = len(self.oscillators)