############################################################# # # # Kremer-Grest's Linear Polymer Melts # # # ############################################################# puts " " puts "===================================================" puts "= Simulating a Kremer-Grest polymer =" puts "===================================================" puts " " puts "Program Information: \n[code_info]\n" ############################################################# # Parameters # ############################################################# # Start of all files written by this script set name "KremerGrest" # System parameters ############################################################# set N_P 16 ;# number of polymer chains set MPC 50 ;# number of monomers in the chain set bond_l 0.97 ;# initial polymer bond length set density 0.85 # Interaction parameters ############################################################# # repulsive Lennard Jones set lj1_eps 1.0 set lj1_sig 1.0 set lj1_cut 1.12246 set lj1_shift [calc_lj_shift $lj1_sig $lj1_cut] set lj1_off 0.0 # attractive FENE bond set fene_k 30.0 set fene_r 1.5 # Integration parameters ############################################################# setmd time_step 0.006 # tuning parameter for frequency of Verlet rebuilds setmd skin 0.4 # coupling to heat bath via dissipation-fluctutation theory thermostat langevin 1.0 0.5 # warmup integration (with capped LJ potential) until particles are at least $min_dist apart # see below for what the values do set warm_step 200 set warm_loop 300 set warm_cap1 10 set warm_incr 10 set min_dist 0.90 # integration (with full LJ potential) for $int_time set int_step 2000 set int_time 10000 # Other parameters ############################################################# set tcl_precision 6 ############################################################# # Setup System # ############################################################# set n_part [expr $N_P*$MPC] set box_l [expr pow($n_part/$density, 1./3.)]; setmd box_l $box_l $box_l $box_l # Interaction setup ############################################################# inter 0 0 lennard-jones $lj1_eps $lj1_sig $lj1_cut $lj1_shift $lj1_off inter 0 FENE $fene_k $fene_r puts "\nSimulate $N_P polymer chains with $MPC monomers of bond-length $bond_l each" puts " in a cubic simulation box of length [setmd box_l] at density $density with gamma [setmd gamma] and temp [setmd temp]." puts "Interactions:\n [inter]" # Particle setup ############################################################# puts -nonewline "Creating LJ-polymers... "; flush stdout polymer $N_P $MPC $bond_l mode SAW 0.2 puts "Done." # Alternatively set up Lennard-Jones particles ############################################################# # #for {set i 0} {$i < $n_part} {incr i} { # # generate here initial random x,y,z-positions of an LJ particle # for random generator use, for example, rand() # # part $i pos $pos_x $pos_y $pos_z type 0 #} # setup topology for polymer analysis ############################################################# analyze set chains 0 $N_P $MPC # Prepare output of VTF file (for VMD) ############################################################# set vtf_file [open "$name.vtf" w] writevsf $vtf_file writevcf $vtf_file ############################################################# # Warm-up Integration (with capped LJ-potential) # ############################################################# puts -nonewline "\nStart warm-up integration (capped LJ-interactions) for maximal [expr $warm_step*$warm_loop] timesteps in $warm_loop loops; " puts "stop if minimal distance is larger than $min_dist." puts " Analysis at t=[setmd time]: mindist=[analyze mindist], re=[lindex [analyze re] 0], rg= [lindex [analyze rg] 0], T=[setmd temp]." setmd time 0 set cap $warm_cap1 set obs_file [open "$name-obs1.dat" "w"] puts $obs_file "# t mindist re rg Temp" puts $obs_file "[setmd time] [analyze mindist] [lindex [analyze re] 0] [lindex [analyze rg] 0] [setmd temp]" for { set j 0 } { $j < $warm_loop } { incr j } { inter ljforcecap $cap integrate $warm_step; set dist [analyze mindist] set temp [expr [analyze energy kin]/$n_part/([degrees_of_freedom]/2.0)] writevcf $vtf_file puts -nonewline " Step [expr ($j+1)*$warm_step]/[expr $warm_step*$warm_loop] (t=[setmd time]): " puts -nonewline "LJ's cap = $cap, Temp = $temp, mindist=[analyze mindist], re=[lindex [analyze re] 0], rg=[lindex [analyze rg] 0]...\r" flush stdout puts $obs_file "[setmd time] [analyze mindist] [lindex [analyze re] 0] [lindex [analyze rg] 0] $temp" if { $dist >= $min_dist } { break } set cap [expr $cap + $warm_incr] } close $obs_file ############################################################# # Integration # ############################################################# set int_loop [expr int($int_time/([setmd time_step]*$int_step))] setmd time 0 set step 0 #putting all averages to zero set re_av 0 set rg_av 0 set p_av 0 puts "\nStart integration (full interactions) with timestep [setmd time_step] until time t>=$int_time (-> $int_loop loops); " puts " Remove capping of LJ-interactions." inter ljforcecap 0 set obs_file [open "$name-obs2.dat" "w"] puts $obs_file "# t mindist re rg Temp pressure p" flush $obs_file for { set j 0 } { $j < $int_loop } { incr j } { integrate $int_step set step [expr ($j+1)*$int_step] puts -nonewline " Step $step/[expr $int_step*$int_loop] (t=[setmd time]): " flush stdout set dist [analyze mindist] set re [lindex [analyze re] 0] set rg [lindex [analyze rg] 0] set ptot [analyze pressure total] set re_av [expr $re_av+$re*$re] set rg_av [expr $rg_av+$rg*$rg] set p_av [expr $p_av+$ptot] set temp [expr [analyze energy kin]/$n_part/([degrees_of_freedom]/2.0)] puts $obs_file "[setmd time] $dist $re $rg $temp $ptot" flush $obs_file puts -nonewline "mindist=$dist, T=$temp, p=$ptot...\r" flush stdout writevcf $vtf_file analyze append } close $obs_file close $vtf_file puts "\nFinished with simulation:" # derive ensemble averages set re_av [expr $re_av/$int_loop] set rg_av [expr $rg_av/$int_loop] set p_av [expr $p_av/$int_loop] set rat2 [expr $re_av/$rg_av] puts " = $re_av" puts " = $rg_av" puts "/ = $rat2 (RW=6)" puts "

= $p_av" # print and sort , , and into .g123-file set g123_file [open "$name-g123.dat" "w"] for {set gx 1} {$gx<=3} {incr gx} { eval set g$gx [list [analyze ]] } for {set gt 0} {$gt<[llength $g1]} {incr gt} { puts $g123_file "[expr $gt*[setmd time_step]] [lindex $g1 $gt] [lindex $g2 $gt] [lindex $g3 $gt]" } close $g123_file puts "Done." exit