#!/usr/bin/perl

 if ($ARGV[0] eq '') {
        print "USAGE: PE_scan.pl atom1 atom2 atom3 atom4 > PEscan.py \nRETURNS: chimera python file for dihedral potential energy scan\n (default 30 degree increment)\nThe atoms must have unique names in the structure file.\nAlso creates clear_reply.py and save_reply.py.\n These are python scripts called from the PE_scan python script used in Chimera.
The name of the output file is reply_log.txt.\n This file is input to pot_energy.pl.\n";
        exit;
 }

#read in atoms
$a1=$ARGV[0];
$a2=$ARGV[1];
$a3=$ARGV[2];
$a4=$ARGV[3];

# generate python script for potential energy scan
print "from chimera import runCommand as rc\n";
#clear reply log (clear_reply.py is created below)
print 'rc("open clear_reply.py")',"\n";
# calculate minimized structure
print 'rc("minimize nogui true")',"\n";
print 'rc("open clear_reply.py")',"\n";
#change $ang if a different angle increment is desired 
$ang=30;
$increments=360/$ang;

# generate potential energy scan for multiple angles 
# default increment 30 degrees starting from minimized structure.
for ($i=1;$i<=$increments;$i++) {
print "\n";
print <<EOF;
rc("rotation 1 \@$a2\@$a3")
rc("rotation 1 $ang")
rc("~rotation 1")
rc("open angle.py")
rc("select \@$a1\@$a2\@$a3\@$a4")
rc("minimize nogui true freeze selected")
rc("~select")
EOF
}
#save in reply log (save_reply.py is created below)
print 'rc("open save_reply.py")';

#create angle.py to obtain appropriate dihedral angle for saving in reply log
open(STDOUT,">angle.py");
print STDOUT "from Midas import angle\n";
$str="\@$a1\@$a2\@$a3\@$a4";
print STDOUT "dihed = angle(\"$str\")";
close(STDOUT);

#create save_reply.py
open(STDOUT,">save_reply.py");
print STDOUT <<EOF;
def save_reply_log(path):
   from chimera import dialogs
   r = dialogs.find('reply')
   text = r.text.get('1.0', 'end')
   f = open(path, 'w')
   f.write(text)
   f.close()

save_reply_log("./reply_log.txt")
EOF
close(STDOUT);

#create clear_reply.py
open(STDOUT,">clear_reply.py");
print STDOUT <<EOF;
from chimera import dialogs
r = dialogs.find('reply')
r.Clear()
EOF
close(STDOUT);
