#レプリカ交換分子動力学法の外部から温度を参照する版です。
#NAMDのディレクトリにあるreplica.namdを一部変更して作成
#使用する場合は自己責任でお願いします
#また、エネルギー最小化は行っていないので別途スクリプトを用意してください
#平衡化計算をレプリカ交換法の計算の前に行っているため、sortreplicasを行う前に、
#catdcd等を用いて最初の平衡化を行ったフレームのみ削除をするひつようがあります。
#温度のファイルは一列のデータで用意してください


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(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 } {
    set replica(loc.a) [expr $r+1]
    set replica(index.a) [expr $r+1]
  }
  if { $r % 2 == 1 && $r > 0 } {
    set replica(loc.a) [expr $r-1]
    set replica(index.a) [expr $r-1]
  }

  if { $r % 2 == 1 && $r+1 < $nr } {
    set replica(loc.b) [expr $r+1]
    set replica(index.b) [expr $r+1]
  }
  if { $r % 2 == 0 && $r > 0 } {
    set replica(loc.b) [expr $r-1]
    set replica(index.b) [expr $r-1]
  }

}

set job_output_root "$output_root.job$i_job"
firsttimestep $i_step

if [catch {open "tdata.dat" r} fd] {
        error "cannot open the file"
}

set I 0
while {[gets $fd line] >= 0} {
        set temp_data($I) [lindex $line 0]
        incr I
}
close $fd

proc replica_temp { i } {
        upvar $t temp_data
        return $temp_data($i)
}

set replica(temperature) [replica_temp temp_data $replica(index)]
set replica(temperature.a) [replica_temp temp_data $replica(index.a)]
set replica(temperature.b) [replica_temp temp_data $replica(index.b)]

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 saved_values saved_array
  foreach label $saved_labels value $saved_values {
    set saved_array($label) $value
  }
}

set NEWTEMP $replica(temperature)
seed [expr int(0*srand(int(100000*rand()) + 100*$replica_id) + 100000*rand())]
langevinTemp $NEWTEMP
langevinPistonTemp $NEWTEMP
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 {
  temperature $NEWTEMP
}

outputEnergies [expr $steps_per_run / 10]
dcdFreq [expr $steps_per_run * $runs_per_frame]

source $namd_config_file

set history_file [open [format "$job_output_root.$replica_id.history" $replica_id] "w"]
fconfigure $history_file -buffering line

if { $i_run == 0 } { 
  run 50000
}

while {$i_run < $num_runs} {

  run $steps_per_run
  save_array
  incr i_step $steps_per_run
  set TEMP $saved_array(TEMP)
  set POTENTIAL [expr $saved_array(TOTAL) - $saved_array(KINETIC)]
  puts $history_file "$i_step $replica(index) $NEWTEMP $TEMP $POTENTIAL"

  if { $i_run % 2 == 0 } {
    set swap a; set other b
  } else {
    set swap b; set other a
  }

  set doswap 0
  if { $replica(index) < $replica(index.$swap) } {
    set temp $replica(temperature)
    set temp2 $replica(temperature.$swap)
    set BOLTZMAN 0.001987191
    set dbeta [expr ((1.0/$temp) - (1.0/$temp2)) / $BOLTZMAN]
    set pot $POTENTIAL
    set pot2 [replicaRecv $replica(loc.$swap)]
    set delta [expr $dbeta * ($pot2 - $pot)]
    set doswap [expr $delta < 0. || exp(-1. * $delta) > rand()]
    replicaSend $doswap $replica(loc.$swap)
    if { $doswap } {
      set rid $replica(index)
      set rid2 $replica(index.$swap)
      puts stderr "EXCHANGE_ACCEPT $rid ($temp) $rid2 ($temp2) RUN $i_run"
      incr replica(exchanges_accepted)
    }
    incr replica(exchanges_attempted)
  }
  if { $replica(index) > $replica(index.$swap) } {
    replicaSend $POTENTIAL $replica(loc.$swap)
    set doswap [replicaRecv $replica(loc.$swap)]
  }

  set newloc $r
  if { $doswap } {
    set newloc $replica(loc.$swap)
    set replica(loc.$swap) $r
  }
  set replica(loc.$other) [replicaSendrecv $newloc $replica(loc.$other) $replica(loc.$other)]
  set oldidx $replica(index)
  if { $doswap } {
    set OLDTEMP $replica(temperature)
    array set replica [replicaSendrecv [array get replica] $newloc $newloc]
    set NEWTEMP $replica(temperature)
    rescalevels [expr sqrt(1.0*$NEWTEMP/$OLDTEMP)]
    langevinTemp $NEWTEMP
    langevinPistonTemp $NEWTEMP
  }

  # puts stderr "iteration $i_run replica $replica(index) now on rank $r"
  # replicaBarrier

  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 temp $replica(temperature)
  set temp2 $replica(temperature.$swap)
  set accepts $replica(exchanges_accepted)
  set ratio [expr 1.0*$accepts/$attempts]
  puts stderr "EXCHANGE_RATIO $temp $temp2 $accepts $attempts $ratio"
}

replicaBarrier