{+ file: shift_solvent.inp +}
{+ directory: xtal_refine +}
{+ description: Move solvent molecules as close as possible to
                a reference molecule +}
{+ authors: Axel T. Brunger, and Paul D. Adams +}
{+ copyright: Yale University +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) are not quoted
   - do not remove any evaluate statements from the file 
   - the selections store1 through store6 are available for general use -}

{- begin block parameter definition -} define(

{======================= molecular structure =========================}

{* molecular topology file *}
{===>} structure_infile="water_pick.mtf";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param";
{===>} parameter_infile_2="CNS_TOPPAR:water_rep.param";
{===>} parameter_infile_3="";
{===>} parameter_infile_4="";
{===>} parameter_infile_5="";

{* molecule coordinate file *}
{===>} coordinate_infile="water_pick.pdb";

{====================== crystallographic data ========================}

{* space group *}
{* use International Table conventions with subscripts substituted
   by parenthesis *}
{===>} sg="P2(1)2(1)2(1)";

{* unit cell parameters in Angstroms and degrees *}
{+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +}
{===>} a=61.76;
{===>} b=40.73;
{===>} c=26.74;
{===>} alpha=90;
{===>} beta=90;
{===>} gamma=90;

{========================== atom selection ===========================}

{* select atoms in reference molecule *}
{===>} atom_select=(known and (not resname TIP));

{* select solvent molecules to be moved *}
{===>} atom_solvent=(resname TIP);

{=========================== output files ============================}

{* output coordinate file *}
{===>} coordinate_outfile="shift_solvent.pdb";

{* format output coordinates for use in o *}
{* if false then the default CNS output coordinate format will be used *}
{+ choice: true false +}
{===>} pdb_o_format=true;

{* output listing file *}
{===>} list_outfile="shift_solvent.list";

{===========================================================================}
{        things below this line do not normally need to be changed          }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates @&coordinate_infile

 parameter
   if ( &BLANK%parameter_infile_1 = false ) then
      @@&parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      @@&parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      @@&parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      @@&parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      @@&parameter_infile_5
   end if
 end

 parameter
   nbonds cutnb=1000 end
 end

 xray

   @CNS_XTALLIB:spacegroup.lib (sg=&sg;
                                sgparam=$sgparam;)

   a=&a b=&b c=&c  alpha=&alpha beta=&beta gamma=&gamma
  
 end

 coord fractionalize end

 show ave(x) ( &atom_select and not ( &atom_solvent ) )
 evaluate ($refx=$result)
 show ave(y) ( &atom_select and not ( &atom_solvent ) )
 evaluate ($refy=$result)
 show ave(z) ( &atom_select and not ( &atom_solvent ) )
 evaluate ($refz=$result)

 coord orthogonalize end

 for $id in id ( tag and ( &atom_solvent ) ) loop id

   evaluate ($min.$id=9999)
   evaluate ($min_symm.$id=1)
   evaluate ($min_dx.$id=0)
   evaluate ($min_dy.$id=0)
   evaluate ($min_dz.$id=0)

 end loop id

 if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 coord copy end

 evaluate ($symm=1)

 while ( $symm <= $symmetry ) loop symm

   coord swap end
   coord copy end

   coord symmetry $symmetry_op_$symm selection=( &atom_solvent )
   end

   coord fractionalize end

   for $id in id ( tag and ( &atom_solvent ) ) loop id1

     show ave(x) ( byresidue ( id $id ) )
     evaluate ($molx=$result)
     show ave(y) ( byresidue ( id $id ) )
     evaluate ($moly=$result)
     show ave(z) ( byresidue ( id $id ) )
     evaluate ($molz=$result)

     for $gid in id ( byresidue ( id $id ) ) loop gid1
       evaluate ($shift_dx.$gid=int($refx-$molx))
       evaluate ($shift_dy.$gid=int($refy-$moly))
       evaluate ($shift_dz.$gid=int($refz-$molz))
     end loop gid1

     do (x=x+$shift_dx.$id) ( byresidue ( id $id ) )
     do (y=y+$shift_dy.$id) ( byresidue ( id $id ) )
     do (z=z+$shift_dz.$id) ( byresidue ( id $id ) )

   end loop id1

   coord orthogonalize end

   do ( store7=x ) ( &atom_solvent )
   do ( store8=y ) ( &atom_solvent )
   do ( store9=z ) ( &atom_solvent )

   for $dx in ( -1 0 1 ) loop dx

     for $dy in ( -1 0 1 ) loop dy

       for $dz in ( -1 0 1 ) loop dz

         do ( x=store7 ) ( &atom_solvent )
         do ( y=store8 ) ( &atom_solvent )
         do ( z=store9 ) ( &atom_solvent )

         coord fractionalize end

         coord translate 
           vector=( $dx $dy $dz ) selection=( &atom_solvent )
         end

         coord orthogonalize end

         for $id in id ( tag and ( &atom_solvent ) ) loop id2

           igroup interaction
             ( byresidue ( id $id ) ) 
             ( &atom_select and not ( &atom_solvent ) )
           end

           distance
             from=( byresidue ( id $id ) )
             to=( &atom_select and not ( &atom_solvent ) )
             cuton=0
             cutoff=1000
             disp=rmsd
           end           

           evaluate ($gmin=9999)
           for $gid in id ( byresidue ( id $id ) ) loop gid2
             if ( $min.$gid < $gmin ) then
               evaluate ($gmin=$min.$gid)
             end if
           end loop gid2

           show min( rmsd ) ( byresidue ( id $id ) )

           if ( $result < $gmin ) then
             for $gid in id ( byresidue ( id $id ) ) loop gid3
               show ( rmsd ) ( id $gid )
               evaluate ($min.$gid=$result)
               do (refx=x) ( id $gid )
               do (refy=y) ( id $gid )
               do (refz=z) ( id $gid )
               evaluate ($min_symm.$gid=$symm)
               evaluate ($min_dx.$gid=$dx+$shift_dx.$gid)
               evaluate ($min_dy.$gid=$dy+$shift_dy.$gid)
               evaluate ($min_dz.$gid=$dz+$shift_dz.$gid)
             end loop gid3
           end if

         end loop id2

       end loop dz

     end loop dy

   end loop dx

   evaluate ($symm=$symm+1)

 end loop symm

 do (x=refx) ( &atom_solvent )
 do (y=refy) ( &atom_solvent )
 do (z=refz) ( &atom_solvent )

 set display=&list_outfile end

 display
 display new position for each molecule is generated by:
 display
 display    r_new = ( R_symm * r_old ) + T_final
 display

 evaluate ($counter=1)

 for $id in id ( tag and &atom_solvent ) loop id

   show (segid) (id $id)
   evaluate ($segid=$result)
   show (resid) (id $id)
   evaluate ($resid=$result)
   show (resname) (id $id)
   evaluate ($resname=$result)

   display ================================================================
   display

   display solvent molecule $counter: $resname $resid $segid
   evaluate ($symm=$min_symm.$id)
   display   minimimum after-> symmetry operator= $symmetry_op_$symm
   display               and-> translation= ( $min_dx.$id[i2] , $min_dy.$id[i2] , $min_dz.$id[i2] )
   display   minimum distance-> $min.$id[f8.3] Angstroms

   display

   evaluate ($counter=$counter+1)

 end loop id

 display ================================================================

 set display=OUTPUT end

 @CNS_XTALMODULE:write_pdb (pdb_o_format=&pdb_o_format;
                            coordinate_outfile=&coordinate_outfile;
                            sgparam=$sgparam;)

 stop