From 433a0d59f8da58c66192acd1b0bd0dea2bf9f5c1 Mon Sep 17 00:00:00 2001 From: Kevin Lewis <55998034+kevindlewis23@users.noreply.github.com> Date: Thu, 15 Aug 2024 10:36:39 -0400 Subject: [PATCH] Unitcell testcases (#702) * Unitcell testcases --- hexrd/material/unitcell.py | 177 ++++++++++++++++---------------- tests/unitcell/test_vec_math.py | 103 +++++++++++++++++++ 2 files changed, 194 insertions(+), 86 deletions(-) create mode 100644 tests/unitcell/test_vec_math.py diff --git a/hexrd/material/unitcell.py b/hexrd/material/unitcell.py index 9f3b3f24..26f3a2fe 100644 --- a/hexrd/material/unitcell.py +++ b/hexrd/material/unitcell.py @@ -42,7 +42,7 @@ def _calcstar(v, sym, mat): if dist < 1E-3: isnew = False break - if(isnew): + if isnew: vp = np.atleast_2d(vp) vsym = np.vstack((vsym, vp)) @@ -158,7 +158,7 @@ def calcmatrices(self): [a*c*cb, b*c*ca, c**2]]) self._vol = np.sqrt(np.linalg.det(self.dmt)) - if(self.vol < 1e-5): + if self.vol < 1e-5: warnings.warn('unitcell volume is suspiciously small') ''' @@ -198,32 +198,34 @@ def calcmatrices(self): choices are 'd' (direct), 'r' (reciprocal) and 'c' (cartesian)''' def TransSpace(self, v_in, inspace, outspace): - if(inspace == 'd'): - if(outspace == 'r'): + if inspace == outspace: + return v_in + if inspace == 'd': + if outspace == 'r': v_out = np.dot(v_in, self.dmt) - elif(outspace == 'c'): + elif outspace == 'c': v_out = np.dot(self.dsm, v_in) else: raise ValueError( - 'inspace in ''d'' but outspace can''t be identified') + 'inspace in "d" but outspace can\'t be identified') - elif(inspace == 'r'): - if(outspace == 'd'): + elif inspace == 'r': + if outspace == 'd': v_out = np.dot(v_in, self.rmt) - elif(outspace == 'c'): + elif outspace == 'c': v_out = np.dot(self.rsm, v_in) else: raise ValueError( - 'inspace in ''r'' but outspace can''t be identified') + 'inspace in "r" but outspace can\'t be identified') - elif(inspace == 'c'): - if(outspace == 'r'): - v_out = np.dot(v_in, self.rsm) - elif(outspace == 'd'): + elif inspace == 'c': + if outspace == 'r': v_out = np.dot(v_in, self.dsm) + elif outspace == 'd': + v_out = np.dot(v_in, self.rsm) else: raise ValueError( - 'inspace in ''c'' but outspace can''t be identified') + 'inspace in "c" but outspace can\'t be identified') else: raise ValueError('incorrect inspace argument') @@ -234,11 +236,11 @@ def TransSpace(self, v_in, inspace, outspace): def CalcDot(self, u, v, space): - if(space == 'd'): + if space == 'd': dot = np.dot(u, np.dot(self.dmt, v)) - elif(space == 'r'): + elif space == 'r': dot = np.dot(u, np.dot(self.rmt, v)) - elif(space == 'c'): + elif space == 'c': dot = np.dot(u, v) else: raise ValueError('space is unidentified') @@ -247,13 +249,13 @@ def CalcDot(self, u, v, space): def CalcLength(self, u, space): - if(space == 'd'): + if space == 'd': mat = self.dmt # vlen = np.sqrt(np.dot(u, np.dot(self.dmt, u))) - elif(space == 'r'): + elif space == 'r': mat = self.rmt # vlen = np.sqrt(np.dot(u, np.dot(self.rmt, u))) - elif(space == 'c'): + elif space == 'c': mat = np.eye(3) # vlen = np.linalg.norm(u) else: @@ -297,7 +299,7 @@ def CalcAngle(self, u, v, space): def CalcCross(self, p, q, inspace, outspace, vol_divide=False): iv = 0 - if(vol_divide): + if vol_divide: vol = self.vol else: vol = 1.0 @@ -306,50 +308,50 @@ def CalcCross(self, p, q, inspace, outspace, vol_divide=False): p[2]*q[0]-p[0]*q[2], p[0]*q[1]-p[1]*q[0]]) - if(inspace == 'd'): + if inspace == 'd': ''' cross product vector is in reciprocal space and can be converted to direct or cartesian space ''' pxq *= vol - if(outspace == 'r'): + if outspace == 'r': pass - elif(outspace == 'd'): + elif outspace == 'd': pxq = self.TransSpace(pxq, 'r', 'd') - elif(outspace == 'c'): + elif outspace == 'c': pxq = self.TransSpace(pxq, 'r', 'c') else: raise ValueError( 'inspace is ''d'' but outspace is unidentified') - elif(inspace == 'r'): + elif inspace == 'r': ''' cross product vector is in direct space and can be converted to any other space ''' pxq /= vol - if(outspace == 'r'): + if outspace == 'r': pxq = self.TransSpace(pxq, 'd', 'r') - elif(outspace == 'd'): + elif outspace == 'd': pass - elif(outspace == 'c'): + elif outspace == 'c': pxq = self.TransSpace(pxq, 'd', 'c') else: raise ValueError( 'inspace is ''r'' but outspace is unidentified') - elif(inspace == 'c'): + elif inspace == 'c': ''' cross product is already in cartesian space so no volume factor is involved. can be converted to any other space too ''' - if(outspace == 'r'): + if outspace == 'r': pxq = self.TransSpace(pxq, 'c', 'r') - elif(outspace == 'd'): + elif outspace == 'd': pxq = self.TransSpace(pxq, 'c', 'd') - elif(outspace == 'c'): + elif outspace == 'c': pass else: raise ValueError( @@ -398,7 +400,7 @@ def GenerateCartesianPGSym(self): self.SYM_PG_c = np.array(self.SYM_PG_c) self.SYM_PG_c[np.abs(self.SYM_PG_c) < eps] = 0. - if(self._pointGroup == self._laueGroup): + if self._pointGroup == self._laueGroup: self.SYM_PG_c_laue = self.SYM_PG_c else: for sop in self.SYM_PG_d_laue: @@ -422,8 +424,7 @@ def GenerateCartesianPGSym(self): supergroup_laue = self._supergroup_laue sym_supergroup_laue = symmetry.GeneratePGSYM(supergroup_laue) - if((self.latticeType == 'monoclinic' or - self.latticeType == 'triclinic')): + if self.latticeType in ('monoclinic', 'triclinic'): ''' for monoclinic groups c2 and c2h, the supergroups are orthorhombic, so no need to convert from direct to @@ -462,7 +463,7 @@ def GenerateCartesianPGSym(self): not be rotated as they have the c* axis already aligned with the z-axis SS 12/10/2020 ''' - if(self.latticeType == 'monoclinic'): + if self.latticeType == 'monoclinic': om = np.array([[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]]) @@ -480,7 +481,7 @@ def GenerateCartesianPGSym(self): triclinic group c1!! SS 12/10/2020 ''' - if(self._pointGroup == 'c1'): + if self._pointGroup == 'c1': om = np.array([[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]]) for i, s in enumerate(self.SYM_PG_supergroup): @@ -536,7 +537,7 @@ def CalcOrbit(self, v, reduceToUC=True): break # if its new add this to the list - if(isnew): + if isnew: asym_pos = np.vstack((asym_pos, rr)) n += 1 @@ -549,21 +550,21 @@ def CalcStar(self, v, space, applyLaue=False): this function calculates the symmetrically equivalent hkls (or uvws) for the reciprocal (or direct) point group symmetry. ''' - if(space == 'd'): + if space == 'd': mat = self.dmt.astype(np.float64) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_d_laue.astype(np.float64) else: sym = self.SYM_PG_d.astype(np.float64) - elif(space == 'r'): + elif space == 'r': mat = self.rmt.astype(np.float64) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_r_laue.astype(np.float64) else: sym = self.SYM_PG_r.astype(np.float64) - elif(space == 'c'): + elif space == 'c': mat = np.eye(3) - if(applyLaue): + if applyLaue: sym = self.SYM_PG_c_laue.astype(np.float64) else: sym = self.SYM_PG_c.astype(np.float64) @@ -750,7 +751,11 @@ def InitializeInterpTable(self): f_anomalous_data = [] self.pe_cs = {} - data = importlib.resources.open_binary(hexrd.resources, 'Anomalous.h5') + data = ( + importlib.resources.files(hexrd.resources) + .joinpath('Anomalous.h5') + .open('rb') + ) with h5py.File(data, 'r') as fid: for i in range(0, self.atom_ntype): @@ -917,7 +922,7 @@ def ChooseSymmetric(self, hkllist, InversionSymmetry=True): laue = InversionSymmetry for i, g in enumerate(hkllist): - if(mask[i]): + if mask[i]: geqv = self.CalcStar(g, 'r', applyLaue=laue) @@ -995,11 +1000,11 @@ def getHKLs(self, dmin): for g in hkl_allowed: # ignore [0 0 0] as it is the direct beam - if(np.sum(np.abs(g)) != 0): + if np.sum(np.abs(g)) != 0: dspace = 1./self.CalcLength(g, 'r') - if(dspace >= dmin): + if dspace >= dmin: hkl_dsp.append(g) ''' @@ -1037,7 +1042,7 @@ def Required_C(self, C): return np.array([C[x] for x in _StiffnessDict[self._laueGroup][0]]) def MakeStiffnessMatrix(self, inp_Cvals): - if(len(inp_Cvals) != len(_StiffnessDict[self._laueGroup][0])): + if len(inp_Cvals) != len(_StiffnessDict[self._laueGroup][0]): x = len(_StiffnessDict[self._laueGroup][0]) msg = (f"number of constants entered is not correct." f" need a total of {x} independent constants.") @@ -1079,16 +1084,16 @@ def inside_spheretriangle(self, conn, dir3, hemisphere, switch): first get vertices of the triangles in the ''' vertex = self.sphere_sector.vertices[switch] - # if(switch == 'pg'): + # if switch == 'pg': # vertex = self.sphere_sector.vertices - # elif(switch == 'laue'): + # elif switch == 'laue': # vertex = self.sphere_sector.vertices_laue - # elif(switch == 'super'): + # elif switch == 'super': # vertex = self.sphere_sector.vertices_supergroup - # elif(switch == 'superlaue'): + # elif switch == 'superlaue': # vertex = self.sphere_sector.vertices_supergroup_laue A = np.atleast_2d(vertex[:, conn[0]]).T @@ -1107,29 +1112,29 @@ def inside_spheretriangle(self, conn, dir3, hemisphere, switch): determinant can be very small positive or negative number ''' - if(np.abs(d1) < eps): + if np.abs(d1) < eps: d1 = 0. - if(np.abs(d2) < eps): + if np.abs(d2) < eps: d2 = 0. - if(np.abs(d3) < eps): + if np.abs(d3) < eps: d3 = 0. ss = np.unique(np.sign([d1, d2, d3])) - if(hemisphere == 'upper'): - if(np.all(ss >= 0.)): + if hemisphere == 'upper': + if np.all(ss >= 0.): mask.append(True) else: mask.append(False) - elif(hemisphere == 'both'): - if(len(ss) == 1): + elif hemisphere == 'both': + if len(ss) == 1: mask.append(True) - elif(len(ss) == 2): - if(0 in ss): + elif len(ss) == 2: + if 0 in ss: mask.append(True) else: mask.append(False) - elif(len(ss) == 3): + elif len(ss) == 3: mask.append(False) mask = np.array(mask) @@ -1161,7 +1166,7 @@ def reduce_dirvector(self, dir3, switch='pg'): ''' idx = np.arange(dir3.shape[0], dtype=np.int32) dir3 = np.ascontiguousarray(np.atleast_2d(dir3)) - if(dir3.ndim != 2): + if dir3.ndim != 2: raise RuntimeError("reduce_dirvector: invalid shape of dir3 array") ''' @@ -1173,10 +1178,10 @@ def reduce_dirvector(self, dir3, switch='pg'): ''' eps = constants.sqrt_epsf - if(np.all(np.abs(np.linalg.norm(dir3, axis=1) - 1.0) < eps)): + if np.all(np.abs(np.linalg.norm(dir3, axis=1) - 1.0) < eps): dir3n = dir3 else: - if(np.all(np.linalg.norm(dir3) > eps)): + if np.all(np.linalg.norm(dir3) > eps): dir3n = dir3/np.tile(np.linalg.norm(dir3, axis=1), [3, 1]).T else: raise RuntimeError( @@ -1200,30 +1205,30 @@ def reduce_dirvector(self, dir3, switch='pg'): ntriangle = self.sphere_sector.ntriangle[switch] connectivity = self.sphere_sector.connectivity[switch] - if(switch == 'pg'): + if switch == 'pg': sym = self.SYM_PG_c - elif(switch == 'super'): + elif switch == 'super': sym = self.SYM_PG_supergroup - elif(switch == 'laue'): + elif switch == 'laue': sym = self.SYM_PG_c_laue - elif(switch == 'superlaue'): + elif switch == 'superlaue': sym = self.SYM_PG_supergroup_laue for sop in sym: - if(dir3_copy.size != 0): + if dir3_copy.size != 0: dir3_sym = np.dot(sop, dir3_copy.T).T mask = np.zeros(dir3_sym.shape[0]).astype(bool) - if(ntriangle == 0): - if(hemisphere == 'both'): + if ntriangle == 0: + if hemisphere == 'both': mask = np.ones(dir3_sym.shape[0], dtype=bool) - elif(hemisphere == 'upper'): + elif hemisphere == 'upper': mask = dir3_sym[:, 2] >= 0. else: for ii in range(ntriangle): @@ -1232,8 +1237,8 @@ def reduce_dirvector(self, dir3, switch='pg'): hemisphere, switch) mask = np.logical_or(mask, tmpmask) - if(np.sum(mask) > 0): - if(dir3_reduced.size != 0): + if np.sum(mask) > 0: + if dir3_reduced.size != 0: dir3_reduced = np.vstack( (dir3_reduced, dir3_sym[mask, :])) idx_red = np.hstack((idx_red, idx[mask])) @@ -1269,7 +1274,7 @@ class which correctly color the orientations for this crystal class. the replaceing barycenter to pi-theta) ''' - if(laueswitch == True): + if laueswitch: ''' this is the case where we color orientations based on the laue group of the crystal. this is always going to be the case with x-ray which @@ -1280,7 +1285,7 @@ class which correctly color the orientations for this crystal class. the dir3, switch='superlaue') switch = 'superlaue' - elif(laueswitch == False): + else: ''' follow the logic in the function description ''' @@ -1317,7 +1322,7 @@ def color_orientations(self, ''' first make sure that the rotation matric is size nx3x3 ''' - if(rmats.ndim == 2): + if rmats.ndim == 2: rmats = np.atleast_3d(rmats).T else: assert rmats.ndim == 3, "rotations matrices need to \ @@ -1358,7 +1363,7 @@ def convert_lp_to_valunits(self, lp): """ lp_valunit = [] for i in range(6): - if(i < 3): + if i < 3: lp_valunit.append( valWUnit('lp', 'length', lp[i], 'nm')) @@ -1415,7 +1420,7 @@ def lparms(self, lp): self.calcmatrices() self.init_max_g_index() self.CalcMaxGIndex() - if(hasattr(self, 'numat')): + if hasattr(self, 'numat'): self.CalcDensity() @property @@ -1541,7 +1546,7 @@ def U(self): def U(self, Uarr): self._U = Uarr self.aniU = False - if(Uarr.ndim > 1): + if Uarr.ndim > 1: self.aniU = True self.calcBetaij() @@ -1569,9 +1574,9 @@ def sgnum(self): @sgnum.setter def sgnum(self, val): - if(not(isinstance(val, int))): + if not(isinstance(val, int)): raise ValueError('space group should be integer') - if(not((val >= 1) and (val <= 230))): + if not((val >= 1) and (val <= 230)): raise ValueError('space group number should be between 1 and 230.') self._sym_sgnum = val diff --git a/tests/unitcell/test_vec_math.py b/tests/unitcell/test_vec_math.py new file mode 100644 index 00000000..794ad4da --- /dev/null +++ b/tests/unitcell/test_vec_math.py @@ -0,0 +1,103 @@ +from pytest import fixture +from hexrd.material import Material, unitcell +import numpy as np + + +@fixture +def cell() -> unitcell.unitcell: + return Material().unitcell + + +def get_space(): + return np.random.choice(['d', 'r', 'c']) + + +def test_calc_dot(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + v2c = cell.TransSpace(v2, space, 'c') + assert np.allclose(np.dot(v1c, v2c), cell.CalcDot(v1, v2, space)) + + +def test_calc_length(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + assert np.allclose(np.linalg.norm(v1c), cell.CalcLength(v1, space)) + + +def test_calc_angle(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + space = get_space() + v1c = cell.TransSpace(v1, space, 'c') + v2c = cell.TransSpace(v2, space, 'c') + norms = np.linalg.norm(v1c) * np.linalg.norm(v2c) + assert np.allclose( + np.arccos(np.dot(v1c, v2c) / norms), cell.CalcAngle(v1, v2, space) + ) + + +def test_calc_cross(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + v2 = np.random.rand(3) * 10 - 5 + inspace = get_space() + outspace = get_space() + v1c = cell.TransSpace(v1, inspace, 'c') + v2c = cell.TransSpace(v2, inspace, 'c') + expected = cell.TransSpace(np.cross(v1c, v2c), 'c', outspace) + assert np.allclose( + expected, cell.CalcCross(v1, v2, inspace, outspace, True) + ) + + +def test_norm_vec(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = get_space() + norm_v1 = cell.NormVec(v1, space) + assert np.allclose(1, cell.CalcLength(norm_v1, space)) + # Make sure we don't change the direction + assert np.allclose( + v1 / np.linalg.norm(v1), norm_v1 / np.linalg.norm(norm_v1) + ) + + +def test_trans_space(cell: unitcell.unitcell): + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + inspace = get_space() + outspace = get_space() + v1_out = cell.TransSpace(v1, inspace, outspace) + assert np.allclose(v1, cell.TransSpace(v1_out, outspace, inspace)) + + +def test_calc_star(cell: unitcell.unitcell): + """ + Just ensuring that the outspace doesn't matter + """ + np.random.seed(0) + for _ in range(100): + v1 = np.random.rand(3) * 10 - 5 + space = np.random.choice(['d', 'r']) + v1c = cell.TransSpace(v1, space, 'c') + assert np.allclose( + cell.CalcStar(v1, space, False), + cell.TransSpace(cell.CalcStar(v1c, 'c', False), 'c', space), + ) + assert np.allclose( + cell.CalcStar(v1, space, True), + cell.TransSpace(cell.CalcStar(v1c, 'c', True), 'c', space), + )