#############################################################
#                                                           #
#            Kremer-Grest's Linear Polymer Melts            #
#                                                           #
#############################################################

puts " "
puts "==================================================="
puts "=     Simulating a Kremer-Grest polymer           ="
puts "==================================================="
puts " "

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

#############################################################
#  Parameters                                               #
#############################################################

# Start of all files written by this script
set name "KremerGrest"

# System parameters
#############################################################

set N_P      16  ;# number of polymer chains
set MPC      50  ;# number of monomers in the chain
set bond_l  0.97 ;# initial polymer bond length
set density 0.85


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

# repulsive Lennard Jones
set lj1_eps     1.0
set lj1_sig     1.0
set lj1_cut     1.12246
set lj1_shift   [calc_lj_shift $lj1_sig $lj1_cut]
set lj1_off     0.0

# attractive FENE bond
set fene_k      30.0
set fene_r       1.5

# Integration parameters
#############################################################

setmd time_step 0.006

# tuning parameter for frequency of Verlet rebuilds
setmd skin      0.4

# coupling to heat bath via dissipation-fluctutation theory
thermostat langevin 1.0 0.5

# warmup integration (with capped LJ potential) until particles are at least $min_dist apart
# see below for what the values do
set warm_step   200
set warm_loop   300
set warm_cap1   10
set warm_incr   10
set min_dist    0.90

# integration (with full LJ potential) for $int_time
set int_step   2000 
set int_time  10000 

# Other parameters
#############################################################

set tcl_precision  6

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

set n_part [expr $N_P*$MPC]
set box_l  [expr pow($n_part/$density, 1./3.)];
setmd box_l $box_l $box_l $box_l

# Interaction setup
#############################################################

inter 0 0 lennard-jones $lj1_eps $lj1_sig $lj1_cut $lj1_shift $lj1_off
inter 0   FENE          $fene_k  $fene_r

puts "\nSimulate $N_P polymer chains with $MPC monomers of bond-length $bond_l each" 
puts "    in a cubic simulation box of length [setmd box_l] at density $density with gamma [setmd gamma] and temp [setmd temp]."
puts "Interactions:\n   [inter]"

# Particle setup
#############################################################

puts -nonewline "Creating LJ-polymers... "; flush stdout

polymer $N_P $MPC $bond_l mode SAW 0.2

puts "Done."

# Alternatively set up Lennard-Jones particles
#############################################################
#
#for {set i 0} {$i < $n_part} {incr i} {
#   
#   generate here initial random x,y,z-positions of an LJ particle
#   for random generator use, for example, rand()
#
#   part $i pos $pos_x $pos_y $pos_z type 0
#}

# setup topology for polymer analysis
#############################################################
analyze set chains 0 $N_P $MPC

# Prepare output of VTF file (for VMD)
#############################################################
set vtf_file [open "$name.vtf" w]
writevsf $vtf_file
writevcf $vtf_file

#############################################################
#  Warm-up Integration (with capped LJ-potential)           #
#############################################################

puts -nonewline "\nStart warm-up integration (capped LJ-interactions) for maximal [expr $warm_step*$warm_loop] timesteps in $warm_loop loops; "
puts "stop if minimal distance is larger than $min_dist."
puts "    Analysis at t=[setmd time]: mindist=[analyze mindist], re=[lindex [analyze re] 0], rg= [lindex [analyze rg] 0], T=[setmd temp]."

setmd time 0
set cap $warm_cap1

set obs_file [open "$name-obs1.dat" "w"]
puts $obs_file "# t mindist re rg Temp"
puts $obs_file "[setmd time] [analyze mindist] [lindex [analyze re] 0] [lindex [analyze rg] 0]  [setmd temp]"

for { set j 0 } { $j < $warm_loop } { incr j } {
    inter ljforcecap $cap

    integrate $warm_step; set dist [analyze mindist]
    set temp [expr [analyze energy kin]/$n_part/([degrees_of_freedom]/2.0)]
    writevcf $vtf_file

    puts -nonewline " Step [expr ($j+1)*$warm_step]/[expr $warm_step*$warm_loop] (t=[setmd time]): " 
    puts -nonewline "LJ's cap = $cap, Temp = $temp, mindist=[analyze mindist], re=[lindex [analyze re] 0], rg=[lindex [analyze rg] 0]...\r"
    flush stdout

    puts $obs_file "[setmd time] [analyze mindist] [lindex [analyze re] 0] [lindex [analyze rg] 0] $temp"

    if { $dist >= $min_dist } { break }

    set cap [expr $cap + $warm_incr]
}

close $obs_file

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

set int_loop [expr int($int_time/([setmd time_step]*$int_step))] 

setmd time 0
set step 0

#putting all averages to zero
set re_av 0
set rg_av 0
set p_av 0

puts "\nStart integration (full interactions) with timestep [setmd time_step] until time t>=$int_time (-> $int_loop loops); "
puts "    Remove capping of LJ-interactions." 
inter ljforcecap 0 

set obs_file [open "$name-obs2.dat" "w"]
puts $obs_file "# t mindist re rg Temp pressure p"
flush $obs_file

for { set j 0 } { $j < $int_loop } { incr j } {
    integrate $int_step
    set step [expr ($j+1)*$int_step]

    puts -nonewline " Step $step/[expr $int_step*$int_loop] (t=[setmd time]): " 
    flush stdout

    set dist [analyze mindist] 
    set re [lindex [analyze re] 0] 
    set rg [lindex [analyze rg] 0]
    set ptot  [analyze pressure total]
    set re_av [expr $re_av+$re*$re]
    set rg_av [expr $rg_av+$rg*$rg]
    set p_av [expr $p_av+$ptot]
    set temp [expr [analyze energy kin]/$n_part/([degrees_of_freedom]/2.0)]

    puts $obs_file "[setmd time] $dist $re $rg $temp $ptot"
    flush $obs_file

    puts -nonewline "mindist=$dist, T=$temp, p=$ptot...\r" 
    flush stdout

    writevcf $vtf_file
    analyze append
    
}
close $obs_file
close $vtf_file

puts "\nFinished with simulation:" 

# derive ensemble averages
set re_av [expr $re_av/$int_loop]
set rg_av [expr $rg_av/$int_loop]
set p_av [expr $p_av/$int_loop]

set rat2 [expr $re_av/$rg_av]

puts "<re^2> = $re_av"
puts "<rg^2> = $rg_av"
puts "<re^2>/<rg^2> = $rat2 (RW=6)"
puts "<p> = $p_av"

# print and sort <g1>, <g2>, and <g3> into .g123-file

set g123_file [open "$name-g123.dat" "w"]
for {set gx 1} {$gx<=3} {incr gx} { eval set g$gx [list [analyze <g$gx>]] }
for {set gt 0} {$gt<[llength $g1]} {incr gt} { 
    puts $g123_file "[expr $gt*[setmd time_step]] [lindex $g1 $gt] [lindex $g2 $gt] [lindex $g3 $gt]"
}
close $g123_file

puts "Done."
exit