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

# System identification:
set name BPAPC
set ident "wall" 
set vmd_output "no"
set tcl_precision  6

set C_PI  3.1415926535897932384626434
set two16 [expr pow(2,(1.0/6.0))]
set three16 [expr pow(3,(1.0/6.0))]
set cshift126 [expr 1.0/4.0]
set cshift93  [expr 2.0/3.0/sqrt(3.0)]
set roffC  0.05
set roffI -0.25
# in between roffC and roffI - only for equilibration
set roffP -0.05 

set verb 50

puts "Restarting from previous calculation\n"
set bfile "|gzip -cd [lindex $argv 0].block.gz"
puts "Reading data from $bfile"
set inp [open $bfile "r"]
while {[blockfile $inp read auto] != "eof"} {}
close $inp

if { [file exists "[lindex $argv 0].chk" ] } {
  puts -nonewline "    Checkpoint found (currently reading it...\n"; flush stdout
  checkpoint_read "[lindex $argv 0]"
  puts " N_P=$N_P BPC=$BPC\n"
  puts -nonewline "    checkpoint read"; flush stdout
} else {
  puts -nonewline "Cannot find checkpoint file [lindex $argv 0].chk\n"; flush stdout
  exit 1
}

set loop [lindex $argv 0]
puts "From checkpoint file:\n"
puts "Particles:\n [part]\n"
puts "Exclusions must be redefined!\n"
puts "Constraints must be redefined!\n [constraint]\n"
puts "setmd:\n [setmd]\n"
puts "integrate:\n  [integrate]\n"
puts "thermostat:\n [thermostat]\n"
set output [open "restart.vtf" "w"]
writevsf $output
#puts $output "atom 0:[expr $nsites-1] resid 0"
for { set nid 0 } { $nid < [expr $xcells*$ycells] } { incr nid } {
    set sitetop1 [expr $nid*10]
    set sitebot1 [expr $nid*10+1]
    puts $output "atom $sitetop1:$sitebot1 resid 1"
    set sitetop2 [expr $nid*10+2]
    set sitebot2 [expr $nid*10+3]
    puts $output "atom $sitetop2:$sitebot2 resid 2"
    set sitetop3 [expr $nid*10+4]
    set sitebot7 [expr $nid*10+9]
    puts $output "atom $sitetop3:$sitebot7 resid 3"
}
for { set cid 0 } { $cid < $N_P } { incr cid } {
    set first_bead [expr $cid * $BPC +$nsites]
    set last_bead [expr $first_bead + $BPC - 1]
    puts $output "atom $first_bead:$last_bead resid [expr $cid+4]"
}
writevcf $output
close $output
writepdb  "restart.pdb"
writepsf  "restart.psf"

# Interaction & Particle setup
#############################################################

# Interaction parameters
#############################################################

constraint wall normal 0.0 0.0  1.0 dist $lowerz type $typeW
constraint wall normal 0.0 0.0 -1.0 dist -$upperz type $typeW
puts "Redefined constraints:\n [constraint]\n"

# EXCLUSIONS
for { set chain 0 } { $chain < $N_P } { incr chain} {
  for { set bead 0 } { $bead < [expr $BPC-1] } { incr bead} {
    set at1 [expr $chain * $BPC + $bead + $nsites]
    set at2 [expr $at1 + 1]
    part $at1 exclude $at2
  }
}
puts "1-2 exclusions added."

for { set chain 0 } { $chain < $N_P } { incr chain} {
  for { set bead 0 } { $bead < [expr $BPC-5] } { set bead [expr $bead+4]} {
    set at1 [expr $chain * $BPC + $bead + $nsites]
    set atP2 [expr $at1 + 2]
    set atP1next [expr $at1 + 4]
    part $atP2 exclude $atP1next
  }
}
puts "1-3 exclusions added.\n"

# SWITCH ON SURFACE SITE POTENTIALS
  puts "\nInteractions:\n[inter]\n"
  puts "Attraction switched on!!!"
  inter $angleIPC angledist [expr 0.1*pow(180.0/$C_PI, 2)] [expr 120.0*$C_PI/180.0] [expr 2.2/$sigma] $C_PI [expr 3.1/$sigma]
  set anglePCW 5
  inter $anglePCW endangledist [expr 0.1*pow(180.0/$C_PI, 2)] [expr 60.0*$C_PI/180.0] [expr 2.2/$sigma] [expr 3.1/$sigma]

  for { set chain 0 } { $chain < $N_P } { incr chain} {
    for { set bead 0 } { $bead < [expr $BPC-1] } { set bead [expr $bead+4]} {
      set atP1 [expr $chain * $BPC + $bead + $nsites]
      set atC [expr $atP1 + 1]
      set atP2 [expr $atP1 + 2]
      if { $bead == 0 } {
#       START P1--C and wall
        part $atP1 bond $anglePCW $atC
      }
      if { $bead == [expr $BPC-3]} {
#       PCW:  END C--P2 and wall
        part $atP2 bond $anglePCW $atC
      }
    }
  }

#### WALL-BEAD INTERACTION ###
#  inter $typeC $typeW lj-gen 1 $C_sigma [expr $three16*$C_sigma] $cshift93 $roffC 9 3 1 1
#  inter $typeI $typeW lj-gen 1 $I_sigma [expr $three16*$I_sigma] $cshift93 $roffI 9 3 1 1
#  inter $typeP $typeW lj-gen 0 $I_sigma [expr $three16*$I_sigma] $cshift93 $roffI 9 3 1 1
  inter  $typeI $typeW tabulated "IW-ljgen.DAT"
  inter  $typeC $typeW tabulated "CW-ljgen.DAT"
  inter  $typeP $typeW tabulated "PW-0.DAT"
  inter $typeP $typeS1 tabulated "PW-ljcos.DAT"
  inter $typeP $typeS2 tabulated "cosine_r1.DAT"
  inter $typeP $typeS3 tabulated "cosine_r4.DAT"

  puts "Particles:\n [part]\n"
  puts "\nInteractions:\n[inter]\n"
  set output [open "restart2.vtf" "w"]
  writevsf $output
#  puts $output "atom 0:[expr $nsites-1] resid 0"
  for { set nid 0 } { $nid < [expr $xcells*$ycells] } { incr nid } {
    set sitetop1 [expr $nid*10]
    set sitebot1 [expr $nid*10+1]
    puts $output "atom $sitetop1:$sitebot1 resid 1"
    set sitetop2 [expr $nid*10+2]
    set sitebot2 [expr $nid*10+3]
    puts $output "atom $sitetop2:$sitebot2 resid 2"
    set sitetop3 [expr $nid*10+4]
    set sitebot7 [expr $nid*10+9]
    puts $output "atom $sitetop3:$sitebot7 resid 3"
  }
  for { set cid 0 } { $cid < $N_P } { incr cid } {
    set first_bead [expr $cid * $BPC]
    set last_bead [expr $first_bead + $BPC - 1]
      puts $output "atom $first_bead:$last_bead resid [expr $cid+4]"
  }
  writevcf $output
  close $output
  writepdb  "restart2.pdb"
  writepsf  "restart2.psf"

#############################################################
#  Setup System                                             #
#############################################################

    setmd time_step 0.001
    setmd skin      0.4
    thermostat set langevin 1.0 0.5
    analyze set chains 0 $N_P $BPC
    puts "Parameters of system now:\n"
    puts "   (setmd):      [setmd]\n"
    puts "   (integrate):  [integrate]\n"
    puts "   (thermostat): [thermostat]\n"
    if { $vmd_output=="yes" } {
      prepare_vmd_connection_special "$name$ident" 
    }

#############################################################
# Integration 
#############################################################

puts "    Integration loop starting... "; flush stdout;
setmd time 0
set time1 [clock seconds]

puts "prevcalc set to 1. Restarting integration."
set int_loop  30
set int_step  1000000
puts "Starting from previous calculation at $loop.\n"
puts -nonewline "I:   Step  time   Temp    Etot     mindist     re     bondlength \n"; flush stdout
set tmp_Temp [expr [analyze energy kin]/$n_part/1.5]
puts -nonewline "I: $counter [setmd time]  "
puts -nonewline "$tmp_Temp   [analyze energy total]  "
puts -nonewline "[analyze mindist $typeC $typeI]   [lindex [analyze re] 0] "
puts -nonewline "[lindex [analyze bond_l] 0] \n"; flush stdout
integrate [expr 0]
writepdb "$loop-2.pdb"
polyBlockWriteAll "$loop-2.block.gz";
checkpoint_set "$loop-2.check.gz" [expr 1] "all" "lowerz upperz" "-"
set output [open "$loop-2.vtf" "w"]
writevsf $output
#puts $output "atom 0:[expr $nsites-1] resid 0"
for { set nid 0 } { $nid < [expr $xcells*$ycells] } { incr nid } {
    set sitetop1 [expr $nid*10]
    set sitebot1 [expr $nid*10+1]
    puts $output "atom $sitetop1:$sitebot1 resid 1"
    set sitetop2 [expr $nid*10+2]
    set sitebot2 [expr $nid*10+3]
    puts $output "atom $sitetop2:$sitebot2 resid 2"
    set sitetop3 [expr $nid*10+4]
    set sitebot7 [expr $nid*10+9]
    puts $output "atom $sitetop3:$sitebot7 resid 3"
}
for { set cid 0 } { $cid < $N_P } { incr cid } {
    set first_bead [expr $cid * $BPC + $nsites]
    set last_bead [expr $first_bead + $BPC - 1]
    puts $output "atom $first_bead:$last_bead resid [expr $cid+4]"
}
writevcf $output
close $output

set counter $loop
for { set j 1 } { $j <= $int_loop } { incr j } {
  integrate [expr $int_step]
  set counter [expr $counter+$int_step]
  set tmp_Temp [expr [analyze energy kin]/$n_part/1.5]
  puts -nonewline "$j: $counter [setmd time]  "
  puts -nonewline "$tmp_Temp   [analyze energy total]  "
  puts -nonewline "[analyze mindist $typeC $typeI]   [lindex [analyze re] 0] "
  puts -nonewline "[lindex [analyze bond_l] 0] \n"; flush stdout
  polyBlockWriteAll "$counter.block.gz";
  checkpoint_set "$counter.check.gz" [expr 1] "all" "lowerz upperz" "-"
  writepdb "$counter.pdb"
  set output [open "$counter.vtf" "w"]
  writevsf $output
#  puts $output "atom 0:[expr $nsites-1] resid 0"
  for { set nid 0 } { $nid < [expr $xcells*$ycells] } { incr nid } {
      set sitetop1 [expr $nid*10]
      set sitebot1 [expr $nid*10+1]
      puts $output "atom $sitetop1:$sitebot1 resid 1"
      set sitetop2 [expr $nid*10+2]
      set sitebot2 [expr $nid*10+3]
      puts $output "atom $sitetop2:$sitebot2 resid 2"
      set sitetop3 [expr $nid*10+4]
      set sitebot7 [expr $nid*10+9]
      puts $output "atom $sitetop3:$sitebot7 resid 3"
  }
  for { set cid 0 } { $cid < $N_P } { incr cid } {
      set first_bead [expr $cid * $BPC + $nsites]
      set last_bead [expr $first_bead + $BPC - 1]
      puts $output "atom $first_bead:$last_bead resid [expr $cid+4]"
  }
  writevcf $output
  close $output
}


set time2 [clock seconds]
puts "I: [expr $int_loop * $int_step] steps took [expr $time2 - $time1] seconds."
puts "I: => [expr ($time2 - $time1) / $int_loop / $int_step] seconds per step."

puts "\nFinished with current system\n"; flush stdout