# 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 monomer_wall_epsilon $lj_epsilon
set monomer_wall_sigma $d_mon
set monomer_wall_cutoff [expr $monomer_wall_sigma * pow(2.0,1./6.)]
set monomer_wall_shift [calc_lj_shift $monomer_wall_sigma $monomer_wall_cutoff]
set monomer_wall_offset 0.0

set ion_wall_epsilon $lj_epsilon
set ion_wall_sigma [expr ($d_mon+$d_ion)/2.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 monomer_part_epsilon $lj_epsilon
set monomer_part_sigma [expr $d_mon*4]
set monomer_part_cutoff [expr $monomer_part_sigma * pow(2.0,1./6.)]
set monomer_part_shift [calc_lj_shift $monomer_part_sigma $monomer_part_cutoff]
set monomer_part_offset 0.0

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


## Pore properties
set x_axis 	1.0
set y_axis 	0.0
set z_axis 	0.0
set rad		[expr 1.5* $lj_sigma]
set leng 	[expr 5. * $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_momomer_epsilon $monomer_momomer_sigma $monomer_momomer_cutoff $monomer_momomer_shift $monomer_momomer_offset
inter 1 1 lennard-jones $ion_ion_epsilon $ion_ion_sigma $ion_ion_cutoff $ion_ion_shift $ion_ion_offset
inter 0 2 lennard-jones $monomer_momomer_epsilon $monomer_momomer_sigma $monomer_momomer_cutoff $monomer_momomer_shift $monomer_momomer_offset
inter 1 2 lennard-jones $ion_wall_epsilon $ion_wall_sigma $ion_wall_cutoff $ion_wall_shift $ion_wall_offset
inter 0 4 lennard-jones $monomer_part_epsilon $monomer_part_sigma $monomer_part_cutoff $monomer_part_shift $monomer_part_offset
inter 1 4 lennard-jones $monomer_part_epsilon $monomer_part_sigma $monomer_part_cutoff $monomer_part_shift $monomer_part_offset
##set up the system
polymer $num_polymers $num_monomers_per_polymer $d_mon pos $x0 $y0 $z0 mode PSAW type 0 
constraint pore center $x_pore $y_pore $z_pore axis $x_axis $y_axis $z_axis radius $rad length [expr 2] voltage 0 type 2 
  constraint wall dist 0. normal 1.0 0. 0. type 2 
# WARMUP INTEGRATION
  puts "\nWarmup:"

  # after warmup, all particle pairs should have at least this distance
  set min_dist [expr 0.9 * $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+5]
      inter ljforcecap $cap
  }
# TURNING ON ELECTROSTATIC FORCE

  if {$electrostatic == "yes"} {
    # Turning on electrostatics
    puts "\nTuning electrostatics:"
    puts [inter coulomb $bjerrum p3m tune accuracy 1e-3]
    puts "\n"
  }