计算化学公社

标题: 求助：怎样从cp2k轨迹文件（xyz）追踪并提取出小分子片段的数量随时间的变化 [打印本页]

作者
Author: 吞木木 时间: 2022-3-27 15:40
标题: 求助：怎样从cp2k轨迹文件（xyz）追踪并提取出小分子片段的数量随时间的变化

接触动力学有一段时间了，通过看文献经常能看到下边的图，就是大家从轨迹文件中提取出各个时间段对应的小分子数量，图片中的是关于甲烷燃烧的过程中各个分子的数量岁时间演化的数量，由于编程基础有点差，我一直不知道怎么才能从轨迹文件中采集这些数据，还希望各位老师和同学们能够指点迷津，有什么能够实现这种功能的脚本或者软件，能够解决我遇到的问题，谢谢啦！

(, 下载次数 Times of downloads: 26)

作者
Author: 丁越 时间: 2022-3-27 22:32
same fragment as {element C} && numbonds=1 CO
same fragment as {element C} && numbonds=4 CH4
same fragment as {element C} && numbonds=2 CO2
same fragment as {element O} && numbonds=2 H2O
same fragment as {element O} && numbonds=1 && not same fragment as {element C} O2
按照上述片段选择就可以定义CO、CH4、CO2等了，然后通过for循环统计每帧的片段数目就行了

作者
Author: sobereva 时间: 2022-3-28 10:35
还要注意在脚本中循环时必须实时更新连接关系，参考下文里的相关命令
VMD初始化文件(vmd.rc)我的推荐设置
http://sobereva.com/545（http://bbs.keinsci.com/thread-16834-1-1.html）

作者
Author: 吞木木 时间: 2022-3-28 11:12

丁越发表于 2022-3-27 22:32
same fragment as {element C} && numbonds=1 CO
same fragment as {element C} && numbonds=4 CH4
same ...

谢谢丁老师的回复

作者
Author: 吞木木 时间: 2022-3-28 11:13

sobereva 发表于 2022-3-28 10:35
还要注意在脚本中循环时必须实时更新连接关系，参考下文里的相关命令
VMD初始化文件(vmd.rc)我的推荐设置
...

谢谢卢天老师的回复

作者
Author: 吞木木 时间: 2022-4-30 14:46
本帖最后由吞木木于 2022-4-30 14:48 编辑

丁越发表于 2022-3-27 22:32
same fragment as {element C} && numbonds=1 CO
same fragment as {element C} && numbonds=4 CH4
same ...

老师你好，我写了几行tcl的代码，但是他报错很离谱，您能抽时间帮忙看一下怎么改吗？

set numframe [molinfo top get numframes]

for {set i 0} {$i <= $numframe} {incr i 1} {
set a [atomselect top "same fragment as {element N} && numbonds=2"]
$a frame $i
$a update
puts $i $a;
}

它报错是提示{set i 0}这一行中的0，报错。我实在是不知道怎么改了，麻烦老师了
(, 下载次数 Times of downloads: 25)

作者
Author: 丁越 时间: 2022-5-11 08:26
本帖最后由丁越于 2022-5-11 12:07 编辑

吞木木发表于 2022-4-30 14:46
老师你好，我写了几行tcl的代码，但是他报错很离谱，您能抽时间帮忙看一下怎么改吗？

set numframe [m ...

你看看TCL语法中puts是怎么用的，显然不是这么接的。
set numframe [molinfo top get numframes]
set a [atomselect top "same fragment as {element N} && numbonds=2"]
set fileid [open num.dat w]
for {set i 0} {$i <= $numframe} {incr i 1} {
$a frame $i
animate goto $i
  mol bondsrecalc all
topo retypebonds
$a update
set N [$a num]
  puts $fileid "$i \t [expr double($N)/3]"
}
close $fileid

[expr double($N)/3]"除以3是我假设你的每个分子片段有三个原子

作者
Author: 吞木木 时间: 2022-5-11 16:17

丁越发表于 2022-5-11 08:26
你看看TCL语法中puts是怎么用的，显然不是这么接的。
set numframe [molinfo top get numframes]
set a ...

谢谢老师的批评指正！
我回去再好好学习一下TCL语法

欢迎光临计算化学公社 (http://bbs.keinsci.com/)