{+ file: generate_extended.inp +}
{+ directory: nmr_calc +}
{+ description: Generates an extended strand with ideal geometry 
                for each connected polymer.  
                The molecular structure file must not contain any 
                closed loops except disulfide bonds which are automatically
                excluded from the generation of the strand conformation.  +}
{+ authors: Axel T. Brunger +}
{+ copyright: Yale University +}

{- begin block parameter definition -} define(

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

{* structure file(s) *}
{===>} structure_file="il8.mtf";

{* parameter file(s) *}
{===>} par_1="CNS_TOPPAR:protein-allhdg.param";
{===>} par_2="";
{===>} par_3="";
{===>} par_4="";
{===>} par_5="";

{======================= input parameters ============================}

{* maximum number of trials to generate an acceptable structure *}
{===>} max_trial=10;

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

{* output coordinates *}
{===>} output_coor="il8_extended.pdb";
                                  
{===========================================================================}
{        things below this line do not normally need to be changed          }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_file end

 parameter
   if (&par_1 # " ") then
      @@&par_1
   end if
   if (&par_2 # " ") then
      @@&par_2
   end if
   if (&par_3 # " ") then
      @@&par_3
   end if
   if (&par_4 # " ") then
      @@&par_4
   end if
   if (&par_5 # " ") then
      @@&par_5
   end if
 end

{ Set force constants for S-S bond lengths and angles to zero  }
parameter
   bonds ( name SG ) ( name SG ) 0. 1. 
end


igroup interaction=(all) (all) end

ident (x) ( all )
do (x=x/5.) ( all )
do (y=random(0.5) ) ( all )
do (z=random(0.5) ) ( all )

flags exclude * include bond angle impr dihe vdw end
parameter
   nbonds
      rcon=50. nbxmod=-3 repel=0.8 cutnb=6. 
      rexp=2 irexp=2 inhibit=0.0 wmin=0.1 tolerance=0.5
   end
end

evaluate ($count=1) 
while ($count < 10 ) loop l1
   do (x=x+gauss(0.1)) ( all ) 
   do (y=y+gauss(0.1)) ( all ) 
   do (z=z+gauss(0.1)) ( all ) 
   minimize powell nstep=200 nprint=10 end
   evaluate ($count=$count+1)
end loop l1


evaluate ($accept=false) 
evaluate ($trial=1) 
while ($accept=false) loop accp

   for $1 in id ( tag ) loop resi
      igroup 
         interaction=( byresidue (id $1 ) and not name SG ) 
                     ( not name SG ) 
      end
      evaluate ($accept=true) 
      print thres=0.1 bonds
      if ($violations > 0) then
         evaluate ($accept=false) 
      end if
      print thres=10. angles 
      evaluate ($angles=$result)
      if ($violations > 0) then
         evaluate ($accept=false) 
      end if
      print thres=10. improper
      if ($violations > 0) then
         evaluate ($accept=false) 
      end if
      if ($accept=false) then
         do (x=x+gauss(0.3)) ( byresidue (id $1 ) ) 
         do (y=y+gauss(0.3)) ( byresidue (id $1 ) ) 
         do (z=z+gauss(0.3)) ( byresidue (id $1 ) ) 
      end if
   end loop resi

   igroup interaction=( all ) ( all ) end

   parameter
      nbonds
         rcon=50. nbxmod=-3 repel=3. cutnb=10. 
      end
   end
   flags exclude angle improper end
   
   minimize powell nstep=200 nprint=10 end

   parameter
      nbonds
         rcon=50. nbxmod=-3 repel=0.8 cutnb=6. 
      end
   end
   flags include angle improper end
   
   evaluate ($count=1) 
   while ($count < 5 ) loop l2
      do (x=x+gauss(0.05)) ( all ) 
      do (y=y+gauss(0.05)) ( all ) 
      do (z=z+gauss(0.05)) ( all ) 
      minimize powell nstep=200 nprint=10 end
      evaluate ($count=$count+1)
   end loop l2
   
   parameter
      nbonds
         rcon=50. nbxmod=3 repel=0.8 cutnb=6. 
      end
   end
   
   minimize powell nstep=300 nprint=10 end   
   minimize powell nstep=300 nprint=10 end

   igroup interaction=( not name SG ) ( not name SG ) end
   energy end

   evaluate ($accept=true) 

   print thres=0.05 bonds
   evaluate ($bonds=$result)
   if ($violations > 0) then
      evaluate ($accept=false) 
   end if

   print thres=10. angles 
   evaluate ($angles=$result)
   if ($violations > 0) then
      evaluate ($accept=false) 
   end if

   print thres=10. improper
   evaluate ($impr=$result)
   if ($violations > 0) then
      evaluate ($accept=false) 
   end if

   print thres=180. dihedral 
   evaluate ($dihe=$result)

   evaluate ($trial=$trial + 1) 
   if ($trial > &max_trial ) then
      exit loop accp
   end if

end loop accp

remarks extended strand(s) generation
remarks input molecular structure file=&structure_file 
remarks final rms deviations (excluding disulfide bonds): 
remarks    bonds=	 $bonds[F8.4] A  
remarks    angles=	 $angles[F8.4] degrees
remarks    impropers= $impr[F8.4] degrees
remarks    dihedrals= $dihe[F8.4] degrees (not used in some parameter sets!)
remarks final van der Waals (repel) energy=$vdw kcal/mole

write coordinates output=&output_coor end  

stop