[ESPResSo-users] Temp

[Muhammad Anwar]

Hello,
Dear Users/Developers,
I am sending you my script, the problem is i am unable to get desired 
temperature. I checked the script with FENE bonding and i got the 
desired temperature but not getting with rigid bonds. Please have a look 
on it and guide me how to sort this problem out. I am using ESPResSo 3.0.1.
Thanks a lot and have a nice day.

Best Regards,
Muhammad Anwar

#!/bin/sh
# tricking... the line after a these comments are interpreted as standard shell 
script \
    exec $ESPRESSO_SOURCE/Espresso $0 $*

puts "[code_info]"

#############################################################
#  System Parameter                                          #
#############################################################

set l_poly  44
set n_poly 45
set T 450
set lj1_epsilon [expr 0.09344/($T * 0.00198721)]    ; # Epsilon in KCal/mol 
divided by KbT in KCal/mol
set lj2_epsilon [expr 0.1455/($T * 0.00198721)]     ; # Epsilon in KCal/mol 
divided by KbT in KCal/mol
set lj3_epsilon [expr 0.2264/($T * 0.00198721)]     ; # Epsilon in KCal/mol 
divided by KbT in KCal/mol

set k1 [expr 120/($T * 0.00198721)]         ; # Stiffness/Bend in KCal/mol 
divided by KbT in KCal/mol
set k2 [expr 1.6/($T * 0.00198721)]         ; # Stiffness/Bend in KCal/mol 
divided by KbT in KCal/mol

set sigma 4.5
set lj1_sig 1
set lj1_cut [expr 2*$lj1_sig]
set lj1_shift [expr - 0.031]
set n_part  [expr $l_poly*$n_poly]
set box_length [expr 39.8886/$sigma]
set bond_l [expr 1.53/$sigma]
set shield 0.1

setmd box_l $box_length $box_length $box_length
setmd time_step 0.0001
setmd skin 0.4

# Thermostat
integrate set nvt
set kT 1.0
thermostat langevin $kT 1


#############################################################
#  Interaction Setup                                          #
#############################################################

inter 1 rigid_bond $bond_l 1e-3 1e-3 
inter 2 angle $k1 [expr [PI]*0.61] ; # Cossquare Bond Angle Potential
inter 3 DIHEDRAL 1 $k2 1           ; # Old Dihedral Potential


inter 0 0 lj-gen $lj1_epsilon $lj1_sig $lj1_cut $lj1_shift 0 12 6 1 2 
inter 0 1 lj-gen $lj2_epsilon $lj1_sig $lj1_cut $lj1_shift 0 12 6 1 2 
inter 1 1 lj-gen $lj3_epsilon $lj1_sig $lj1_cut $lj1_shift 0 12 6 1 2 

puts "Interactions:\n   [inter]"

#############################################################
#  Polymer & Particle setup                                    #
#############################################################

polymer $n_poly $l_poly $bond_l mode SAW $shield bond 1

for { set i 0 } { $i < $n_part } { incr i } {
if {$i % 44 == 0} {part $i type 1} elseif {$i % 44 == 43} {part $i type 1} else 
{part $i type 0}
}


part auto_exclusions 1

for {set i 0} { $i < $n_part -2} {incr i} {
if {$i % 44 != 42 && $i % 44 != 43} {part [expr $i+1] bond 2 [expr $i] [expr 
$i+2]}
}

for {set i 0} { $i < $n_part - 3} {incr i} {
if {$i % 44 != 41 && $i % 44 != 42 && $i % 44 != 43} {part [expr $i+1] bond 3 
[expr $i] [expr $i+2] [expr $i+3]}
}


set struct [open "|gzip -c - >phy.gz" w]
blockfile $struct write variable box_l
blockfile $struct write particles "id pos type " all
blockfile $struct write bonds all
close $struct


#############################################################
#  VMD Connection                                          #
#############################################################

set vmd "no"
if { $vmd == "yes" } {
prepare_vmd_connection tutorial 3000
exec sleep 4
imd positions
}

#############################################################
#  Warm up Integration                                             #
#############################################################

set min 0
set cap 500
while { $min < 0.32} {
    inter ljforcecap $cap
    integrate 10000
    set min [analyze mindist]
    incr cap 500
}
puts "Warmup finished. Minimal distance now $min"


#############################################################
#  Main Integration                                             #
#############################################################

set int_loop 3000
set int_step 1000

set obs [open "rg.dat" "w"]
setmd time 0; inter ljforcecap 0 
puts "Done."

set i 0 
while { $i<$int_loop } {
set Temp [expr [analyze energy kinetic]/$n_part/([degrees_of_freedom]/2.0)]
set p1 [analyze pressure total]
analyze set chains 0 $n_poly $l_poly

set rg [lindex [analyze rg] 0]
set re [lindex [analyze re] 0]
set Cn [expr $re * 4.5*$re * 4.5/(43 * 1.53 * 1.53)] ;# Characteristic Ratio
puts $obs "[setmd time] $rg   $Cn"
puts -nonewline  "t= [expr [setmd time]] P = $p1 Rg = $rg  T=$Temp  C= 
$Cn....\r";  flush stdout
integrate $int_step
if { $vmd == "yes" } { imd positions }
incr i
}
close $obs

puts "\nFinished"
exit