#########################################################################################
#
# Cylindrical cell model : MD mimic Deserno et al Macromolecules 2000,33, 199-206
#
# Generate integrated charge distribution around the rod
#
#########################################################################################

puts " "
puts "======================================================="
puts "=                                                     ="
puts "=    MD simulation of a charged rod                   ="
puts "=                                                     ="
puts "======================================================="
puts " "
puts "[code_info]"


####################
# PARAMETERS
####################

# System parameters
set L 25.0;				# box_length
set r0 1.0;				# rod radius
set bjerrum_length 1.0;
set lambda 1.0;				# line charge density of the rod 
set ion_diameter 1.0;			# ion diameter
set valency_ci 1.0;			# valency of the counterions
set num_rod_beads 50;			# number of beads of the rod

# Output parameters
set energy_file "rod-energy.dat"
set dist_file "rod-dist.dat"
set vtf_file "rod.vtf"
set checkpoint_filename "rod_%.0f.chk"

# Run parameters
set measure_distribution 0;		# measure the distribution?
set max_frames 100;			# MD frames to go
set steps_per_frame 100;		# number of timesteps per frame
set timestep 0.1;			# timestep in simulation units
set skin 0.4;				# Verlet skin
set accuracy 1e-3;			# accuracy of p3m algorithm

# P(r) distribution 
set r_min 1.0;				# starting distance for P(r)
set r_max [expr $L/2];			# end distance for P(r)
set r_bins 100;				# Number of bins
set log_flag 1;				# Logarithmic value if 1
set int_flag 1;				# Integrated distribution if 1

####################
# CONSTANTS
####################
# particle types 
set rod_id 1
set ci_id 2

####################
# DERIVED PARAMETERS
####################
# total charge of the rod
set total_rod_charge [expr $lambda*$L]
# charge per rod bead
set rod_charge [expr -$total_rod_charge/$num_rod_beads]
# distance of beads on the rod
set rod_distance [expr $L/$num_rod_beads]
# number of counterions
set num_ci [expr $total_rod_charge/$valency_ci]
    
####################
# SETTING UP THE SYSTEM
####################

# GLOBAL PARAMETERS
setmd box_l $L $L $L
setmd periodic 1 1 1
setmd time_step $timestep
setmd skin $skin
thermostat langevin 1.0 0.5

#
# SETTING UP THE INTERACTIONS
# 
# ion-ion interaction
set lj0_eps 1.0
set lj0_sig $ion_diameter
set lj0_cut [expr pow(2,1/6)*$ion_diameter]
set lj0_off 0.0
set lj0_shift 0.25
inter $ci_id $ci_id lennard-jones $lj0_eps $lj0_sig $lj0_cut $lj0_shift $lj0_off

# ion-rod interaction
set lj1_eps 1.0
set lj1_sig $ion_diameter
set lj1_cut [expr pow(2,1/6)*$ion_diameter]
set lj1_off [expr $r0-$ion_diameter]
set lj1_shift 0.25
puts "Offset for the ion-rod interaction = $lj1_off"
inter $ci_id $rod_id lennard-jones $lj1_eps $lj1_sig $lj1_cut $lj1_shift $lj1_off

if {[llength $argv] == 1} then {
    # start from checkpoint
    set checkpointfile [lindex $argv 0]

    puts "Starting from checkpoint $checkpointfile..."
    set checkpoint [open "$checkpointfile" "r"]
    while { [blockfile $checkpoint read auto] != "eof" } {}
    close $checkpoint

    # open files for appending data
    set efile [open $energy_file "a"]
    set vtffile [open $vtf_file "a"]
} else {
    # CREATE ROD
    # set rod along the center of the box
    set px [expr $L/2]
    set py [expr $L/2]
    set pz 0.0
    for {set i 0} {$i < $num_rod_beads} {incr i} {
	part $i \
	    pos $px $py $pz \
	    type $rod_id \
	    q $rod_charge \
	    fix 1 1 1
	set pz [expr $pz+$rod_distance]
    }
    
    #
    # CREATE COUNTERIONS
    # 
    set start_pid [expr [setmd max_part] + 1]
    for { set i 0 } { $i < $num_ci } { incr i } {
	part [expr $start_pid + $i] \
	    pos [expr [t_random]*$L] [expr [t_random]*$L] [expr [t_random]*$L] \
	    q $valency_ci type $ci_id
    }
    
    ####################
    # WARMING UP
    ####################
    puts "Warming up...\n"
    set min_dist 0.90
    set act_min_dist [analyze mindist]
    set i 0
    set cap 1.0
    inter ljforcecap $cap
    for { set i 0 } { $i < 4000 && $act_min_dist < $min_dist } { incr i } {
	integrate 200
	# Warmup criterion
	set act_min_dist [analyze mindist [list $ci_id] [list $ci_id]]
	# Increase LJ cap
	set cap [expr $cap+1.0]
	inter ljforcecap $cap
	puts "Warming step: $i min_dist=$act_min_dist cap=$cap\r"
    }
    inter ljforcecap 0.0
    puts "Warming up finished!"

    setmd time 0

    # init data files
    set efile [open $energy_file "w"]
    puts $efile "\#step\ttotal\tkin\tcoulomb\tlj"

    set vtffile [open $vtf_file "w"]
    writevsf $vtffile
}

####################
# SETTING UP THE COULOMB INTERACTION
####################
puts "Tuning P3M..."
inter coulomb $bjerrum_length p3m tunev2 accuracy $accuracy
set params [inter coulomb]
puts "Tuning p3m has finished: $params"

puts "[thermostat]"
puts "[inter]"

####################
# STARTING SIMULATION
####################
puts "Starting simulation..."
if { $measure_distribution } then {
    # create a list of length $r_bins
    for {set i 0} { $i < $r_bins } { incr i } {
	lappend p_ci 0
    }
    set dist {}
}


set frame_count 0
set i 0
for { set i 0 } { $i < $max_frames } { incr i } {
    puts -nonewline "frame: $i/$max_frames\r"
    flush stdout
    integrate $steps_per_frame
    puts $efile [join [list \
			   [setmd time] \
			   [analyze energy total] \
			   [analyze energy kinetic] \
			   [analyze energy coulomb] \
			  ] "\t"]
    flush $efile

    # off-line visualisation
    writevcf $vtffile 

    # measure distribution function
    if {$measure_distribution} then {
	set dist \
	    [analyze distribution \
		 $ci_id $rod_id \
		 $r_min $r_max $r_bins \
		 $log_flag $int_flag]
	# add result to list
	set j 0
	foreach data [lindex $dist 1] old_p $p_ci {
	    lset p_ci $j [expr $old_p + [lindex $data 1]]
	    incr j
	}
	incr frame_count
    }
}

if { $measure_distribution } then {
    # now get the radii
    set rlist {}
    set j 0
    foreach data [lindex $dist 1] old_p $p_ci {
	lappend rlist [lindex $data 0]
	incr j
    }
}
puts "\nSimulation finished!"

####################
# ANALYZE
####################
# compute the average distribution
if { $measure_distribution } then {
    set p_ci [vecscale [expr 1.0/$frame_count] $p_ci] 
    set distfile [open $dist_file "w"]
    puts $distfile "#"
    puts $distfile "# Integrated charge distribution around the rod"
    puts $distfile "#r\tP(r)"
    foreach r $rlist value $p_ci { puts $distfile "$r\t$value" }
    
    puts "Wrote integrated charge to $dist_file."
    close $distfile
}

puts "Wrote energies to $energy_file."
puts "Wrote trajectory to $vtf_file."
close $efile
close $vtffile

####################
# WRITE CHECKPOINT
####################

set cp_file [format $checkpoint_filename [setmd time]]
puts "Writing checkpoint to $cp_file..."
set cpfile [open $cp_file "w"]
blockfile $cpfile write particles "id type pos q v f fix" all
blockfile $cpfile write variable {time}
close $cpfile
puts "Checkpoint written."

exit