This is a use case for the cry
and rgl
packages. This package provides tools for selected area electron
diffraction (SAED) pattern and crystal structure visualization using the
rgl
and cry
packages. In particular, the
cry_demo()
and dp_demo()
function read files
in CIF (Crystallographic Information Framework) format and display SAED
patterns and crystal structures. The dp_demo()
function
also performs simple simulations of powder X-ray diffraction (PXRD)
patterns, and the results can be saved to a file in the working
directory.
The package has been tested on several platforms, including Linux on Crostini with a Core™ m3-8100Y Chromebook, I found that even on this low-powered platform, the performance was acceptable.
install.packages("rgl.cry")
A CIF file is read, and a reciprocal lattice map with a cell widget is drawn. In this example, a file is not specified, so the system default is used.
dp_demo()
A CIF file is read, and a crystal structure with a axis widget is drawn.
cry_demo()
The following is a working snapshot. Note that while this WebGL graphics is interactive, the behavior is simulated using JavaScript instead of being processed in actual R.
If an error occurs when specifying a CIF file with the message “Error in sites[i, 2] : subscript out of bounds,” it may be that the CIF file has not been read correctly. In this case, please try the following troubleshooting steps:
> remove.packages("cry")
> library(devtools)
> devtools::install_github("SaitouToshihide/cry", type="source", ref="2df5dcc")
The crystal and diffraction pattern are aligned and displayed.
align("a")
align("c*")
align("30 30") # x, y (deg)
Select one or more atoms or reciprocal lattice points in the window. The labels and Miller indices of the selected atoms or lattice points will be displayed.
> dp_demo()
[1] 1
> select()
To select points, drag the left mouse button.
To finish, press ESC.
.
[1] "1 1 -3" "1 0 -2" "2 0 -2" "1 1 -1" "1 0 0" "2 0 0" "1 1 1" "1 0 2"
> cry_demo()
[1] 2
> select()
To select points, drag the left mouse button.
To finish, press ESC.
.
[1] "Ti1" "Ti1"
dp_demo()
can perform PXRD pattern simulation. The
result is saved as a file in the current directory by specifying options
like this:
## Output the simulation results of the PXDR pattern as a file.
dp_demo(xrd = TRUE)
The file looks like this:
% sort -n +7 rgl.cry.dp.demo.2024-02-26_000000.dat
h k l d absF lp twotheta
137 0 0 0 Inf 308.117522 Inf 0.00000
136 -1 0 0 5.583222 17.389409 101.955123 15.87322
138 1 0 0 5.583222 17.389409 101.955123 15.87322
128 0 -2 0 4.453238 40.628676 63.820476 19.93782
146 0 2 0 4.453238 40.628676 63.820476 19.93782
129 1 -2 0 4.142661 17.444203 54.849532 21.44963
145 -1 2 0 4.142661 17.444203 54.849532 21.44963
You can plot the data as follows:
df <- read.table(
file = "data/rgl.cry.dp.demo.2024-02-26_000000.dat",
header = TRUE,
sep = "",
comment.char = "#",
blank.lines.skip = TRUE,
skipNul = TRUE)
df <- df[-which(df$absF == 0 | df$d == Inf | df$twotheta == 0), ]
df2 <- data.frame(twotheta = df$twotheta, I = (df$lp * df$absF^2))
df2 <- aggregate(df2$I, by = list(df2$twotheta), sum)
names(df2) <- c("twotheta", "I")
df2 <- data.frame(twotheta = df2$twotheta, I = 100 * df2$I / max(df2$I))
plot(df2,
xlim = c(10, 80), ylim = c(0, 105), type = "p", cex = 0.02, xaxs = "i", yaxs = "r",
xlab = "", ylab = ""
)
mtext(~ italic("2θ"), side = 1, line = 3)
mtext(~ italic("I"), side = 2, line = 3)
segments(x0 = df2$twotheta, y0 = 0, x1 = df2$twotheta, y1 = df2$I)
In the above plot, a line is drawn perpendicular to the expected intensity at the angle at which diffraction occurs. Even if there is multiple diffraction data at the same position, they are not added together. On the other hand, the following plot is a distribution curve rather than a straight line so that it has values even around the diffraction angle, and all the diffraction intensities are added together and displayed.
f1 <- function(x) {
w <- 0.02 # hwhm (= fwhm/2)
r <- 0.01 # ratio of Lorentzian and Gaussian
sum(apply(df2, 1, function(v) { # sum of pseudo-Voigt functions at each position
h <- v["I"]
t <- v["twotheta"]
(h * ((1 - r) * exp(-log(2) * ((x - t) / w)**2) + r / (1 + ((x - t) / w)**2)))
}))
}
x <- do.call(c, sapply(df2$twotheta, function(v) {
as.list(seq(v - 3, v + 3, 0.2)) # A smaller list length result in shorter calculation time.
}))
x <- unique(sort(x))
df3 <- data.frame(twotheta = x, I = sapply(x, f1))
df3 <- data.frame(twotheta = df3$twotheta, I = 100 * df3$I / max(df3$I))
plot(df3,
xlim = c(10, 80), ylim = c(0, 105), type = "l", xaxs = "i", yaxs = "r",
xlab = "", ylab = ""
)
mtext(~ italic("2θ"), side = 1, line = 3)
mtext(~ italic("I"), side = 2, line = 3)
I would like to express my gratitude to the rgl, cry packages, and R for making this work possible. The data for the scattering factors is based on the data from the KEK Report Hanashima. For the coloring of the atoms, I used Helmenstine (2019), and for the atomic radii, I used Wikipedia contributors (2023).
I’m considering implementing Kikuchi line drawing for my studies at https://github.com/SaitouToshihide/rgl.cry/tree/DiffractionCone and my study zone.