set init_seed [clock clicks]
append seeds $init_seed
t_random seed [expr {int(rand()*10000)}] [expr {int(rand()*10000)}] [expr {int(rand()*10000)}]
puts "$init_seed"
puts " "
puts "============================================="
puts "==========ion in electric field=============="
puts "============================================="
puts " "
puts "program Information:\n[code_info]\n"
###MD parameters
setmd time_step 0.01
setmd skin 0.4
set temperature 1.0
set gamma 15.34
#interaction parameters

set lj_sigma 1.0 
set d_ion 1.0
set ion_wall_epsilon 1.0
set ion_wall_sigma 1.0
set ion_wall_cutoff [expr $ion_wall_sigma*pow(2.0,1./6)]
set ion_wall_shift [calc_lj_shift $ion_wall_sigma $ion_wall_cutoff]
set ion_wall_offset 0.0
set init_num_step 7000
####pore properties
set x_axis 1.0
set y_axis 0.0
set z_axis 0.0
set rad [expr $lj_sigma+$d_ion/2]
set leng [expr 5*$lj_sigma]

#########set the box size
set box_n 15.0
###box size
set box_lx [expr $box_n]
set box_ly $box_n
set box_lz $box_n
puts "box_lx=$box_lx,box_ly=$box_ly,box_lz=$box_lz"
setmd box_l $box_lx $box_ly $box_lz
##interaction
inter 1 2 lennard-jones $ion_wall_epsilon $ion_wall_sigma  $ion_wall_cutoff $ion_wall_shift $ion_wall_offset
##initial condition
set x0 [expr $box_lx/2]
puts "x0=$x0"
set y0 [expr $box_ly/2]
puts "y0=$y0"
set z0 [expr $box_lz/2]
puts "z0=$z0"
part 0 pos $x0 $y0 $z0 type 1 
#constraint
constraint pore center [expr $box_lx/2] [expr $box_ly/2] [expr $box_ly/2] axis 1 0 0 radius $rad length [expr $leng/2] type 2 

####open file for output
set posfile [open "pos.dat" "w"]
set chekfile [open "check.dat" "w"]
####open vtf file
set pos_vtf_file [open "pos.vtf" "w"]
set pos_file [open "pos.dat" "w"]
set energy_file [open "energy.dat" "w"]
#####seting up the thermostat
thermostat langevin $temperature $gamma
##warmup integration
#puts "don,t need warmup"
writevsf $pos_vtf_file
writevcf $pos_vtf_file
flush $pos_vtf_file
  
for { set i 0 } { $i < $init_num_step } { incr i } {
     set x [lindex [part 0 print pos] 0] 
    set y [lindex [part 0 print pos] 1] 
    set z [lindex [part 0 print pos] 2]
part 0 ext_force +3 0 0
puts $chekfile "run=$i,mindist=[constraint mindist_position $x $y $z],e_nonbond=[analyze energy nonbonded 1 2]"
puts "[constraint mindist_position $x $y $z]"
    integrate 1
      writevcf $pos_vtf_file
      flush $pos_vtf_file
puts $pos_file "$i,time=[setmd time],pos=[part 0 prin pos]"
puts $energy_file "$i,time=[setmd time], e_nonbond= [analyze energy nonbonded 1 2],e_tot=[analyze energy total],e_kin=[analyze energy kinetic]"
puts "run$i,[ part 0 print folded_pos ]"
}