remarks file  nmr/sa.inp
remarks  Simulated annealing protocol for NMR structure determination.
remarks  The starting structure for this protocol can be any structure with
remarks  a reasonable geometry, such as randomly assigned torsion angles or 
remarks  extended strands.  
remarks
remarks  Author: Michael Nilges

!
! Please cite the following papers when using this protocol:
!
! 1. Nilges, M., Clore, G.M., and Gronenborn, A.M.  (1988). 
!    Determination of three-dimensional structures of proteins from 
!    interproton distance data by dynamical simulated annealing
!    from a random array of atoms,  FEBS Lett. 239, 129-136. 
! 2. Nilges, M., Kuszewski, J., and Brunger, A.T. (1991). 
!    In: Computational Aspects of the Study of Biological
!    Macromolecules by NMR, (J.C. Hoch, ed.).  New York: Plenum Press. 
!


{====>} 
evaluate ($init_t = 1000 )       {*Initial simulated annealing temperature.*}
{====>} 
evaluate ($high_steps= 6000 )         {*Total number of steps at high temp.*}
{====>} 
evaluate ($cool_steps = 3000 )      {*Total number of steps during cooling.*}

parameter                                        {*Read the parameter file.*}
{====>} 
   @TOPPAR:parallhdg.pro
end

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

{====>} 
coordinates @generate_template.pdb                  {*Read the coordinates.*}

    
noe
{====>}                                                         
   nres=3000             {*Estimate greater than the actual number of NOEs.*}
   class all 
{====>} 
   @g_protein_noe.tbl                           {*Read NOE distance ranges.*}
end

{====>} 
@g_protein_dihe.tbl                       {*Read dihedral angle restraints.*}


{====>}
               {*If protein contains S-S bridges, appropriately modify and *}
               {*then uncomment the following lines. The S-S covalent bonds*}
               {*will be represented as fake NOE distances.                *}
noe  
!   assign (resid  5 and name sg) (resid 55 and name sg)  2.02 0.1 0.1 
!   assign (resid  14 and name sg) (resid 38 and name sg)  2.02 0.1 0.1 
!   assign (resid 30 and name sg) (resid 51 and name sg)  2.02 0.1 0.1 
end


flags exclude * include bonds angle impr vdw noe cdih end

                        {*Friction coefficient for MD heatbath, in 1/ps.   *}
vector do (fbeta=10) (all)
                        {*Uniform heavy masses to speed molecular dynamics.*}
vector do (mass=100) (all)

noe                             {*Parameters for NOE effective energy term.*}
  ceiling=1000                     
  averaging  * cent
  potential  * soft
  scale      * 50.                           
  sqoffset   * 0.0  
  sqconstant * 1.0
  sqexponent * 2
  soexponent * 1
  asymptote  * 0.1                         {*Initial value--modified later.*}
  rswitch    * 0.5
end

parameter                       {*Parameters for the repulsive energy term.*}
    nbonds
      repel=1.                   {*Initial value for repel--modified later.*}
      rexp=2 irexp=2 rcon=1.
      nbxmod=3                
      wmin=0.01 
      cutnb=4.5 ctonnb=2.99 ctofnb=3. 
      tolerance=0.5 
   end
end

restraints dihedral 
      scale=5.                 
end

{====>} 
evaluate ($end_count=10)          {*Loop through a family of 10 structures.*}

coor copy end

evaluate ($count = 0)
while ($count < $end_count ) loop main 
    
   evaluate ($count=$count+1)

   coor swap end
   coor copy end

     {* ============================================= Initial minimization.*}
   restraints dihedral   scale=5.   end 
   noe asymptote * 0.1  end
   parameter  nbonds repel=1.   end end 
   constraints interaction 
           (not name SG) (all) weights * 1  vdw 0.002 end end
   minimize powell nstep=50 drop=10.  nprint=25 end


     {* ======================================== High-temperature dynamics.*}
   constraints interaction (not name SG) (all)  
               weights * 1  angl 0.4  impr 0.1 vdw 0.002 end end

   evaluate ($nstep1=int($high_steps * 2. / 3. ) ) 
   evaluate ($nstep2=int($high_steps * 1. / 3. ) ) 

   dynamics  verlet
      nstep=$nstep1   timestep=0.005   iasvel=maxwell   firstt=$init_t 
      tcoupling=true  tbath=$init_t  nprint=50  iprfrq=0
   end


     {* ============== Tilt the asymptote and increase weights on geometry.*}
   noe asymptote * 1.0  end

   constraints interaction 
              (not name SG) (all) weights * 1  vdw 0.002  end end
   dynamics  verlet
      nstep=$nstep2   timestep=0.005    iasvel=current   tcoupling=true 
      tbath=$init_t  nprint=50  iprfrq=0 
   end

    {* ==================================================  Cool the system.*}

   restraints dihedral   scale=200.   end  

   evaluate ($final_t = 100)      { K }
   evaluate ($tempstep = 50)      { K }
   
   evaluate ($ncycle = ($init_t-$final_t)/$tempstep)
   evaluate ($nstep = int($cool_steps/$ncycle))
  
   evaluate ($ini_rad  = 0.9)        evaluate ($fin_rad  = 0.75)
   evaluate ($ini_con=  0.003)       evaluate ($fin_con=  4.0)

   evaluate ($bath  = $init_t)
   evaluate ($k_vdw = $ini_con)
   evaluate ($k_vdwfact = ($fin_con/$ini_con)^(1/$ncycle))
   evaluate ($radius=    $ini_rad)
   evaluate ($radfact = ($fin_rad/$ini_rad)^(1/$ncycle))

   evaluate ($i_cool = 0)
   while ($i_cool < $ncycle) loop cool
      evaluate ($i_cool=$i_cool+1)

      evaluate ($bath  = $bath  - $tempstep)       
      evaluate ($k_vdw=min($fin_con,$k_vdw*$k_vdwfact))
      evaluate ($radius=max($fin_rad,$radius*$radfact))

      parameter  nbonds repel=$radius   end end 
      constraints interaction (not name SG) (all) 
                     weights * 1. vdw $k_vdw end end

      dynamics  verlet
         nstep=$nstep time=0.005 iasvel=current firstt=$bath 
         tcoup=true tbath=$bath nprint=$nstep iprfrq=0  
      end


{====>}                                                  {*Abort condition.*}
      evaluate ($critical=$temp/$bath)
      if ($critical >  10. ) then
         display  ****&&&& rerun job with smaller timestep (i.e., 0.003) 
         stop
      end if

   end loop cool

   {* ================================================= Final minimization.*}

   constraints interaction (all) (all) weights * 1. vdw 1. end end
   parameter                  
      nbonds  
         repel=0.80  
         rexp=2 irexp=2 rcon=1. 
         nbxmod=3  
         wmin=0.01  
         cutnb=6.0 ctonnb=2.99 ctofnb=3.  
         tolerance=1.5  
      end  
   end

   minimize powell nstep=1000 drop=10.0 nprint=25 end

   {* =================================== Write out the final structure(s).*}
   print threshold=0.5 noe 
   evaluate ($rms_noe=$result)
   evaluate ($violations_noe=$violations)
   print threshold=5. cdih
   evaluate ($rms_cdih=$result)
   evaluate ($violations_cdih=$violations)
   print thres=0.05 bonds          
   evaluate ($rms_bonds=$result)
   print thres=5. angles
   evaluate ($rms_angles=$result)
   print thres=5. impropers
   evaluate ($rms_impropers=$result)
   remarks =============================================================== 
   remarks            overall,bonds,angles,improper,vdw,noe,cdih
   remarks energies: $ener, $bond, $angl, $impr, $vdw, $noe, $cdih 
   remarks =============================================================== 
   remarks            bonds,angles,impropers,noe,cdih
   remarks rms-d: $rms_bonds,$rms_angles,$rms_impropers,$rms_noe,$rms_cdih
   remarks =============================================================== 
   remarks               noe,  cdih
   remarks violations.: $violations_noe, $violations_cdih
   remarks =============================================================== 

{====>}                        {*Name(s) of the family of final structures.*}
   evaluate ($filename="sa_"+encode($count)+".pdb")

   write coordinates output =$filename end
      
end loop main

stop