如何在 R 中的能量图中绘制曲线？

Question

我写了以下 R 脚本：

#energy diagram
x <- c(0.1, 0.3, 0.5, 0.7, 0.9 )    #chosen randomly, reaction axis
y <- c(-5.057920, -5.057859, -5.057887,-5.057674, -5.057919 ) #energy of the educt, intermediate, transtition states and product 
plot(x,y, type="p",
     xlim=c(0,1),
     ylim=c(-5.058,-5.0575),
     xlab="reaction axis",
     ylab=expression(paste(E[el] ," / ",10^6," ",kJ/mol)),
     xaxt="n"        #hide x-axis
     )
#h- and v-lines, so i can draw curves by hand
abline(v=seq(0,1,0.1),h=seq(-5.0600,-5.0500,0.00005),col="black",lty=1,lwd=1)
abline(h=c(-5.057920, -5.057859, -5.057887,-5.057674), col="blue", lty=1,lwd=0.7)

是否可以通过看起来像能量图的点绘制曲线。能量图示例如下：

Answer 1

您可以使用样条拟合。沿能量图定义一些点，然后使用样条函数拟合它们。您提供的分数越多，您的合身性就越好。您可以查看 stats 包中的 smooth.splines 函数，了解样条拟合的一种实现方式。

Answer 2

可以做很多工作来简化/矢量化此代码，但对于一个较小的图表来说，它工作得很好：

# get that data
x <- c(0.1, 0.3, 0.5, 0.7, 0.9 )    # reaction axis
y <- c(-5.057920, -5.057859, -5.057887,-5.057674, -5.057919 ) # energies

我将制作一条小贝塞尔曲线将每个点连接到下一个点---这样我们就可以确保平滑的线穿过数据，而不仅仅是靠近它。我会给每个点一个 'control point' 来定义斜率。通过对一个点及其控制点使用相同的 y 值，该点的斜率将为 0。我将该点和控制点之间的偏移量称为 delta。我们将从一对点开始：

library(Hmisc)
delta = 0.15
bezx = c(0.1, 0.1 + delta, 0.3 - delta, 0.3)
bezy = rep(y[1:2], each = 2)
plot(bezx, bezy, type = 'b', col = "gray80")
lines(bezier(bezx, bezy), lwd = 2, col = "firebrick4")

这里我用灰色绘制了点和控制点，用红色绘制了平滑线，这样我们就可以看到发生了什么。它看起来很有希望，让我们把它变成一个我们可以应用于每对点的函数：

bezf = function(x1, x2, y1, y2, delta = 0.15) {
    bezier(x = c(x1, x1 + delta, x2 - delta, x2), y = c(y1, y1, y2, y2))
}

你可以玩玩delta这个参数，我觉得0.1挺好看的

plot(x, y, xlab = "Reaction coordinate", ylab = "E", axes = F)
box(bty = "L")
axis(side = 2)
for(i in 1:(length(x) - 1)) {
    lines(bezf(x1 = x[i], x2 = x[i + 1], y1 = y[i], y2 = y[i + 1], delta = 0.1))
}

您当然可以调整情节、添加标签和 ablines，就像您的原作一样。（使用我的 for 循环和 lines 命令只绘制平滑线。）我留下点来表明我们正在通过它们，而不仅仅是接近它们。

我更喜欢在 ggplot2 中绘图，如果你也这样做，你需要将数据提取到 data.frame:

bezlist = list()
for (i in 1:(length(x) - 1)) {
    bezlist[[i]] = bezf(x1 = x[i], x2 = x[i + 1], y1 = y[i], y2 = y[i + 1], delta = 0.1)
}


xx = unlist(lapply(bezlist, FUN = '[', 'y'))
yy = unlist(lapply(bezlist, FUN = '[', 'y'))

bezdat = data.frame(react = xx, E = yy)
library(ggplot2)
ggplot(bezdat, aes(x = react, y = E)) + 
    geom_line() +
    labs(x = "Reaction coordinate")

如何在 R 中的能量图中绘制曲线？

How to draw the curves in an energy diagram in R?

r

diagram

curve

energy