subroutine int 2FOFC int FOFC char SEGMENTS char PROT_SEG

! set cutoffs
    set vari TOO_CLOSE local real = 2.3
    set vari TOO_FAR local real = 4.5
    set vari PAIR_CUT local real = 2.3
    set vari PAIR_LIST local real = 2.5

    set vari 2FO-FC_CUT local real = 1.1
    set vari FO-FC_CUT local real = 2.3
    set vari FO-FC_ION local real = 1.2
    set vari TEMP_CUT local real = 70.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 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      

! 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.

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

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 xxx

make atom from point reindex
key xxx select .not old end  
key 3sigma sele xxx .a map FOFC 3.1 1000. end

! show key xxx

! 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 xxx 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 xxx ) .and. xxx end
delete atom select .not. -
    ( around key PROTEIN dist TOO_FAR from xxx .or 3sigma ) .and. xxx end
delete atom select around key dummy dist TOO_FAR from xxx .and. xxx end
delete atom select around key orig dist TOO_CLOSE from xxx end
!delete atom select around key orig dist TOO_CLOSE from xxx .a .not 3sigma end
!delete atom select around key orig dist TOO_CLOSE * 0.6 from 3sigma 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 xxx end atom OH2
calc bond distance PEAK_SIZE select xxx end
make segm select xxx end
make atom merge ( resi select xxx 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 xxx end resi H2O
set temp select xxx end = 30.0
set weigh select xxx end = 1.0
rena select xxx end seg WAT2


! 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 xxx .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. xxx 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. ( xxx .or. SYM_ATOM ) end
del atom sele by sequence out .and. ( xxx .or. SYM_ATOM ) 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.

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

write over file pair.list pair

! 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. xxx end auto

! clean keys etc...

delete atom sele segm name #* end
key xxx drop
key dw drop
key dummy drop
key nnn drop
key NNN drop
key out drop
key PROTEIN drop
key SYM_ATOM drop
key 3sigma drop
key orig drop


! 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
return