## Initialization of random generator
set init_seed [clock clicks]
append seeds $init_seed
t_random seed [expr {int(rand()*10000)}] [expr {int(rand()*10000)}] [expr {int(rand()*10000)}] [expr {int(rand()*10000)}]
puts "$init_seed"

puts "Program Information: \n[code_info]\n"

## MD parameters
setmd time_step 0.01
### these tune the speed
#setmd max_num_cells 1000
setmd skin 0.4

set grid_size 1.0

# particle parameters
set num_polymers 1
set num_monomers_per_polymer 45
puts "number of monomers is $num_monomers_per_polymer"

# interaction parameters
set lj_sigma 1.0
set lj_epsilon 1.0
set bjerrum [expr 0.7 * $lj_sigma]

set q_mon 3.0
set q_ion 1.0

set d_mon $lj_sigma
set d_ion [expr 0.425*$lj_sigma]

set monomer_momomer_epsilon $lj_epsilon
set monomer_momomer_sigma $d_mon
set monomer_momomer_cutoff [expr $monomer_momomer_sigma * pow(2.0,1./6.)]
set monomer_momomer_shift [calc_lj_shift $monomer_momomer_sigma $monomer_momomer_cutoff]
set monomer_momomer_offset 0.0

set monomer_ion_epsilon $lj_epsilon
set monomer_ion_sigma [expr ($d_mon+$d_ion)/2.0]
set monomer_ion_cutoff [expr $monomer_ion_sigma * pow(2.0,1./6.)]
set monomer_ion_shift [calc_lj_shift $monomer_ion_sigma $monomer_ion_cutoff]
set monomer_ion_offset 0.0

set ion_ion_epsilon $lj_epsilon
set ion_ion_sigma $d_ion
set ion_ion_cutoff [expr $ion_ion_sigma * pow(2.0,1./6.)]
set ion_ion_shift [calc_lj_shift $ion_ion_sigma $ion_ion_cutoff]
set ion_ion_offset 0.0

set bond_K 30.
set bond_rcut [expr 1.2 * $lj_sigma]
set temperature 1.0

set gamma 15.34
## Information of steps
set init_num_step [expr int(1500*sqrt($num_monomers_per_polymer/10))]
set init_internal_step 500

## Set the box size

set box_l [expr $lj_sigma * int(pow($num_monomers_per_polymer / 10. * pow(15.0,3.),1./3.))]
set Rc [expr 1.7*$lj_sigma*sqrt($num_monomers_per_polymer)*10/12.*9/(0.7*5*$lj_sigma)]
## Box size
set box_lx [expr int(($box_l/2. + 4*$lj_sigma + $Rc + 5 + $box_l/2.)/$grid_size)*$grid_size]
set box_ly [expr int($box_l/$grid_size)*$grid_size]
set box_lz [expr int($box_l/$grid_size)*$grid_size]

setmd box_l $box_lx $box_ly $box_lz

## position for centermass of polymer
set x_0 	[expr $box_l/2. + 4*$lj_sigma + $Rc + 5/2.]
set y_0	[expr $box_ly / 2.]
set z_0	[expr $box_lz / 2.]
##Boundary conditions for set up conterion
set wall1 [expr $x_0 - 5/2 - $lj_sigma]
set wall2 [expr $x_0 + 5/2 + $lj_sigma]

## SET UP THE INTERACTIONS

inter 0 fene $bond_K $bond_rcut
inter 0 0 lennard-jones $monomer_momomer_epsilon $monomer_momomer_sigma $monomer_momomer_cutoff $monomer_momomer_shift $monomer_momomer_offset
inter 0 1 lennard-jones $monomer_ion_epsilon $monomer_ion_sigma $monomer_ion_cutoff $monomer_ion_shift $monomer_ion_offset
inter 1 1 lennard-jones $ion_ion_epsilon $ion_ion_sigma $ion_ion_cutoff $ion_ion_shift $ion_ion_offset

set vtf_output "yes"

  #############################################################
  # SET UP THE SYSTEM

  setmd time 0.0

  set x0 [expr $x_0 - 5/2.0 - $Rc]
  set y0 [expr $y_0]
  set z0 [expr $z_0]

  polymer $num_polymers $num_monomers_per_polymer $d_mon pos $x0 $y0 $z0 type 0 bond 0

  set px [lindex [part [expr $num_monomers_per_polymer/2] print pos] 0]
  set py [lindex [part [expr $num_monomers_per_polymer/2] print pos] 1]
  set pz [lindex [part [expr $num_monomers_per_polymer/2] print pos] 2]

  set distx [expr $x0-$px]
  set disty [expr $y0-$py]
  set distz [expr $z0-$pz]

  for { set i 0 } { $i < [setmd n_part]} { incr i } {
    set posx [lindex [part $i print pos] 0]
    set posy [lindex [part $i print pos] 1]
    set posz [lindex [part $i print pos] 2]
    set typ [part $i print type] 
    part $i pos [expr $posx + $distx] [expr $posy + $disty] [expr $posz + $distz] type $typ q 0
    puts "$i [part $i print pos type q]"
  }
 
    for { set i 0 } { $i < $num_monomers_per_polymer } {incr i} {
      part $i q [expr -$q_mon]
  }
    set num_counterion [expr int($num_polymers*$num_monomers_per_polymer*$q_mon/$q_ion)]
    set min_init_dist [analyze mindist]
    for {set i $num_monomers_per_polymer} { $i < [expr $num_counterion + $num_monomers_per_polymer] } {incr i} {
      set flag 0
      while {$flag != 1} {
        set posx [expr $wall1*[t_random]]
        set posy [expr $box_ly*[t_random]]
        set posz [expr $box_lz*[t_random]]
        set flag 1
        set dist 1000
        part $i pos $posx $posy $posz type 1 q $q_ion
        set dist [analyze mindist]
        if ($dist<$min_init_dist) {
          set flag 0
        }
      }
      puts "$i [part $i print pos type q] [setmd n_part] [analyze mindist] $dist"
    }
  ##############################################################

 
  # open VTF file
  set pre_vtf_file [open "pos.vtf" "w"]
  writevsf $pre_vtf_file
  writevcf $pre_vtf_file
  flush $pre_vtf_file

  #############################################################
  # SETTING UP THE THERMOSTAT


  thermostat langevin $temperature $gamma
   
  
  part [expr $num_monomers_per_polymer/2] fix

  #############################################################
  # WARMUP INTEGRATION
  puts "\nWarmup:"

  # after warmup, all particle pairs should have at least this distance
  set min_dist [expr 0.85 * $lj_sigma]

  # compute the minimal distance between two particles
  set act_min_dist [analyze mindist]
  puts "$act_min_dist"

  # set LJ capping
  set cap 20
  inter ljforcecap $cap

  for {set i 0} { $act_min_dist < $min_dist } {incr i} {
      integrate 100

      # Warmup criterion
      set act_min_dist [analyze mindist]
      puts "\tstep=$i, cap=$cap, act_min_dist = $act_min_dist"

      # Increase LJ cap
      set cap [expr $cap+10]
      inter ljforcecap $cap
  }

  # turn off capping
  inter ljforcecap 0

  #############################################################
  # TURNING ON ELECTROSTATIC FORCE

    puts [inter coulomb $bjerrum p3m tune accuracy 1e-3]
 

  #############################################################
  # EQUILLIBRATION INTEGRATION
set start [clock clicks -milliseconds]
  for { set i 0 } { $i < $init_num_step } { incr i } {

    integrate $init_internal_step

    if { $i%50==0 } {
      writevcf $pre_vtf_file
      flush $pre_vtf_file
    }
   puts "$i"
    }
    
  set last [clock clicks -milliseconds]
 set time_taken [expr $last-$start]
puts "$time_taken"