#
# Copyright (C) 2010,2012 The ESPResSo project
# Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 
#   Max-Planck-Institute for Polymer Research, Theory Group
#  
# Simulates a simple Lennard-Jones liquid and calculates the
# radial distribution function.
#  
# This script is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#  
# ESPResSo is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#  
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>. 
# 

# set to yes to visualize your simulation on the fly
set vmd_output "yes"

# number of particles and density
set n_part  200
set density   0.7

# periodic box of given density
set box_l [expr pow($n_part/$density,1./3.)]
setmd box_l $box_l $box_l $box_l
setmd periodic 1 1 1

# add particles at random positions
for {set i 0} {$i < $n_part} {incr i} {
    set posx [expr $box_l*[t_random]]
    set posy [expr $box_l*[t_random]]
    set posz [expr $box_l*[t_random]]
    
    part $i pos $posx $posy $posz
}

# integrator settings
setmd time_step 0.01
# how many steps to perform per cycle 
set integ_steps 100
# Verlet list skin, a tuning parameter
setmd skin 0.3

# configure thermostat
set temperature 1.0
set gamma       1.0
thermostat langevin $temperature $gamma

# interaction parameters
# change this to activate your new potential!
set sig 1.0
set eps 1.0
set cut [expr 1.12246*$sig]
inter 0 0 lennard-jones $eps $sig $cut auto 0

# prepare vmd connection
if { $vmd_output=="yes" } {
    prepare_vmd_connection "vmd"
    imd listen 10000
}

# warmup with capped forces
set maxcap 500
for {set cap 20} {$cap <= $maxcap} {incr cap 20} {
    puts "eq cap = $cap/$maxcap t=[format %#.4g [setmd time]] E=[analyze energy total]"
    # if switched on, write positions to vmd
    if {$vmd_output=="yes"} {imd positions}

    inter forcecap $cap
    integrate $integ_steps
}
inter forcecap 0

# we hope here it is more or less equilibrated...

# initialize the list for the average rdf
for {set i 0} {$i < 100} {incr i} { lappend sum_rdf 0 }

# main integration loop for sampling
set samples 1000
for {set cnt 1} {$cnt <= $samples} { incr cnt} {
    puts "step $cnt of $samples: t=[format %#.4g [setmd time]] E=[analyze energy total]"
    # if switched on, write positions to vmd
    if { $vmd_output=="yes" } {imd positions}

    integrate $integ_steps

    # calculate rdf
    set rdf [analyze rdf 0 0 0.9 [expr $box_l/2] 100]

    # unzip the r,rdf-pairs for averaging
    set rlist ""
    set rdflist ""
    foreach value [lindex $rdf 1] {
	lappend rlist   [lindex $value 0]
	lappend rdflist [lindex $value 1]
    }
    
    # average the rdf
    set sum_rdf [vecadd $sum_rdf $rdflist]
}

# normalize the rdf by the number of samples
set avg_rdf [vecscale [expr 1.0/$samples] $sum_rdf]

# write averaged rdfs to file rdf.data
set plot [open "liquid_rdf.data" "w"]
puts $plot "\# r rdf(r)"
foreach r $rlist rdf $avg_rdf {puts $plot "[format %#.4g $r] $rdf"}
close $plot

exit 0