Nothing
#
# This file is part of the cry package
#
#########
### Functions for calculations related to the crystal lattice
#########
#' Calculation of useful lattice parameters
#'
#' @param a A real number. One of the unit cell's side lengths, in angstroms.
#' @param b A real number. One of the unit cell's side lengths, in angstroms.
#' @param c A real number. One of the unit cell's side lengths, in angstroms.
#' @param aa A real number. One of the unit cell's angles, in degrees.
#' @param bb A real number. One of the unit cell's angles, in degrees.
#' @param cc A real number. One of the unit cell's angles, in degrees.
#' @return A named vector of real numbers and length 16. The names are:
#' \itemize{
#' \item{sa. Sine(aa)}
#' \item{sb. Sine(bb)}
#' \item{sc. Sine(cc)}
#' \item{ca. Cosine(aa)}
#' \item{cb. Cosine(bb)}
#' \item{cc. Cosine(cc)}
#' \item{ar. a* (reciprocal cell side length)}
#' \item{br. b* (reciprocal cell side length)}
#' \item{cr. c* (reciprocal cell side length)}
#' \item{sar. Sine of a reciprocal cell angle}
#' \item{sbr. Sine of a reciprocal cell angle}
#' \item{scr. Sine of a reciprocal cell angle}
#' \item{car. Cosine of a reciprocal cell angle}
#' \item{cbr. Cosine of a reciprocal cell angle}
#' \item{ccr. Cosine of a reciprocal cell angle}
#' \item{V. Volume of the unit cell in cubic angstroms}
#' }
#' @examples
#' datadir <- system.file("extdata",package="cry")
#' fname <- file.path(datadir,"1dei_phases.mtz")
#' hdr <- readMTZHeader(fname,message=FALSE)
#' ucell <- hdr$CELL
#' vtmp <- lattice_stuff(ucell[1],ucell[2],ucell[3],ucell[4],ucell[5],ucell[6])
#' vtmp[1:3]
#' vtmp[4:6]
#' vtmp[7:9]
#' vtmp[10:12]
#' vtmp[13:15]
#' vtmp[16]
#' @export
lattice_stuff <- function(a,b,c,aa,bb,cc)
{
# Save cc in another variable, as this is used later
gg <- cc
aa <- aa*pi/180
bb <- bb*pi/180
gg <- gg*pi/180
sa <- sin(aa)
sb <- sin(bb)
sc <- sin(gg)
ca <- cos(aa)
cb <- cos(bb)
cc <- cos(gg)
# To avoid NaN generated by rounding off errors, use for cell-derived quantities formulas
# derived previously by computationa crystallographers
sang <- 0.5*(aa+bb+gg)
V2 <- sqrt(sin(sang-aa)*sin(sang-bb)*sin(sang-gg)*sin(sang))
V <- 2*a*b*c*V2
ar <- b*c*sa/V
br <- a*c*sb/V
cr <- a*b*sc/V
car <- (cb*cc-ca)/(sb*sc)
cbr <- (ca*cc-cb)/(sa*sc)
ccr <- (ca*cb-cc)/(sa*sb)
sar <- sqrt(1-car*car)
sbr <- sqrt(1-cbr*cbr)
scr <- sqrt(1-ccr*ccr)
ll <- c(sa,sb,sc,ca,cb,cc,ar,br,cr,sar,sbr,scr,car,cbr,ccr,V)
names(ll) <- c("SIN_ALPHA","SIN_BETA","SIN_GAMMA",
"COS_ALPHA","COS_BETA","COS_GAMMA",
"A*","B*","C*",
"SIN_ALPHA*","SIN_BETA*","SIN_GAMMA*",
"COS_ALPHA*","COS_BETA*","COS_GAMMA*","V")
return(ll)
}
#' Calculates resolution, given the Miller indices
#'
#' @param h An integer, A Miller index.
#' @param k An integer, A Miller index.
#' @param l An integer, A Miller index.
#' @param a A real number. One of the unit cell's side lengths, in angstroms.
#' @param b A real number. One of the unit cell's side lengths, in angstroms.
#' @param c A real number. One of the unit cell's side lengths, in angstroms.
#' @param aa A real number. One of the unit cell's angles, in degrees.
#' @param bb A real number. One of the unit cell's angles, in degrees.
#' @param cc A real number. One of the unit cell's angles, in degrees.
#' @return A positive, real number. The resolution associated with (h,k,l), in angstroms.
#' @examples
#' datadir <- system.file("extdata",package="cry")
#' fname <- file.path(datadir,"1dei_phases.mtz")
#' hdr <- readMTZHeader(fname,message=FALSE)
#' ucell <- hdr$CELL
#' reso1 <- hkl_to_reso(1,0,0,ucell[1],ucell[2],ucell[3],ucell[4],ucell[5],ucell[6])
#' print(reso1) # Low resolution
#' reso2 <- hkl_to_reso(20,20,20,ucell[1],ucell[2],ucell[3],ucell[4],ucell[5],ucell[6])
#' reso2 # High resolution
#' @export
hkl_to_reso <- function(h,k,l,a,b,c,aa,bb,cc)
{
aa <- aa*pi/180
bb <- bb*pi/180
cc <- cc*pi/180
top <- 1-(cos(aa))^2-(cos(bb))^2-(cos(cc))^2+2*cos(aa)*cos(bb)*cos(cc)
b1 <- h^2*(sin(aa))^2/a^2
b2 <- k^2*(sin(bb))^2/b^2
b3 <- l^2*(sin(cc))^2/c^2
b4 <- 2*h*k*(cos(aa)*cos(bb)-cos(cc))/(a*b)
b5 <- 2*h*l*(cos(aa)*cos(cc)-cos(bb))/(a*c)
b6 <- 2*k*l*(cos(bb)*cos(cc)-cos(aa))/(b*c)
d2 <- top/(b1+b2+b3+b4+b5+b6)
return(sqrt(d2))
}
#' Crystal system corresponding to given space group.
#'
#' @param gn A natural integer (1,2,3,...). the space group number.
#' @return A character string, the name of the crystal system associated
#' with the given space group. If the input integer does not correspond
#' any existing space group, the function returns NULL and throws a
#' warning.
#' @examples
#' # P1 is part of the TRICLINIC system
#' crystal_system(1)
#'
#' # The object returned is a string
#' csys <- crystal_system(1)
#' class(csys)
#' @export
crystal_system <- function(gn) {
# Given the space group number (gn), returns the crystal system
# bs = 1 TRICLINIC
# bs = 2 MONOCLINIC
# bs = 3 ORTHOROMBIC
# bs = 4 TETRAGONAL
# bs = 5 CUBIC
# bs = 6 HEXAGONAL
# bs = 7 TRIGONAL
if (gn >= 1 & gn <= 2) bs <- 1
if (gn >= 3 & gn <= 15) bs <- 2
if (gn >= 16 & gn <= 74) bs <- 3
if (gn >= 75 & gn <= 142) bs <- 4
if (gn >= 143 & gn <= 167) bs <- 7
if (gn >= 168 & gn <= 194) bs <- 6
if (gn >= 195 & gn <= 230) bs <- 5
if (gn < 1 | gn > 230) {
warning("No space group associated with this number.")
return(NULL)
}
bs_name <- c("TRICLINIC","MONOCLINIC","ORTHOROMBIC","TETRAGONAL","CUBIC","HEXAGONAL","TRIGONAL")
return(bs_name[bs])
}
#' Crystal family corresponding to given space group.
#'
#' @param gn A natural integer (1,2,3,...). the space group number.
#' @return A character string, the name of the crystal family associated
#' with the given space group. If the input integer does not correspond
#' any existing space group, the function returns NULL and throws a
#' warning.
#' @examples
#' # P1 is part of the TRICLINIC family
#' crystal_family(1)
#'
#' # The object returned is a string
#' cfam <- crystal_family(1)
#' class(cfam)
#' @export
crystal_family <- function(gn) {
# Given the space group number (gn), returns the crystal family
# bs = 1 TRICLINIC
# bs = 2 MONOCLINIC
# bs = 3 ORTHOROMBIC
# bs = 4 TETRAGONAL
# bs = 5 CUBIC
# bs = 6 HEXAGONAL
if (gn >= 1 & gn <= 2) bs <- 1
if (gn >= 3 & gn <= 15) bs <- 2
if (gn >= 16 & gn <= 74) bs <- 3
if (gn >= 75 & gn <= 142) bs <- 4
if (gn >= 143 & gn <= 194) bs <- 6
if (gn >= 195 & gn <= 230) bs <- 5
if (gn < 1 | gn > 230) {
warning("No space group associated with this number.")
return(NULL)
}
bs_name <- c("TRICLINIC","MONOCLINIC","ORTHOROMBIC","TETRAGONAL","CUBIC","HEXAGONAL")
return(bs_name[bs])
}
#' Matrix for cell orthogonalisation (first choice)
#'
#' Given the cell parameters, this function returns a matrix for
#' transforming fractional to orthogonal coordinates, corresponding
#' to the first choice in Giacovazzo's book.
#'
#' @param a A real number. One of the unit cell's side lengths, in angstroms.
#' @param b A real number. One of the unit cell's side lengths, in angstroms.
#' @param c A real number. One of the unit cell's side lengths, in angstroms.
#' @param aa A real number. One of the unit cell's angles, in degrees.
#' @param bb A real number. One of the unit cell's angles, in degrees.
#' @param cc A real number. One of the unit cell's angles, in degrees.
#' @return A \eqn{3\times }$ matrix \eqn{M} that transforms a \eqn{3\times 1} vector of fractional
#' coordinates into a \eqn{3\times 1} vector of orthogonal coordinates.
#' @examples
#' # Fractional coordinates
#' Xf = c(0.1,0.4,0.8)
#'
#' # Orthorombic unit cell
#' M = xtal_mat01(10,40,20,90,90,90)
#'
#' # Cartesian coordinates
#' Xc = M%*%Xf
#' @export
xtal_mat01 <- function(a,b,c,aa,bb,cc) {
lp <- lattice_stuff(a,b,c,aa,bb,cc)
m_1 <- matrix(c(a,0,0,b*lp[6],b*lp[3],0,c*lp[5],-c*lp[2]*lp[13],1/lp[9]),
nrow=3,ncol=3,byrow=TRUE)
return(m_1)
}
#' Matrix for cell orthogonalisation (second choice)
#'
#' Given the cell parameters, this function returns a matrix for
#' transforming fractional to orthogonal coordinates, corresponding
#' to the second choice in Giacovazzo's book.
#'
#' @param a A real number. One of the unit cell's side lengths, in angstroms.
#' @param b A real number. One of the unit cell's side lengths, in angstroms.
#' @param c A real number. One of the unit cell's side lengths, in angstroms.
#' @param aa A real number. One of the unit cell's angles, in degrees.
#' @param bb A real number. One of the unit cell's angles, in degrees.
#' @param cc A real number. One of the unit cell's angles, in degrees.
#' @return A \eqn{3\times }$ matrix \eqn{M} that transforms a \eqn{3\times 1} vector of fractional
#' coordinates into a \eqn{3\times 1} vector of orthogonal coordinates.
#' @examples
#' # Fractional coordinates
#' Xf = c(0.1,0.4,0.8)
#'
#' # Orthorombic unit cell
#' M = xtal_mat02(10,40,20,90,90,90)
#'
#' # Cartesian coordinates
#' Xc = M%*%Xf
#' @export
xtal_mat02 <- function(a,b,c,aa,bb,cc) {
lp <- lattice_stuff(a,b,c,aa,bb,cc)
m_1 <- matrix(c(1/lp[7],-lp[15]/(lp[7]*lp[12]),a*lp[5],0,1/(lp[8]*lp[12]),
b*lp[4],0,0,c),nrow=3,ncol=3,byrow=TRUE)
return(m_1)
}
#' From fractional to orthogonal coordinates
#'
#' This function transforms any number of fractional coordinates \eqn{(x_f,y_f,z_f)},
#' arranged as a vector or in a matrix or data frame, into the corresponding number
#' of orthogonal coordinates \eqn{(x,y,z)}, arranged in the same format.
#' \enumerate{
#' \item ochoice = 1: X axis along a; Y axis normal to a, in the (a,b) plane;
#' Z axis normal to X and Y (and therefore parallel to
#' c*).
#' \item ochoice = 2: this is also called "Cambridge setting". The X axis is
#' along a*; the Y axis lies in the (a*,b*) plane; the Z
#' axis is, consequently, along c.
#' }
#'
#' @param xyzf A vector or \eqn{n\times 3} matrix or data frame of fractional crystal coordinates.
#' @param a A real number. One of the unit cell's side lengths, in angstroms.
#' @param b A real number. One of the unit cell's side lengths, in angstroms.
#' @param c A real number. One of the unit cell's side lengths, in angstroms.
#' @param aa A real number. One of the unit cell's angles, in degrees.
#' @param bb A real number. One of the unit cell's angles, in degrees.
#' @param cc A real number. One of the unit cell's angles, in degrees.
#' @param ochoice A natural integer indicating the choice of orthogonal transformation.
#' 1 corresponds to the first choice and 2 to the second choice in
#' Giacovazzo's book (see \code{\link{xtal_mat01}} and \code{\link{xtal_mat02}}).
#' @return A \eqn{n\times 3} matrix or data frame of orthogonal coordinates corresponding
#' to the fractional coordinates provided in the input.
#' @examples
#' # Matrix containing 3 fractional coordinates
#' xyzf <- matrix(c(0.1,0.2,0.3,0.2,0.6,0.7,0.15,0.28,0.55),ncol=3,byrow=TRUE)
#'
#' # Cartesian coordinates
#' xyz <- frac_to_orth(xyzf,10,30,20,90,90,90,1)
#' @export
frac_to_orth <- function(xyzf, a, b, c, aa, bb, cc, ochoice=1) {
# Check xyzf object is either vector, matrix or data.frame
if (!is.vector(xyzf) & !is.matrix(xyzf) & !is.data.frame(xyzf))
stop("Input object is not of a valid type (vector, matrix or data.frame)")
# If xyzf is a vector, turn it into a matrix
if (is.vector(xyzf)) dim(xyzf) <- c(1,length(xyzf))
# If xyzf has less or more than 3 columns stop
if (dim(xyzf)[2] != 3) stop("Input object has more or less than 3 columns")
# If xyzf is a data.frame turn it into a matrix
if (is.data.frame(xyzf)) xyzf <- as.matrix(xyzf)
# Orthogonalization matrix
if (ochoice == 1) M_1 <- xtal_mat01(a,b,c,aa,bb,cc)
if (ochoice == 2) M_1 <- xtal_mat02(a,b,c,aa,bb,cc)
# Transformed coordinates
#x <- av[1]*xf+bv[1]*yf+cv[1]*zf
#y <- av[2]*xf+bv[2]*yf+cv[2]*zf
#z <- av[3]*xf+bv[3]*yf+cv[3]*zf
xyz <- xyzf%*%M_1
# Turn matrix back into a data frame
xyz <- as.data.frame(xyz)
colnames(xyz) <- c("x","y","z")
return(xyz)
}
#' From orthogonal to fractional coordinates
#'
#' This function transforms any number of orthogonal coordinates \eqn{(x,y,z)},
#' arranged as a vector or in a matrix or data frame, into the corresponding number
#' of fractional coordinates \eqn{(x_f,y_f,z_f)}, arranged in the same format.
#' \enumerate{
#' \item ochoice = 1: X axis along a; Y axis normal to a, in the (a,b) plane;
#' Z axis normal to X and Y (and therefore parallel to
#' c*).
#' \item ochoice = 2: this is also called "Cambridge setting". The X axis is
#' along a*; the Y axis lies in the (a*,b*) plane; the Z
#' axis is, consequently, along c.
#' }
#'
#' @param xyz A vector or \eqn{n\times 3} matrix or data frame of orthogonal crystal coordinates.
#' @param a A real number. One of the unit cell's side lengths, in angstroms.
#' @param b A real number. One of the unit cell's side lengths, in angstroms.
#' @param c A real number. One of the unit cell's side lengths, in angstroms.
#' @param aa A real number. One of the unit cell's angles, in degrees.
#' @param bb A real number. One of the unit cell's angles, in degrees.
#' @param cc A real number. One of the unit cell's angles, in degrees.
#' @param ochoice A natural integer indicating the choice of orthogonal transformation.
#' 1 corresponds to the first choice and 2 to the second choice in
#' Giacovazzo's book (see \code{\link{xtal_mat01}} and \code{\link{xtal_mat02}}).
#' @return A \eqn{n\times 3} matrix or data frame of fractional coordinates corresponding
#' to the orthogonal coordinates provided in the input.
#' @examples
#' # Matrix containing 3 orthogonal coordinates
#' xyz <- matrix(c(5, 15, 10, 2, 10, 8, 1, 1, 1),ncol=3,byrow=TRUE)
#'
#' # Fractional coordinates
#' xyzf <- orth_to_frac(xyz,10,30,20,90,90,90,1)
#' @export
orth_to_frac <- function(xyz,a,b,c,aa,bb,cc,ochoice=1) {
# Check xyz object is either vector, matrix or data.frame
if (!is.vector(xyz) & !is.matrix(xyz) & !is.data.frame(xyz))
stop("Input object is not of a valid type (vector, matrix or data.frame)")
# If xyz is a vector, turn it into a matrix
if (is.vector(xyz)) dim(xyz) <- c(1,length(xyz))
# If xyz has less or more than 3 columns stop
if (dim(xyz)[2] != 3) stop("Input object has more or less than 3 columns")
# If xyz is a data.frame turn it into a matrix
if (is.data.frame(xyz)) xyz <- as.matrix(xyz)
# Orthogonalization matrix
if (ochoice == 1) M_1 <- xtal_mat01(a,b,c,aa,bb,cc)
if (ochoice == 2) M_1 <- xtal_mat02(a,b,c,aa,bb,cc)
# Inverse matrix
M <- solve(M_1)
# Transform coordinates
xyzf <- xyz%*%M
# Turn matrix back into a data frame
xyzf <- as.data.frame(xyzf)
colnames(xyzf) <- c("xf","yf","zf")
return(xyzf)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.