VMD脚本用于分析CP2K AIMD轨迹文件中分子种类和数量 - 计算化学公社

# Define your analysis procedure
proc my_analysis { frame } {
global all myfile timestep startframe pbc_parameters mol vmd_frame
pbc set $pbc_parameters
$all frame $frame
mol bondsrecalc top
mol reanalyze top
pbc wrap -compound fragment
pbc join fragment -bondlist
mol bondsrecalc top
mol reanalyze top
puts "Doing frame $frame"
set nfrag [llength [lsort -unique -integer [$all get fragment]]]
# This is where you set up your molecule analysis
set molecules {H2 H2O CO CO2 CH4 NH3 CHN C2H4}
array set numc {}
foreach molecule $molecules {
set numc($molecule) 0
}
for {set ifrag 0} {$ifrag < $nfrag} {incr ifrag} {
# 基本原子计数
set nhtmp [[atomselect top "fragment $ifrag and name H"] num]
set nctmp [[atomselect top "fragment $ifrag and name C"] num]
set notmp [[atomselect top "fragment $ifrag and name O"] num]
set nntmp [[atomselect top "fragment $ifrag and name N"] num]
# 键合数计数
set H_nbtmp_1 [[atomselect top "fragment $ifrag and {name H and numbonds 1}"] num]
set O_nbtmp_2 [[atomselect top "fragment $ifrag and {name O and numbonds 2}"] num]
set C_nbtmp_2 [[atomselect top "fragment $ifrag and {name C and numbonds 2}"] num]
set C_nbtmp_3 [[atomselect top "fragment $ifrag and {name C and numbonds 3}"] num]
set C_nbtmp_4 [[atomselect top "fragment $ifrag and {name C and numbonds 4}"] num]
set N_nbtmp_3 [[atomselect top "fragment $ifrag and {name N and numbonds 3}"] num]
# H2: 2个H原子，每个H有1个键合
if {$nhtmp == 2 && $nctmp == 0 && $notmp == 0 && $nntmp == 0 && $H_nbtmp_1 == 2} {
incr numc(H2)
# H2O: 2个H原子，1个O原子，O有2个键合
} elseif {$nhtmp == 2 && $nctmp == 0 && $notmp == 1 && $nntmp == 0 && $O_nbtmp_2 == 1} {
incr numc(H2O)
# CO: 1个C原子，1个O原子，C有2个键合
} elseif {$nhtmp == 0 && $nctmp == 1 && $notmp == 1 && $nntmp == 0 && $C_nbtmp_2 == 1} {
incr numc(CO)
# CO2: 1个C原子，2个O原子，C有4个键合
} elseif {$nhtmp == 0 && $nctmp == 1 && $notmp == 2 && $nntmp == 0 && $C_nbtmp_4 == 1} {
incr numc(CO2)
# CH4: 1个C原子，4个H原子，C有4个键合
} elseif {$nhtmp == 4 && $nctmp == 1 && $notmp == 0 && $nntmp == 0 && $C_nbtmp_4 == 1} {
incr numc(CH4)
# NH3: 3个H原子，1个N原子，N有3个键合
} elseif {$nhtmp == 3 && $nctmp == 0 && $notmp == 0 && $nntmp == 1 && $N_nbtmp_3 == 1} {
incr numc(NH3)
# CHN: 1个C原子，1个H原子，1个N原子，C和N都有2个键合
} elseif {$nhtmp == 1 && $nctmp == 1 && $notmp == 0 && $nntmp == 1 && $C_nbtmp_2 == 1} {
incr numc(CHN)
# C2H4: 2个C原子，4个H原子，C有3个键合
} elseif {$nctmp == 2 && $nhtmp == 4 && $notmp == 0 && $nntmp == 0 && $C_nbtmp_3 == 2} {
incr numc(C2H4)
}
# 清理临时变量
unset nhtmp nctmp notmp nntmp H_nbtmp_1 O_nbtmp_2 C_nbtmp_2 C_nbtmp_3 C_nbtmp_4 N_nbtmp_3
}
set frame_real [expr {$frame + $startframe - 1}]
puts -nonewline $myfile [expr $timestep * $frame_real]\
foreach molecule $molecules {
puts -nonewline $myfile $numc($molecule)\
}
puts $myfile " "
unset nfrag molecules numc ifrag frame_real
}
# Call bigdcd command to manage memory and call your analysis procedure for each frame
proc bigdcd { script type args } {
global bigdcd_frame bigdcd_proc bigdcd_lastframe vmd_frame bigdcd_running startframe endframe
set bigdcd_running 1
set bigdcd_frame 0
set bigdcd_lastframe [molinfo top get numframes]
set bigdcd_proc $script
# backwards "compatibility". type flag is omitted.
if {[file exists $type]} {
set args [linsert $args 0 $type]
set type auto
}
uplevel #0 trace variable vmd_frame w bigdcd_callback
foreach dcd $args {
if { $type == "auto" } {
mol addfile $dcd first $startframe last $endframe waitfor 0
} else {
mol addfile $dcd type $type first $startframe last $endframe waitfor 0
}
}
after idle bigdcd_wait
}
proc bigdcd_callback { tracedvar mol op } {
global bigdcd_frame bigdcd_proc bigdcd_lastframe vmd_frame
set msg {}
set thisframe $vmd_frame($mol)
incr bigdcd_frame
if { [catch {uplevel #0 $bigdcd_proc $bigdcd_frame} msg] } {
puts stderr "bigdcd aborting at frame $bigdcd_frame\n$msg"
bigdcd_done
return
}
animate delete beg $thisframe end $thisframe $mol
return $msg
}
proc bigdcd_done { } {
global bigdcd_running
if {$bigdcd_running > 0} {
uplevel #0 trace vdelete vmd_frame w bigdcd_callback
puts "End of frame"
puts "Bigdcd_done"
set bigdcd_running 0
}
}
proc bigdcd_wait { } {
global bigdcd_running bigdcd_frame bigdcd_lastframe
while {$bigdcd_running > 0} {
global bigdcd_oldframe
set bigdcd_oldframe $bigdcd_frame
# run global processing hooks (including loading of scheduled frames)
display update ui
# Stop when reaching maximum frame
if { $bigdcd_frame >= $bigdcd_lastframe } {
bigdcd_done
}
}
}
# Set up periodic boundary condition
set pbc_parameters {11.117 11.117 11.117 90.0000 90.0000 90.0000}
# Set up timestep (fs)
set timestep 0.5
# Set up frame index (start from 0)
set startframe 0
set endframe 1604
# Set up file path
set inputfile "C:/Users/Lenovo/Desktop/1/qingniaosuan-pos-1.xyz"
set outputfile "C:/Users/Lenovo/Desktop/1/qingniaosuan-pos-1.txt"
# Check if input file exists
if {![file exists $inputfile]} {
puts "Error: Input file not found: $inputfile"
return
}
# Job starts
set myfile [open $outputfile a]
set mol [mol new $inputfile first $startframe last $endframe type xyz waitfor all]
set all [atomselect top all]
# Call bigdcd command to manage memory and call your analysis procedure for each frame
bigdcd my_analysis auto $inputfile
# Wait for all frames to be processed
bigdcd_wait
# Close the result file
close $myfile

复制代码