标题: 想要统计离子与膜之间的静电作用，但是运行下面的代码会发生“vecsub: two vectors... [打印本页] 作者Author: JYT 时间: 2023-6-6 18:53 标题: 想要统计离子与膜之间的静电作用，但是运行下面的代码会发生“vecsub: two vectors... 想要统计离子与膜之间的静电作用，但是运行下面的代码会发生“vecsub: two vectors don't have the same size”报错，请问原因是？
# Load the necessary files
mol load psf XXXXX.psf
mol addfile XXXXXXX.dcd first 1000 last 4000 waitfor all
# Add ions to the system
set ion [atomselect top "resname XXX"]
set num_ions [$ion num]
set wat [atomselect top "resname XXX"]
set num_wats [$wat num]
set sel [atomselect top "all"]
$sel set beta 0.0
set num_atoms [expr {[molinfo top get numatoms] - $num_ions - $num_wats}]
$sel set beta 1.0
# Set up output file for energy vs. time data
set outfile [open "energy_vs_time.dat" w]
puts $outfile "Time (ps) Energy (kJ/mol)"
# Calculate the Coulombic interaction between ions and XXX
set sel_ion [atomselect top "resname XXX and name XXX"]
set sel_cof [atomselect top "segname XXXXX"]
set ion_charge 1.0
set dielectric 80.0
set rij_cutoff 12.0
set k_e 1389.3548
set start_frame 0
set end_frame [molinfo top get numframes]
for {set frame $start_frame} {$frame < $end_frame} {incr frame} {
animate goto $frame
set time [expr {[molinfo top get timesteps] * $frame / 1000.0}]
set total_energy 0.0
for {set i 0} {$i < $num_ions} {incr i} {
set ion_pos [$sel_ion get {x y z}]
set ion_q [expr {$ion_charge * $k_e}]
for {set j 0} {$j < $num_atoms} {incr j} {
set cof_pos [$sel_cof get {x y z}]
set rij [vecdist $ion_pos $cof_pos]
if {$rij < $rij_cutoff} {
set q_cof [$sel_cof get charge]
set q_ion $ion_q
set energy [expr {$k_e * $q_ion * $q_cof / ($dielectric * $rij)}]
set total_energy [expr {$total_energy + $energy}]
}
}
}
puts $outfile "$time $total_energy"
}
close $outfile