Re: [ESPResSo-users] Langevin thermostat can not keep temperature as the

On Fri, Dec 9, 2016 at 5:57 AM, Georg Rempfer <address@hidden> wrote:

Hello Qi,

I tried to reproduce your problem but failed to compile the current master branch with your feature set (see myconfig.hpp attached).

We don't want to spend much time supporting the TCL version any more as we focus our development methods on the new Espresso version with a Python interface instead of the old TCL interface. If you could rewrite this script using the Python interface and the Espresso python branch (https://github.com/espressomd/espresso/tree/python), we would be very interested in your results.

Greetings,
Georg

On Wed, Dec 7, 2016 at 6:59 AM, RunningQ <address@hidden> wrote:
Hello,

I forgot to include the function file functions_dev.tcl. It has the following 2 functions.

proc save_sim {cfile range } {
# write all available sim information to channel cfile
# in block format

blockfile $cfile write variable all
# blockfile $cfile write tclvariable all
blockfile $cfile write particles "id pos v f quat omega_lab torque_lab" $range
blockfile $cfile write particles "id pos dip type" $range
blockfile $cfile write interactions
# blockfile $cfile write bonds
blockfile $cfile write random
blockfile $cfile write seed
# blockfile $cfile write bitrandom
# blockfile $cfile write bitseed
}

proc saveDipoles {cfile n_part} {
for {set i 0} {$i < $n_part} {incr i} {
puts $cfile [part $i print id dip]
}
puts $cfile " "
}

Regards,
Qi

On Wed, Dec 7, 2016 at 12:24 AM, RunningQ <address@hidden> wrote:
Hello, everyone,

I'm currently using the development version with TCL interface and found Langevin thermostat gives me the wrong temperature.

My simulation is about magnetic dipoles with LJ potential and dipolar p3m. All parameters are non-dimensional. I used the released version before and had no problem. I know Langevin thermostat has got some changes last year, so I actually copied new files from this commit https://github.com/espressomd/espresso/pull/358 to the released version 3.3.1 and recompiled the software. This gives me the correct result.

Now I need to use lees-edwards feature which is only in the development version. However, in the development version, when I set the temperature as 1, the result is fluctuating between 1.7 to 2.1. In some cases it even goes up to 4.

When I turn off the thermostat, the old and new version give the same energies and temperature. I thought the problem might be the thermostat, but development version with thermostat files from commit 358 also gives wrong temperature.

I don't know whether this is a bug or I made some mistake in the setup. If someone can help me with this, I will be very grateful!

Regards,
Qi

My features I turned on are { Compilation status { ENGINE } { CONSTRAINTS } { ROTATIONAL_INERTIA } { COMFIXED } { NPT } { GHMC } { COMFORCE } { GHOST_FLAG } { LENNARD_JONES } { MODES } { EXTERNAL_FORCES } { ROTATION } { LANGEVIN_PER_PARTICLE } { LENNARD_JONES_GENERIC } { DP3M } { ROTATION_PER_PARTICLE } { GHOSTS_HAVE_BONDS } { MOLFORCES } { FFTW } { DPD } { MPI_CORE } { NEMD } { LEES_EDWARDS } { CONFIGTEMP } { FORCE_CORE } { DIPOLES } { COLLISION_DETECTION } { MASS } { EXCLUSIONS } }

The following is my code. ( I omitted the sampling iteration)

################################################################
# Source (call) functions to be used
source functions_dev.tcl

puts " "
puts "======================================================="
puts "= lj_liquid_tutorial.tcl ="
puts "======================================================="
puts " "
puts "Espresso Code Base : \n[code_info]\n"
puts " "
puts " "

########################################################################

set dir "MvsAlpha_dev/"

if { [file exists $dir]==0 } {
exec mkdir $dir
puts "The directory $dir is created."
}

#############################################################
# Parameters
#############################################################

set n_part 1000
set phi 0.05
set lambda 2
set alpha 5

set pi [expr {4*atan(1)}]
puts "pi = $pi"

#############################################################
# Integration parameters

set method "p3m"
set skin 3
set temp 1; set gammat 10; set gammar 3
#thermostat off
thermostat langevin $temp $gammat $gammar $gammar $gammar
puts "thermostat = [thermostat]"

set dipm [expr {pow($lambda,0.5)}]
set Hz [expr {$alpha*$temp/$dipm}]
constraint ext_magn_field 0 0 $Hz
puts "dipm = $dipm"
puts "Hz = $Hz"
puts [constraint]
puts "constraint H is set up.\r"

set warm_dt 0.0015
set warm_interval 10
set warm_iterations 140
setmd time_step $warm_dt

set sp_dt 0.0015
set sampling_interval 15000
set sampling_iterations 134

setmd skin $skin
#############################################################
# System Setup
#############################################################
# Particle setup

set tcl_precision 16
#setting a seed for the random number generator
expr srand([pid])

########################################################################
# initialize particles
#######################################################################

set box_length [expr {pow($n_part/6.0*$pi/$phi,1.0/3.0)}]
puts "box_length = $box_length"
setmd box_l $box_length $box_length $box_length

set partFile [open "$dir/particles.dat" "w"]
set dipFile [open "$dir/dipoles.dat" "w"]
set sampDipFile [open "$dir/sampDipoles.dat" "w"]
set trajFile [open "$dir/traj.vft" "w"]
set sampTrajFileF [open "$dir/sampTrajF.vft" "w"]
set sampTrajFileA [open "$dir/sampTrajA.vft" "w"]
set energyFile [open "$dir/energy.dat" "w"]

for {set i 0} {$i < $n_part} {incr i} {

set aa [expr ($i+1)/100]
set bb [expr ($i+1)/10-$aa*10]
set cc [expr $i+1-$aa*100-$bb*10]
set pos_x [expr {$aa/10.0*$box_length}]
set pos_y [expr {$bb/10.0*$box_length}]
set pos_z [expr {$cc/10.0*$box_length}]
set dip_x $aa
set dip_y $bb
set dip_z $cc
set norm [expr {pow($dip_x*$dip_x+$dip_y*$dip_y+$dip_z*$dip_z,0.5)}]
set dip_x [expr {$dip_x/$norm*$dipm}]
set dip_y [expr {$dip_y/$norm*$dipm}]
set dip_z [expr {$dip_z/$norm*$dipm}]

part $i pos $pos_x $pos_y $pos_z mass 1 rinertia 0.1 0.1 0.1 dip $dip_x $dip_y $dip_z type 0
}
set act_min_dist [analyze mindist]
puts "act_min_dist = $act_min_dist.\n"

for {set i 0} {$i < $n_part} {incr i} {
puts $partFile [part $i print id pos dip]
}
sort_particles

saveDipoles $dipFile $n_part
writevsf $trajFile
writevsf $sampTrajFileF
writevsf $sampTrajFileA
puts $energyFile "#"
puts $energyFile "# Total Magnetic Kinetic Temperature"
puts $energyFile "# "
writevcf $trajFile

set wid_start 1
set sid_start 1
puts "wid_start = $wid_start.\r"
puts "sid_start = $sid_start.\r"

#############################################################
# Interaction setup

set lj_eps 1.0
set lj_sig 1.0
set lj_cut [expr {pow(2,1.0/6.0)*$lj_sig}]
puts "lj_cut = $lj_cut"
set lj_shift 0.25
set lj_offset 0.0

inter 0 0 lennard-jones $lj_eps $lj_sig $lj_cut $lj_shift $lj_offset
puts "lennard-jones is set up."

###############################################################
# p3m tuning

set lb [expr {1.0/$temp}]
puts "lb = $lb"

if { $method == "p3m" } {
# cellsystem domain_decomposition -no_verlet_list
setmd periodic 1 1 1
puts [inter magnetic $lb p3m tune accuracy 1.6e-5]
} elseif {$method == "dawaanr"} {
cellsystem domain_decomposition -no_verlet_list
setmd periodic 1 1 1
puts [inter magnetic 1.0 dawaanr]
} else {
setmd periodic 1 1 1
cellsystem domain_decomposition -no_verlet_list
}
puts [inter magnetic]

if { [regexp "ROTATION" [code_info]] } {
set deg_free 6
} else {
set deg_free 3
}
puts "there are $deg_free degrees of freedom."

#############################################################
# preparation
#############################################################

set act_min_dist [analyze mindist]
set energies [analyze energy]
puts "energies = $energies"
set total [expr {[analyze energy total]/$n_part}]
set mag [expr {[analyze energy magnetic]/$n_part}]
set kinetic [expr {[analyze energy kinetic]/$n_part}]
set kinetic_temp [expr {$kinetic/($deg_free/2.0)}]
puts $energyFile " "
puts $energyFile "[expr {$wid_start-1}] $total $mag $kinetic $kinetic_temp"


#############################################################
# Warmup Integration
#############################################################

puts "\nWarm-up integration started"
puts "Maximum time steps $warm_interval times $warm_iterations"
puts "Start with minimal distance $act_min_dist and temperature $kinetic_temp"

set flag 1
for {set wid $wid_start} {$wid <= $warm_iterations} {incr wid} {

if {$flag == 1} {
# set rCap [expr {$wid*5000}]
set rCap 5000
} else {
set rCap 0
}

sort_particles
inter forcecap $rCap
integrate $warm_interval

writevcf $trajFile folded
saveDipoles $dipFile $n_part

set act_min_dist [analyze mindist]
set energies [analyze energy]
puts ""
puts "energies = $energies"
set total [expr {[analyze energy total]/$n_part}]
set mag [expr {[analyze energy magnetic]/$n_part}]
set kinetic [expr {[analyze energy kinetic]/$n_part}]
set kinetic_temp [expr {$kinetic/($deg_free/2.0)}]
puts $energyFile "$wid $total $mag $kinetic $kinetic_temp"

if {$kinetic_temp < $temp} {
set flag 0
}

set simFile [open [format "$dir/simulation_%i_%i.dat" $wid 0] "w"]
save_sim $simFile "all"
close $simFile

set binFile [open [format "$dir/binary_%i_%i.dat" $wid 0] "w"]
writemd $binFile "posx" "posy" "posz" "vx" "vy" "vz" "fx" "fy" "fz" "mx" "my" "mz"
close $binFile


set finishFile [open "$dir/finished_iter.dat" "w"]
puts $finishFile "$wid"
puts $finishFile 0
close $finishFile

puts "Warm-up iteration $wid of $warm_iterations \
at LJ rCap=$rCap, min dist = $act_min_dist, temperature = $kinetic_temp\r"
flush stdout

}
inter forcecap 0

puts "\nWarm-up finished with minimal distance $act_min_dist and temperature $kinetic_temp\n"

From:	Qi Liu
Subject:	Re: [ESPResSo-users] Langevin thermostat can not keep temperature as the set value in the development version
Date:	Fri, 9 Dec 2016 16:20:01 -0500