计算化学公社

标题: 如何用origin分段拟合求Tg？ [打印本页]
作者
Author: 古德猫宁341    时间: 2023-5-25 09:35
标题: 如何用origin分段拟合求Tg？
各位老师好，我在NPT的条件下对聚乙烯进行降温，然后用体积和时间进行作图。文献上说在Tg附近体积会发生突变，可以用分段拟合求交点的方式求Tg。我实际操作下来发现，前后拟合区间不同，得到的结果上下波动极大，请问有什么好的方法对数据进行拟合吗？或者如何合理地选择拟合的区间呢？（因为模拟的PE分子量很大，最终的Tg应该在300K-400K之间）
我刚刚接触模拟，请各位老师指教，谢谢。


(, 下载次数 Times of downloads: 27) (, 下载次数 Times of downloads: 31) (, 下载次数 Times of downloads: 27)

作者
Author: pal    时间: 2023-5-25 10:06
感觉最前面取一些点，后面取一些点分别拟合合理点，毕竟中间的都弯曲了
作者
Author: 王二葛    时间: 2023-5-25 12:34
如果按“突变”去拟合，那么确实没有一个良好定义的点

只能是退而求其次，找二阶导数最大的点。对于有噪声的数据是有点麻烦的

建议把数据发上来
作者
Author: 古德猫宁341    时间: 2023-5-25 12:57
王二葛 发表于 2023-5-25 12:34
如果按“突变”去拟合，那么确实没有一个良好定义的点

只能是退而求其次，找二阶导数最大的点。对于有噪 ...

谢谢老师您的回答。我用origin求体积对温度的二阶导数，得到的图像波动也很大。由于lammps模拟得到的温度不是线性变化的，无法使用savitzky-golay平滑对导数图像进行平滑处理，这种情况该怎么办呢？我已上传数据，请老师指教。 (, 下载次数 Times of downloads: 2)
作者
Author: iuhgnor    时间: 2023-5-25 13:03
可以参考这篇文献使用双曲线函数拟合确定Tg点：
High-Throughput Molecular Dynamics Simulations and Validation of Thermophysical Properties of Polymers for Various Applications
作者
Author: 王二葛    时间: 2023-5-25 13:31
数据已收到

为了便于拟合，令 x = T / 100, y = V / 100,000，得到 y = 1.06 exp( 0.18*(x-1.66) ) + 5.53。此时 R^2=0.999，说明拟合很好

既然此数据符合指数增长，那么从数学上就没法确定一个突变的点，因为 e^x 在左右两侧没有渐进行为，楼主这种方法可能不太行

至于怎么确定突变点，对于此函数我没有经验，抱歉
作者
Author: hdhxx123    时间: 2023-5-25 14:40
本帖最后由 hdhxx123 于 2023-5-25 14:45 编辑

我这儿有个简单脚本是椭偏数据中温度和厚度变化找Tg，找朋友帮忙写的，你可以参考一下。主要是调控m和addition控制选区，实验上一般在Tg正负20-30K温度外取点
作者
Author: 古德猫宁341    时间: 2023-5-26 17:27
hdhxx123 发表于 2023-5-25 14:40
我这儿有个简单脚本是椭偏数据中温度和厚度变化找Tg，找朋友帮忙写的，你可以参考一下。主要是调控m和addit ...

谢谢老师您的回答。我对python不是很熟，我不太明白脚本里面的m和addtion是如何控制选区的，能请您讲解一下吗？
作者
Author: hdhxx123    时间: 2023-5-28 18:19
古德猫宁341 发表于 2023-5-26 17:27
谢谢老师您的回答。我对python不是很熟，我不太明白脚本里面的m和addtion是如何控制选区的，能请您讲解一 ...

这里的m是控制最少几个点成线，addition是排除Tg附近addition个采样点的干扰，就是满足正常作玻璃化转变温度的基本要求
作者
Author: 王二葛    时间: 2024-2-12 02:05
王二葛 发表于 2023-5-25 13:31
数据已收到

为了便于拟合，令 x = T / 100, y = V / 100,000，得到 y = 1.06 exp( 0.18*(x-1.66) ) + 5. ...

最近看到知乎王赟提到了一个观点，曲率半径最小值可认为是函数的拐点
https://www.zhihu.com/question/376116261

把这个函数带入到曲率公式，可得到曲率半径的最小值在 9，也就是拐点在 900 K
作者
Author: dayu8278    时间: 2024-2-17 22:20
计算Tg一般有几种方法。
第一种是先计算体系的tau_alpha，然后绘制tau_alpha和1/T的曲线，对该曲线进行拟合，外推至tau_alpha = 100s所对应的T就是Tg。
第二种就是你这个直接用rho和T进行拟合，也有的是用E_coul以及E_vdw和T进行拟合，得到体系的Tg以及局部Tg‘。
你试试用150-250拟合一段，用400-500拟合另外一段，尽量避开300-400这个区间？确实到Tg附近，rho有个转折呢
作者
Author: Daniel_Arndt    时间: 2024-2-19 11:10
也许可以试试看piecewise linear regression，既然你已经知道应该分成两段了。我当时是用R语言的一个名叫“segmented”的包来做的，相应的网址为 https://cran.r-project.org/web/packages/segmented/index.html ，需要注意的是如果在正式发表的工作中用到“segmented”这个包的话，要引用网址开头提到的几篇文献。安装过程我简单说一下，得先装好R语言和RStudio，然后打开RStudio，在左边输入“install.packages(c("nlme", "MASS", "segmented"))”。

然后在Excel中从第二行开始写好两列数据，在A1单元格写“x”，在B1单元格写“y”，然后单击“File”，再单击“Save”，把“Save as type:”改成“CSV (Comma delimited) (*.csv)”。把csv文件放在桌面上，为方便起见，文件名就是file.csv。我当时是在英文版的Windows上操作的，中文版上有些许不同。

然后打开RStudio，输入以下代码：
  1. #pay attention to the format of the directory
  2. setwd("C:\\Users\\[USERNAME]\\Desktop")
  3. df <- read.csv("file.csv", sep = ",")
  4. View(df)

  6. library(segmented)

  8. #fit simple linear regression model
  9. fit <- lm(y ~ x, data=df)

  11. #fit piecewise regression model to original model, estimating a breakpoint at x=12
  12. segmented.fit <- segmented(fit, seg.Z = ~x, psi=12)

  14. #view summary of segmented model
  15. summary(segmented.fit)
复制代码
代码的第二行就是桌面的路径，里面的“[USERNAME]”得换成用户名。我从未在Windows操作系统上用过中文用户名，如果你是中文用户名的话，建议换个目录存储csv文件，相应的，代码的第二行需要改动。Windows上R语言里面路径的格式有点繁琐，需要留心。代码的第12行需要对所谓的“变点”提供一个初始猜测，这个实际操作中用的就是肉眼观察到的斜率发生变化的区域的自变量的值。

然后，输出如下（我直接用的我的一次操作中的输出）：
  1.     ***Regression Model with Segmented Relationship(s)***

  3. Call:
  4. segmented.lm(obj = fit, seg.Z = ~x, psi = 12)

  6. Estimated Break-Point(s):
  7.          Est. St.Err
  8. psi1.x 13.42  1.179

  10. Meaningful coefficients of the linear terms:
  11.             Estimate Std. Error t value Pr(>|t|)   
  12. (Intercept)  38.1363     0.9392  40.605  < 2e-16 ***
  13. x            -1.0005     0.1233  -8.112 5.73e-11 ***
  14. U1.x          1.0832     0.1439   7.529       NA   
  15. ---
  16. Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  18. Residual standard error: 2.451 on 55 degrees of freedom
  19. Multiple R-Squared: 0.7354,  Adjusted R-squared: 0.7209

  21. Boot restarting based on 7 samples. Last fit:
  22. Convergence attained in 2 iterations (rel. change 1.1647e-08)
复制代码
输出的解读略有点繁琐，如果x小于13.42的话，y=-1.0005x+38.1363。如果x大于等于13.42的话，y=-1.0005*13.42+38.1363+(-1.0005+1.0832)(x-13.42)=0.0827x+23.5998。




欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3