计算化学公社

 找回密码 Forget password
 注册 Register
Views: 1707|回复 Reply: 11
打印 Print 上一主题 Last thread 下一主题 Next thread

[综合交流] 如何用origin分段拟合求Tg?

[复制链接 Copy URL]

21

帖子

0

威望

689

eV
积分
710

Level 4 (黑子)

跳转到指定楼层 Go to specific reply
楼主
各位老师好,我在NPT的条件下对聚乙烯进行降温,然后用体积和时间进行作图。文献上说在Tg附近体积会发生突变,可以用分段拟合求交点的方式求Tg。我实际操作下来发现,前后拟合区间不同,得到的结果上下波动极大,请问有什么好的方法对数据进行拟合吗?或者如何合理地选择拟合的区间呢?(因为模拟的PE分子量很大,最终的Tg应该在300K-400K之间)
我刚刚接触模拟,请各位老师指教,谢谢。




517

帖子

1

威望

2414

eV
积分
2951

Level 5 (御坂)

12#
发表于 Post on 2024-2-19 11:10:52 | 只看该作者 Only view this author
也许可以试试看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)

  5. library(segmented)

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

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

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

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

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

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

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

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

  16. Boot restarting based on 7 samples. Last fit:
  17. 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。

53

帖子

0

威望

313

eV
积分
366

Level 3 能力者

11#
发表于 Post on 2024-2-17 22:20:53 | 只看该作者 Only view this author
计算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有个转折呢

135

帖子

1

威望

4135

eV
积分
4290

Level 6 (一方通行)

10#
发表于 Post on 2024-2-12 02:05:23 | 只看该作者 Only view this author
王二葛 发表于 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
十八介姑娘一蕾花呀,白白介牙齿、红红介嘴唇,得人惜

444

帖子

0

威望

2580

eV
积分
3024

Level 5 (御坂)

娃娃儿鱼

9#
发表于 Post on 2023-5-28 18:19:35 | 只看该作者 Only view this author
古德猫宁341 发表于 2023-5-26 17:27
谢谢老师您的回答。我对python不是很熟,我不太明白脚本里面的m和addtion是如何控制选区的,能请您讲解一 ...

这里的m是控制最少几个点成线,addition是排除Tg附近addition个采样点的干扰,就是满足正常作玻璃化转变温度的基本要求
真·探

21

帖子

0

威望

689

eV
积分
710

Level 4 (黑子)

8#
 楼主 Author| 发表于 Post on 2023-5-26 17:27:22 | 只看该作者 Only view this author
hdhxx123 发表于 2023-5-25 14:40
我这儿有个简单脚本是椭偏数据中温度和厚度变化找Tg,找朋友帮忙写的,你可以参考一下。主要是调控m和addit ...

谢谢老师您的回答。我对python不是很熟,我不太明白脚本里面的m和addtion是如何控制选区的,能请您讲解一下吗?

444

帖子

0

威望

2580

eV
积分
3024

Level 5 (御坂)

娃娃儿鱼

7#
发表于 Post on 2023-5-25 14:40:37 | 只看该作者 Only view this author
本帖最后由 hdhxx123 于 2023-5-25 14:45 编辑

我这儿有个简单脚本是椭偏数据中温度和厚度变化找Tg,找朋友帮忙写的,你可以参考一下。主要是调控m和addition控制选区,实验上一般在Tg正负20-30K温度外取点

Thickness_Temperature.py

5.07 KB, 下载次数 Times of downloads: 17

真·探

135

帖子

1

威望

4135

eV
积分
4290

Level 6 (一方通行)

6#
发表于 Post on 2023-5-25 13:31:16 | 只看该作者 Only view this author
数据已收到

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

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

至于怎么确定突变点,对于此函数我没有经验,抱歉
十八介姑娘一蕾花呀,白白介牙齿、红红介嘴唇,得人惜

48

帖子

0

威望

4930

eV
积分
4978

Level 6 (一方通行)

5#
发表于 Post on 2023-5-25 13:03:48 | 只看该作者 Only view this author

21

帖子

0

威望

689

eV
积分
710

Level 4 (黑子)

4#
 楼主 Author| 发表于 Post on 2023-5-25 12:57:55 | 只看该作者 Only view this author
王二葛 发表于 2023-5-25 12:34
如果按“突变”去拟合,那么确实没有一个良好定义的点

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

谢谢老师您的回答。我用origin求体积对温度的二阶导数,得到的图像波动也很大。由于lammps模拟得到的温度不是线性变化的,无法使用savitzky-golay平滑对导数图像进行平滑处理,这种情况该怎么办呢?我已上传数据,请老师指教。 温度-体积.xlsx (10.67 KB, 下载次数 Times of downloads: 2)

135

帖子

1

威望

4135

eV
积分
4290

Level 6 (一方通行)

3#
发表于 Post on 2023-5-25 12:34:45 | 只看该作者 Only view this author
如果按“突变”去拟合,那么确实没有一个良好定义的点

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

建议把数据发上来
十八介姑娘一蕾花呀,白白介牙齿、红红介嘴唇,得人惜

323

帖子

0

威望

1986

eV
积分
2309

Level 5 (御坂)

2#
发表于 Post on 2023-5-25 10:06:58 | 只看该作者 Only view this author
感觉最前面取一些点,后面取一些点分别拟合合理点,毕竟中间的都弯曲了

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2026-2-26 14:53 , Processed in 0.195241 second(s), 29 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list