MentDB WeakMQL Recipes
← innov-ai.com All profiles Portal
Profiles / Sciences
CM

Crystallographer & materials physicist

Miller indices, metric tensors, elastic tensors and structure factors.

10 recipes on this page · open in the interactive view

Recipes
  1. Miller indices, d-spacing and Bragg angle in a monoclinic cell
  2. Elastic tensor of silicon: stability and moduli
  3. Structure factors and powder pattern of NaCl
  4. Symmetry operations: rotation angle, order and the crystallographic restriction
  5. Indexing the powder pattern of a cubic metal
  6. Debye-Waller factor of copper
  7. Crystallite size and microstrain: Williamson-Hall plot
  8. Twinning in titanium: the 10-12 twin
  9. X-ray crystal density and the Avogadro constant
  10. Lattice parameters from an orientation matrix by QR
executed on MentDB

Miller indices, d-spacing and Bragg angle in a monoclinic cell

A plane cuts the axes at 2a, 3b, 6c. Get its Miller indices, build the metric tensor of a monoclinic cell, the cell volume, the reciprocal metric tensor and the d-spacing, then the Bragg angle 2 theta for Cu K-alpha1.

Pipeline
01 · MillerReciprocals times the lcmarith lcm→02 · metricG and V = sqrt(det G)matrix det→03 · reciprocalG* = G^-1matrix inverse→04 · d-spacing1/d^2 = h G* hvector dot→05 · Bragg2 theta for 1.5406 Amath calc
Script
admin
#1. Plane cutting the axes at 2a, 3b, 6c: reciprocals times the lcm give the Miller indices;
-> "[m]" (arith lcm "[2, 3, 6]");
-> "[h]" (arith frac_mul (arith frac_div 1 2) [m]);
-> "[k]" (arith frac_mul (arith frac_div 1 3) [m]);
-> "[l]" (arith frac_mul (arith frac_div 1 6) [m]);
-> "[hkl]" (concat "[" [h] ", " [k] ", " [l] "]");
#2. Metric tensor of a monoclinic cell a 5.20, b 9.00, c 7.40 A, beta 104.5 deg;
-> "[G]" (math calc "5.20*7.40*cos(104.5*pi/180)");
-> "[G]" (concat "[[27.04, 0, " [G] "], [0, 81, 0], [" [G] ", 0, 54.76]]");
-> "[V]" (math calc (concat "sqrt(" (matrix det [G]) ")"));
#3. Reciprocal metric tensor and d-spacing: 1/d^2 = h G* h;
-> "[Gs]" (matrix inverse [G]);
-> "[d]" (math calc (concat "1/sqrt(" (vector dot [hkl] (matrix mul [Gs] [hkl])) ")"));
#4. Bragg angle for Cu K-alpha1 1.5406 A;
-> "[tt]" (math calc (concat "2*asin(1.5406/(2*" [d] "))*180/pi"));
concat "Miller (" [h] " " [k] " " [l] "), gcd " (arith gcd (concat "[" [h] ", " [k] ", " [l] "]")) "; V = " (math round [V] 3) " A3; d = " (math round [d] 4) " A = " (math round (unit convert [d] "angstrom" "nm") 5) " nm; 2theta = " (math round [tt] 3) " deg"
Result
mentdb
Miller (3 2 1), gcd 1; V = 335.289 A3; d = 1.4655 A = 0.14655 nm; 2theta = 63.421 deg
Go further
  • Compute the angle between two plane normals with the reciprocal metric tensor
  • Index a whole powder pattern with a for loop over hkl
executed on MentDB

Elastic tensor of silicon: stability and moduli

From the cubic stiffness constants of silicon, check the Born mechanical stability (Cholesky factor, positive eigenvalues), invert to the compliance tensor and compute the Young moduli along 100 and 111, the Poisson ratio, the bulk modulus and the Zener anisotropy.

Pipeline
01 · stabilityCholesky factormatrix cholesky→02 · eigenvaluesC11+2C12, C11-C12, C44matrix eigen→03 · complianceS = C^-1matrix inverse→04 · moduliE100, E111, numatrix get
Script
admin
#1. Stiffness of silicon in Voigt notation, GPa: C11 165.7, C12 63.9, C44 79.6;
-> "[C]" "[[165.7, 63.9, 63.9, 0, 0, 0], [63.9, 165.7, 63.9, 0, 0, 0], [63.9, 63.9, 165.7, 0, 0, 0], [0, 0, 0, 79.6, 0, 0], [0, 0, 0, 0, 79.6, 0], [0, 0, 0, 0, 0, 79.6]]";
#2. Born stability: C positive definite, the Cholesky factor exists and all eigenvalues are > 0;
-> "[L]" (json load select "rx_chol" (matrix cholesky [C]) /L);
-> "[ev]" (json load select "rx_eig" (matrix eigen [C]) /values);
#3. Compliance S = C^-1: Young modulus and Poisson ratio along 100;
-> "[S]" (matrix inverse [C]);
-> "[E100]" (math calc (concat "1/" (matrix get [S] 1 1)));
-> "[nu]" (math calc (concat "-" (matrix get [S] 1 2) "/" (matrix get [S] 1 1)));
#4. Young modulus along 111: 1/E = S11 - 2 S11 - S12 - S44/2 x 1/3;
-> "[E111]" (math calc (concat "1/(" (matrix get [S] 1 1) "-2*(" (matrix get [S] 1 1) "-" (matrix get [S] 1 2) "-" (matrix get [S] 4 4) "/2)/3)"));
#5. Bulk modulus and Zener anisotropy;
-> "[B]" (math calc "(165.7+2*63.9)/3");
-> "[Az]" (math calc "2*79.6/(165.7-63.9)");
concat "eigenvalues " (string replace (string replace [ev] [_n_] "") " " "") " = C11+2C12, C11-C12, C44; L11 = " (math round (matrix get [L] 1 1) 4) "; E100 = " (math round [E100] 1) " GPa, E111 = " (math round [E111] 1) " GPa, nu = " (math round [nu] 4) ", B = " (math round [B] 1) " GPa, Zener A = " (math round [Az] 3) ", rank " (matrix rank [C]) ", trace " (matrix trace [C]) ", spectral norm " (matrix norm [C]) ", cond " (math round (matrix cond [C]) 3)
Result
mentdb
eigenvalues [293.5,101.8,101.8,79.6,79.6,79.6] = C11+2C12, C11-C12, C44; L11 = 12.8725; E100 = 130.1 GPa, E111 = 187.9 GPa, nu = 0.2783, B = 97.8 GPa, Zener A = 1.564, rank 6, trace 735.9, spectral norm 293.5, cond 3.687
Go further
  • Try an unstable set C12 > C11: matrix cholesky then fails
  • Average the moduli of a polycrystal (Voigt and Reuss) from C and S
executed on MentDB

Structure factors and powder pattern of NaCl

Sum the structure factor F = sum f exp(2 pi i h.x) over the 8 ions of the NaCl cell for the first reflections, then plot the powder pattern with multiplicity, Lorentz-polarisation and a simple fall-off: the odd reflections are weak (|F| = 4 (fCl - fNa)).

Pipeline
01 · phase2 pi (hx + ky + lz)math calc→02 · sumf exp(i phase)complex exp→03 · intensityF times its conjugatecomplex conj→04 · 2 thetaBragg anglemath calc→05 · patternBar chartplot chart
Script
admin
#1. NaCl, a = 5.640 A: Na+ 10 electrons on the fcc sites, Cl- 18 electrons shifted by a/2;
-> "[atoms]" "10,0,0,0;10,0.5,0.5,0;10,0.5,0,0.5;10,0,0.5,0.5;18,0.5,0,0;18,0,0.5,0;18,0,0,0.5;18,0.5,0.5,0.5";
-> "[refl]" "1,1,1,8;2,0,0,6;2,2,0,12;3,1,1,24;2,2,2,8;4,0,0,6;3,3,1,24;4,2,0,24";
-> "[tt]" ""; -> "[I]" ""; -> "[txt]" "";
for (-> "[r]" 1) (<= [r] (atom size [refl] ";")) (++ "[r]") {
	-> "[hkl]" (atom get [refl] [r] ";");
	-> "[h]" (atom get [hkl] 1 ","); -> "[k]" (atom get [hkl] 2 ","); -> "[l]" (atom get [hkl] 3 ",");
	#2. Structure factor F = sum f exp 2 pi i h.x;
	-> "[F]" "0";
	for (-> "[j]" 1) (<= [j] (atom size [atoms] ";")) (++ "[j]") {
		-> "[at]" (atom get [atoms] [j] ";");
		-> "[ph]" (math calc (concat "2*pi*(" [h] "*" (atom get [at] 2 ",") "+" [k] "*" (atom get [at] 3 ",") "+" [l] "*" (atom get [at] 4 ",") ")"));
		-> "[F]" (complex add [F] (complex mul (atom get [at] 1 ",") (complex exp (concat [ph] "i"))));
	};
	#3. Intensity = multiplicity x Lorentz-polarisation x F F* x fall-off B = 3 A2;
	-> "[F2]" (complex re (complex mul [F] (complex conj [F])));
	-> "[th]" (math calc (concat "asin(1.5406*sqrt(" [h] "^2+" [k] "^2+" [l] "^2)/(2*5.640))"));
	-> "[lp]" (math calc (concat "(1+cos(2*" [th] ")^2)/(sin(" [th] ")^2*cos(" [th] "))"));
	-> "[tt]" (concat [tt] "," (math round (math calc (concat "2*" [th] "*180/pi")) 2));
	-> "[I]" (concat [I] "," (math round (math calc (concat (atom get [hkl] 4 ",") "*" [lp] "*" [F2] "*exp(-6*(sin(" [th] ")/1.5406)^2)")) 0));
	-> "[txt]" (concat [txt] " " [h] [k] [l] " " (math round (complex abs [F]) 0) "/" (math round (complex arg [F]) 2));
};
#4. Bar chart of the normalised intensities;
-> "[Iv]" (string replace (concat "[" [I] "]") "[," "[");
-> "[In]" (vector scale [Iv] (math calc (concat "100/" (stat max [Iv]))));
plot chart (concat "{\"title\": \"NaCl, Cu Ka, |F|/phase:" [txt] "\", \"x_label\": \"2 theta (deg)\", \"y_label\": \"I (%)\", \"series\": [{\"name\": \"NaCl\", \"type\": \"bars\", \"x\": " (string replace (concat "[" [tt] "]") "[," "[") ", \"y\": " [In] "}]}")
Result
mentdb
The chart is shown in the Plot tab of the editor.
Go further
  • Replace Na by K (KCl): fK+ = fCl- and the odd reflections vanish
  • Use form factors depending on sin(theta)/lambda with num interpolate in a table
executed on MentDB

Symmetry operations: rotation angle, order and the crystallographic restriction

Three symmetry operations of a hexagonal crystal are given as integer matrices in the a, b, c basis. Their determinant tells proper from improper, the trace gives the rotation angle (1 + 2cos theta), and the smallest power equal to the identity gives the order; a 5-fold axis is impossible.

Pipeline
01 · proper?Determinant +1 or -1matrix det→02 · angleTrace = 1 + 2 cos thetamatrix trace→03 · orderSmallest n with W^n = Imatrix pow→04 · restriction2cos(2pi/5) not an integermath calc
Script
admin
#1. Symmetry operations given as matrices in the hexagonal basis a, b, c;
-> "[R3]" "[[0, -1, 0], [1, -1, 0], [0, 0, 1]]";
-> "[R6]" "[[1, -1, 0], [1, 0, 0], [0, 0, 1]]";
-> "[M]" "[[0, -1, 0], [-1, 0, 0], [0, 0, 1]]";
#2. Determinant +1 proper rotation, -1 improper, trace = 1 + 2 cos theta for a rotation;
-> "[out]" "";
-> "[ops]" "R3 R6 M";
for (-> "[i]" 0) (< [i] 3) (++ "[i]") {
	-> "[name]" (string substring [ops] (* [i] 3) (+ (* [i] 3) 2));
	-> "[W]" [R3]; if (== [i] 1) {-> "[W]" [R6];}; if (== [i] 2) {-> "[W]" [M];};
	-> "[d]" (matrix det [W]);
	-> "[th]" (math round (math rad_to_deg (math acos (math calc (concat "(" [d] "*" (matrix trace [W]) "-1)/2")))) 6);
	#Order n: smallest power giving the identity, a finite-order matrix with trace 3;
	-> "[n]" 1;
	while (!= (matrix trace (matrix pow [W] [n])) 3) {++ "[n]";};
	-> "[out]" (concat [out] [name] ": det " [d] ", angle " [th] " deg, order " [n] "; ");
};
#3. Crystallographic restriction: 2 cos theta must be an integer, only n = 1, 2, 3, 4, 6;
concat [out] "for n = 5, 2cos(2pi/5) = " (math round (math calc "2*cos(2*pi/5)") 4) " is not an integer: no 5-fold axis in a lattice"
Result
mentdb
R3: det 1, angle 120 deg, order 3; R6: det 1, angle 60 deg, order 6; M: det -1, angle 180 deg, order 2; for n = 5, 2cos(2pi/5) = 0.618 is not an integer: no 5-fold axis in a lattice
Go further
executed on MentDB

Indexing the powder pattern of a cubic metal

Five peaks of a copper powder pattern measured with Cu K-alpha1: the ratios of sin^2 theta give the integers h^2+k^2+l^2 (3, 4, 8, 11, 12, so FCC), each peak gives a lattice parameter, and the cell with 4 atoms gives the X-ray density.

Pipeline
01 · ratios3.sin^2 theta / sin^2 theta1math calc→02 · latticea = lambda.sqrt(N)/2 sin theta per peakjson select→03 · meanMean a and its spreadstat pstddev→04 · densityrho = Z.M / NA.a^3const get
Script
admin
#1. Powder pattern of a cubic metal with Cu K-alpha1, 2theta of the first five peaks in degrees;
-> "[lambda]" 1.5406;
-> "[tt]" "[43.30, 50.43, 74.13, 89.93, 95.14]";
json load "rx_tu1_pk" [tt];
#2. sin^2 theta of each peak divided by the first one, times 3: FCC gives 3 4 8 11 12, BCC 2 4 6 8 10;
-> "[s1]" (math calc (concat "sin(" (json select "rx_tu1_pk" "/[0]") "*pi/360)^2"));
-> "[N]" "["; -> "[a]" "[";
for (-> "[i]" 0) (< [i] 5) (++ "[i]") {
	-> "[th]" (math calc (concat (json select "rx_tu1_pk" (concat "/[" [i] "]")) "*pi/360"));
	-> "[hkl2]" (math round (math calc "3*sin([th])^2/[s1]") 0);
	-> "[N]" (concat [N] [hkl2] ",");
	#Lattice parameter a = lambda.sqrt of h2+k2+l2 over 2 sin theta;
	-> "[a]" (concat [a] (math calc "[lambda]*sqrt([hkl2])/(2*sin([th]))") ",");
};
-> "[a]" (concat (string substring [a] 0 (- (string length [a]) 1)) "]");
#3. Mean lattice parameter, then X-ray density with Z = 4 atoms of Cu per cell;
-> "[am]" (stat mean [a]);
-> "[rho]" (math calc (concat "4*63.546/(" (const get "NA") "*([am]*1e-8)^3)"));
-> "[N]" (concat (string substring [N] 0 (- (string length [N]) 1)) "]");
concat "h2+k2+l2 = " [N] ": 111 200 220 311 222, h k l unmixed, FCC; a = " (math round [am] 4) " A, spread " (math round_sig (stat pstddev [a]) 2) " A; density = " (math round [rho] 3) " g/cm3"
Result
mentdb
h2+k2+l2 = [3,4,8,11,12]: 111 200 220 311 222, h k l unmixed, FCC; a = 3.6156 A, spread 0.00062 A; density = 8.93 g/cm3
Go further
  • Correct the zero shift by fitting a against cos^2 theta / sin theta with stat fit (Nelson-Riley)
  • Get the crystallite size from the peak widths with the Scherrer equation
executed on MentDB

Debye-Waller factor of copper

From the Debye temperature of copper, the script computes the Debye function, the Debye-Waller factor B with the Planck, Boltzmann and atomic mass constants, the mean square displacement and the intensity kept by the 111, 200 and 400 reflections.

Pipeline
01 · constantsh, kB, uconstant physics_planck→02 · boltzmannBoltzmann constantconstant physics_k→03 · massAtomic mass unitconstant physics_amu→04 · debyeDebye functionnum integrate→05 · intensityexp(-2 B s^2)math calc
Script
admin
#1. Copper at 293 K, Debye temperature 343 K, M = 63.546 u: constants h, kB and u;
-> "[h]" (string replace (atom get (constant physics_planck) 1 " ") "(" "");
-> "[kB]" (string replace (atom get (constant physics_k) 1 " ") "(" "");
-> "[m]" (math calc (concat "63.546*" (string replace (atom get (constant physics_amu) 1 " ") "(" "")));
-> "[T]" 293; -> "[TD]" 343; -> "[x]" (math calc "[TD]/[T]");
#2. Debye function phi x = 1/x integral of t / e^t - 1 from 0 to x;
-> "[phi]" (math calc (concat (num integrate "t/(exp(t)-1)" "t" 0 [x]) "/[x]"));
#3. Debye-Waller factor B = 6 h^2 T / m kB TD^2 x phi x / x + 1/4, in angstrom^2;
-> "[B]" (math calc "6*[h]^2*[T]/([m]*[kB]*[TD]^2)*([phi]/[x]+0.25)*1e20");
#4. Intensity lost by thermal vibration exp -2B s^2, s = 1/2d, a = 3.615 A;
-> "[out]" "";
for (-> "[i]" 1) (<= [i] 3) (++ "[i]") {
	-> "[hkl]" (atom get "111,200,400" [i] ",");
	-> "[d]" (math calc (concat "3.615/sqrt(" (atom get "3,4,16" [i] ",") ")"));
	-> "[out]" (concat [out] " " [hkl] ": " (math round (math calc (concat "exp(-2*[B]/(4*" [d] "^2))*100")) 1) " %");
};
concat "phi(x) = " (math round [phi] 4) ", B = " (math round [B] 3) " A2, mean square displacement u2 = B/8 pi^2 = " (math round (math calc "[B]/(8*pi^2)") 5) " A2; intensity kept:" [out]
Result
mentdb
phi(x) = 0.7449, B = 0.399 A2, mean square displacement u2 = B/8 pi^2 = 0.00505 A2; intensity kept: 111: 95.5 % 200: 94.1 % 400: 78.3 %
Go further
  • Go back from measured intensities at two temperatures to the Debye temperature with num solve
  • Add the zero-point term at 0 K
executed on MentDB

Crystallite size and microstrain: Williamson-Hall plot

Six peaks of a CeO2 nanopowder are corrected for the instrument width. The script compares the Scherrer size of the first peak with a Williamson-Hall least-squares line, which separates size and microstrain.

Pipeline
01 · correctsqrt(b^2 - bi^2)math calc→02 · cosbeta cos theta = beta / sec thetamath sec→03 · fitLeast squares linematrix solve→04 · strain4 eps tan thetamath tan
Script
admin
#1. CeO2 nanopowder, Cu Ka 1.5406 A: 2 theta and FWHM in degrees, instrument FWHM 0.07 deg;
-> "[tt]" "28.55,33.08,47.48,56.33,69.40,76.70";
-> "[fw]" "0.452,0.478,0.556,0.613,0.715,0.786";
-> "[X]" ""; -> "[Y]" "";
for (-> "[i]" 1) (<= [i] 6) (++ "[i]") {
	-> "[th]" (math calc (concat (atom get [tt] [i] ",") "/2*pi/180"));
	#2. Sample broadening in radians, Gaussian correction b^2 - bi^2;
	-> "[b]" (math calc (concat "sqrt(" (atom get [fw] [i] ",") "^2-0.07^2)*pi/180"));
	#3. Williamson-Hall: b cos theta = K lambda / D + 4 eps sin theta, with cos = 1/sec;
	-> "[Y]" (concat [Y] "," (math calc (concat [b] "/" (math sec [th]))));
	-> "[X]" (concat [X] ",[1," (math calc (concat "4*" (math sin [th]))) "]");
	if (== [i] 1) { -> "[D1]" (math calc (concat "0.9*1.5406/(" [b] "/" (math sec [th]) ")/10")); -> "[tan1]" (math tan [th]); -> "[b1]" [b]; };
};
#4. Least squares line: intercept K lambda / D, slope = microstrain;
-> "[fit]" (matrix solve (concat "[" (string substr [X] 1) "]") (concat "[" (string substr [Y] 1) "]"));
-> "[D]" (math calc (concat "0.9*1.5406/" (matrix get (concat "[" [fit] "]") 1 1) "/10"));
-> "[eps]" (matrix get (concat "[" [fit] "]") 1 2);
concat "Scherrer on 111 alone: D = " (math round [D1] 1) " nm; Williamson-Hall: D = " (math round [D] 1) " nm, strain = " (math round (math calc "[eps]*100") 3) " %, strain part of the 111 width b_strain = 4 eps tan theta = " (math round (math calc "4*[eps]*[tan1]/[b1]*100") 1) " %"
Result
mentdb
Scherrer on 111 alone: D = 18.4 nm; Williamson-Hall: D = 25.2 nm, strain = 0.208 %, strain part of the 111 width b_strain = 4 eps tan theta = 27.1 %
Go further
  • Use Lorentzian profiles: subtract the widths linearly
  • Draw the Williamson-Hall plot with plot fit
executed on MentDB

Twinning in titanium: the 10-12 twin

In hexagonal alpha-Ti, the script builds the Cartesian basis, the normal of the (10-12) twin plane from the reciprocal basis, the angle with the c axis using atan2, the 85 degree c-axis rotation, the twinning shear and the d-spacing.

Pipeline
01 · reciprocalRows of B^-1matrix inverse→02 · normalh a* + k b* + l c*vector add→03 · angleatan2(|n x c|, n.c)math atan2→04 · shears = (g^2 - 3)/sqrt 3 gmath signum→05 · checkd from hypotmath hypot
Script
admin
#1. alpha-Ti, hexagonal a = 2.951 A, c = 4.684 A: Cartesian basis a1, a2, c;
-> "[a]" 2.951; -> "[c]" 4.684; -> "[g]" (math calc "[c]/[a]");
-> "[Bm]" (concat "[[" [a] ", " (math calc "-[a]/2") ", 0], [0, " (math calc "[a]*sqrt(3)/2") ", 0], [0, 0, " [c] "]]");
#2. Normal of the twin plane 10-12 = 1 0 2 in 3 indices: h a* + k b* + l c*, rows of B^-1;
-> "[Bi]" (matrix inverse [Bm]);
-> "[n]" (vector add (vector add (matrix row [Bi] 1) (vector scale (matrix row [Bi] 2) 0)) (vector scale (matrix row [Bi] 3) 2));
-> "[cz]" "[0, 0, 1]";
#3. Angle plane normal / c axis with atan2 of |n x c| and n.c, robust near 0 and 90 deg;
-> "[ang]" (math atan2 (vector norm (vector cross [n] [cz])) (vector dot [n] [cz]));
#4. Twinning shear of 10-12 s = gamma^2 - 3 / sqrt 3 gamma, d-spacing 1/|n|;
-> "[s]" (math calc "([g]^2-3)/(sqrt(3)*[g])");
concat "c/a = " (math round [g] 4) "; angle between (10-12) normal and c = " (math round (math calc "[ang]*180/pi") 2) " deg, so the twin rotates c by " (math round (math calc "2*[ang]*180/pi") 2) " deg; twin shear s = " (math round [s] 4) " (c/a < sqrt 3: sign " (math signum [s]) "); d(10-12) = " (math round (math calc (concat "1/" (vector norm [n]))) 4) " A, check " (math round (math calc (concat "1/" (math hypot (math calc "2/(sqrt(3)*[a])") (math calc "2/[c]")))) 4) " A"
Result
mentdb
c/a = 1.5873; angle between (10-12) normal and c = 42.5 deg, so the twin rotates c by 85 deg; twin shear s = -0.1748 (c/a < sqrt 3: sign -1); d(10-12) = 1.7266 A, check 1.7266 A
Go further
  • Compare with magnesium, c/a = 1.624
  • Compute the angle between [2-1-10] and [10-10] directions
executed on MentDB

X-ray crystal density and the Avogadro constant

The XRCD method counts the atoms of a 28Si sphere: N_A = 8 M / (rho a^3). The script compares the result with the exact SI 2019 value (the same in 'constant' and 'const get'), shows the sensitivity to a, and solves the reverse problem for NaCl and iron with cube roots.

Pipeline
01 · avogadroN_A = 8 M / rho a^3math calc→02 · constantN_A from the constantsconstant physics_n→03 · exactExact SI 2019 valueconst get→04 · cube roota of NaClmath cbrt→05 · n-th roota of bcc ironmath nroot
Script
admin
#1. XRCD route to the Avogadro constant: silicon sphere density and molar mass of enriched 28Si;
-> "[rho]" 2320.070; -> "[M]" 27.97697e-3; -> "[a]" 543.0996e-12;
#2. N_A = 8 M / rho a^3, 8 atoms in the diamond cell, against the exact SI 2019 value;
-> "[NA]" (math calc "8*[M]/([rho]*[a]^3)");
-> "[cst]" (string regexp_extract (constant physics_n) "[-0-9.]+E?[-0-9]*" 0);
-> "[si]" (const get "NA");
#3. Reverse problem: lattice parameter of NaCl from its density 2165 kg/m3, Z = 4, cube root;
-> "[aN]" (math cbrt (math calc (concat "4*58.443e-3/(2165*" [si] ")")));
#4. Same with nroot for a hypothetical bcc iron cell of density 7874, Z = 2;
-> "[aFe]" (math nroot (math calc (concat "2*55.845e-3/(7874*" [si] ")")) 3);
-> "[da]" (math calc "-3*0.1e-15/[a]*1e9");
concat "N_A = " (math round_sig (math calc "[NA]/1e23") 9) "E23 /mol, SI 2019 " [cst] " = " [si] " (gap " (math round (math calc "([NA]-[si])/[si]*1e9") 1) " ppb, an error of only 0.1 fm on a shifts N_A by " (math round [da] 0) " ppb); a(NaCl) = " (math round (math calc "[aN]*1e10") 4) " A, a(Fe) = " (math round (math calc "[aFe]*1e10") 4) " A"
Result
mentdb
N_A = 6.02214387E23 /mol, SI 2019 6.02214076E23 = 6.02214076E23 (gap 517.1 ppb, an error of only 0.1 fm on a shifts N_A by -552 ppb); a(NaCl) = 5.6389 A, a(Fe) = 2.8665 A
Go further
  • Propagate the uncertainties of rho, M and a with num uncertainty
  • Get the density of copper from a = 3.615 A
executed on MentDB

Lattice parameters from an orientation matrix by QR

A diffractometer gives the orientation matrix UB of a monoclinic crystal in the lab frame. A QR decomposition separates the rotation Q from the upper-triangular cell matrix R, from which the script reads a, b, c, beta and the volume.

Pipeline
01 · orientationUB = Rz Rx B0matrix mul→02 · decomposeUB = Q Rmatrix qr→03 · columnsCell vectorsmatrix col→04 · anglebeta from the dot productvector dot→05 · volumeabs det Rmatrix det
Script
admin
#1. Orientation matrix of a monoclinic crystal on a diffractometer: cell 5.20, 7.10, 9.40 A, beta 104 deg, rotated in the lab frame;
-> "[be]" (math calc "104*pi/180");
-> "[B0]" (concat "[[5.20, 0, " (math calc (concat "9.40*" (math cos [be]))) "], [0, 7.10, 0], [0, 0, " (math calc (concat "9.40*" (math sin [be]))) "]]");
-> "[Rz]" (concat "[[" (math cos 0.5) ", " (math calc "-sin(0.5)") ", 0], [" (math sin 0.5) ", " (math cos 0.5) ", 0], [0, 0, 1]]");
-> "[Rx]" (concat "[[1, 0, 0], [0, " (math cos 0.3) ", " (math calc "-sin(0.3)") "], [0, " (math sin 0.3) ", " (math cos 0.3) "]]");
-> "[UB]" (matrix mul (matrix mul [Rz] [Rx]) [B0]);
#2. QR splits UB into a rotation Q and an upper triangular cell matrix R, the standard setting;
-> "[R]" (json load select "rx_p4b_ub" (matrix qr [UB]) /R);
-> "[Q]" (json select "rx_p4b_ub" /Q);
json unload "rx_p4b_ub";
#3. Cell parameters back from the columns of R, angle beta from the dot product;
-> "[ra]" (matrix col [R] 1); -> "[rb]" (matrix col [R] 2); -> "[rc]" (matrix col [R] 3);
-> "[beta]" (math calc (concat "acos(" (vector dot [ra] [rc]) "/(" (vector norm [ra]) "*" (vector norm [rc]) "))*180/pi"));
#4. Volume = abs det R = abs det UB, and Q is orthogonal: det = +1 or -1, a Householder reflection;
concat "a = " (math round (vector norm [ra]) 4) ", b = " (math round (vector norm [rb]) 4) ", c = " (math round (vector norm [rc]) 4) " A, beta = " (math round [beta] 3) " deg; V = " (math round (math abs (matrix det [R])) 3) " A3 = " (math round (math abs (matrix det [UB])) 3) "; det Q = " (math round (matrix det [Q]) 6)
Result
mentdb
a = 5.2, b = 7.1, c = 9.4 A, beta = 104 deg; V = 336.739 A3 = 336.739; det Q = -1
Go further
  • Flip the sign of a row of R to get a proper rotation with det Q = +1
  • Build the reciprocal cell from R^-1 transposed

Also in Sciences

MA
Mathematician
Calculus, linear algebra, number theory and complex numbers: exact and numerical in the same script.
10 recipes
MT
Math teacher & student
Equations, tables of values, fractions and probabilities: exact answers and their decimal check.
10 recipes
PH
Physicist
Mechanics with units, constants and differential equations: trajectories, pendulums, oscillators.
10 recipes
OP
Optics & photonics
Refraction, lenses, photons and diffraction with SI units and CODATA constants.
10 recipes
TH
Thermodynamics & heat transfer
Conduction, heat pumps, ideal gases and cooling laws with units, integrals and ODEs.
10 recipes
NR
Nuclear & radiation physicist
Decay, half-lives, generators and shielding: dates, fits and solvers for radiation work.
10 recipes
AQ
Atomic & quantum physicist
Spectral lines, matter waves, blackbody radiation and CODATA constants with their uncertainties.
10 recipes
AE
Acoustics engineer
Decibels, Doppler, spectra and reverberation: the daily formulas of noise and room acoustics.
10 recipes
CN
Cryptographer & number theorist
RSA and Diffie-Hellman by hand, primes, modular arithmetic and exact big integers.
10 recipes
OR
Operations researcher
Optimisation with constraints, Markov chains for maintenance, newsvendor inventory: solve and plot.
10 recipes
ME
Metrologist
Uncertainty budgets following the GUM, exact decimals and calibration curves.
10 recipes
MP
Mathematical physicist
Special functions at work: Bessel modes, Hermite and Laguerre states, gamma and zeta, with the CAS.
10 recipes