replicaBarrier
set nr [numReplicas]
if { $num_replicas != $nr } {
error "restart with wrong number of replicas"
}
set r [myReplica]
set replica_id $r
if {[info exists restart_root]} { #restart
set restart_root [format $restart_root $replica_id]
source $restart_root.$replica_id.tcl
} else {
set i_job 0
set i_run 0
set i_step 0
if {[info exists first_timestep]} {
set i_step $first_timestep
}
set replica(index) $r
set replica(ParamID) $r
set replica(loc.a) $r
set replica(index.a) $r
set replica(loc.b) $r
set replica(index.b) $r
set replica(exchanges_attempted) 0
set replica(exchanges_accepted) 0
if { $r % 2 == 0 && $r+1 < $nr } {
if { $drude_model != 0 } {
if { $r < $blockerb } {
set replica(loc.a) [expr $r+1]
set replica(index.a) [expr $r+1]
}
} else {
set replica(loc.a) [expr $r+1]
set replica(index.a) [expr $r+1]
}
}
if { $r % 2 == 1 && $r > 0 } {
if { $drude_model != 0 } {
if { $r <= $blockerb } {
set replica(loc.a) [expr $r-1]
set replica(index.a) [expr $r-1]
}
} else {
set replica(loc.a) [expr $r-1]
set replica(index.a) [expr $r-1]
}
}
if { $r % 2 == 1 && $r+1 < $nr } {
if { $drude_model != 0 } {
if { $r < $blockerb } {
set replica(loc.b) [expr $r+1]
set replica(index.b) [expr $r+1]
}
} else {
set replica(loc.b) [expr $r+1]
set replica(index.b) [expr $r+1]
}
}
if { $r % 2 == 0 && $r > 0 } {
if { $drude_model != 0 } {
if { $r <= $blockerb } {
set replica(loc.b) [expr $r-1]
set replica(index.b) [expr $r-1]
}
} else {
set replica(loc.b) [expr $r-1]
set replica(index.b) [expr $r-1]
}
}
}
array set lambda_values { 0 0.0000 1 0.045 2 0.090 3 0.14546 4 0.22425 5 0.30303 6 0.38182 7 0.46061 8 0.5394 9 0.61819 10 0.6970 11 0.77576 12 0.85455 13 0.910 14 0.955 15 1.0000 }
set job_output_root "$output_root.job$i_job"
firsttimestep $i_step
proc replica_lambda { i } {
global lambda_values num_replicas num_replicasa num_replicasb num_replicasc
if { $i < $num_replicas } {
return [ expr $lambda_values($i) ]
} else {
return 1.0
}
}
proc replica_sptscale { i } {
global num_replicas min_temp max_temp
set half_replicas [expr $num_replicas/2]
if { $i < $half_replicas } {
set temp [expr ($min_temp * exp( log(1.0*$max_temp/$min_temp)*(1.0*$i/($half_replicas-1)) ) )]
return [ expr $min_temp/$temp ]
} else {
set temp [expr ($max_temp * exp( log(1.0*$min_temp/$max_temp)*(1.0*($i-$half_replicas+1)/($half_replicas)) ) )]
return [ expr $min_temp/$temp ]
}
}
proc setup_parameters { ID } {
global num_replicas restart_root
set IDN [expr ($ID + 1)]
if { $IDN < $num_replicas } {
set Lambda [replica_lambda $ID]
set Lambda2 [replica_lambda $IDN]
alchLambda $Lambda
alchLambda2 $Lambda2
} else {
alchLambda 1.0
alchLambda2 1.0
}
}
proc save_callback {labels values} {
global saved_labels saved_values
set saved_labels $labels
set saved_values $values
}
callback save_callback
proc save_array {} {
global saved_labels
global saved_values
global saved_array
foreach label $saved_labels value $saved_values {
set saved_array($label) $value
}
}
seed [expr int(0*srand(int(100000*rand()) + 100*$replica_id) + 100000*rand())]
outputname [format $job_output_root.$replica_id $replica_id]
if {$i_run} { #restart
bincoordinates $restart_root.$replica_id.coor
binvelocities $restart_root.$replica_id.vel
extendedSystem $restart_root.$replica_id.xsc
} else {
binCoordinates equ_site_6.coor
extendedSystem equ_site_6.xsc
temperature $TEMP
}
outputEnergies [expr $steps_per_run]
dcdFreq [expr $steps_per_run * $runs_per_frame]
source $namd_config_file
alchOutFile [format $job_output_root.$replica_id.wham.fepout $replica_id]
setup_parameters $replica(ParamID)
set history_file [open [format "$job_output_root.$replica_id.history" $replica_id] "w"]
fconfigure $history_file -buffering line
set history_bonds_file [open [format "$job_output_root.$replica_id.bonds_history" $replica_id] "w"]
fconfigure $history_bonds_file -buffering line
while {$i_run < $num_runs} {
if { $i_run == 0 } {
minimize 200
}
run $steps_per_run
save_array
incr i_step $steps_per_run
set POTENTIAL [expr $saved_array(TOTAL) - $saved_array(KINETIC)]
set BONDS [expr $saved_array(BOND) + $saved_array(ANGLE) + $saved_array(DIHED) + $saved_array(IMPRP)]
set doswap 0
if { $i_run % 2 == 0 } {
set swap a; set other b
} else {
set swap b; set other a
}
if { $replica(index) < $replica(index.$swap) } {
set POTENTIAL2 [replicaRecv $replica(loc.$swap)]
}
if { $replica(index) > $replica(index.$swap) } {
replicaSend $POTENTIAL $replica(loc.$swap)
}
if { $replica(index) > $replica(index.$swap) } {
set POTENTIAL2 [replicaRecv $replica(loc.$swap)]
}
if { $replica(index) < $replica(index.$swap) } {
replicaSend $POTENTIAL $replica(loc.$swap)
}
if { $replica(index) != $replica(index.$swap) } {
set replica(ParamID) $replica(index.$swap)
setup_parameters $replica(ParamID)
run 0
save_array
set POTENTIAL_NEW [expr $saved_array(TOTAL) - $saved_array(KINETIC)]
set BONDS_NEW [expr $saved_array(BOND) + $saved_array(ANGLE) + $saved_array(DIHED) + $saved_array(IMPRP)]
if { $replica(index) < $replica(index.$swap) } {
set POTENTIAL_NEW2 [replicaRecv $replica(loc.$swap)]
}
if { $replica(index) > $replica(index.$swap) } {
replicaSend $POTENTIAL_NEW $replica(loc.$swap)
}
if { $replica(index) > $replica(index.$swap) } {
set POTENTIAL_NEW2 [replicaRecv $replica(loc.$swap)]
}
if { $replica(index) < $replica(index.$swap) } {
replicaSend $POTENTIAL_NEW $replica(loc.$swap)
}
if { $replica(index) < $replica(index.$swap) } {
set BOLTZMAN 0.001987191
set delta [expr ($POTENTIAL_NEW + $POTENTIAL_NEW2 - $POTENTIAL - $POTENTIAL2)/($BOLTZMAN * $TEMP)]
set doswap [expr $delta < 0. || exp(-1. * $delta) > rand()]
replicaSend $doswap $replica(loc.$swap)
puts $history_file "$i_step $replica(index) $replica(index.$swap) $TEMP $POTENTIAL $POTENTIAL_NEW $POTENTIAL2 $POTENTIAL_NEW2 $doswap"
puts $history_bonds_file "$i_step $replica(index) $replica(index.$swap) $TEMP $BONDS $BONDS_NEW $doswap"
if { $doswap } {
puts stderr "EXCHANGE_ACCEPT $replica(index) $replica(index.$swap) RUN $i_run"
incr replica(exchanges_accepted)
}
incr replica(exchanges_attempted)
}
if { $replica(index) > $replica(index.$swap) } {
set doswap [replicaRecv $replica(loc.$swap)]
puts $history_file "$i_step $replica(index) $replica(index.$swap) $TEMP $POTENTIAL $POTENTIAL_NEW $POTENTIAL2 $POTENTIAL_NEW2 $doswap"
puts $history_bonds_file "$i_step $replica(index) $replica(index.$swap) $TEMP $BONDS $BONDS_NEW $doswap"
}
}
set newloc $r
if { $doswap } {
set newloc $replica(loc.$swap)
set replica(loc.$swap) $r
} else {
set replica(ParamID) $replica(index)
setup_parameters $replica(ParamID)
}
set replica(loc.$other) [replicaSendrecv $newloc $replica(loc.$other) $replica(loc.$other)]
if { $doswap } {
array set replica [replicaSendrecv [array get replica] $newloc $newloc]
if { [expr $num_runs - $i_run] == 1 } {
set replica(ParamID) $replica(index)
}
}
incr i_run
if { $i_run % ($runs_per_frame * $frames_per_restart) == 0 ||
$i_run == $num_runs } { # restart
set restart_root "$job_output_root.restart$i_run"
output [format $restart_root.$replica_id $replica_id]
set rfile [open [format "$restart_root.$replica_id.tcl" $replica_id] "w"]
puts $rfile [list array set replica [array get replica]]
close $rfile
replicaBarrier
if { $replica_id == 0 } {
set rfile [open [format "$restart_root.tcl" ""] "w"]
puts $rfile [list set i_job [expr $i_job + 1]]
puts $rfile [list set i_run $i_run]
puts $rfile [list set i_step $i_step]
puts $rfile [list set restart_root $restart_root]
close $rfile
if [info exists old_restart_root] {
set oldroot [format $old_restart_root ""]
file delete $oldroot.tcl
}
}
replicaBarrier
if [info exists old_restart_root] {
set oldroot [format $old_restart_root $replica_id]
file delete $oldroot.$replica_id.tcl
file delete $oldroot.$replica_id.coor
file delete $oldroot.$replica_id.vel
file delete $oldroot.$replica_id.xsc
}
set old_restart_root $restart_root
}
}
set attempts $replica(exchanges_attempted)
if $attempts {
set i $replica(index)
if { $replica(index.a) > $i } {
set swap a
} else {
set swap b
}
set accepts $replica(exchanges_accepted)
set ratio [expr 1.0*$accepts/$attempts]
puts stderr "EXCHANGE_RATIO $replica(index) $replica(index.$swap) $accepts $attempts $ratio"
}
replicaBarrier