#!/usr/bin/perl -w
use strict;
use lib $ENV{ 'MAIN_CONF'}  ;
use main_config ;
my %commands = () ;

if ( @ARGV <= 0 ) { push ( @{$commands{'help'}}, ' ' ) ; } 
else { (%commands) = main_config::get_commands ( @ARGV ) ; }

my $conf_file ;
if ( -f "main_config.data" )
{    $conf_file = "main_config.data" ; }
else
{    $conf_file = join ('', $ENV{"MAIN_CONF"}, "main_config.data" ) ; }

my %my_main = main_config::read_main_config_file ( $conf_file ) ; 

#####################################################################

if ( !exists($my_main{"BREF_OVER_STEP"}  ) )  { push ( @{$my_main{"BREF_OVER_STEP"}}, -1 ) ; }
if ( !exists($my_main{"BREF_IND_STEP"}  ) )  { push ( @{$my_main{"BREF_IND_STEP"}}, 20 ) ; }
if ( !exists($my_main{"BREF_GROUP"}  ) )  { push ( @{$my_main{"BREF_GROUP"}}, "SINGLE" ) ; }
if ( !exists($my_main{"BREF_FORCE_GROUP"}  ) )  { push ( @{$my_main{"BREF_FORCE_GROUP"}}, -1.0 ) ; }
if ( !exists($my_main{"BREF_FORCE_BOND"}  ) )  { push ( @{$my_main{"BREF_FORCE_BOND"}}, 0.004 ) ; }
if ( !exists($my_main{"BREF_FORCE_ANGLE"}  ) )  { push ( @{$my_main{"BREF_FORCE_ANGLE"}}, 0.003 ) ; }
if ( !exists($my_main{"BREF_KICK_B0"}  ) )  { push ( @{$my_main{"BREF_KICK_B0"}}, "0" ) ; }
if ( !exists($my_main{"BREF_KICK_B1"}  ) )  { push ( @{$my_main{"BREF_KICK_B1"}}, "0" ) ; }
if ( !exists($my_main{"BREF_STEP_COOL"}  ) )  { push ( @{$my_main{"BREF_STEP_COOL"}}, 0 ) ; }
if ( !exists($my_main{"PREF_BVAL_MIN"}  ) )  { push ( @{$my_main{"PREF_BVAL_MIN"}}, 5 ) ; }
if ( !exists($my_main{"PREF_BVAL_MAX"}  ) )  { push ( @{$my_main{"PREF_BVAL_MAX"}}, 60 ) ; }
if ( !exists($my_main{"BREF_NCS_FLAG"}  ) )  { push ( @{$my_main{"BREF_NCS_FLAG"}}, "OFF" ) ; }
if ( !exists($my_main{"BREF_NCS_FORCE"} ) )  { push ( @{$my_main{"BREF_NCS_FORCE"}}, 0.004 ) ; }
if ( !exists($my_main{"PREF_TARGET"}  ) )  { push ( @{$my_main{"PREF_TARGET"}}, "ML" ) ; }
if ( !exists($my_main{"PREF_ML_MISSIN"}  ) )  { push ( @{$my_main{"PREF_ML_MISSIN"}}, 100 ) ; }
if ( !exists($my_main{"BREF_SCAL_XTARG"} ) )  { push ( @{$my_main{"BREF_SCAL_XTARG"}}, 2.0 ) ; }
if ( !exists($my_main{"PREF_MAP_XTARG"}  ) )  { push ( @{$my_main{"PREF_MAP_XTARG"}}, 3 ) ; }
if ( !exists( $my_main{"BULK_SOLV_METHOD"} ) ) { push ( @{$my_main{"BULK_SOLV_METHOD"}}, "NONE" ) } ; 
if ( !exists( $my_main{"DFR_RFREE"} ) ) { push ( @{$my_main{"DFR_RFREE"}}, "OFF" ) } ; 

######################################################

my $MACRO_FILE = "refine_b.cmds" ;
my $PROC_ERROR = "FALSE" ;
my $doit = "FALSE" ;

my $keyword ; 
foreach $keyword ( sort keys (%commands) ) 
{
    if ( $keyword eq 'h' || substr($keyword, 0,4) eq "help" ) { goto HELP }
    if ( $keyword eq 'v' || substr($keyword, 0,4) eq "over" ) { goto OVER }
    if ( $keyword eq 's' || substr($keyword, 0,4) eq "step" ) { goto STEP }
    if ( $keyword eq 'b' || substr($keyword, 0,6) eq "bramp" ) { goto BRAMP }
    if ( $keyword eq 'k' || substr($keyword, 0,4) eq "kick" ) { goto KICK }
    if ( $keyword eq 'c' || substr($keyword, 0,4) eq "cool" ) { goto COOL }

    if ( $keyword eq 'g' || substr($keyword, 0,5) eq "group" ) { goto GROUP }
    if ( $keyword eq 'd' || substr($keyword, 0,6) eq "define" ) { goto DEFINE }
    if ( $keyword eq 'l' || substr($keyword, 0,5) eq "angle" ) { goto ANGLE }
    if ( $keyword eq 'p' || substr($keyword, 0,4) eq "bond" ) { goto BOND }

#    if ( $keyword eq 'r' || substr($keyword, 0,4) eq "diff" ) { goto DIFF }
#    if ( $keyword eq 'l' || substr($keyword, 0,3) eq "lsq" ) { goto LSQ }
#    if ( $keyword eq 'x' || substr($keyword, 0,4) eq "xmap" ) { goto XMAP }
#    if ( $keyword eq 'q' || substr($keyword, 0,4) eq "qmap" ) { goto QMAP }

    if ( $keyword eq 't' || substr($keyword, 0,6) eq "target" ) { goto TARGET }
#    if ( $keyword eq 'a' || substr($keyword, 0,5) eq "atoms" ) { goto ATOMS }
    if ( $keyword eq 'd' || substr($keyword, 0,4) eq "dens" ) { goto DENS }
#    if ( $keyword eq 'i' || substr($keyword, 0,4) eq "diff" ) { goto DIFF }
    if ( $keyword eq 'r' || substr($keyword, 0,5) eq "ratio"  ) { goto RATIO }
    if ( $keyword eq 'm' || substr($keyword, 0,3) eq "map"  ) { goto MAP }
    if ( $keyword eq 'x' || substr($keyword, 0,4) eq "xmap" ) { goto XMAP }
    if ( $keyword eq 'q' || substr($keyword, 0,4) eq "qmap" ) { goto QMAP }

    if ( $keyword eq 'n' || substr($keyword, 0,3) eq "ncs" ) { goto NCS }
    if ( $keyword eq 'f' || substr($keyword, 0,4) eq "fncs" ) { goto FNCS }

    if ( $keyword eq 'o' || substr($keyword, 0,6) eq "output" ) { goto OUTPUT }
    if ( substr($keyword, 0,4) eq "doit" ) { $doit = "TRUE" ; next ; }
    print "unknown option >$keyword<\n" ; $PROC_ERROR = "TRUE" ;

#  -a|--atoms)   missing atoms for estimation of ML [@{$my_main{"PREF_ML_MISSIN"}}]

  HELP:
    print << "END_OF_TEXT" ;
  -h|--help)    prints this message with available options and current status 
  -v|--over)    overall B-value refinement  steps [-1 off; >0 on] = [$my_main{"BREF_OVER_STEP"}[0]]
  -s|--step)    ind. B-value refinement  steps [-1 off; >0 on ] = [$my_main{"BREF_IND_STEP"}[0]]

  -k|--kick)    sets starting and final maximal B-kick [0 - 30] current = [$my_main{"BREF_KICK_B0"}[0] - $my_main{"BREF_KICK_B1"}[0]]
  -c|--cool)    sets number of cooling cycles = [$my_main{"BREF_STEP_COOL"}[0]]

  -d|--define)  define group organisation [SINGLE/ATOM/RESIDUE/SEGMENT] = [$my_main{"BREF_GROUP"}[0]]
  -l|--angle)   B-value angle inter-group constrain force [-1.0 off; > 0.0 on] = [$my_main{"BREF_FORCE_ANGLE"}[0]]
  -p|--bond)    B-value bond inter-group constrain force [-1.0 off; > 0.0 on] = [$my_main{"BREF_FORCE_BOND"}[0]]
  -g|--group)   B-value intra-group constrain force [-1.0 hard; > 0.0 soft] = [$my_main{"BREF_FORCE_GROUP"}[0]]
  -b|--bramp)   bramp interval of valid B-values = [$my_main{"PREF_BVAL_MIN"}[0] - $my_main{"PREF_BVAL_MAX"}[0]]

  -t|--target)  target function for refinement [LSQ/ML/ML_EST/ML_APR]
                current [$my_main{"PREF_TARGET"}[0]]
  -r|--ratio)   approximate ratio target for standard deviations of boned B-values [1. - 10.]
                current = [@{$my_main{"BREF_SCAL_XTARG"}}]
  -q|--qmap)    MAP for rec. space refinement [@{$my_main{"PREF_MAP_XTARG"}}]

  -n|--ncs)     Non-Cryst-Similarity NCS B-value flag [ON/OFF] = [$my_main{"BREF_NCS_FLAG"}[0]]
  -f|--fncs)    B-value NCS force  = [$my_main{"BREF_NCS_FORCE"}[0]]

  -o|--output)  set file name of the created macro [$MACRO_FILE]
     --doit)    create the macro
END_OF_TEXT
    exit ;
OUTPUT: $MACRO_FILE = shift ( @{$commands{$keyword}} ) ; next ;
STEP:    @{$my_main{"BREF_IND_STEP"}} = shift ( @{$commands{$keyword}} ) ; next ;
OVER:    @{$my_main{"BREF_OVER_STEP"}} = shift ( @{$commands{$keyword}} ) ; next ;
BRAMP: 
    @{$my_main{"BREF_MIN"}} = shift ( @{$commands{$keyword}} ) ; 
    @{$my_main{"BREF_MAX"}} = shift ( @{$commands{$keyword}} ) ; 
    next ;
COOL:  @{$my_main{"BREF_STEP_COOL"}} = shift ( @{$commands{$keyword}} ) ; next ;
KICK:  
    @{$my_main{"BREF_KICK_B0"}} = shift ( @{$commands{$keyword}} ) ; 
    @{$my_main{"BREF_KICK_B1"}} = shift ( @{$commands{$keyword}} ) ; 
    next ;

DEFINE: @{$my_main{"BREF_GROUP"}} = shift ( @{$commands{$keyword}} ) ; next ;
ANGLE:  @{$my_main{"BREF_FORCE_ANGLE"}} = shift ( @{$commands{$keyword}} ) ; next ;
BOND:   @{$my_main{"BREF_FORCE_BOND"}} = shift ( @{$commands{$keyword}} ) ; next ;
GROUP:  @{$my_main{"BREF_FORCE_GROUP"}} = shift ( @{$commands{$keyword}} ) ; next ;

TARGET: @{$my_main{"PREF_TARGET"}} = shift ( @{$commands{$keyword}} ) ; next ;
#ATOMS:  @{$my_main{"PREF_ML_MISSIN"}} = shift ( @{$commands{$keyword}} ) ; next ;
QMAP:  @{$my_main{"PREF_MAP_XTARG"}} = shift ( @{$commands{$keyword}} ) ; next ;
RATIO: @{$my_main{"BREF_SCAL_XTARG"}} = shift ( @{$commands{$keyword}} ) ; next ;

#LSQ:    @{$my_main{"PREF_SCAL_XLSQ"}} = shift ( @{$commands{$keyword}} ) ; next ;
NCS:    @{$my_main{"BREF_NCS_FLAG"}} = shift ( @{$commands{$keyword}} ) ; next ;
FNCS:  @{$my_main{"BREF_NCS_FORCE"}} = shift ( @{$commands{$keyword}} ) ; next ;
}

if ( $PROC_ERROR eq "TRUE" ) { exit ; }

#######################################


$conf_file = "main_config.data" ; 
main_config::write_main_config_file ( $conf_file, %my_main ) ; 

###  . ${MAIN_CONF}main_config.data.save.sh

if ( $doit eq "FALSE" )  { exit ; }

#####################################

print "creating \"$MACRO_FILE\" \n" ;
if ( -f $MACRO_FILE ) 
{ 
    my $back = join ('', $MACRO_FILE, ".old" ) ; 
    system ( "mv  $MACRO_FILE $back" ) ; 
}
open ( MACRO,  "> $MACRO_FILE" ) || die "FILE OPEN ERROR: >$MACRO_FILE< \n" ;

my $I = $my_main{"BREF_KICK_B0"}[0] ;
#print "if test 137 ", $I, $my_main{"BREF_STEP_COOL"}[0], "\n" ;
# if ( $I <= 0 && $my_main{"BREF_STEP_COOL"}[0] <= 0 )
# {
#     print << "END_OF_TEXT" ;
#     create_refine_b.pl> warning:: COLL STEPS [$my_main{"BREF_STEP_COOL"}[0]] and/or 
# 				    KICK VALUE [$my_main{"BREF_KICK_B0"}[0]] unspecified
# END_OF_TEXT
# }
my $BB_ADJ = $my_main{"BREF_SCAL_XTARG"}[0] ;
print MACRO << "END_OF_TEXT" ;
subroutine char SEGMENTS int THE_MAP 

! renames new waters to prevent deletions of already accepted solvents
rename segm WAT sele segm name WAT2 end  
set vari II int segm sele segm name WAT end
if ( II .gt. 0 ) then
set vari WORK_SEGM char global - 
    segm sele segm name WORK_SEGM .o segm name WAT end
end_if

set weight sele atom name H* .a segm name SEGMENTS end = 0.0 
set temp sele atom name H* .a segm name SEGMENTS end = 0.0 
set temp sele temp -0.1 0.1 .a .not atom name H* .a \\ 
    segm name SEGMENTS end = 25.0 ! resets Bs for rebuilt atoms


mini adjust ( bval = $BB_ADJ step = 20 )

END_OF_TEXT

######## NCS region definition

if ( @{$my_main{"SG_NMOL"}}[0] > @{$my_main{"SG_NGROUP"}}[0] )
{
    print MACRO
	"! include all the non equivalent residues here\n", 
	" key out sele .not all end\n",
	" <define_key_out.com\n\n",
	"! defining groups \n" ;

    my ( $IGR, $GROUP );
    my $NCS_START = "FALSE" ;
    for ( $IGR = 1; $IGR <= @{$my_main{"SG_NGROUP"}}[0] ; $IGR++ )
    {
	$GROUP= join ( '', "SG_GROUP" , $IGR ) ;
#	print " groups: $IGR, <$GROUP> \n" ;
	if ( $NCS_START eq "FALSE" )
	{ 
		print MACRO " define init constr ncs \n" ;
		$NCS_START = "TRUE" ;
	    }
	    else
	    { 
		print MACRO " define constraint ncs next \n" ; 
	    }
	
	    my $MEMBER ;
	    foreach $MEMBER ( @{$my_main{$GROUP}} )
	    {
#	print " member: $MEMBER> \n" ;
		print MACRO " define constr ncs sele atom name CA N C O H CB %G* %D* \\ \n",
		"          .a ( segm name $MEMBER .a .not out ) end \n",
	    }
	    print MACRO " define constr ncs force @{$my_main{\"BREF_NCS_FORCE\"}}[0] \n", 
	    " define constr ncs b-force @{$my_main{\"BREF_NCS_FORCE\"}}[0] \n\n" ;
	}
    print MACRO " ener ncs @{$my_main{\"BREF_NCS_FLAG\"}}[0] \n" ;
}
if ( $my_main{"BREF_GROUP"}[0] eq "SINGLE" )
{ print MACRO "define sele segm name SEGMENTS end init group single\n" ; }
elsif ( $my_main{"BREF_GROUP"}[0] eq "ATOM" )
{ print MACRO "define sele segm name SEGMENTS end init group atom\n" ; }
elsif ( $my_main{"BREF_GROUP"}[0] eq "RESIDUE" )
{ print MACRO "define sele segm name SEGMENTS end init group residue\n" ; } 
elsif ( $my_main{"BREF_GROUP"}[0] eq "SEGMENT" )
{ print MACRO "define sele segm name SEGMENTS end init group segment\n" ; }
else 
{ 
    print MACRO "%ERROR: BREF_GROUP $my_main{\"BREF_GROUP\"}[0]: bad value\n" ; 
    print "%ERROR: BREF_GROUP $my_main{\"BREF_GROUP\"}[0]: bad value\n" ; 
}

if ( $my_main{"PREF_MAP_XTARG"}[0] !~ m"[0-9]" )
{
print MACRO << "END_OF_TEXT" ;

 set vari MAP_TARG = nmaps + 1
END_OF_TEXT
}
else
{
print MACRO << "END_OF_TEXT" ;

 set vari MAP_TARG = @{$my_main{"PREF_MAP_XTARG"}}[0] 
END_OF_TEXT
}

print MACRO << "END_OF_TEXT" ;
 make map MAP_TARG auto_grid 4.0 init cell real 
  if ( MAP_TARG .ne. nmaps ) return

ener b-val group = $my_main{"BREF_FORCE_GROUP"}[0] \\ 
     b-val bond = $my_main{"BREF_FORCE_BOND"}[0] \\ 
     b-val angle = $my_main{"BREF_FORCE_ANGLE"}[0] \\ 
     b-val mini = $my_main{"PREF_BVAL_MIN"}[0] \\ 
     b-val maxi = $my_main{"PREF_BVAL_MAX"}[0]

 ener x-targ on \\ 
      x-targ map MAP_TARG \\
      x-targ scale = 100

END_OF_TEXT


    if ( $my_main{"BREF_OVER_STEP"}[0] > 0 ) 
{
print MACRO << "END_OF_TEXT" ;

! Over all B-value Refinement based on Wilson plot differences
 ener dens map MAP_TARG

 mini b-val over \\ 
      sele segm name SEGMENTS end \\ 
      step = $my_main{"BREF_OVER_STEP"}[0] \\
      phase \( sele WORK_REFL end \\
END_OF_TEXT

    if ( $my_main{"BULK_SOLV_METHOD"}[0] ne "NONE" )
{    print MACRO " 	     partial \\ \n" ; }

    print MACRO << "END_OF_TEXT" ;
            ) 

END_OF_TEXT
}

my ( $STEP_COOL , $NCOOL ) ;
if ( $my_main{"BREF_IND_STEP"}[0] <= 0 ) { goto FINISH }

    if ( $my_main{"BREF_STEP_COOL"}[0] > 1 )
    {
	$STEP_COOL = $my_main{"BREF_STEP_COOL"}[0] - 1 ;
	$NCOOL = $STEP_COOL ;
    }
    else { $NCOOL=0 ; }

    for ( $I = 0; $I <= $NCOOL; $I++ )
    {
	if ( $my_main{"BREF_STEP_COOL"}[0] == 1 )
	{
	    print MACRO << "END_OF_TEXT" ;

! kick : 4. the usual one 
set temp kick $my_main{"BREF_KICK_B0"}[0] sele segm name SEGMENTS end 
END_OF_TEXT
        }
	elsif ( $my_main{"BREF_STEP_COOL"}[0] > 1 )
	{
#    K=$I
#    FRACTION=`expr $my_main{"BREF_STEP_COOL"}[0] - $K`
#    FRACTION=5
#    echo "ffffffffffffrrrrrrrr $my_main{"BREF_STEP_COOL"}[0] - $I" 
	    my $KICK_VAL = ( $STEP_COOL - $I) * ( $my_main{"BREF_KICK_B0"}[0] - $my_main{"BREF_KICK_B1"}[0] ) / $STEP_COOL ;
	    print MACRO "set temp kick $KICK_VAL sele segm name SEGMENTS end \n" ;
	}
	else { print MACRO "! no kick applied \n" ; }

my $REF_TARGET = $my_main{"PREF_TARGET"}[0] ;

if ( $REF_TARGET eq "ML_EST" ) 
{
print MACRO << "END_OF_TEXT" ;

< ?MAIN_CMDS:re_phase.cmds WORK_SEGM MAP_WORK
refl 
move fobs fsave
exit

END_OF_TEXT
}

print MACRO << "END_OF_TEXT" ;

 mini b-val individua  \\ 
      sele segm name SEGMENTS end \\ 
      step = $my_main{"BREF_IND_STEP"}[0] \\
      phase \( sele WORK_REFL end \\
END_OF_TEXT

if ( $my_main{"DFR_RFREE"}[0] eq "ON" )
{    print MACRO "             free sele TEST end \\ \n" ; }
    if ( $my_main{"BULK_SOLV_METHOD"}[0] ne "NONE" )
{    print MACRO " 	     partial \\ \n" ; }

if ( $REF_TARGET eq "LSQ" )
{    print MACRO " 	     target ( least_square ) \\ \n" ; }
elsif ( $REF_TARGET eq "ML" )
{    print MACRO " 	     target ( max_like ) \\ \n" ; }
elsif ( $REF_TARGET eq "ML_EST" )
{    print MACRO " 	     target ( estim_ml missing @{$my_main{\"PREF_ML_MISSIN\"}}[0] ) \\ \n" ; }
elsif ( $REF_TARGET eq "ML_APR" )
{    print MACRO " 	     target ( aprox_ml ) \\ \n" ; }

    print MACRO << "END_OF_TEXT" ;
            ) \\ 
     writ = 1 
END_OF_TEXT

}

FINISH:
## echo after DONE ===================

if ( $my_main{"SG_NMOL"}[0] > $my_main{"SG_NGROUP"}[0] )
{ print MACRO "ener ncs off" } ;

print MACRO << "END_OF_TEXT" ;

delete map MAP_TARG

 ener  dens ( on map MAP_2FOFC poly 2 ) \\ 
       ncs off
mini posit writ 10 step 80
return
END_OF_TEXT
close (MACRO ) ;
