#
# 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 1:1-salt.
#  
# 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"

set ident "freeze_salt"

# number of particles and density
set n_part  100
set density   0.35

# 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
puts "set up periodic box of size [setmd box_l]."

# add particles at random positions
set q    1
set type 0
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 q $q type $type 
}
set q    -0.5
set type 1
for {set i 0} {$i < 2*$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 [expr $n_part + $i] pos $posx $posy $posz q $q type $type 
}
puts "created [setmd n_part] particles."

# 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 temp  1.0
set gamma 1.0
thermostat langevin $temp $gamma
puts "set up thermostat [thermostat]."

# interaction parameters
# short range Lennard-Jones/Weeks-Chandler-Andersen
set sig   1.0
set eps   1.0
set cut   [expr 1.12246*$sig]
inter 0 0 lennard-jones $eps $sig $cut auto 0
inter 1 0 lennard-jones $eps [expr 0.75*$sig] [expr 0.75*$cut] auto 0
inter 1 1 lennard-jones $eps [expr 0.5*$sig] [expr 0.5*$cut] auto 0

######### turn on Coulomb interaction ###########
#
#

set bjerrum [expr 5/$temp]
set method "p3m"

switch $method {
    "memd" {
	# change the cellsystem according to MEMDs needs
	cellsystem domain_decomposition -no_verlet_list

	# calculate MEMD parameters
	# minima and maxima for the mesh. A good idea to prevent
	# inaccuracy or huge, RAM-exceeding grids.
	set min_memd_mesh 16
	set max_memd_mesh 64
	# calculate mesh like in the user guide. rounded up to the
	# next higher 2^n.
	set memd_mesh [ expr int( pow( 2, ceil( ( log10($box_l/$sig)/log10(2) ) ) ) )]
	if {$memd_mesh > $max_memd_mesh} {
	    set memd_mesh $max_memd_mesh
	}
	if {$memd_mesh < $min_memd_mesh} {
	    set memd_mesh $min_memd_mesh
	}
	set f_mass [expr 100.0*pow( ([setmd time_step]*$memd_mesh/$box_l) , 2.0 ) ]
	puts "memd parameters: mesh=$memd_mesh, f_mass=$f_mass"
	puts [inter coulomb $bjerrum maggs $f_mass $memd_mesh]
    }
    "p3m" {
	set mesh 16
	set accuracy 1e-2
	puts "p3m tuning with mesh $mesh and accuracy $accuracy"
	puts [inter coulomb $bjerrum p3m tunev2 accuracy $accuracy mesh $mesh]
    }
}
#
#
#####################################################

# 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

for {set temp 1} {$temp > 0.1} {set temp [expr $temp*0.98]} {
    thermostat langevin $temp $gamma
    puts "set up thermostat [thermostat]."

    set bjerrum [expr 5/$temp]

    switch $method {
	"memd" {
	    puts [inter coulomb $bjerrum maggs $f_mass $memd_mesh]
	}
	"p3m" {
	    puts [inter coulomb $bjerrum p3m tunev2 accuracy $accuracy mesh $mesh]
	}
    }

    puts "simulated annealing 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
}

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

# open energy file
set ene_file [open "${ident}_energy.data" "w"]
puts $ene_file "\#t\tE_tot\tE_kin\tE_ele"

# open vtf file
set vtf_file [open "${ident}.vtf" "w"]
writevsf $vtf_file
writevcf $vtf_file

# main integration loop for sampling
set samples 100
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 
    
    # output detailed energies
    set e_tot [analyze energy total]
    set e_kin [analyze energy kinetic]
    set e_ele [analyze energy coulomb]
    puts $ene_file "[setmd time]\t$e_tot\t$e_kin\t$e_ele"
    flush $ene_file

    # output configuration
    writevcf $vtf_file

    analyze append
}
close $ene_file
close $vtf_file

# compute the radial distribution functions

puts "compute RDFs..."

set rdf_file [open "${ident}_rdf.data" "w"]
puts $rdf_file "#r\tg+-(r)\tg++(r)\tg--(r)"

# The '<>' around <rdf> mean to analyze the stored configurations and
# compute the mean rdf
set rdf01 [lindex [analyze <rdf> 0 1 0 [expr 0.5*$box_l] 100] 1]
set rdf00 [lindex [analyze <rdf> 0 0 0 [expr 0.5*$box_l] 100] 1]
set rdf11 [lindex [analyze <rdf> 1 1 0 [expr 0.5*$box_l] 100] 1]
foreach v01 $rdf01 v00 $rdf00 v11 $rdf11 {
    set x   [lindex $v01 0]
    set g01 [lindex $v01 1]
    set g00 [lindex $v00 1]
    set g11 [lindex $v11 1]
    puts $rdf_file "$x\t$g01\t$g00\t$g11"
}
close $rdf_file

puts "done!"

exit 0