remarks file:  xtalmr/search.inp 
remarks "Direct" rotation search using the constant domain only.  

{===>} structure @fab2hfl.psf end                    {*Read structure file.*}

{===>} coor @fab2hfl.pdb                                {*Read coordinates.*}

evaluate ($wa=10000.)                {* Use this weight if the only active *}
                                     {* energy term is XREF, otherwise, use*}
                                     {* the weight obtained from the       *}
                                     {* "check.inp" protocol.              *}


xrefin
{===>}                                             {*Unit cell for crystal.*}
   a=44.144 b=164.69 c=70.17 alpha=90. beta=108.50 gamma=90.

{===>}
   symmetry=(x,y,z)                 {*Operators for crystal symmetry P2(1).*}
   symmetry=(-x,y+1/2,-z) 


   SCATter ( chemical C* ) 
     2.31000 20.8439 1.02000 10.2075 1.58860 .568700 .865000 51.6512 .215600
   SCATter ( chemical N* )
     12.2126 .005700 3.13220 9.89330 2.01250 28.9975 1.16630 .582600 -11.529
   SCATter ( chemical O* )
      3.04850 13.2771 2.28680 5.70110 1.54630 .323900 .867000 32.9089 .250800
   SCATter ( chemical S* )
      6.90530 1.46790 5.20340 22.2151 1.43790 .253600 1.58630 56.1720 .866900
   SCATter ( chemical P* )
      6.43450 1.90670 4.17910 27.1570 1.78000 0.52600 1.49080 68.1645 1.11490
   SCATter ( chemical FE* )
      11.1764 4.61470 7.38630 0.30050 3.39480 11.6729 0.07240 38.5566 0.97070

{===>} 
   nreflections=50000 
   reflection @3r9a_det.fob end                         {*Read reflections.*}

{===>}   
   resolution 10.0 4.                                   {*Resolution range.*}

   reduce   
   do amplitude ( fobs = fobs * step(fobs - 2.0*sigma))     {*Sigma cutoff.*}
   fwind=0.1=100000 
      

{===>}    
   method=fft       
   fft 
       memory=1000000                   {*FFT method with memory statement.*}
   end

   wa=$wa 
   target=E2E2                 {*Specify target used for both PC-refinement*}
                               {*and translation search.                   *}
   mbins=20                        {*Number of bins used for E calculation.*}

   tolerance=0. lookup=false              {*This makes the minimizer happy.*}

                                          {*Expand data to a P1 hemisphere.*}
   expand 

   selection=(( resid 107:212 and segid L ) or ( resid 117:213 and segid H )) 
end

flags exclude * include xref end              {*Use only XREF energy terms.*}

{===>} set display=search.dat end                            {*Output file.*}

coor copy end

evaluate ($pcmax=-9999.)
evaluate ($pcave=0.)
evaluate ($pc2ave=0.)
evaluate ($count=0.)

set message=off end
set echo=off end

evaluate ($t1=0. )                                    {*Theta1 start value.*}
evaluate ($t1_int=5.)                                    {*Theta1 interval.*}
evaluate ($t1_end=360.01)                               {*Theta1 end value.*}
while ($t1 < $t1_end ) loop m1
  evaluate ($t2=0. )                                  {*Theta2 start value.*}
  evaluate ($t2_int=5.)                                  {*Theta2 interval.*}
  evaluate ($t2_end=90.01)                           {*Theta2 ending value.*}
  while ($t2 < $t2_end ) loop m2
     evaluate ($t3=0. )                               {*Theta3 start value.*}
     evaluate ($t3_int=5.)                               {*Theta3 interval.*}
     evaluate ($t3_end=360.01)                          {*Theta3 end value.*}
     while ($t3 < $t3_end ) loop m3
        coor swap end
        coor copy end       
        coor rotate euler=( $t1 $t2 $t3 ) end
        xrefin
          update 
          print target
        end
        display   $t1[F12.5]  $t2[F12.5]  $t3[F12.5]  $corr[F12.5]

        if ($corr > $pcmax) then
           evaluate ($pcmax=$corr) 
           evaluate ($t1max=$t1)
           evaluate ($t2max=$t2)
           evaluate ($t3max=$t3)
        end if
        evaluate ($pcave=$pcave+$corr)
        evaluate ($pc2ave=$pc2ave+$corr^2)
        evaluate ($count=$count+1)

        evaluate ($t3=$t3+$t3_int)
     end loop m3
  evaluate ($t2=$t2+$t2_int)
  end loop m2
evaluate ($t1=$t1+$t1_int)
end loop m1

close search.dat end

{====>} set display=search.summary end         {* Output file for summary. *}

evaluate ($pcave=$pcave/$count)
evaluate ($pc2ave=sqrt($pc2ave/$count-$pcave^2))

display   max. peak  (theta1,       theta2,       theta3,        value): 
display          $t1max[F12.5]  $t2max[F12.5]  $t3max[F12.5]  $pcmax[F12.5]
display   mean of rotation function: $pcave[F12.5]
display   standard deviation around mean: $pc2ave[F12.5]


      {* To sort the rotation function values we can use system utilities. *}

      {* UNIX operating systems: *}
system  sort -o search.sort +3 -n -r search.dat

      {* VMS operating systems: *}
!system sort/key=(position:43,size:14,decimal,desc) search.dat search.sort


stop