VMD脚本自动判断成键并添加键 - 计算化学公社

# 格式: 符号 -> {共价半径电负性是否金属是否半金属}
array set element_data {
Lp {0.30 0 0 0}
H {0.31 2.20 0 0}
He {0.28 0 0 0}
Li {1.28 0.98 1 0}
Be {0.96 1.57 1 0}
B {0.84 2.04 0 1}
C {0.76 2.55 0 0}
N {0.71 3.04 0 0}
O {0.66 3.44 0 0}
F {0.57 3.90 0 0}
Ne {0.58 0 0 0}
Na {1.66 0.93 1 0}
Mg {1.41 1.31 1 0}
Al {1.21 1.61 1 0}
Si {1.11 1.90 0 1}
P {1.07 2.19 0 0}
S {1.05 2.58 0 0}
Cl {1.02 3.16 0 0}
Ar {1.06 0 0 0}
K {2.03 0.82 1 0}
Ca {1.76 1.00 1 0}
Sc {1.70 1.36 1 0}
Ti {1.60 1.54 1 0}
V {1.53 1.63 1 0}
Cr {1.39 1.66 1 0}
Mn {1.39 1.55 1 0}
Fe {1.32 1.83 1 0}
Co {1.26 1.88 1 0}
Ni {1.24 1.91 1 0}
Cu {1.32 1.90 1 0}
Zn {1.22 1.65 1 0}
Ga {1.22 1.81 1 0}
Ge {1.20 2.01 0 1}
As {1.19 2.18 0 1}
Se {1.20 2.55 0 0}
Br {1.20 2.96 0 0}
Kr {1.16 0 0 0}
Rb {2.20 0.82 1 0}
Sr {1.95 0.95 1 0}
Y {1.90 1.22 1 0}
Zr {1.75 1.33 1 0}
Nb {1.64 1.60 1 0}
Mo {1.54 2.16 1 0}
Tc {1.47 2.10 1 0}
Ru {1.46 2.20 1 0}
Rh {1.42 2.28 1 0}
Pd {1.39 2.20 1 0}
Ag {1.45 1.93 1 0}
Cd {1.44 1.69 1 0}
In {1.42 1.78 1 0}
Sn {1.39 1.96 1 0}
Sb {1.39 2.05 0 1}
Te {1.38 2.10 0 1}
I {1.39 2.66 0 0}
Xe {1.40 0 0 0}
Cs {2.44 0.79 1 0}
Ba {2.15 0.89 1 0}
La {2.07 1.10 1 0}
Ce {2.04 1.12 1 0}
Pr {2.03 1.13 1 0}
Nd {2.01 1.14 1 0}
Pm {1.99 1.15 1 0}
Sm {1.98 1.17 1 0}
Eu {1.98 1.18 1 0}
Gd {1.96 1.20 1 0}
Tb {1.94 1.21 1 0}
Dy {1.92 1.22 1 0}
Ho {1.92 1.23 1 0}
Er {1.89 1.24 1 0}
Tm {1.90 1.25 1 0}
Yb {1.87 1.26 1 0}
Lu {1.87 1.00 1 0}
Hf {1.75 1.30 1 0}
Ta {1.70 1.50 1 0}
W {1.62 1.70 1 0}
Re {1.51 1.90 1 0}
Os {1.44 2.20 1 0}
Ir {1.41 2.20 1 0}
Pt {1.36 2.20 1 0}
Au {1.36 2.40 1 0}
Hg {1.32 1.90 1 0}
Tl {1.45 1.80 1 0}
Pb {1.46 1.80 1 0}
Bi {1.48 1.90 1 0}
Po {1.40 2.00 1 0}
At {1.50 2.20 1 0}
Rn {1.50 0 0 0}
Fr {2.60 0.70 1 0}
Ra {2.21 0.90 1 0}
Ac {2.15 1.10 1 0}
Th {2.06 1.30 1 0}
Pa {2.00 1.50 1 0}
U {1.96 1.38 1 0}
Np {1.90 1.36 1 0}
Pu {1.87 1.28 1 0}
Am {1.80 1.13 1 0}
Cm {1.69 1.28 1 0}
Bk {1.50 1.30 1 0}
Cf {1.50 1.30 1 0}
Es {1.50 1.30 1 0}
Fm {1.50 1.30 1 0}
Md {1.50 1.30 1 0}
No {1.50 1.30 1 0}
Lr {1.50 1.29 1 0}
Rf {1.60 0 1 0}
Db {1.60 0 1 0}
Sg {1.60 0 1 0}
Bh {1.60 0 1 0}
Hs {1.60 0 1 0}
Mt {1.70 0 1 0}
Ds {1.70 0 1 0}
Rg {1.70 0 1 0}
Cn {1.70 0 1 0}
Nh {1.70 0 1 0}
Fl {1.70 0 1 0}
Mc {1.70 0 1 0}
Lv {1.70 0 1 0}
Ts {1.70 0 1 0}
Og {1.70 0 1 0}
}
# 辅助函数：检查元素是否存在
proc element_exists {sym} {
global element_data
return [info exists element_data($sym)]
}
# 辅助函数：获取元素信息
proc get_element_info {sym} {
global element_data
if {[element_exists $sym]} {
return $element_data($sym)
}
return ""
}
# 从原子名称提取元素符号
proc extract_element_from_name {name} {
# 去除数字和特殊字符
set clean_name [regsub -all {[0-9_\.\+\-]} $name ""]
# 尝试提取两个字符的元素符号
if {[string length $clean_name] >= 2} {
set two_char [string range $clean_name 0 1]
# 首字母大写，第二个字母小写（标准元素符号格式）
set standard_two [string toupper [string index $two_char 0]][string tolower [string index $two_char 1]]
if {[element_exists $standard_two]} {
return $standard_two
}
}
# 尝试提取一个字符的元素符号
if {[string length $clean_name] >= 1} {
set one_char [string toupper [string index $clean_name 0]]
if {[element_exists $one_char]} {
return $one_char
}
}
# 如果都失败，返回原始名称的前两个字符（大写）
if {[string length $clean_name] >= 2} {
return [string toupper [string range $clean_name 0 1]]
} elseif {[string length $clean_name] >= 1} {
return [string toupper [string index $clean_name 0]]
}
return "UNKNOWN"
}
# 获取原子的元素符号（优先使用element属性，其次从name推断）
proc get_atom_element {atom_name element_attribute} {
# 如果element属性存在且不为空，使用它
if {$element_attribute != "" && [element_exists $element_attribute]} {
return $element_attribute
}
# 否则从name中提取
return [extract_element_from_name $atom_name]
}
# 获取元素的共价半径
proc get_covalent_radius {sym} {
set info [get_element_info $sym]
if {$info ne ""} {
return [lindex $info 0]
}
# 如果元素未知，返回默认值
return 1.2
}
# 获取元素的电负性
proc get_electronegativity {sym} {
set info [get_element_info $sym]
if {$info ne ""} {
return [lindex $info 1]
}
return 0.0
}
# 判断元素是否为金属
proc is_metal {sym} {
set info [get_element_info $sym]
if {$info ne ""} {
return [lindex $info 2]
}
return 0
}
# 判断元素是否为半金属
proc is_metalloid {sym} {
set info [get_element_info $sym]
if {$info ne ""} {
return [lindex $info 3]
}
return 0
}
# 预计算成键判断所需的元素属性
proc precompute_bonding_params {elem_i elem_j tolerance} {
global bonding_cache
# 创建缓存键
set cache_key "${elem_i}_${elem_j}_${tolerance}"
# 检查缓存
if {[info exists bonding_cache($cache_key)]} {
return $bonding_cache($cache_key)
}
# 初始化返回值：{can_bond max_dist_sq needs_eneg_check eneg_threshold}
set result [list 0 0.0 0 0.0]
# 规则1: 金属-金属不成键
if {[is_metal $elem_i] && [is_metal $elem_j]} {
set bonding_cache($cache_key) $result
return $result
}
# 规则3: 金属-H不成键
if {([is_metal $elem_i] && ($elem_j == "H" || $elem_j == "D")) || \
(($elem_i == "H" || $elem_i == "D") && [is_metal $elem_j])} {
set bonding_cache($cache_key) $result
return $result
}
# 计算最大成键距离的平方（避免sqrt计算）
set rad_i [get_covalent_radius $elem_i]
set rad_j [get_covalent_radius $elem_j]
set max_dist [expr {$tolerance * ($rad_i + $rad_j)}]
set max_dist_sq [expr {$max_dist * $max_dist}]
# 规则2: 金属-半金属 (需要电负性差 > 0.3)
if {([is_metal $elem_i] && [is_metalloid $elem_j]) || \
([is_metalloid $elem_i] && [is_metal $elem_j])} {
set eneg_i [get_electronegativity $elem_i]
set eneg_j [get_electronegativity $elem_j]
# 电负性存在且差值大于0.3
if {$eneg_i > 0 && $eneg_j > 0} {
set result [list 1 $max_dist_sq 1 0.3]
}
set bonding_cache($cache_key) $result
return $result
}
# 规则4: 默认规则 - 基于共价半径
set result [list 1 $max_dist_sq 0 0.0]
set bonding_cache($cache_key) $result
return $result
}
# 快速成键判断函数（使用预计算参数和距离平方）
proc fast_bonding_check {dist_sq bonding_params eneg_diff} {
lassign $bonding_params can_bond max_dist_sq needs_eneg_check eneg_threshold
if {!$can_bond} {
return 0
}
if {$dist_sq >= $max_dist_sq} {
return 0
}
if {$needs_eneg_check && $eneg_diff <= $eneg_threshold} {
return 0
}
return 1
}
# 添加键的函数
proc addbond {molid atom1 atom2} {
if {$atom1 == $atom2} {
return
}
if {$atom1 > $atom2} {
set tmp $atom1
set atom1 $atom2
set atom2 $tmp
}
set sel [atomselect $molid "index $atom1 $atom2"]
lassign [$sel getbonds] bond1 bond2
set id [lsearch -exact $bond1 $atom2]
if {$id == -1} {
lappend bond1 $atom2
}
set id [lsearch -exact $bond2 $atom1]
if {$id == -1} {
lappend bond2 $atom1
}
$sel setbonds [list $bond1 $bond2]
$sel delete
}
# 主要函数：自动添加键
proc auto_addbonds {args} {
global bonding_cache
# 清空成键参数缓存
array unset bonding_cache
array set bonding_cache {}
# 默认参数
set molid top
set tolerance 1.1
set cutoff 5.0
set selection ""
set quiet 0
# 解析命令行参数
set arg_count [llength $args]
for {set i 0} {$i < $arg_count} {incr i} {
set arg [lindex $args $i]
switch $arg {
"-molid" {
incr i
if {$i < $arg_count} {
set molid [lindex $args $i]
} else {
puts "错误: -molid 需要参数"
return
}
}
"-tolerance" {
incr i
if {$i < $arg_count} {
set tolerance [lindex $args $i]
} else {
puts "错误: -tolerance 需要参数"
return
}
}
"-cutoff" {
incr i
if {$i < $arg_count} {
set cutoff [lindex $args $i]
} else {
puts "错误: -cutoff 需要参数"
return
}
}
"-sel" {
incr i
if {$i < $arg_count} {
set selection [lindex $args $i]
} else {
puts "错误: -sel 需要参数"
return
}
}
"-quiet" {
set quiet 1
}
default {
puts "错误: 未知参数 '$arg'"
puts "用法: auto_addbonds \[-molid <id>\] \[-tolerance <val>\] \[-cutoff <val>\] \[-sel <selection>\] \[-quiet\]"
return
}
}
}
if {!$quiet} {
puts "开始自动添加键 (分子ID: $molid, 容差: $tolerance, 截断距离: $cutoff Ang)"
if {$selection ne ""} {
puts "原子选择: $selection"
} else {
puts "原子选择: 所有原子"
}
}
# 创建原子选择
if {$selection ne ""} {
set selected_atoms [atomselect $molid $selection]
} else {
set selected_atoms [atomselect $molid all]
}
set num_atoms [$selected_atoms num]
if {$num_atoms == 0} {
if {$selection ne ""} {
puts "错误: 选择语法 '$selection' 没有选中任何原子!"
} else {
puts "错误: 没有找到原子!"
}
$selected_atoms delete
return
}
if {!$quiet} {
puts "选中的原子数: $num_atoms"
}
# 批量获取原子信息
set indices [$selected_atoms get index]
set names [$selected_atoms get name]
set elements [$selected_atoms get element]
set x_coords [$selected_atoms get x]
set y_coords [$selected_atoms get y]
set z_coords [$selected_atoms get z]
# 预处理：确定元素符号和电负性
set atom_elements [list]
set atom_enegs [list]
for {set i 0} {$i < $num_atoms} {incr i} {
set element_attr [lindex $elements $i]
set name_attr [lindex $names $i]
set elem [get_atom_element $name_attr $element_attr]
lappend atom_elements $elem
lappend atom_enegs [get_electronegativity $elem]
}
# 计算截断距离的平方
set cutoff_sq [expr {$cutoff * $cutoff}]
# 统计变量
set bonds_added 0
set processed_pairs 0
# 使用measure contacts进行快速距离筛选
set contact_sel [atomselect $molid $selection]
set contacts [measure contacts $cutoff $contact_sel $contact_sel]
$contact_sel delete
# 提取接触对
set contact_list1 [lindex $contacts 0]
set contact_list2 [lindex $contacts 1]
set num_contacts [llength $contact_list1]
if {!$quiet} {
puts "距离筛选后的原子对数: $num_contacts"
}
# 遍历所有接触对
for {set k 0} {$k < $num_contacts} {incr k} {
set idx_i [lindex $contact_list1 $k]
set idx_j [lindex $contact_list2 $k]
# 跳过自身和重复对
if {$idx_i >= $idx_j} {
continue
}
# 在indices列表中查找位置
set pos_i [lsearch -exact $indices $idx_i]
set pos_j [lsearch -exact $indices $idx_j]
if {$pos_i == -1 || $pos_j == -1} {
continue
}
# 获取元素信息
set elem_i [lindex $atom_elements $pos_i]
set elem_j [lindex $atom_elements $pos_j]
# 预计算成键参数
set bonding_params [precompute_bonding_params $elem_i $elem_j $tolerance]
# 快速检查是否可能成键
lassign $bonding_params can_bond max_dist_sq needs_eneg_check eneg_threshold
if {!$can_bond} {
continue
}
# 计算距离平方（避免sqrt）
set dx [expr {[lindex $x_coords $pos_i] - [lindex $x_coords $pos_j]}]
set dy [expr {[lindex $y_coords $pos_i] - [lindex $y_coords $pos_j]}]
set dz [expr {[lindex $z_coords $pos_i] - [lindex $z_coords $pos_j]}]
set dist_sq [expr {$dx*$dx + $dy*$dy + $dz*$dz}]
incr processed_pairs
# 计算电负性差（如果需要）
set eneg_diff 0.0
if {$needs_eneg_check} {
set eneg_i [lindex $atom_enegs $pos_i]
set eneg_j [lindex $atom_enegs $pos_j]
set eneg_diff [expr {abs($eneg_j - $eneg_i)}]
}
# 快速成键判断
if {[fast_bonding_check $dist_sq $bonding_params $eneg_diff]} {
addbond $molid $idx_i $idx_j
incr bonds_added
if {!$quiet} {
set dist [expr {sqrt($dist_sq)}]
puts " 添加键: $elem_i (原子$idx_i) - $elem_j (原子$idx_j) 距离: [format "%.3f" $dist] Ang"
}
}
}
$selected_atoms delete
if {!$quiet} {
puts "完成!"
puts " 检查原子对: $processed_pairs"
puts " 添加键数: $bonds_added"
} else {
puts "添加了 $bonds_added 个键"
}
# 更新分子显示
mol reanalyze $molid
display update
return $bonds_added
}
# 显示帮助
proc show_bonding_help {} {
puts "自动添加键函数使用说明:"
puts ""
puts "用法:"
puts " auto_addbonds \[options\]"
puts ""
puts "选项:"
puts " -molid <id> 分子ID (默认: top，当前分子)"
puts " -tolerance <val> 容差系数，控制成键距离的宽松程度 (默认: 1.1)"
puts " -cutoff <val> 最大考虑距离，大于此值的原子对不检查 (默认: 5.0 Ang)"
puts " -sel <selection> VMD选择语法，选择要处理的原子 (默认: 所有原子)"
puts " -quiet 安静模式，减少输出信息"
puts ""
puts "示例:"
puts " auto_addbonds"
puts " auto_addbonds -molid 0 -tolerance 1.1 -cutoff 4.0"
puts " auto_addbonds -sel "resname His""
puts " auto_addbonds -sel "resname UiO and within 10 of resname LIG""
puts " auto_addbonds -sel "element C O N" -tolerance 1.2"
puts " auto_addbonds -sel "protein" -cutoff 4.0 -quiet"
}
# 快捷函数
proc addbonds_sel {selection {molid top} {tolerance 1.1} {cutoff 5.0}} {
return [auto_addbonds -molid $molid -tolerance $tolerance -cutoff $cutoff -sel $selection]
}
# 快捷函数 - 安静模式
proc addbonds_sel_quiet {selection {molid top} {tolerance 1.1} {cutoff 5.0}} {
return [auto_addbonds -molid $molid -tolerance $tolerance -cutoff $cutoff -sel $selection -quiet]
}

复制代码