set mypdb 2D8R
set logstr ""

proc getboxV {mysel} {
	set tmp [measure minmax $mysel]
	set L_x [expr [lindex [lindex $tmp 1] 0] - [lindex [lindex $tmp 0] 0]]
	set L_y [expr [lindex [lindex $tmp 1] 1] - [lindex [lindex $tmp 0] 1]]
	set L_z [expr [lindex [lindex $tmp 1] 2] - [lindex [lindex $tmp 0] 2]]
	return [expr $L_x*$L_y*$L_z]
}

mol load webpdb $mypdb
set allsel [atomselect top "not water"]
set logstr "${logstr}box volume (initial): [getboxV $allsel]\n\n"
$allsel writepdb $mypdb.pdb
$allsel delete

package require solvate 1.5
set mylog [open "solvate_rotate.log" w]
array set bounds { -x 0 +x 0 -y 0 +y 0 -z 0 +z 0 }
::Solvate::rotate_save_water $mypdb.pdb "not water" 36 $mylog
close $mylog
set allsel  [atomselect top all]
set logstr "${logstr}box volume (solvate,rotate_z_x): [getboxV $allsel]\n\n"
set V0zx [getboxV $allsel]
$allsel delete

mol delete all
mol load pdb $mypdb.pdb
set allsel  [atomselect top all]
$allsel moveby [vecsub {0.0 0.0 0.0} [measure center $allsel]]
set logstr "${logstr}box volume (initial): [getboxV $allsel]\n\n"

set N_rot 36
set d_phi [expr 360.0/$N_rot]
set d_theta [expr 360.0/$N_rot]
set d_gamma [expr 360.0/$N_rot]
set V0 [getboxV $allsel]
set bufcoord [$allsel get {x y z}]
set kV1 0
set kV2 0
set kV3 0
for {set k1 1} {$k1 < [expr $N_rot + 1]} {incr k1} {
  $allsel move [trans axis z $d_phi deg]
  for {set k2 1} {$k2 < [expr $N_rot + 1]} {incr k2} {
    $allsel move [trans axis x $d_theta deg]
    for {set k3 1} {$k3 < [expr $N_rot + 1]} {incr k3} {
      $allsel move [trans axis y $d_gamma deg]
      set V [getboxV $allsel]
      if {$V < $V0} {
        set V0 $V
        set kV1 $k1
        set kV2 $k2
        set kV3 $k3
        set bufcoord [$allsel get {x y z}]
      }
    }
  }
}
puts $V0  
$allsel set {x y z} $bufcoord
$allsel writepdb $mypdb.rotate_zxy.pdb
set logstr "${logstr}box volume (rotate_z_x_y): [getboxV $allsel]\n\n"
set V0zxy [getboxV $allsel]

puts $logstr
puts "The system was rotated by [expr $kV1*$d_phi] degrees around Z axis and [expr $kV2*$d_theta] degrees around X axis and [expr $kV3*$d_gamma] degrees around Y axis."
puts "[expr ($V0zx-$V0zxy)/$V0zx]"