VMD脚本用于分析CP2K AIMD轨迹文件中分子种类和数量

ma455173220

该脚本通过VMD对CP2K AIMD产生的轨迹文件进行统计分析，遍历轨迹文件并记录随时间变化的不同产物数量，以探究相关反应机理和路径。

本脚本参考了sob老师在CP2K培训班中第一性原理分子动力学部分的实例2：高温下烷烃裂解的模拟的产物各类片段数目的变化分析脚本MD\decane_pyrolysis\product_time.tcl，并做了相应修改。为了优化内存占用，结合了bigdcd脚本(ks.uiuc.edu/Research/vmd/script_library/scripts/bigdcd/bigdcd.tcl)，以防止统计大量帧数时内存占用过高的问题。

使用方法：

• 打开VMD，但不要不要不要加载轨迹文件。
• 点击菜单中的Extensions -> Tk Console。
• 在打开的VMD TkConsole窗口中点击File -> Load File，加载脚本。在加载脚本前，需要根据自己的需求进行一些修改。
  

需要修改的内容：

• 第130行：设定 pbc_parameters，以符合你的系统设置。
• 第132行：设定 timestep，表示加载的轨迹文件中每一帧的时间步长。
• 第134-135行：设定要分析的起始帧和终止帧。
• 第137-138行：设定文件路径。
• 第17-45行：用来设定要分析的分子。根据自己的需求修改这部分内容。脚本中内容是我用来分析C2H5OH高温裂解的C2H2,C2H4,C2H6分子数量的，里面个人觉得写的蛮清晰的，所以这部分需要你自己来编译成你想分析的分子，不算特别难。
  

内存占用情况：

• 下图是未结合bigdcd脚本的内存占用，我只分析了500帧，内存已经顶到了18GB，所以我就没有继续跑，之前测试我的32GB内存电脑大概跑个1500帧左右就炸了（当然了，应该跟你分析的系统大小有关）：
  
  

• 下图是结合了bigdcd脚本的内存占用，该例我分析了2000帧，可以看到内存只占用了500MB左右，继续算也基本不会有太大变化，分析个几万帧应该没啥压力：
  
  

补充：

• 若分析的帧数过多，或者分析的帧数并不是从第一帧开始（比如分析的是5000帧到10000帧；换句话说就是脚本中的startframe并不是0），那么开始运行脚本时可能会卡顿一会儿才正式开始运行，因为那时脚本是正在加载轨迹文件，所以不要乱动。
• 刚出炉的原数据曲线可能会非常凌乱，因为模拟体系较小，分子数目波动过于剧烈，所以需要手动后续用Origin进行平滑处理：Analysis -> Signal processing -> Smooth，Points of Windows按需设定。附一张示例图如下：
  
  

lywbanner

wildon

您好，使用您的脚本时，VMD会提示下面的错误，请问是不是我哪里没有设置正确。我是通过CP2K进行CP2K计算后，获取了xyz轨迹文件，然后按照帖子中的内容对脚本进行了修改。运行后会出现下面的错误，生成的txt为空文件。
bigdcd aborting at frame 1
wrong # args: should be "set varName ?newValue?"

wildon

wildon 发表于 2024-12-16 15:54
您好，使用您的脚本时，VMD会提示下面的错误，请问是不是我哪里没有设置正确。我是通过CP2K进行CP2K计算后 ...

感谢，已解决。

1. # Define your analysis procedure
   
2. proc my_analysis { frame } {
   
3.     global all myfile timestep startframe pbc_parameters mol vmd_frame
   
4.    
   
5.     pbc set $pbc_parameters
   
6.         
   
7.     $all frame $frame
   
8.     mol bondsrecalc top
   
9.     mol reanalyze top
   
10.     pbc wrap -compound fragment
    
11.     pbc join fragment -bondlist
    
12.     mol bondsrecalc top
    
13.     mol reanalyze top   
    
14.     puts "Doing frame $frame"
    
15.    
    
16.     set nfrag [llength [lsort -unique -integer [$all get fragment]]]
    
17.    
    
18. # This is where you set up your molecule analysis
    
19.     set molecules {H2 H2O CO CO2 CH4 NH3 CHN C2H4}
    
20.    
    
21.     array set numc {}
    
22.     foreach molecule $molecules {
    
23.         set numc($molecule) 0
    
24.     }
    
25.    
    
26.     for {set ifrag 0} {$ifrag < $nfrag} {incr ifrag} {
    
27.         # 基本原子计数
    
28.         set nhtmp [[atomselect top "fragment $ifrag and name H"] num]
    
29.         set nctmp [[atomselect top "fragment $ifrag and name C"] num]
    
30.         set notmp [[atomselect top "fragment $ifrag and name O"] num]
    
31.         set nntmp [[atomselect top "fragment $ifrag and name N"] num]
    
32.         
    
33.         # 键合数计数
    
34.         set H_nbtmp_1 [[atomselect top "fragment $ifrag and {name H and numbonds 1}"] num]
    
35.         set O_nbtmp_2 [[atomselect top "fragment $ifrag and {name O and numbonds 2}"] num]
    
36.         set C_nbtmp_2 [[atomselect top "fragment $ifrag and {name C and numbonds 2}"] num]
    
37.         set C_nbtmp_3 [[atomselect top "fragment $ifrag and {name C and numbonds 3}"] num]
    
38.         set C_nbtmp_4 [[atomselect top "fragment $ifrag and {name C and numbonds 4}"] num]
    
39.         set N_nbtmp_3 [[atomselect top "fragment $ifrag and {name N and numbonds 3}"] num]
    
40.         
    
41.         # H2: 2个H原子，每个H有1个键合
    
42.         if {$nhtmp == 2 && $nctmp == 0 && $notmp == 0 && $nntmp == 0 && $H_nbtmp_1 == 2} {
    
43.             incr numc(H2)
    
44.         # H2O: 2个H原子，1个O原子，O有2个键合
    
45.         } elseif {$nhtmp == 2 && $nctmp == 0 && $notmp == 1 && $nntmp == 0 && $O_nbtmp_2 == 1} {
    
46.             incr numc(H2O)
    
47.         # CO: 1个C原子，1个O原子，C有2个键合
    
48.         } elseif {$nhtmp == 0 && $nctmp == 1 && $notmp == 1 && $nntmp == 0 && $C_nbtmp_2 == 1} {
    
49.             incr numc(CO)
    
50.         # CO2: 1个C原子，2个O原子，C有4个键合
    
51.         } elseif {$nhtmp == 0 && $nctmp == 1 && $notmp == 2 && $nntmp == 0 && $C_nbtmp_4 == 1} {
    
52.             incr numc(CO2)
    
53.         # CH4: 1个C原子，4个H原子，C有4个键合
    
54.         } elseif {$nhtmp == 4 && $nctmp == 1 && $notmp == 0 && $nntmp == 0 && $C_nbtmp_4 == 1} {
    
55.             incr numc(CH4)
    
56.         # NH3: 3个H原子，1个N原子，N有3个键合
    
57.         } elseif {$nhtmp == 3 && $nctmp == 0 && $notmp == 0 && $nntmp == 1 && $N_nbtmp_3 == 1} {
    
58.             incr numc(NH3)
    
59.         # CHN: 1个C原子，1个H原子，1个N原子，C和N都有2个键合
    
60.         } elseif {$nhtmp == 1 && $nctmp == 1 && $notmp == 0 && $nntmp == 1 && $C_nbtmp_2 == 1} {
    
61.             incr numc(CHN)
    
62.         # C2H4: 2个C原子，4个H原子，C有3个键合
    
63.         } elseif {$nctmp == 2 && $nhtmp == 4 && $notmp == 0 && $nntmp == 0 && $C_nbtmp_3 == 2} {
    
64.             incr numc(C2H4)
    
65.         }
    
66.         
    
67.         # 清理临时变量
    
68.         unset nhtmp nctmp notmp nntmp H_nbtmp_1 O_nbtmp_2 C_nbtmp_2 C_nbtmp_3 C_nbtmp_4 N_nbtmp_3
    
69.     }
    
70.    
    
71.     set frame_real [expr {$frame + $startframe - 1}]
    
72.     puts -nonewline $myfile [expr $timestep * $frame_real]\
    
73.     foreach molecule $molecules {
    
74.         puts -nonewline $myfile $numc($molecule)\
    
75.     }
    
76.     puts $myfile " "
    
77. 
    
78.     unset nfrag molecules numc ifrag frame_real
    
79. }
    
80. 
    
81. # Call bigdcd command to manage memory and call your analysis procedure for each frame
    
82. proc bigdcd { script type args } {
    
83.     global bigdcd_frame bigdcd_proc bigdcd_lastframe vmd_frame bigdcd_running startframe endframe
    
84.    
    
85.     set bigdcd_running 1
    
86.     set bigdcd_frame 0
    
87.     set bigdcd_lastframe [molinfo top get numframes]
    
88.     set bigdcd_proc $script
    
89.    
    
90.     # backwards "compatibility". type flag is omitted.
    
91.     if {[file exists $type]} {
    
92.         set args [linsert $args 0 $type]
    
93.         set type auto
    
94.     }
    
95.    
    
96.     uplevel #0 trace variable vmd_frame w bigdcd_callback
    
97.    
    
98.     foreach dcd $args {
    
99.         if { $type == "auto" } {
    
100.             mol addfile $dcd first $startframe last $endframe waitfor 0
     
101.         } else {
     
102.             mol addfile $dcd type $type first $startframe last $endframe waitfor 0
     
103.         }
     
104.     }
     
105.    
     
106.     after idle bigdcd_wait
     
107. }
     
108. 
     
109. proc bigdcd_callback { tracedvar mol op } {
     
110.     global bigdcd_frame bigdcd_proc bigdcd_lastframe vmd_frame
     
111.     set msg {}
     
112.    
     
113.     set thisframe $vmd_frame($mol)
     
114.    
     
115.     incr bigdcd_frame
     
116.     if { [catch {uplevel #0 $bigdcd_proc $bigdcd_frame} msg] } {
     
117.         puts stderr "bigdcd aborting at frame $bigdcd_frame\n$msg"
     
118.         bigdcd_done
     
119.         return
     
120.     }
     
121.    
     
122.     animate delete beg $thisframe end $thisframe $mol
     
123.     return $msg
     
124. }
     
125. 
     
126. proc bigdcd_done { } {
     
127.     global bigdcd_running
     
128.    
     
129.     if {$bigdcd_running > 0} {
     
130.         uplevel #0 trace vdelete vmd_frame w bigdcd_callback
     
131.         puts "End of frame"
     
132.         puts "Bigdcd_done"
     
133.         set bigdcd_running 0
     
134.     }
     
135. }
     
136. 
     
137. proc bigdcd_wait { } {
     
138.     global bigdcd_running bigdcd_frame bigdcd_lastframe
     
139.    
     
140.     while {$bigdcd_running > 0} {
     
141.         global bigdcd_oldframe
     
142.         set bigdcd_oldframe $bigdcd_frame
     
143.         
     
144.         # run global processing hooks (including loading of scheduled frames)
     
145.         display update ui
     
146.         
     
147.         # Stop when reaching maximum frame
     
148.         if { $bigdcd_frame >= $bigdcd_lastframe } {
     
149.             bigdcd_done
     
150.         }
     
151.     }
     
152. }
     
153. 
     
154. # Set up periodic boundary condition
     
155. set pbc_parameters {11.117 11.117 11.117  90.0000  90.0000  90.0000}
     
156. # Set up timestep (fs)
     
157. set timestep 0.5
     
158. # Set up frame index (start from 0)
     
159. set startframe 0
     
160. set endframe 1604
     
161. # Set up file path
     
162. set inputfile "C:/Users/Lenovo/Desktop/1/qingniaosuan-pos-1.xyz"
     
163. set outputfile "C:/Users/Lenovo/Desktop/1/qingniaosuan-pos-1.txt"
     
164. 
     
165. # Check if input file exists
     
166. if {![file exists $inputfile]} {
     
167.     puts "Error: Input file not found: $inputfile"
     
168.     return
     
169. }
     
170. 
     
171. # Job starts
     
172. set myfile [open $outputfile a]
     
173. set mol [mol new $inputfile first $startframe last $endframe type xyz waitfor all]
     
174. set all [atomselect top all]
     
175. 
     
176. # Call bigdcd command to manage memory and call your analysis procedure for each frame
     
177. bigdcd my_analysis auto $inputfile
     
178. 
     
179. # Wait for all frames to be processed
     
180. bigdcd_wait
     
181. 
     
182. # Close the result file
     
183. close $myfile

复制代码