I'm currently using the development version with TCL interface and found Langevin thermostat gives me the wrong temperature.
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!
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"