Rocking curve computation of Germanium

[Didou09]

Motivations:

• Thanks to @adriendaros, we now have a routine to compute the rocking curve of Quartz crystals, for any cut, taking into account a small miscut angle and ambient temperature changes
• It would be nice to extend that routine so it can also compute rocking curves for Germanium

Proposition:

• @cjperks7, from MIT, proposes to handle that task.
• He should be able to do it under the supervision of @adriendaros and myself @Didou09

Getting started:

• Create a branch called Issue654_Germanium, from devel, that will be dedicated to this particular task
• The current routine is in tofu/spectro/_rockingcurve.py
• The main routine is compute_rockingcurve()
• It currently takes Miller indices, assuming the crystal is Quartz

Proposed todo list:

• For @adriendaros (preparations):
  • docstr: add all the references you used and advise to use at the end of the docstring, so Conor can find them easily
  • Tutorial: write here (in the comments below) a tutorial on how to use the routine
  • input arguments: add an extra keyword argument dcryst, it should be a dictionnary containing at least:
    • name: the name of the crystal material (ex: Quartz, Germanium)
    • symbol: a short version of the name of the of the crystal material (ex: Quartz, Ge)
    • all the geometrical / atomic parameters useful for the computation,
  • I advise to store all the default crystal dictionnaries in a dedicated new module tofu/spectro/_rockingcurve_def.py as a global variable:

_DCRYST = {                                                                     
    'Quartz': {                                                                 
        'name': 'Quartz',                                                       
        'symbol': 'qz',                                                         
        'atoms': ['Si', 'O'],                                                   
        'meshtype': 'hexagonal',                                                
        '',                                                                     
    },                                                                          
                                                                                
    'Ge': {                                                                     
        'name': 'Germanium',                                                    
        'symbol': 'Ge',                                                         
    },                                                                          
}                 

• For @cjperks7:
  • Fill in the crystal dict for Ge, following Adrien's guidance and template
  • Test the routine with the Ge crystal, report here any bug / new feature needed

[adriendaros]

Motivation:

The compute_rockingcurve() routine available at tofu/spectro/_rockingcurve.py makes it possible to compute the diffraction pattern of specific crystals (material, Miller indices, etc). In order to symplify the code, all the main characteristics of the crystals implemented have been relocated into tofu/spectro/_rockingcurve_def.py, creating specific dictionnaries of crystals.

Example:

See below an example for the ArXVII crystal:

How it is written:

All the geometrical and atomic parameters useful for the computation are stored into sub-dictionnaries such as:

• the Miller indices 'Miller indices';
• the type of mesh 'Meshtype';
• the inter-atomic distances within the elementary unit volume 'Inter-atomic distances';
• etc.

All values which cannot be directly written into the dictionnary (see None values in front of each sub-dictionnary) are computed in the following parts of the file:

• The positions of atoms in the elementary unit volume to fill the sub-dict ['mesh positions']['atoms']:

• The elementary box volume and the associated inter-reticular spacing, filling the sub-dict ['Volume'] & ['Inter-reticular spacing']:

• The whole crystal structure factor computation depending on the atoms componing it:

According to the last point, only the phases are stored inside the dictionnary. The formulas of the linear absorption coefficients and the components of the atomic scattering factor are called directly from the compute_rockingcurve() routine because of the dependency on the wavelength.

[adriendaros]

The use of the _compute_rockingcurve() routine works interactively with the tofu/spectro/_rockingcurve_def.py file. The uploading of a CrystalClass is still needed. The line code launching the routine is showed below:

The arguments showed here are sufficient:

• the name of the sub-dict related to the wanted crystal (here 102-Quartz but also available 110-Quartz)
• the wavelength targeted lamb (already in Angstroms)
• the lattice modifications such as the miscut angle use_non_parallelism and its range in radians and also the temperature changes causing thermal expansion of the reticular planes therm_exp and its range in Celsius degrees.

For more information about the entry arguements, everything is explained at the beginning of the tofu/spectro/_rockingcurve.py routine:

Results:

• Quartz crystal, (1, 1, 0), supposed ideal:

• Quartz crystal, (1, 0, 2), supposed ideal:

[Didou09]

Well done @adriendaros !

If I may, I would suggest the following minor changes in the dictionnary:

• When comments include quantities, include instead that quantity as a number (i.e.: can later be loaded as a parameter)
  • e.g.: temperature in comment => temperature as float
• Try, when possible, to use shorter keys, and to group items into sub-dict of common interest, like:

'mesh' : {
    'type': 'hexagonal',
    'positions': {                                                    
        'Si': {                                                             
            'u': np.r_[0.465],                                              
            'x': None,                                                      
            'y': None,                                                      
            'z': None,                                                      
            'N': None,                                                      
        },                                                                  
        'O': {                                                              
            'u': np.r_[0.415, 0.272, 0.120],                                
            'x': None,                                                      
            'y': None,                                                      
            'z': None,                                                      
            'N': None,                                                      
        },                                                                  
    },                                                          
    'source': 'R.W.G. Wyckoff, Crystal Structures (1963)',
}  

and

'Inter-atomic': {
    'distances': {
        'a0': 4.91304,                                                      
        'c0': 5.40463,                                                      
       },        
    'units': 'A',
    'comments': 'within the unit cell',
    ' Tref': {                         # 'refetence temperature', but shorter
        'data': 25.,
        'units': 'C', 
    },

and

'Thermal expansion': {
    'coefs': {                                            
        'alpha_a': 1.337e-5,                                                
        'alpha_c': 7.97e-6,                                                 
    },
    'units': '1/C',                          
    'comments':   'in parallel directions to a0 and c0',               
    'source': 'R.W.G. Wyckoff, Crystal Structures',
}

and

'atomic scattering': { 
    'factor': {                                          
        'Si': np.r_[                                                        
            12., 11., 9.5, 8.8, 8.3, 7.7, 7.27, 6.25, 5.3,                  
            4.45, 3.75, 3.15, 2.7, 2.35, 2.07, 1.87, 1.71, 1.6,             
        ],                                                                  
        'O': np.r_[                                                         
            9., 7.836, 5.756, 4.068, 2.968, 2.313, 1.934, 1.710, 1.566,     
            1.462, 1.373, 1.294,                                            
        ],                                                                  
    },                                                                      
    'source': 'Intern. Tables for X-Ray Crystallography, Vol.I,II,III,IV (1985)',
}

=> it will allow to re-use some key words like 'source', 'units', ... which makes them more generic

[Didou09]

@cjperks7 , once @adriendaros will have implemented those last minor changes, you're good to go

[Didou09]

@adriendaros, @cjperks7 , just to let you know:

• Minor changes: I'm changing the names of some keys of the _DCRYST dict in _rockingcurve_def.py to:
  • make sure all keys are valid python variable identifier (can be checked using str.isidentifier(): no special characters except _, in particular no - or /, (orspaces` in the names + no upper case for first letter + no starting with a number). This makes the library more robust in case we need to use these keys as input arguments of functions later.
  • try to harmonize and remove ambiguities (general rule: all keys that can be group in a theme shall start with a common scheme)
    • atoms, atomic number and number of atoms become atoms, atoms_Z and atoms_nb

I am also adding a basic unit test to make sure the computation always works. It just runs it for akk known crystals at the same wavelength. It simply makes sure it is not broken and there is no bug that would stop the computation, does not check the value of the result and does not check plots.

[Didou09]

@adriendaros , I have a question:

The key sin(theta)/lambda (now called sin_theta_lambda, see above for reasons), seems to never be used in the routine compute_rockingcurve(). => why is it never used ? => is it useful in the computation ? => if it is not useful, should we remove it ?

[adriendaros]

@adriendaros , I have a question:

The key sin(theta)/lambda (now called sin_theta_lambda, see above for reasons), seems to never be used in the routine compute_rockingcurve(). => why is it never used ? => is it useful in the computation ? => if it is not useful, should we remove it ?

@Didou09 don't do that you fool, it's very useful to get the values of the atomic scattering factor, real and imaginary parts, of each kind of atoms within the crystal.

Into tofu/spectro/_rocking_def.py, there is data corresponding to the variation of the atomic scattering[factors] with respect to this variable sin_theta_lambda for each atom.

In order to interpolate the value of the scattering factor in our crystals, the routine CrystBragg_comp_lattice_spacing() called by compute_rockingcurve() is computing the appropriated value of sin_theta_lambda for the crystal used, depending on the Miller indices and the wavelength targetted.

Then, you'll find into the compute_rockingcurve() routine line 214 part Calculation of the structure factor that this constant is called in the computation of each real and imaginary scattering factor for each atom.

So no, we can't do without it because of its great importance in the determination of the crystal structure factor.

[Didou09]

OK, good reasons, keeping it (I did not delete it, was waiting for your answer)

[cjperks7]

Documentation as of the last commit

Changes implemented --- (i) One can now calculate the rocking curve for Quartz with flexible choice on Miller indice by either:

• (1) Select a pre-defined crystal name as appropriate for WEST within ['Quartz_110', 'Quartz_102']

# Computes rocking curve
dout1 = rc.compute_rockingcurve(
    crystal = 'Quartz_102', 
    lamb = 3.96, 
    plot_power_ratio = False
    )

• (2) Or one can define a suitable crystal material dictionary, i.e. as appropriate for the CMOD He-like Ar Johann spectrometer

# C-Mod's crystal
dcry2 = {
  'material': 'Quartz',
  'name': 'HIREX',
  'symbol': 'Qz102',
  'miller': np.r_[1., 0., 2.,],
  'target': {
      'ion': 'Ar16+',
      'lamb': 3.96, # e-10
      'units': 'm',
      },
  'd_hkl': 2.281124272650091,
  }


# Computes rocking curve
dout2 = rc.compute_rockingcurve(
  crystal = dcry2['name'], 
  din=dcry2, 
  lamb = 3.96, 
  plot_power_ratio = False
  )

• Note that the two methods shown above give the same rocking curve since the target wavelength in the pre-defined crystals get over-written
• Benchmarking against the x0h code for two crystals of potential interest as been conducted to good agreement between the two models
• 

(ii) It is now simple to loop over instances of tofu.spectro._rockingcurve.compute_rockingcurve() for Quartz with a target wavelength and arbitrary Miller indices to scope Bragg angle v. integrated reflectivity. All one needs to do is alter the array in dcry2['miller'] as given above

• 
• Note that this looping will be inserted into the tofu_sparc workflow under Issue015

(iii) These changes have been propagated up into tofu.data._class5_check() so that one can initiate a crystal using coll.add_crystal() choosing the dmat input to either be a string or dictionary as discussed in point (i). For example, building CMOD's Johann He-like Ar spectrometer:

dcry ={
    # Defines crystal geometry
    'dgeom': {
        'cent': np.array([1,1,1]), # [m], center in global axis
        'nin': np.array([1,0,0]), # inward normal vec
        'e0': np.array([0,1,0]), # horizontal vec
        'e1': np.array([0,0,1]), # vertical vec
        'curve_r': np.array([1385e-3, 1385e-3]), # [m], crystal radius of curvature
        'extenthalf': np.array([64e-3/2, 27e-3/2]), # [m], disp. of corners wrt cent in (hor. ver.)
    },
    # Defines crystal material
    'dmat': {
        'material': 'Quartz',
        'name': 'HIREX',
        'symbol': 'Qz102',
        'miller': np.r_[1., 0., 2.],
        'target':{
            'ion': 'Ar16+',
            'lamb': 3.96, # [1e-10 m], photon wavelength
            'units': 'm',
            }, 
        'd_hkl': 2.281124272650091, # [m], crystal lattice spacing
        },
    }

# Adds diffraction crystal object to diagnostic
coll.add_crystal(
    key='cry',
    dgeom=dcry['dgeom'],
    dmat=dcry['dmat'],
    #dmat='Quartz_102'
    )

• It has been benchmarked that the expected rocking curve is stored in the Collection object
• 

[cjperks7]

@Didou09 @adriendaros Germanium rocking curves work now 🎉

[cjperks7]

Will do a PR in ToFu and tofu_sparc tomorrow