#
# 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 .
#
# set to yes to visualize your simulation on the fly
set vmd_output "yes"
set ident "salt"
# 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
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]]
set q [expr -$q]
set type [expr 1-$type]
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 $sig $cut auto 0
inter 1 1 lennard-jones $eps $sig $cut auto 0
######### turn on Coulomb interaction ###########
#
#
set bjerrum 1.0
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
# 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 500
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
# store the last 20 configurations
analyze push 20
}
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)"
# The '<>' around mean to analyze the stored configurations and
# compute the mean rdf
set rdf01 [lindex [analyze 0 1 0 [expr 0.5*$box_l] 100] 1]
foreach v01 $rdf01 {
set x [lindex $v01 0]
set g01 [lindex $v01 1]
puts $rdf_file "$x\t$g01"
}
close $rdf_file
puts "done!"
exit 0