<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Germanium telluride ferroelectric switching contrast
This is an example script to show how to use pyasf to calculate the switching contrast of several Bragg reflections as a function of energy, e.g., near an absorption edge.
Step1:  Let's look at the positions in the asymmetric unit cell...
Step2:  Comparing with the high-temperature cubic unit cell one can see that the atoms are displaced away from the high symmetry positions (1/4 and 3/4) by the same amount ($\delta = 0.0097$ in z direction).
In Practice, this displacement may change significantly (piezo- or pyroelectricity?) or inverted (ferroelectric switching).
Therefore, it is useful to define a variable which describes this displacement
Step3:  Yes, 1/4 is more precise than 0.25. There is a numerical limit to the precision so 0.2499999999999 = 0.25.
Using sympy Symbols is more precise and in crystallography it is important especially when describing special (high symmetry) Wyckoff positions.
The subs dictionary contains the numerical values used as substitute for final evaluation.
Let's calculate the site symmetry and build the unit cell
Step4:  Now we have a positions dictionary containing the whole unit cell and some other stuff
Step5:  Now we define an energy range in eV an play around a bit...
Step6:  So we see already there is a significant contrast when switching the ferroelectric domains (inverting the displaecment). In practice, a strong relative change in intensity is easier to measure which is why weak reflections are often better...
Let's now do a proper calculation of switching contrast for all reflections in a certain range of Q (momentum transfter). To increase speed, we don't recalculate the structure factor each time (is done in .DAFS). The crystal is not modified except for the numerical value of $\delta$, so no refreshing is needed...
<ASSISTANT_TASK:>
Python Code:
%matplotlib notebook
import pylab as pl
import pyasf
import sympy as sp # symbolic computing
from IPython.display import display, Math
print_latex = lambda x: display(Math(sp.latex(x)))
# sp.init_printing()
ls
GeTe = pyasf.unit_cell("ICSD_188458.cif") # initialize crystal structure object
Explanation: Germanium telluride ferroelectric switching contrast
This is an example script to show how to use pyasf to calculate the switching contrast of several Bragg reflections as a function of energy, e.g., near an absorption edge.
End of explanation
print_latex(GeTe.AU_positions)
Explanation: Let's look at the positions in the asymmetric unit cell...
End of explanation
delta = sp.Symbol("delta") # Just a symbol
delta_val = 0.0097 # The initial value of the symbol
GeTe.subs[delta] = delta_val # needed later
GeTe.AU_positions["Ge1"][2] = sp.S("1/4") - delta # redefine the positions
GeTe.AU_positions["Te1"][2] = sp.S("3/4") + delta # redefine the positions
Explanation: Comparing with the high-temperature cubic unit cell one can see that the atoms are displaced away from the high symmetry positions (1/4 and 3/4) by the same amount ($\delta = 0.0097$ in z direction).
In Practice, this displacement may change significantly (piezo- or pyroelectricity?) or inverted (ferroelectric switching).
Therefore, it is useful to define a variable which describes this displacement:
End of explanation
GeTe.get_tensor_symmetry()
GeTe.build_unit_cell()
Explanation: Yes, 1/4 is more precise than 0.25. There is a numerical limit to the precision so 0.2499999999999 = 0.25.
Using sympy Symbols is more precise and in crystallography it is important especially when describing special (high symmetry) Wyckoff positions.
The subs dictionary contains the numerical values used as substitute for final evaluation.
Let's calculate the site symmetry and build the unit cell:
End of explanation
print_latex(GeTe.positions)
GeTe.multiplicity
Explanation: Now we have a positions dictionary containing the whole unit cell and some other stuff:
End of explanation
E = pl.linspace(7000, 30000, 1001)
fig, ax = pl.subplots(1,2, figsize=(9.5,3))
GeTe.subs[delta] = delta_val 
F = GeTe.DAFS(E, (0,0,3), Uaniso=False) # because no anisotropic Uij is defined...
ax[0].plot(E, abs(F)**2)
GeTe.subs[delta] = -delta_val 
F = GeTe.DAFS(E, (0,0,3), Uaniso=False)
ax[0].plot(E, abs(F)**2)
GeTe.subs[delta] = delta_val 
F = GeTe.DAFS(E, (3,3,6), Uaniso=False)
ax[1].plot(E, abs(F)**2)
GeTe.subs[delta] = -delta_val 
F = GeTe.DAFS(E, (3,3,6), Uaniso=False)
ax[1].plot(E, abs(F)**2)
print_latex(GeTe.F_0)
Explanation: Now we define an energy range in eV an play around a bit...
End of explanation
results = [] # some container to store results
GeTe.calc_structure_factor(Temp = False, DD=False, DQ=False) # calculate (`cache`) all structure factors
qmax = 1.5 # q=2*sin(theta)/lambda 
for R in GeTe.iter_rec_space(qmax): # iterate over all reflections in Ewald sphere
    #print R,
    GeTe.subs[delta] = delta_val
    F1 = GeTe.DAFS(E, R,  force_refresh=False, fwhm_ev=1)
    GeTe.subs[delta] = -delta_val
    F2 = GeTe.DAFS(E, R,  force_refresh=False, fwhm_ev=1)
    # minimum energy necessary to reach reflection:
    ## .theta is a symbolic expression
    ## .energy is just a symbol as well
    Emin = float(sp.sin(GeTe.theta.subs(GeTe.subs)) * GeTe.energy)
    
    if Emin>20000.:
        continue # skip if energy too high
    ind = E>Emin
    I1 = abs(F1)**2 # calculate intensity
    I2 = abs(F2)**2
    I1 = I1[ind]
    I2 = I2[ind]
    
    if I1.max()<1e-2 and I2.max()<1e-2: # Intensity too small?
        continue
    
    c = abs(pl.log10(I1/I2)) # relative contrast
    idx_max = c.argmax()
    
    cmax = 10**c[idx_max] # get the maximum contrast
    Imax = 0.5 * (I1+I2)[idx_max] # ... and its intensity
    Emax = E[ind][idx_max] # ... and its energy
    theta_max = GeTe.theta_degrees(Emax) # ... as well as the theta value
    h,k,l = R
    
    results.append([h, k, l, cmax, Imax, Emax, theta_max])
    print(results[-1])
    
    
results = pl.array(results)
isort = results[:,3].argsort() 
results = results[isort[::-1]] # sort after contrast
columns=("h", "k", "l", "contrast", "intensity", "energy", "theta")
pl.savetxt("switch_contrast.txt", results,
           fmt="  %5i %5i %5i %10.3f %10.1f %10.1f %10.2f",
           header = "%5s %5s %5s %10s %10s %10s %10s"%columns)
results
Explanation: So we see already there is a significant contrast when switching the ferroelectric domains (inverting the displaecment). In practice, a strong relative change in intensity is easier to measure which is why weak reflections are often better...
Let's now do a proper calculation of switching contrast for all reflections in a certain range of Q (momentum transfter). To increase speed, we don't recalculate the structure factor each time (is done in .DAFS). The crystal is not modified except for the numerical value of $\delta$, so no refreshing is needed...
End of explanation
<END_TASK>