|
From: | RunningQ |
Subject: | Re: [ESPResSo-users] Langevin thermostat can not keep temperature as the set value in the development version |
Date: | Wed, 7 Dec 2016 00:59:05 -0500 |
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,QiMy 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 usedsource functions_dev.tclputs " "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 $dirputs "The directory $dir is created."}############################################################ # # Parameters############################################################ # set n_part 1000set phi 0.05set lambda 2set alpha 5set pi [expr {4*atan(1)}]puts "pi = $pi"############################################################ # # Integration parametersset method "p3m"set skin 3set temp 1; set gammat 10; set gammar 3#thermostat offthermostat langevin $temp $gammat $gammar $gammar $gammarputs "thermostat = [thermostat]"set dipm [expr {pow($lambda,0.5)}]set Hz [expr {$alpha*$temp/$dipm}]constraint ext_magn_field 0 0 $Hzputs "dipm = $dipm"puts "Hz = $Hz"puts [constraint]puts "constraint H is set up.\r"set warm_dt 0.0015set warm_interval 10set warm_iterations 140setmd time_step $warm_dtset sp_dt 0.0015set sampling_interval 15000set sampling_iterations 134setmd skin $skin############################################################ # # System Setup############################################################ # # Particle setupset tcl_precision 16#setting a seed for the random number generatorexpr 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_lengthset 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 $aaset dip_y $bbset dip_z $ccset 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_particlessaveDipoles $dipFile $n_partwritevsf $trajFilewritevsf $sampTrajFileFwritevsf $sampTrajFileAputs $energyFile "#"puts $energyFile "# Total Magnetic Kinetic Temperature"puts $energyFile "# "writevcf $trajFileset wid_start 1set sid_start 1puts "wid_start = $wid_start.\r"puts "sid_start = $sid_start.\r"########################################################### ## # Interaction setupset lj_eps 1.0set lj_sig 1.0set lj_cut [expr {pow(2,1.0/6.0)*$lj_sig}]puts "lj_cut = $lj_cut"set lj_shift 0.25set lj_offset 0.0inter 0 0 lennard-jones $lj_eps $lj_sig $lj_cut $lj_shift $lj_offsetputs "lennard-jones is set up."############################################################ ### # p3m tuningset lb [expr {1.0/$temp}]puts "lb = $lb"if { $method == "p3m" } {# cellsystem domain_decomposition -no_verlet_listsetmd periodic 1 1 1puts [inter magnetic $lb p3m tune accuracy 1.6e-5]} elseif {$method == "dawaanr"} {cellsystem domain_decomposition -no_verlet_listsetmd periodic 1 1 1puts [inter magnetic 1.0 dawaanr]} else {setmd periodic 1 1 1cellsystem 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 1for {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_particlesinter forcecap $rCapintegrate $warm_intervalwritevcf $trajFile foldedsaveDipoles $dipFile $n_partset 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 $simFileset 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 $binFileset finishFile [open "$dir/finished_iter.dat" "w"]puts $finishFile "$wid"puts $finishFile 0close $finishFileputs "Warm-up iteration $wid of $warm_iterations \at LJ rCap=$rCap, min dist = $act_min_dist, temperature = $kinetic_temp\r"flush stdout}inter forcecap 0puts "\nWarm-up finished with minimal distance $act_min_dist and temperature $kinetic_temp\n"
[Prev in Thread] | Current Thread | [Next in Thread] |