from __future__ import division, print_function
from haslib.science_basics import A, m_He4
from multiprocessing import Pool
from tabulate import tabulate
import itertools
import numpy as np
from haslib.icc_2d_2 import iCC_2D_2 as icc2
[docs]def f(p):
"""
Calculate elastic intensity for a given peak
Args:
p: Parameters(a, h, D, xi, Ei, theta, g1, g2, gopen,
gclose, sta_z, end_z, st_p_wl, cc, maxrad,
m, peak) for CC simulation
Returns:
scalar: intensity for given peak
"""
# unpack parameters
(a, h, D, xi, Ei, theta, g1, g2,
gopen, gclose, sta_z, end_z, st_p_wl,
cc, maxrad, m, peak) = p
# use elastic simulation only
floq = 0
AP = 0
phon = [0, 0, 0]
# simulate
iccResult = icc2(h, AP, a, D, xi, Ei,
theta, g1, g2,
gopen, gclose, phon,
floq, sta_z, end_z, st_p_wl,
CC_border=cc, maxrad=maxrad, verbose=False)
# get result for specified peak
spec_0 = iccResult[(iccResult[:, 0] == peak[0]) &
(iccResult[:, 1] == peak[1])]
# return peak intensity
return spec_0[:, 6]
[docs]def auto_conv(D, xi, h, beam, theta, crystal, maxrad=20,
cpus=1, m=m_He4, error=1e-7, gclosemax=70,
st_p_wl_max=120, zrange=(20, 50), peak=(0, 0),
verbose=False, rel=False):
"""
Find convergent parameter z, gclose and steps per wavelength
Args:
D: potentialdepth in meV
xi: stiffness in 1/A
h: corrugation
beam: beam energy in meV
theta: in degrees
crystal: Crystal class see haslib.crystal_class
maxrad: max radius(Default: 20)
cpus: cpus to use in worker Pool(Default: 1)
m: beam mass(Default: m_He4)
error: minimum accuracy(Default: 1e-7)
gclosemax: maximum closed channels to check(Default: 70)
st_p_wl_max: maximum steps per wavelength to check(Default: 120)
zrange: tuple of (start_z, end_z) to check in A(Default: (20, 50))
peak: tuple of two indices for g1 and g2(Default: (0, 0))
verbose: show results for each iteration(Default: False)
rel: use error as relative error(Default: False)
Returns:
array: array(peak_intensity, gclose, end_z, start_z, st_p_wl)
"""
# get g vectors for given crystal
g1 = crystal.b1[0:2]
g2 = crystal.b2[0:2]
# use all open channels
gopen = -1
# define parameter range to check
gclose = np.arange(0, gclosemax, 2).astype(int)
end_z = np.arange(0, zrange[1], 1)*A
sta_z = -np.arange(0, zrange[0], 1)*A
st_p_wl = np.arange(10, st_p_wl_max, 5)
# save all parameterranges as an array
params = [gclose, end_z, sta_z, st_p_wl]
# set the last values of the defines ranges as default/starting parameters
default = [np.array([param[-1]]) for param in params]
# loop over every type of parameter
for ind in range(len(params)):
# take the default parameters and replace the current
# parameter type with its range
p = default[:]
p[ind] = params[ind]
# create search space out of 3 default parameters
# and one parameter range
search = list(itertools.product(*p))
# create parameter list for each item in the search space
# and store in payload
# these parameters will be given to the cc-simulation
cc = 10e-5
payload = []
for p in search:
gclose, end_z, sta_z, st_p_wl = p
payload.append([crystal.a1[0], h, D, xi, beam, theta, g1,
g2, int(gopen), int(gclose), sta_z,
end_z, int(st_p_wl), cc, maxrad, m, peak])
# check if multiple threads should be used and and run simulations
if cpus >= 1:
p = Pool(cpus)
res = p.map(f, payload)
p.close()
else:
res = map(f, payload)
# create array from search space list
search = np.vstack(search)
# plot current results if verbose is defined
if verbose:
print(tabulate(np.hstack([res, search]),
headers=["gclose", "end", "start", "step"]))
# find the results which are converged
if rel:
# use a relative error
good = abs(res/res[-1]-1) <= error
else:
# use an absolute error
good = abs(res-res[-1]) <= error
# search for position where the result converged
ind2 = len(res)-1
for ind3 in range(len(res)):
# check if all results are good after index ind3
if good[ind3:].all():
ind2 = ind3
break
# if there is no convergence raise an exception
if ind2 == len(res)-1:
d = {0: "gclose", 1: "end_z", 2: "start_z", 3: "steps"}
raise Warning("Boundries too low - {}({})".format(d[ind], ind))
# set the converged parameter as the new default parameter
default[ind] = np.array([search[ind2, ind]])
# print default parameters if verbose=True
if verbose:
print(tabulate(np.array(default)[None, :],
headers=["gclose", "end", "start", "step"]))
# return default parameters after every parameter has converged
return np.array(default)