#############################################################
#                                                           #
#  Sample System 1: Lennard Jones Liquid                    #
#                                                           #
#############################################################
puts " "
puts "======================================================="
puts "=          Sample Script 1: lj.tcl                    ="
puts "======================================================="
puts " "

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

require_feature LENNARD_JONES

# Global Tcl variables
# Throughout this script we use LJ reduced units.
set box_size 15
set density 0.316
# try T=0.3, T=1.0 and T=2.0
set temperature 2.0

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

# define the system size
setmd box_l $box_size $box_size $box_size

# set up the integrator time step
setmd time_step 0.01

# the skin has no effect on the result, only on the speed
setmd skin 0.4

# uncomment the following to output the box size
#puts [setmd box_l]

# set up the thermostat
set langevin_gamma 1.0
thermostat langevin $temperature $langevin_gamma

#############################################################
# SET UP THE PARTICLES
set volume [expr $box_size*$box_size*$box_size]
set n_part [expr floor($volume*$density)]

# generate $n_part particles at random positions
for {set i 0} { $i < $n_part } {incr i} {
    set x [expr $box_size*[t_random]]
    set y [expr $box_size*[t_random]]
    set z [expr $box_size*[t_random]]
 
    part $i pos $x $y $z type 0
}

puts "Simulate $n_part particles in a cubic simulation box "
puts "[setmd box_l] at density $density and temperature $temperature"

#############################################################
# SET UP THE INTERACTION
set lj_epsilon 1.0
set lj_sigma 1.0
set lj_cutoff 2.5

inter 0 0 lennard-jones $lj_epsilon $lj_sigma $lj_cutoff

puts "Interactions:\n[inter]"

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

# after warmup, all particle pairs should have at least this distance
set min_dist 0.85

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

# 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

#############################################################
# MAIN INTEGRATION
puts "\nMain integration:"

# open energy file
set dat_file [open "lj-$temperature.dat" "w"]
puts $dat_file "\#t\tE_tot\tE_kin\tE_LJ"

# open VTF file
set vtf_file [open "lj-$temperature.vtf" "w"]
writevsf $vtf_file

for { set i 0 } { $i < 100 } { incr i } {
    integrate 1000

    # write out observables
    set e_tot [analyze energy total]
    set e_kin [analyze energy kinetic]
    set e_lj [analyze energy nonbonded 0 0]

    # output energies
    puts "\tstep=$i, time=[setmd time], e_tot=$e_tot, e_kin=$e_kin, e_lj=$e_lj"
    puts $dat_file "[setmd time]\t$e_tot\t$e_kin\t$e_lj"

    # output configuration
    writevcf $vtf_file

    # store the last 20 configurations
    analyze push 20
}

close $dat_file
close $vtf_file

#############################################################
# ANALYSIS

# Compute the radial distribution function
set rdf_file [open "lj-$temperature.rdf" "w"]
puts $rdf_file "#r\tg"

set rdf [lindex [analyze <rdf> 0 0] 1]
foreach v $rdf {
    foreach {x g} $v {}
    puts $rdf_file "$x\t$g"
}
close $rdf_file

# terminate program
puts "\n\nFinished"