Merge pull request #31 from DanPorter/supercell_fixes

Supercell scale fix

Author: DanPorter

Dans_Diffraction/__init__.py

@@ -31,8 +31,8 @@
 Diamond
 2017-2025
 
-Version 3.3.3
-Last updated: 06/02/2025
+Version 3.3.4
+Last updated: 12/04/2025
 
 Version History:
 02/03/18 1.0    Version History started.
@@ -84,6 +84,7 @@
 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
 06/02/25 3.3.3  Added scattering options for polarised neutron and x-ray scattering. Thanks dragonyanglong!
+12/04/25 3.3.4  Improved superstructure calculations by fixing scale parameter
 
 Acknoledgements:
     2018        Thanks to Hepesu for help with Python3 support and ideas about breaking up calculations
@@ -113,7 +114,7 @@
     Dec 2024    Thanks to dragonyanglong for pointing out the error with magnetic neutron scattering
 
 -----------------------------------------------------------------------------
-   Copyright 2024 Diamond Light Source Ltd.
+   Copyright 2018-2025 Diamond Light Source Ltd.
 
    Licensed under the Apache License, Version 2.0 (the "License");
    you may not use this file except in compliance with the License.

Dans_Diffraction/classes_crystal.py

@@ -24,8 +24,8 @@
 Diamond
 2017
 
-Version 3.3.0
-Last updated: 06/02/25
+Version 3.3.1
+Last updated: 06/04/25
 
 Version History:
 27/07/17 1.0    Version History started.
@@ -49,7 +49,8 @@
 15/11/21 3.2.2  Added Cell.orientation, updated Cell.UV()
 12/01/21 3.2.3  Added Symmetry.axial_vector
 22/05/23 3.2.4  Added Symmetry.wyckoff_label(), Symmetry.spacegroup_dict
-06/05/25 3.3.0  Symmetry.from_cif now loads operations from find_spacegroup if not already loaded
+06/05/24 3.3.0  Symmetry.from_cif now loads operations from find_spacegroup if not already loaded
+06/04/25 3.3.1  scale parameter of superlattice improved
 
 @author: DGPorter
 """
@@ -68,7 +69,7 @@
 from .classes_multicrystal import MultiCrystal
 from .classes_plotting import Plotting, PlottingSuperstructure
 
-__version__ = '3.3.0'
+__version__ = '3.3.1'
 
 
 class Crystal:
@@ -2386,7 +2387,8 @@ def __init__(self, Parent, P):
         self.Parent = Parent
         newUV = Parent.Cell.calculateR(P)
         self.new_cell(fl.basis2latpar(newUV))
-        self.scale = Parent.scale * np.prod(self.P)
+        parent_cells_in_supercell = fc.calc_vol(P)
+        self.scale = Parent.scale * parent_cells_in_supercell
 
         # Add exta functions
         self.Plot = PlottingSuperstructure(self)

Dans_Diffraction/classes_plotting.py

@@ -1500,12 +1500,12 @@ def parent_generate_intensity_cut(self, x_axis=(1, 0, 0), y_axis=(0, 1, 0), cent
         c_cart = self.xtl.Parent.Cell.calculateQ(centre)
 
         # Generate lattice of reciprocal space points
-        maxq = np.sqrt(q_max**2 + q_max**2 + np.sum(c_cart**2)) # generate all reflections
-        print('Max Q distance: %4.2f A-1'%maxq)
+        maxq = np.sqrt(q_max**2 + q_max**2 + np.sum(c_cart**2))  # generate all reflections
+        print('Max Q distance: %4.2f A-1' % maxq)
         hmax, kmax, lmax = fc.maxHKL(maxq, self.xtl.Cell.UVstar())
         HKL = fc.genHKL([hmax, -hmax], [kmax, -kmax], [lmax, -lmax])
-        HKL = HKL #+ centre  # reflection about central reflection
-        print('Number of reflections in sphere: %1.0f'%len(HKL))
+        HKL = HKL  # + centre  # reflection about central reflection
+        print('Number of reflections in sphere: %1.0f' % len(HKL))
 
         # Determine the directions in cartesian space
         x_cart = self.xtl.calculateQ_parent(x_axis)
@@ -1535,10 +1535,12 @@ def parent_generate_intensity_cut(self, x_axis=(1, 0, 0), y_axis=(0, 1, 0), cent
             pHKL = self.xtl.superhkl2parent(HKLinbox)
             pHKL, inten = self.xtl.Parent.Symmetry.symmetric_intensity(pHKL, inten)
             HKL = self.xtl.parenthkl2super(pHKL)
+            print('Adding parent symmetry domains, adding %d reflections' % (len(HKL) - len(pHKL)))
         else:
             HKL = HKLinbox
         q = self.xtl.calculateQ_parent(HKL)
 
+        # remove reflections not in plot
         box_coord = fg.index_coordinates(q - c_cart, CELL)
         incell = np.all(np.abs(box_coord) <= 0.5, axis=1)
         plane_coord = 2 * q_max * box_coord[incell, :]
@@ -1661,21 +1663,26 @@ def simulate_intensity_cut(self, x_axis=(1, 0, 0), y_axis=(0, 1, 0), centre=(0,
         mesh_vec_b = fg.index_coordinates(Q_vec_b, CELL) * 2 * q_max
 
         # Vector arrows and lattice point labels
-        cen_lab = '(%1.3g,%1.3g,%1.3g)' % (centre[0], centre[1], centre[2])
-        vec_a_lab = '(%1.3g,%1.3g,%1.3g)' % (vec_a[0] + centre[0], vec_a[1] + centre[1], vec_a[2] + centre[2])
-        vec_b_lab = '(%1.3g,%1.3g,%1.3g)' % (vec_b[0] + centre[0], vec_b[1] + centre[1], vec_b[2] + centre[2])
+        supx, supy, supz = 0.0 + np.around(self.xtl.parenthkl2super(centre)[0], 3)
+        svax, svay, svaz = 0.0 + np.around(self.xtl.parenthkl2super(vec_a + centre)[0], 3)
+        svbx, svby, svbz = 0.0 + np.around(self.xtl.parenthkl2super(vec_b + centre)[0], 3)
+        cen_lab = '(%1.3g,%1.3g,%1.3g)$_{p}$' % (centre[0], centre[1], centre[2])
+        # cen_lab += '\n(%1.3g,%1.3g,%1.3g)$_{s}$' % (supx, supy, supz)
+        vec_a_lab = '(%1.3g,%1.3g,%1.3g)$_{p}$' % (vec_a[0] + centre[0], vec_a[1] + centre[1], vec_a[2] + centre[2])
+        vec_a_lab += '\n(%1.3g,%1.3g,%1.3g)$_{s}$' % (svax, svay, svaz)
+        vec_b_lab = '(%1.3g,%1.3g,%1.3g)$_{p}$' % (vec_b[0] + centre[0], vec_b[1] + centre[1], vec_b[2] + centre[2])
+        vec_b_lab += '\n(%1.3g,%1.3g,%1.3g)$_{s}$' % (svbx, svby, svbz)
 
         lattQ = fp.axis_lattice_points(mesh_vec_a, mesh_vec_b, plt.axis())
         fp.plot_lattice_lines(lattQ, mesh_vec_a, mesh_vec_b, lw=0.5, c='grey')
-        fp.plot_vector_arrows(mesh_vec_a, mesh_vec_b, vec_a_lab, vec_b_lab)
+        fp.plot_vector_arrows(mesh_vec_a, mesh_vec_b, vec_a_lab, vec_b_lab, color='k')
         plt.text(0 - (0.2 * q_max), 0 - (0.1 * q_max), cen_lab, fontname=fp.DEFAULT_FONT, weight='bold', size=18)
 
         # Plot labels
         xlab = r'Q || (%1.3g,%1.3g,%1.3g) [$\AA^{-1}$]' % (x_axis[0], x_axis[1], x_axis[2])
         ylab = r'Q || (%1.3g,%1.3g,%1.3g) [$\AA^{-1}$]' % (y_axis[0], y_axis[1], y_axis[2])
-        supercentre = self.xtl.parenthkl2super(centre)[0]
         ttl = '%s\n(%1.3g,%1.3g,%1.3g)$_{p}$ = (%1.3g,%1.3g,%1.3g)$_{s}$' \
-              % (self.xtl.name, centre[0], centre[1], centre[2], supercentre[0], supercentre[1], supercentre[2])
+              % (self.xtl.name, centre[0], centre[1], centre[2], supx, supy, supz)
         fp.labels(ttl, xlab, ylab)
 
 

Dans_Diffraction/functions_crystallography.py

@@ -72,6 +72,7 @@
 PENGFILE = os.path.join(datadir, 'peng.dat')
 NSLFILE = os.path.join(datadir, 'neutron_isotope_scattering_lengths.dat')
 NSLFILE_SEARS = os.path.join(datadir, 'neutron_isotope_scattering_lengths_sears.dat')
+ASFFILE = os.path.join(datadir, 'atomic_scattering_factors.npy')
 
 # List of Elements in order sorted by length of name
 ELEMENT_LIST = [
@@ -1236,8 +1237,7 @@ def atomic_scattering_factor(element, energy_kev=None):
     :param energy_kev: float or list energy in keV (None to return original, uninterpolated list)
     :return: f1, f2, shape dependent on shapes of element and energy_kev:  float, or [ene] or [ele, ene]
     """
-    asf_file = os.path.join(datadir, 'atomic_scattering_factors.npy')
-    asf = np.load(asf_file, allow_pickle=True)
+    asf = np.load(ASFFILE, allow_pickle=True)
     asf = asf.item()
 
     element = np.asarray(element, dtype=str).reshape(-1)
@@ -1295,8 +1295,7 @@ def xray_dispersion_corrections(elements, energy_kev=None):
     :param energy_kev: float or list energy in keV (None to return original, uninterpolated list)
     :return: f', f" with shape (len(energy), len(elements))
     """
-    asf_file = os.path.join(datadir, 'atomic_scattering_factors.npy')
-    asf = np.load(asf_file, allow_pickle=True)
+    asf = np.load(ASFFILE, allow_pickle=True)
     asf = asf.item()
 
     energy_kev = np.asarray(energy_kev, dtype=float).reshape(-1)
@@ -3446,7 +3445,7 @@ def group_intensities(q_values, intensity, min_overlap=0.01):
 def calc_vol(UV):
     """Calculate volume in Angstrom^3 from unit vectors"""
     a, b, c = UV
-    return np.dot(a, np.cross(b, c))
+    return np.abs(np.dot(a, np.cross(b, c)))
 
 
 def cif2table(cif):

Examples/example_supercell.py

@@ -3,11 +3,11 @@
 Build a supercell from multiple unit cells and simulate the diffraction pattern
 """
 
-import sys,os
+import sys, os
 import numpy as np
-import matplotlib.pyplot as plt # Plotting
+import matplotlib.pyplot as plt  # Plotting
 cf = os.path.dirname(__file__)
-sys.path.insert(0,os.path.join(cf,'..'))
+sys.path.insert(0, os.path.join(cf, '..'))
 import Dans_Diffraction as dif
 
 
@@ -17,14 +17,14 @@
 
 #P = [[2,-1,0],[1,3,0],[0,0,1]] # 1/7th Supercell
 #P = [[-1,3,0],[4,3,0],[0,0,1]] # Square Supercell
-P = [[3,0,0],[4,5,0],[0,0,1]] # Stripe Supercell
+P = [[3, 0, 0], [4, 5, 0], [0, 0, 1]]  # Stripe Supercell
 #P = [[1,3,0],[3,-1,0],[0,0,1]]
 #P = [[2,0,0],[0,2,0],[0,0,1]] # Double supercell
 #sup = xtl.generate_superstructure(P)
 
 # Set discrete occupancies for average structure
-xtl.Atoms.occupancy[2]=0
-xtl.Atoms.occupancy[3]=1
+xtl.Atoms.occupancy[2] = 0
+xtl.Atoms.occupancy[3] = 1
 xtl.generate_structure()
 
 # Generate the superstructure, repeating the parent/average structure to fill the supercell
@@ -34,8 +34,8 @@
 # Stripe Cell
 # Na1 6,7 16,17 26,27 36,37 46,47 56,57 66,67 76,77 86,87 96,97 106,107 116,117 126,127 136,137 146,147
 # Na2 8,9 18,19 28,29 38,39 48,49 58,59 68,69 78,79 88,89 98,99 108,109 118,119 128,129 138,139 148,149
-sup.Structure.occupancy[[6 ,16,26,77,87,107]] = 1 # Na1
-sup.Structure.occupancy[[8 ,18,38,28,48,58, 139, 119, 149, 109, 89,79]] = 0 # Na2
+sup.Structure.occupancy[[6, 16, 26, 77, 87, 107]] = 1  # Na1
+sup.Structure.occupancy[[8, 18, 38, 28, 48, 58, 139, 119, 149, 109, 89, 79]] = 0  # Na2
 
 # Plot the Na layers showing the ordering
 sup.Plot.plot_layers(layers=[0.25, 0.75], layer_width=0.01, show_labels=True)