#!/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{"WAT_TEMP"} ) ) { push ( @{$my_main{"WAT_TEMP"}}, "70.0" ) } ; 
if ( !exists( $my_main{"WAT_DENS"} ) ) { push ( @{$my_main{"WAT_DENS"}}, "1.1" ) } ; 
if ( !exists( $my_main{"WAT_PEAK"} ) ) { push ( @{$my_main{"WAT_PEAK"}}, "2.3" ) } ; 
if ( !exists( $my_main{"WAT_CLOSE"} ) ) { push ( @{$my_main{"WAT_CLOSE"}}, "2.0" ) } ; 
if ( !exists( $my_main{"WAT_FAR"} ) ) { push ( @{$my_main{"WAT_FAR"}}, "4.5" ) } ; 
if ( !exists( $my_main{"WAT_PAIR_CUT"} ) ) { push ( @{$my_main{"WAT_PAIR_CUT"}}, $my_main{"WAT_CLOSE"}[0] ) } ; 
if ( !exists( $my_main{"WAT_PAIR_LIST"} ) ) { push ( @{$my_main{"WAT_PAIR_LIST"}}, "2.5" ) } ; 
if ( !exists( $my_main{"WAT_ION"} ) ) { push ( @{$my_main{"WAT_ION"}}, "1.3" ) } ; 
if ( !exists( $my_main{"WAT_WHAT"} ) ) { push ( @{$my_main{"WAT_WHAT"}}, "BOTH" ) } ; 


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

my $MACRO_FILE = "gen_solvent.com" ;
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 'o' || substr($keyword, 0,6) eq "output" ) { goto OUTPUT }
    if ( $keyword eq 'b' || substr($keyword, 0,7) eq "b_value" ) { goto BVAL }
    if ( $keyword eq 'p' || substr($keyword, 0,4) eq "peak" ) { goto PEAK }
    if ( $keyword eq 'd' || substr($keyword, 0,7) eq "density" ) { goto DENSITY }
    if ( $keyword eq 'c' || substr($keyword, 0,5) eq "close" ) { goto CLOSE }
    if ( $keyword eq 'f' || substr($keyword, 0,3) eq "far" ) { goto FAR }
    if ( $keyword eq 'm' || substr($keyword, 0,5) eq "merge" ) { goto MERGE }
    if ( $keyword eq 'l' || substr($keyword, 0,4) eq "list" ) { goto LIST }
    if ( $keyword eq 'i' || substr($keyword, 0,3) eq "ion" ) { goto ION }
    if ( $keyword eq 'w' || substr($keyword, 0,4) eq "what" ) { goto WHAT }

    if ( substr($keyword, 0,4) eq "doit" ) { $doit = "TRUE" ; next ; }
    print "unknown option >$keyword<\n" ; $PROC_ERROR = "TRUE" ;

HELP:
    print << "END_OF_TEXT" ;
  -h|--help)      prints this message with available options and current status 
  -b|--b_value)   B-value cut for previous rounds in WAT [$my_main{"WAT_TEMP"}[0]]
  -p|--peak)      peak threshold in diff density map (fo-fc) [$my_main{"WAT_PEAK"}[0]]
  -d|--density)   density map (2fo-fc) cutoff [$my_main{"WAT_DENS"}[0]]
  -c|--close)     the closest density pick existing atom distance [$my_main{"WAT_CLOSE"}[0]]
  -f|--far)       the farthest density pick distance to an existing atom [$my_main{"WAT_FAR"}[0]]
  -m|--merge)     remove symmetry related peaks if they get closer than [$my_main{"WAT_PAIR_CUT"}[0]]
  -l|--list)      key "check_list" contains residues closer than [$my_main{"WAT_PAIR_LIST"}[0]] to peaks
  -i|--ion)       key "ion_list" conatins WAT atoms with fo-fc dens higher than [$my_main{"WAT_ION"}[0]]
  -w|--what)      what to do: [CHECK/NEW/BOTH] [$my_main{"WAT_WHAT"}[0]]
     --doit)      create the macro
END_OF_TEXT
    exit ;
OUTPUT: $MACRO_FILE = shift ( @{$commands{$keyword}} ) ; next ;
BVAL: @{$my_main{"WAT_TEMP"}} = shift ( @{$commands{$keyword}} ) ; next ;
PEAK: @{$my_main{"WAT_PEAK"}} = shift ( @{$commands{$keyword}} ) ; next ;
DENSITY: @{$my_main{"WAT_DENS"}} = shift ( @{$commands{$keyword}} ) ; next ;
CLOSE: @{$my_main{"WAT_CLOSE"}} = shift ( @{$commands{$keyword}} ) ; next ;
FAR: @{$my_main{"WAT_FAR"}} = shift ( @{$commands{$keyword}} ) ; next ;
MERGE: @{$my_main{"WAT_PAIR_CUT"}} = shift ( @{$commands{$keyword}} ) ; next ;
LIST: @{$my_main{"WAT_PAIR_LIST"}} = shift ( @{$commands{$keyword}} ) ; next ;
ION: @{$my_main{"WAT_ION"}} = shift ( @{$commands{$keyword}} ) ; next ;
WHAT: @{$my_main{"WAT_WHAT"}} = 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 ) ; 

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" ;

print MACRO << "END_OF_TEXT" ; 
subroutine int 2FOFC int FOFC char SEGMENTS char PROT_SEG

! set cutoffs
    set vari TOO_CLOSE local real = $my_main{"WAT_CLOSE"}[0]
    set vari TOO_FAR local real = $my_main{"WAT_FAR"}[0]
    set vari PAIR_CUT local real = $my_main{"WAT_PAIR_CUT"}[0]
    set vari PAIR_LIST local real = $my_main{"WAT_PAIR_LIST"}[0]

    set vari 2FO-FC_CUT local real = $my_main{"WAT_DENS"}[0]
    set vari FO-FC_CUT local real = $my_main{"WAT_PEAK"}[0]
    set vari FO-FC_ION local real = $my_main{"WAT_ION"}[0]
    set vari TEMP_CUT local real = $my_main{"WAT_TEMP"}[0]

! 2 peaks are one if closer than PEAK_SIZE
set vari PEAK_SIZE local real = 1.5

!
! below this line we usually do not change anything !!!!!!
!

set vari MAP_DIFF  local int = FOFC + 1

if ( 2FOFC .gt. FOFC ) set vari MAP_DIFF local int = 2FOFC + 1
make map MAP_DIFF from FOFC init 9999 - 
         around 8 sele segm name SEGMENTS end copy

delete atom sele segment name WAT2 \#* end

! assignes atoms to keys
! old water molecules
!DEFINE KEY
key wat select .not segment name PROT_SEG end   
! PROTEIN nonhydrogen atoms
key PROTEIN select segment name PROT_SEG .and. .not. atom name H* end  
! all nonhydrogen atoms
key orig select segment name SEGMENTS .and. .not. atom name H* end   
! all dummy heavy PROTEIN atoms
key dummy select weight -0.1 0.1 .and. PROTEIN end
key old select all end   

! removes all water molecules with oxygens placed at positions
! with lower electron density than 2FO-FC_CUT

key DW select map 2FOFC -20.0 2FO-FC_CUT - 0.3 from wat .and. atom name O* end
delete atom select by resi DW end      

! removes all water molecules with oxygen temperature factors higher 
! than TEMP_CUT

key DW select temp TEMP_CUT 1000.0 .and. atom name O* .and. wat end
delete atom select by resi DW end      
key DW drop

! Writes the solvent molecules with still positive difference electron
! density to the file ion.list. Their electron number seems to be too
! small for their position.

key ion_list select .not. atom name H* .and. -  
   map MAP_DIFF FO-FC_ION 100. from wat end 
write over file ion.list select ion_list end coordinates  xplor
END_OF_TEXT


    if ( $my_main{"WAT_WHAT"}[0] eq "BOTH" || $my_main{"WAT_WHAT"}[0] eq "NEW" )
{
print MACRO << "END_OF_TEXT" ;

make point init from map MAP_DIFF extrem FO-FC_CUT 200.0 
show point dens

! makes atom from points (peaks in the Fo-Fc map)
! and assignes the newly created peaks to the key XXXX

make atom from ( point )
key XXXX select .not old end  
key 3sigma sele XXXX .a map FOFC 3.1 1000. end

! show key XXXX

! deletes all the peaks that are placed at positions where the 
! electron density map (1) ,2Fobs-Fcalc, is lower
! than 2FO-FC_CUT sigma. 

 delete atom select map 2FOFC -300. 2FO-FC_CUT from XXXX end

! deletes peaks closer than TOO_CLOSE to any "orig" atom
! deletes peaks further than TOO_FAR to any "PROTEIN" atom
! deletes peaks closer than TOO_FAR to any dummy PROTEIN atom

!key tmp sele PROTEIN .a .not atom name CA N C CB O H end
!set weight sele around key 3sigma dist TOO_CLOSE * 0.6 from tmp end = 0.0

 delete atom select .not. -
     ( around key PROTEIN dist TOO_FAR from XXXX ) .and. XXXX end
 delete atom select .not. -
     ( around key PROTEIN dist TOO_FAR from XXXX .or 3sigma ) .and. XXXX end
 delete atom select around key dummy dist TOO_FAR from XXXX .and. XXXX end
 delete atom select around key orig dist TOO_CLOSE from XXXX end

! Since the newly created atoms (peaks) were created as dummy atoms (X)
! they had to be renamed to oxygen(OH2) further treatment.
! the peaks are merged according to their cluster formation.

rename select XXXX end atom OH2
calc bond distance PEAK_SIZE select XXXX end
make segm select XXXX end
make atom merge ( resi select XXXX end )

! renames the clusters into a single segment WAT2 and residues to H2O
! set their temperature factors to 30.0 and 
! crystallographical weights to 1.0

rena select XXXX end resi H2O
set temp select XXXX end = 30.0
set weigh select XXXX end = 1.0
rena select XXXX end seg WAT2

END_OF_TEXT
}

if ( $my_main{"WAT_WHAT"}[0] eq "CHECK" )
{
print MACRO "\nkey XXXX sele .not all end\n" ;
}

print MACRO << "END_OF_TEXT" ; 

! Prepare symmetry overlap check.
! initializes the variable ii and sets its value to the current number
! of segments (the variable is later applied in the command file
! REMOVE_SYMM.COM where the symmetry equivalent solvent atoms are removed
! selects all current atoms
! generates the symmetry related atoms from all solvent oxygen atoms
! that are closer than 4.0A to any atom in nnn (dist 4.0 cut) and does 
! not apply the identity operation (delete on)
! assigns key "symm" to all not "nnn" atoms


key NNN select XXXX .or. ( wat .and. atom name O* ) end

symmetry select NNN end \\
   around select PROTEIN end dist TOO_FAR cut delete on
key SYM_ATOM select segm name \#* end

! calculates the pair list between the solvent oxygen atoms
! and their symmetry related ones.
! a pair of atoms (peaks) is assumed to be identical when they are closer
! than PAIR_CUT

calc pair ( select NNN end select SYM_ATOM end rang 0.0 1.5 initialize )
key out sele .not all end
 <MAIN_UTILS:gen_solvent_remove_symm.com  
 dele atom sele out end

calc pair ( select SYM_ATOM .or. XXXX end -
    select orig .and. .not. atom name H* end -
    rang 0.0 PAIR_CUT initialize )

! removes the symmetry related peaks by checking the PROTEIN atoms

key out sele by pair orig .and. ( XXXX .or. SYM_ATOM ) end
 del atom sele by sequence out .and. ( XXXX .or. SYM_ATOM ) end

make pair init
calc pair ( select PROTEIN end initial -
          select ( SYM_ATOM .or. XXXX .or. wat ) .and. atom name O* end -
          rang 0.0 TOO_CLOSE )

key kk sele by pair atom name * end
dele atom sele by seq kk .a .n by resi PROTEIN  end 

! Calculates the remaining peaks closer than PAIR_LIST to any PROTEIN atom.
! The positions of the paired PROTEIN residues should be checked 
! at the display.

make pair init
calc pair ( select PROTEIN end initial -
          select ( SYM_ATOM .or. XXXX .or. wat ) .and. atom name O* end -
          rang 0.0 PAIR_LIST )

! write over file pair.list pair
key check_list sele -
  ( by resi pair site 2 .or by resi pair site 1 ) .and PROTEIN end 
! key kk sele by pair atom name * end
! dele atom sele by seq kk .a .n by resi PROTEIN  end 

! Perform hydrogen bond tests: all solvent atoms that have no hydrogen donnor or 
! acceptor atom (O*, N*) in their
!vicinity (3.4A) are deleted with the command file ATOM_NEIGH.COM.

<MAIN_UTILS:gen_solvent_atom_neigh  3.5

rename seq W sele wat .or. XXXX end auto

! clean keys etc...

delete atom sele segm name \#* end
key XXXX drop
key dummy drop
key NNN drop
key out drop
key PROTEIN drop
key SYM_ATOM drop
key 3sigma drop
key orig drop

END_OF_TEXT

if ( $my_main{"WAT_WHAT"}[0] ne "CHECK" )
{
print MACRO << "END_OF_TEXT" ;

! It remains to generate the missing hydrogens and update the WORK_SEGM
! variable.

build
sele segm name WAT2 end
fill
exit
calc coor append sele segm name WAT2 end
set vari WORK_SEGM global char segm sele segm name WORK_SEGM WAT2 end

! optimize the positions of the newly generated hydrogen atoms
! exploiting non-bonding interactions with the environment

<DEF_ALL WORK_SEGM
mini sele atom name H* .a segm name WAT2 end sele segm name WORK_SEGM end

image sele check_list end col 111 bond stick
set rad sele ion_list end = 0.5
image sele ion_list end col 111 atom ball

image sele segm name WAT2  end col 140 bond stick

return
END_OF_TEXT
}
close ( MACRO ) ;

