rgl.cry

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.

Installation

install.packages("rgl.cry")

Example

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.

Troubleshooting

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:

  1. Refer to GitHub pull request #6: https://github.com/jfoadi/cry/pull/6 for common solutions.
  2. You can either modify the CIF file or try the patched source here:
> remove.packages("cry")
> library(devtools)
> devtools::install_github("SaitouToshihide/cry", type="source", ref="2df5dcc")

Utility functions

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"

Extras

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)

References

Hanashima, T. “Contents of Atomic Scattering Factors’ Table in S. Sasaki (1987), KEK Report 87-3.” Constructed in web page in 2001. https://www2.kek.jp/imss/pf/tools/sasaki/sinram/sinram.html.
Helmenstine, Todd. 2019. “Molecule Atom Colors – CPK Colors.” https://sciencenotes.org/molecule-atom-colors-cpk-colors/.
Wikipedia contributors. 2023. “Atomic Radius — Wikipedia, the Free Encyclopedia.” https://en.wikipedia.org/w/index.php?title=Atomic_radius&oldid=1179864711.

Acknowledgment

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).

Features under consideration

I’m considering implementing Kikuchi line drawing for my studies at https://github.com/SaitouToshihide/rgl.cry/tree/DiffractionCone and my study zone.

BCC structure along the [0 0 1] axis.

A zoomed-in view.