remarks  file  sa_omit_map.inp
remarks  Make an annealed, sigmaa-weighted 2fo-fc omit map


{===>}                                                           {* Resolution limits. *}
eval ($low_res=6.)                                                  {* For refinement. *}
eval ($low_res_map=50.)                                        {* For map calculation. *}
eval ($high_res=2.3)

{===>} parameter @TOPPAR:parhcsdx.pro end                            {*Read parameters.*}
              
{===>} structure @../../tutorial/generate/generate.psf end      {* Read structure file.*}

{===>} coor @../../tutorial/xtalrefine/slowcool.pdb                {* Read coordinates.*}

{===>} 
 evaluate ($init_temp=1000. )                                   {*Starting temperature.*}  
       
{===>}                                                          {* Minimization steps. *}
eval ($nstep_mini=20)  !20


{===>} 
 vector do (charge=0) ( all )                                 {* Turn off all charges. *}


 flags                                       
    include pele pvdw xref                   
    ?                                        
 end                                         
                                             
 xrefine 
                                  
{===>}
   a=61.76 b=40.73 c=26.74 alpha=90.0 beta=90.0 gamma=90.0               {* Unit cell.*}

{===>}
   eval ($SG="P2(1)2(1)2(1)") 
   
   @SYMMETRY:spacegroup.lib

   @SYMMETRY:scatter.lib

{===>} 
   nreflections=15000
   reflection @../../tutorial/xtalrefine/amy.fob end               {* Read reflections.*}


  resolution $low_res $high_res 

{===>}
   do ( fobs=0)  ( ampl(fobs) <= 2.0*sigma)                 {*Sigma cutoff.*}
   fwind=0.1=100000

   method=FFT          
   
   fft 
      memory=1000000   
   end                 
      

   tolerance=0.2                                  {*Tolerance for dynamics.*}

{===>}
   wa=26000                                  {*Weight from job "check.inp".*}
 end                  
   

{===>}                                             {*Define omitted region.*}
vector ident ( store1 ) ( byresidue ( residue 20 around 8. ) )

                         {*Define immediate surroundings of omitted region.*}
vector ident ( store2 )  ( byres ( ( store1 ) around 3.0 and not store1 ) )   

                       {*Exclude omitted region for empirical energy terms.*}
constraints interaction=( not store1  )=( not store1 ) end 

              {*Keep omitted region fixed during dynamics and minimization.*} 
constraints fix=( store1 ) end  
             
                                   {*Exclude omitted region for X-ray term.*}
xrefin
   selection=( not store1 and not hydrogen)
end 


                     {*Apply harmonic restraints to immediate surroundings.*}
vector do ( harmonic=10. ) ( store2 ) 
vector do ( refx=x ) ( all )
vector do ( refy=y ) ( all )
vector do ( refz=z ) ( all )
flags include harm end

     
           {*The following dynamics and minimization sections represent the*}
           {*slow-cooling SA protocol.                                     *}     

 set seed=432324368 end {*Set the initial random seed for the v-assignment.*}
     
 vector do (vx=maxwell($init_temp)) ( all )
 vector do (vy=maxwell($init_temp)) ( all )
 vector do (vz=maxwell($init_temp)) ( all )

 vector do (fbeta=100.) ( all )

 evaluate ($1=$init_temp)
 while ($1 > 300.0) loop main     
          
    dynamics  verlet 
      timestep=0.0005            
      nstep=50
      iasvel=current 
      nprint=5  iprfrq=0
      tcoupling=true tbath=$1   
    end                                         
    evaluate ($1=$1-25)
 end loop main

 xrefin     
   tolerance=0.0    lookup=false          {*This makes the minimizer happy.*}             
 end
       
 minimize powell                                      {*Final minimization.*}
   nstep=$nstep_mini                                         
   drop=10.0
 end

                                  {*Compute the 2 fobs - fcalc SA omit map.*}
 xrefin

     resolution  $low_res_map 2.0                      {*Resolution range. *}
                                   {*We use a lower limit than was used in *}
                                   {*refinement to improve connectivity.   *}
     update




{===>}                                                                    {* Scaling. *}
   do (fcalc=scale[bins=1](fobs,fcalc)*fcalc)  
                  ($low_res >= d >= $high_res and ampl(fobs) > 0 )

                                                                 {* Sigmaa weighting. *}

   do (OA=norm(amplitude(FOBS)))    ( $low_res >= d >= $high_res and ampl(fobs) > 0 ) 
   do (OB=norm(amplitude(FCALC)))   ( $low_res >= d >= $high_res and ampl(fobs) > 0 )
   do (OC=SIGA(OA,OB))              ( $low_res >= d >= $high_res and ampl(fobs) > 0 )
   do (OD= OC sqrt( save(amplitude(fobs)^2)  / save(amplitude(fcalc)^2) ) )
                                    ( $low_res >= d >= $high_res and ampl(fobs) > 0 )
   do (fom= 2 OC OA OB / ( max(0.0001,1-OC^2) ))
                                    ( $low_res >= d >= $high_res and ampl(fobs) > 0 )
   do (fom= i1overi0(fom))                            ( acentric and ampl(fobs) > 0 )
   do (fom= tanh(fom/2))                              ( centric  and ampl(fobs) > 0 )


                                                {* Compute FT of (fom 2fo- D fc) map. *}
   do (deri1 = combine(2 fom ampl(fobs) - OD ampl(fcalc), phase(fcalc)) ) 
               ( $low_res >= d >= $high_res and acentric and ampl(fobs) > 0 ) 
   do (deri1 = combine(fom ampl(fobs), phase(fcalc)) ) 
               ( $low_res >= d >= $high_res and centric and ampl(fobs) > 0 ) 

   do (map1=ft(deri1)) ( $low_res >= d >= $high_res and ampl(fobs) > 0 ) 
 

    remarks  sigmaa-weighted 2fofc map

    write map                                        

      from=map1
 
{===>}                                       {* Selection of map region to be written. *}
      extend=molecule                               
      cushion=6.0                                    
      selection=( store1 ) 
                                                    

{===>}  output=sa_omit_map.map                                       {*Output file.*}

    end                                         

 stop