#############################################################
#                                                           #
#  Sample System 2: Stretched Polymer                       #
#                                                           #
#############################################################
puts ""
puts "======================================================="
puts "=         Sample script 2: stretched_polymer.tcl      ="
puts "======================================================="
puts "Creating polymer configurations"
puts ""
puts "[code_info]\n"

require_feature LENNARD_JONES
require_feature EXTERNAL_FORCES
require_feature EXCLUSIONS

# particle ids
set monomer_id 0
set pbond_id 0
set self_interaction_id 0

# global parameters
set num_monomers 100
set L_x 15.32
set L_y $L_x
set L_z 20.0

# simulation parameters
set temperature 1.0
set gamma 0.5
set skin 3.0
set timestep 0.005
set num_frames 10000
set steps_per_frame 10
set frames_to_cnf 10

# warmup parameters
set warmup_maxenergy [expr 2.0*$num_monomers]
set warmup_cap_min 10.0
set warmup_cap_step 10.0
set warmup_steps 100

# INTERACTION PARAMETERS
# self-interaction
set self_interaction 1
set self_interaction_exclusion_range 2
set si_sigma 1.0
set si_epsilon 1.0
# bonds
set pbond_k  10.0
set pbond_r0 0.5

# other parameters
set tcl_precision 16
set lj_factor [expr pow(2., 1./6.)]

# create filenames
set vtffilename "stretched_polymer.vtf"
set datfilename "stretched_polymer.dat"

#############################################################
# Initializing the system
#############################################################
puts "Initializing the system..."

# SET UP GLOBAL PARAMETERS
setmd box_l $L_x $L_y $L_z
puts "box_length=[setmd box_l]"
setmd time_step $timestep
setmd skin $skin
#t_random seed [expr abs([clock clicks])]

# SET UP INTERACTIONS
if { $self_interaction } then {
    puts "Creating self-interaction."
    inter $monomer_id $monomer_id \
	lennard-jones $si_epsilon $si_sigma \
	[expr $lj_factor*$si_sigma] [expr 0.25*$si_epsilon] 0

}

inter $pbond_id harmonic $pbond_k $pbond_r0
puts "inter=[inter]"

set lastid [expr $num_monomers-1]

set vtffile [open $vtffilename w]

# SETTING UP PARTICLES AND WARMING UP
puts "Setting up polymer..."

part deleteall

set pid 0
set x [expr 0.5*$L_x]
set y [expr 0.5*$L_y]
set z 0

puts "Creating half a polymer..."
# place the first two beads
part $pid \
    pos $x $y $z \
    type $monomer_id
set last_pid $pid
incr pid

set z [expr $z + $pbond_r0]
part $pid \
    pos $x $y $z \
    type $monomer_id
part $pid bond $pbond_id $last_pid
set last_pid $pid
incr pid

set last_x $x
set last_y $y
set last_z $z

# now make a random walk in xy and walk linear in z
# thus creating half the polymer
set lift_z [expr ($L_z - (4-$num_monomers%2)*$pbond_r0) / ($num_monomers-(4-$num_monomers%2))]
set step_xy [expr sqrt($pbond_r0*$pbond_r0 - $lift_z*$lift_z)]
puts "lift_z=$lift_z step_xy=$step_xy"
for {} { $pid < $num_monomers/2+$num_monomers%2 } { incr pid } {
    set z [expr $last_z + $lift_z]
    set r [expr [t_random] * 2.0 * [PI]]
    set x [expr $last_x + cos($r)*$step_xy]
    set y [expr $last_y + sin($r)*$step_xy]
    part $pid \
	pos $x $y $z \
	type $monomer_id
    # create the bond to the previous bead
    part $pid bond $pbond_id $last_pid
    # save the position of the last particle
    set last_x $x
    set last_y $y
    set last_z $z
    set last_pid $pid
}
puts "Middle bead ($pid): ($x $y $z)"

# mirror the half-polymer
puts "Mirroring half a polymer..."
# put the central bead into the middle
if {$num_monomers % 2 == 0} then { 
    part $pid pos $x $y [expr 0.5*$L_z]
    part $pid bond $pbond_id $last_pid
    set last_pid $pid
    incr pid
}
set mpid [expr $pid - 2 + $num_monomers%2]
for {} { $pid < $num_monomers } { incr pid } {
    set mpos [part $mpid print pos]
    set x [lindex $mpos 0]
    set y [lindex $mpos 1]
    set z [expr $L_z - [lindex $mpos 2]]
    part $pid \
	pos $x $y $z \
	type $monomer_id
    # create the bond to the previous bead
    part $pid bond $pbond_id $last_pid
    set last_pid $pid
    incr mpid -1
}
puts "Last bead ($last_pid): ($x, $y, $z)"
# create the cyclic bond
part 0 bond $pbond_id $lastid

# turn on exclusions
if { $self_interaction } then {
    part auto_exclusions $self_interaction_exclusion_range
}

puts "Polymer is complete."

# write the initial configuration
if { ! $self_interaction } then {
    writevsf $vtffile radius [list 0 [expr 0.5*$pbond_r0]]
} else {
    writevsf $vtffile
}
puts $vtffile "bond 0::$lastid"
writevcf $vtffile

thermostat langevin $gamma $temperature
puts "thermostat=[thermostat]"

if { $self_interaction } then {
    #############################################################
    # Warming up
    #############################################################
    puts "Warming up..."
    
    set wcap $warmup_cap_min
    set e_tot [analyze energy total]
    puts "e_tot=$e_tot"
    while {$e_tot > $warmup_maxenergy} {
	puts -nonewline "forcecap=$wcap"
	inter ljforcecap $wcap
	integrate $warmup_steps
	writevcf $vtffile
	set e_tot [analyze energy total]
	puts " e_tot=$e_tot"
	set wcap [expr $wcap + $warmup_cap_step]
    }
    
    puts "Warmup finished."

    inter ljforcecap 0
}

#############################################################
# Simulating the system
#############################################################

puts "Starting simulation..."

set datfile [open $datfilename w]
puts $datfile "\#t\te_tot\te_kin\te_bond\tp_tot\tp_bond\tbl_mean\tbl_stddev\tbl_max\tbl_min"

# reset the simulation time
setmd time 0.0
for { set frame 1 } { $frame <= $num_frames } { incr frame } {
    puts -nonewline "$frame/$num_frames\r"
    flush stdout
    integrate $steps_per_frame

    set bl [analyze bond_l 0 1 $num_monomers]
    set bl_mean [lindex $bl 0]
    set bl_stddev [lindex $bl 1]
    set bl_max [lindex $bl 2]
    set bl_min [lindex $bl 3]
    set e_bond [expr [analyze energy bonded $pbond_id] / $num_monomers]
    # write observables
    puts $datfile \
	[format "%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f" \
	     [setmd time] \
	     [analyze energy total] \
	     [analyze energy kinetic] \
	     [analyze energy bonded $pbond_id] \
	     [analyze pressure total] \
	     $e_bond \
	     $bl_mean $bl_stddev $bl_max $bl_min
	]
    flush $datfile

    if { $frame % $frames_to_cnf == 0 } then {
	# write VTF
	writevcf $vtffile
	flush $vtffile
    }
}
puts "$frame/$num_frames"

close $vtffile
close $datfile

puts "Simulation complete."