Merge pull request #26 from DanPorter/new_cell_basis

New cell basis

Author: DanPorter

Dans_Diffraction/Structures/README.md

@@ -36,4 +36,5 @@ The files are downloaded or generated from a range of data sources, for example:
 * [MAGNDATA](http://webbdcrista1.ehu.es/magndata/)
 * [ISODISTORT](https://stokes.byu.edu/iso/isodistort.php)
 * [Jana2006](http://www-xray.fzu.cz/jana/jana.html)
-* [Vesta](http://jp-minerals.org/vesta/en/)
+* [Vesta](http://jp-minerals.org/vesta/en/)
+* [CrystalMaker](https://crystalmaker.com/index.html)

Dans_Diffraction/Structures/Rutile.cif

@@ -0,0 +1,40 @@
+data_Rutile
+_audit_creation_method         'generated by CrystalDiffract 6 for OS X'
+_cell_length_a                   4.5870
+_cell_length_b                   4.5870
+_cell_length_c                   2.9540
+_cell_angle_alpha               90.0000
+_cell_angle_beta                90.0000
+_cell_angle_gamma               90.0000
+
+_symmetry_space_group_name_H-M     'P 42/m n m'
+
+loop_
+_symmetry_equiv_pos_as_xyz
+'+x,+y,+z'
+'-x,-y,+z'
+'1/2-y,1/2+x,1/2+z'
+'1/2+y,1/2-x,1/2+z'
+'1/2-x,1/2+y,1/2-z'
+'1/2+x,1/2-y,1/2-z'
+'+y,+x,-z'
+'-y,-x,+z'
+'-x,-y,-z'
+'+x,+y,-z'
+'1/2+y,1/2-x,1/2-z'
+'1/2-y,1/2+x,1/2-z'
+'1/2+x,1/2-y,1/2+z'
+'1/2-x,1/2+y,1/2+z'
+'-y,-x,-z'
+'+y,+x,+z'
+
+loop_
+_atom_site_label
+_atom_site_type_symbol
+_atom_site_occupancy
+_atom_site_fract_x
+_atom_site_fract_y
+_atom_site_fract_z
+ Label  Ti 1.0000  0.0000  0.0000  0.0000
+ Label  O 1.0000  0.3051  0.3051  0.0000
+

Dans_Diffraction/Structures/triclinic.cif

@@ -0,0 +1,31 @@
+
+#======================================================================
+# CRYSTAL DATA
+#----------------------------------------------------------------------
+data_VESTA_phase_1
+
+_chemical_name_common                  'Test_Triclinic'
+_cell_length_a                         6.867000
+_cell_length_b                         16.705999
+_cell_length_c                         5.399000
+_cell_angle_alpha                      62.980000
+_cell_angle_beta                       42.250000
+_cell_angle_gamma                      52.160000
+_cell_volume                           328.878918
+_space_group_name_H-M_alt              'P 1'
+_space_group_IT_number                 1
+
+loop_
+_space_group_symop_operation_xyz
+   'x, y, z'
+
+loop_
+   _atom_site_label
+   _atom_site_occupancy
+   _atom_site_fract_x
+   _atom_site_fract_y
+   _atom_site_fract_z
+   _atom_site_adp_type
+   _atom_site_U_iso_or_equiv
+   _atom_site_type_symbol
+   Fe1        1.0     0.000000     0.000000     0.000000    Uiso  0.050000 Fe

Dans_Diffraction/__init__.py

@@ -31,8 +31,8 @@
 Diamond
 2017
 
-Version 3.2.4
-Last updated: 25/09/24
+Version 3.3.2
+Last updated: 20/11/2024
 
 Version History:
 02/03/18 1.0    Version History started.
@@ -81,6 +81,8 @@
 04/09/24 3.2.3  Updated method of generating charge states in classes_orbtials.py, Thanks Seonghun!
 25/09/24 3.2.4  Fixed error of missing rotation matrices after load_spacegroup. Thanks asteppke!
 26/09/24 3.3.0  Added complex neutron scattering lengths for isotopes from package periodictable. Thanks thamnos!
+06/11/24 3.3.1  Fixed incorrect cell basis for triclinic cells. Added functions_lattice.py and tests. Thanks LeeRichter!
+20/11/24 3.3.2  Added alternate option for neutron scattering lengths
 
 Acknoledgements:
     2018        Thanks to Hepesu for help with Python3 support and ideas about breaking up calculations
@@ -105,6 +107,8 @@
     May 2024    Thanks to paul-cares pointing out a silly spelling error in the title!
     Aug 2024    Thanks to Seonghun for pointing out the error with charge states in Hf
     Aug 2024    Thanks to MaxPelly for spell checks in examples
+    Sep 2024    Thanks to thamnos for suggestion to add complex neutron scattering lengths
+    Oct 2024    Thanks to Lee Richter for pointing out the error in triclinic basis definition
 
 -----------------------------------------------------------------------------
    Copyright 2024 Diamond Light Source Ltd.
@@ -146,6 +150,7 @@
 # Dans Diffraction
 from . import functions_general as fg
 from . import functions_plotting as fp
+from . import functions_lattice as fl
 from . import functions_crystallography as fc
 from .classes_crystal import Crystal
 from .classes_multicrystal import MultiCrystal
@@ -154,8 +159,12 @@
 from .functions_crystallography import readcif
 
 
-__version__ = '3.2.4'
-__date__ = '25/09/24'
+__all__ = ['fg', 'fp', 'fl', 'fc', 'Crystal', 'MultiCrystal', 'readcif',
+           'Structures', 'Fdmnes', 'FdmnesAnalysis']
+
+
+__version__ = '3.3.2'
+__date__ = '2024/11/20'
 
 
 # Build

Dans_Diffraction/classes_crystal.py

@@ -57,8 +57,9 @@
 from warnings import warn
 
 # Internal functions
-from . import functions_general as fg, functions_crystallography as fc
+from . import functions_general as fg
 from . import functions_crystallography as fc
+from . import functions_lattice as fl
 from .classes_orientation import Orientation
 from .classes_properties import Properties
 # from .classes_orbitals import CrystalOrbitals
@@ -514,9 +515,10 @@ def __init__(self, a=1.0, b=1.0, c=1.0, alpha=90., beta=90.0, gamma=90.0):
         self.alpha = float(alpha)
         self.beta = float(beta)
         self.gamma = float(gamma)
+        self._basis_function = fl.basis_3
         self.orientation = Orientation()
 
-    def latt(self, lattice_parameters=(), *args, **kwargs):
+    def latt(self, *lattice_parameters, **kwargs):
         """ 
         Generate lattice parameters with list
           latt(1) -> a=b=c=1,alpha=beta=gamma=90
@@ -527,7 +529,7 @@ def latt(self, lattice_parameters=(), *args, **kwargs):
           latt(a=1,b=2,c=3,alpha=10,beta=20,gamma=30]) -> a=1,b=2,c=3,alpha=10,beta=20,gamma=30
         """
 
-        lp = fc.gen_lattice_parameters(lattice_parameters, *args, **kwargs)
+        lp = fc.gen_lattice_parameters(*lattice_parameters, **kwargs)
         self.a = lp[0]
         self.b = lp[1]
         self.c = lp[2]
@@ -585,35 +587,46 @@ def lp(self):
         """
         return self.a, self.b, self.c, self.alpha, self.beta, self.gamma
 
+    def choose_basis(self, option='default'):
+        """
+        Choose the basis function
+        Options:
+            1. c || z, b* || y - basis choice of Materials Project
+            2. a || x, c* || z - basis choice of Vesta
+            3. c || z, a* || x - basis choice of Busing & Levy (Default)
+        :param option: name or number of basis
+        """
+        self._basis_function = fl.choose_basis(option)
+
     def UV(self):
         """
          Returns the unit cell as a [3x3] array, [A,B,C]
          The vector A is directed along the x-axis
         """
-        return self.orientation(fc.latpar2uv_rot(*self.lp()))
+        return self.orientation(self._basis_function(*self.lp()))
 
     def UVstar(self):
         """
         Returns the reciprocal unit cell as a [3x3] array, [A*,B*,C*]
         :return: [a*;b*;c*]
         """
-        return fc.RcSp(self.UV())
+        return fl.reciprocal_basis(self.UV())
 
     def volume(self):
         """
         Returns the volume of the unit cell, in A^3
         :return: volume
         """
-        return fc.calc_vol(self.UV())
+        return fl.lattice_volume(*self.lp())
 
     def Bmatrix(self):
         """
         Calculate the Busing and Levy B matrix from a real space UV
          "choose the x-axis parallel to a*, the y-axis in the plane of a* and b*, and the z-axis perpendicular to
          that plane"
-         W. R. Busing & H. A. Levy, Acta  Cryst.  (1967). 22,  457
+        W. R. Busing & H. A. Levy, Acta  Cryst.  (1967). 22,  457
         """
-        return fc.Bmatrix(self.UV())
+        return fl.busingandlevy(*self.lp())
 
     def info(self):
         """
@@ -660,7 +673,7 @@ def indexQ(self, Q):
             > HKL = [1,0,0]
         """
         Q = np.reshape(np.asarray(Q, dtype=float), [-1, 3])
-        return fc.indx(Q, self.UVstar())
+        return fl.index_lattice(Q, self.UVstar())
 
     def calculateR(self, UVW):
         """
@@ -685,7 +698,7 @@ def indexR(self, R):
             > UVW = [0.1,0,0]
         """
         R = np.reshape(np.asarray(R, dtype=float), [-1, 3])
-        return fc.indx(R, self.UV())
+        return fl.index_lattice(R, self.UV())
 
     def moment(self, mxmymz):
         """Calcualte moment from value stored in cif"""
@@ -892,7 +905,7 @@ def find_close_reflections(self, hkl, energy_kev, max_twotheta=2, max_angle=10):
         selected = (tth_dif < max_twotheta) * (all_angles < max_angle)
         return all_hkl[selected, :]
 
-    def reciprocal_space_plane(self, x_axis=[1, 0, 0], y_axis=[0, 1, 0], centre=[0, 0, 0], q_max=4.0, cut_width=0.05):
+    def reciprocal_space_plane(self, x_axis=(1, 0, 0), y_axis=(0, 1, 0), centre=(0, 0, 0), q_max=4.0, cut_width=0.05):
         """
         Returns positions within a reciprocal space plane
           x_axis = direction along x, in units of the reciprocal lattice (hkl)
@@ -1481,7 +1494,7 @@ def check(self):
     def fitincell(self):
         """Adjust all atom positions to fit within unit cell"""
         uvw = fc.fitincell(self.uvw())
-
+        self.u, self.v, self.w = uvw.T
 
     def uvw(self):
         """
@@ -2366,7 +2379,7 @@ def __init__(self, Parent, P):
         self.P = np.asarray(P, dtype=float)
         self.Parent = Parent
         newUV = Parent.Cell.calculateR(P)
-        self.new_cell(fc.UV2latpar(newUV))
+        self.new_cell(fl.basis2latpar(newUV))
         self.scale = Parent.scale * np.prod(self.P)
 
         # Add exta functions

Dans_Diffraction/classes_orientation.py

@@ -15,11 +15,9 @@
 @author: DGPorter
 """
 
-import sys, os, re
 import numpy as np
-from warnings import warn
 
-from . import functions_general as fg
+from . import functions_lattice as fl
 from . import functions_crystallography as fc
 
 __version__ = '1.0.0'
@@ -174,10 +172,9 @@ class CrystalOrientation:
       x-axis : vector normal to phi axis where phi=0 (toward ceiling (+y) in lab frame)
       y-axis : vector normal to x,z axes (parallel to beam (+z) in lab frame)
     """
-    def __init__(self, lattice_parameters=(), *args, **kwargs):
+    def __init__(self, *lattice_parameters, **kwargs):
         self.a, self.b, self.c, self.alpha, self.beta, self.gamma = fc.gen_lattice_parameters(
-            lattice_parameters,
-            *args,
+            *lattice_parameters,
             **kwargs
         )
         self.orientation = Orientation()
@@ -218,17 +215,15 @@ def lp(self):
 
     def lp_star(self):
         """Return tuple of reciprocal lattice parameters in inverse-angstroms and degrees"""
-        uv = self._uv()
-        return fc.latpar_reciprocal(uv)
+        return fl.reciprocal_lattice_parameters(*self.lp())
 
     def _uv(self):
-        """Return unit vectors [a,b,c] in default frame (b along X)"""
-        return fc.latpar2uv_rot(self.a, self.b, self.c, self.alpha, self.beta, self.gamma)
+        """Return unit vectors [a,b,c] in default frame (a* along X)"""
+        return fl.basis_3(*self.lp())
 
     def _uvstar(self):
         """Return unit vectors in recirpocal space in default frame (a* along X)"""
-        uv = self._uv()
-        return fc.RcSp(uv)
+        return fl.reciprocal_basis(self._uv())
 
     def unit_vectors(self):
         """Return real space unit vectors [a, b, c]"""
@@ -240,7 +235,7 @@ def reciprocal_unit_vectors(self):
 
     def bmatrix(self):
         """Return the B matrix from Busing & Levy in the diffractometer frame"""
-        return fc.Bmatrix(self._uv())
+        return fl.busingandlevy(*self.lp())
 
     def ubmatrix(self):
         """Return UB matrix from Busing & Levy in the diffractometer frame"""

Dans_Diffraction/classes_scattering.py

@@ -104,6 +104,9 @@ class Scattering:
     # Uses the coefficients for analytical approximation to the scattering factors from:
     #        "Waasmaier and Kirfel, Acta Cryst. (1995) A51, 416-431"
     _use_waaskirf_scattering_factor = False
+
+    # Use the neutron scattering factors from the Internation tables (Sears 1995)
+    _use_sears_scattering_lengths = False
 
     # Thermal Factors
     _use_isotropic_thermal_factor = True
@@ -161,6 +164,8 @@ def __str__(self):
         out += '      max twotheta: %s\n' % self._scattering_max_twotheta
         out += ' ---X-Ray Settings---\n'
         out += ' Waasmaier scattering factor: %s\n' % self._use_waaskirf_scattering_factor
+        out += ' ---Neutron Settings---\n'
+        out += ' Sears (ITC) scattering lengths: %s\n' % self._use_sears_scattering_lengths
         out += ' ---Magnetic Settings---\n'
         out += '   Mag. scattering: %s\n' % self._calclate_magnetic_component
         out += '  Mag. form factor: %s\n' % self._use_magnetic_form_factor
@@ -185,7 +190,7 @@ def setup_scatter(self, scattering_type=None, energy_kev=None, wavelength_a=None
                       powder_units=None, powder_pixels=None, powder_lorentz=None, powder_overlap=None,
                       specular=None, parallel=None, theta_offset=None,
                       min_theta=None, max_theta=None, min_twotheta=None, max_twotheta=None,
-                      output=True, scattering_factors=None, magnetic_formfactor=None,
+                      output=True, scattering_factors=None, scattering_lengths=None, magnetic_formfactor=None,
                       polarisation=None, polarisation_vector=None, azimuthal_reference=None, azimuth=None, flm=None):
         """
         Simple way to set scattering parameters, each parameter is internal to xtl (self)
@@ -205,6 +210,7 @@ def setup_scatter(self, scattering_type=None, energy_kev=None, wavelength_a=None
         specular    : self._scattering_specular_direction : [h,k,l] : reflections normal to sample surface
         parallel    : self._scattering_parallel_direction : [h,k,l] : reflections normal to sample surface
         scattering_factors: self._use_waaskirf_scattering_factor : xray scattering factor ['waaskirf', 'itc']
+        scattering_lengths: self._use_sears_scattering_lengths : neutron scattering lengths ['sears', 'default']
         magnetic_formfactor: self._use_magnetic_form_factor: True/False magnetic form factor for magnetic SF
         polarisation: self._polarisation                  : beam polarisation setting ['ss', 'sp'*, 'sp', 'pp']
         polarisation_vector: _polarisation_vector_incident: [x,y,z] incident polarisation vector
@@ -263,6 +269,14 @@ def setup_scatter(self, scattering_type=None, energy_kev=None, wavelength_a=None
                 print('Using scattering factors from: International Tables of Crystallography Vol. C, Table 6.1.1.4')
                 self._use_waaskirf_scattering_factor = False
 
+        if scattering_lengths is not None:
+            if scattering_lengths.lower() in ['sears', 'itc', 'alternate', 'alt']:
+                print('Using scattering lengths from International Tables of Crystallography Vol. C, Table 4.4.4.1')
+                self._use_sears_scattering_lengths = True
+            else:
+                print('Using scattering lengths from Neutron Data Booklet')
+                self._use_sears_scattering_lengths = False
+
         if magnetic_formfactor is not None:
             self._use_magnetic_form_factor = magnetic_formfactor
 
@@ -470,8 +484,11 @@ def structure_factor(self, hkl=None, scattering_type=None, int_hkl=True, **kwarg
                     qmag = fg.mag(_q)      # q magnitude
                     # Scattering factors
                     if scattering_type in fs.SCATTERING_TYPES['neutron']:
-                        # ff = fc.atom_properties(atom_type, 'Coh_b')
-                        ff = fc.neutron_scattering_length(atom_type)
+                        if self._use_sears_scattering_lengths:
+                            ff = fc.neutron_scattering_length(atom_type, 'Sears')
+                        else:
+                            # ff = fc.atom_properties(atom_type, 'Coh_b')
+                            ff = fc.neutron_scattering_length(atom_type)
                     elif scattering_type in fs.SCATTERING_TYPES['electron']:
                         ff = fc.electron_scattering_factor(atom_type, qmag)
                     elif scattering_type in fs.SCATTERING_TYPES['xray fast']: