define(
{* selection of atoms other than hydrogens for which coordinates
   will be generated *}
{* to generate coordinates for all unknown atoms use: (not(known)) *}
{===>} atom_build=(not(known));
{============================= output files ================================}

{* input coordinate file *}
{===>} coordinate_infile="1xxx_9.pdb";

{* output coordinate file *}
{===>} coordinate_outfile="1xxx_9_cns.pdb";
)

 checkversion 1.1

 evaluate ($log_level=verbose)
 evaluate ($par_nonbonded="PROLSQ")

 !@generateProtonsJFD.inp

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

remarks changed Tue Sep  2 09:48:07 CDT 2003
remarks by jfd to include a fes residue

!@generate_tmoc.inp

topology reset end
structure reset end

topology 
{===>} @TOPPAR:topallhdg5.3.pro
end

! no chain id.
segment  name="    " 
    chain 
        @TOPPAR:topallhdg5.3.pep
        sequence
        	@1xxx.seq
	end
    end
end

! removes the hydrogen on the cys-en
! and changes the atom type from SH1E to S.
patch DISU reference=1=( resid 50 )
           reference=2=( resid 53 ) 
           end

! turn peptide from trans to cis
patch CISP reference=nil=( resid 92 ) end

coor @&coordinate_infile
parameter @TOPPAR:parallhdg5.3.pro end

!!!!!!!!!!!!!!!! Start of atom generation.

show sum(1) ( not(hydrogen) and not(known) )
 if ( $select = 0 ) then
   display  %INFO: There are no coordinates missing for non-hydrogen atoms
 end if
 set echo=on end

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


 identity (store1) (none)

 identity (store1) (&atom_build) 
 identity (store1) (store1 or hydrogen)
 
 show sum(1) (store1)
 evaluate ($tobuild=$result)

 !evaluate ($tobuild=0)

 if ( $tobuild > 0 ) then

   fix selection=(not(store1)) end

   show sum(1) (store1)
   evaluate ($moving=$result)

   if ( $moving > 0 ) then
     for $id in id (tag and byres(store1)) loop avco

       show ave(x) (byres(id $id) and known)
       evaluate ($ave_x=$result)
       show ave(y) (byres(id $id) and known)
       evaluate ($ave_y=$result)
       show ave(z) (byres(id $id) and known)
       evaluate ($ave_z=$result)

       do (x=$ave_x) (byres(id $id) and store1)
       do (y=$ave_y) (byres(id $id) and store1)
       do (z=$ave_z) (byres(id $id) and store1)
 
     end loop avco 

     do (x=x+random(2.0)) (store1)
     do (y=y+random(2.0)) (store1)
     do (z=z+random(2.0)) (store1)

     {- start parameter for the side chain building -}
     parameter
       nbonds
         rcon=20. nbxmod=-2 repel=0.9  wmin=0.1 tolerance=1.
         rexp=2 irexp=2 inhibit=0.25
       end
     end

     {- Friction coefficient, in 1/ps. -}
     do (fbeta=100) (store1)

     evaluate ($bath=300.0)
     evaluate ($nstep=500)
     evaluate ($timestep=0.0005)

     do (refy=mass) (store1)

     do (mass=20) (store1)

     igroup interaction 
       (store1) (store1 or known)
     end

     {- turn on initial energy terms -}
     flags exclude * include bond angle vdw end
 
     minimize powell nstep=50  nprint=10 end

     do (vx=maxwell($bath)) (store1)
     do (vy=maxwell($bath)) (store1)
     do (vz=maxwell($bath)) (store1)

     flags exclude vdw include impr end

     dynamics cartesian
       nstep=50
       timestep=$timestep
       tcoupling=true temperature=$bath
       nprint=$nstep
       cmremove=false
     end

     flags include vdw end

     minimize powell nstep=50 nprint=10 end

     do (vx=maxwell($bath)) (store1)
     do (vy=maxwell($bath)) (store1)
     do (vz=maxwell($bath)) (store1)

     dynamics cartesian
       nstep=50
       timestep=$timestep
       tcoupling=true temperature=$bath
       nprint=$nstep
       cmremove=false
     end

     parameter
       nbonds
         rcon=2. nbxmod=-3 repel=0.75
       end
     end

     minimize powell nstep=100 nprint=25 end

     do (vx=maxwell($bath)) (store1)
     do (vy=maxwell($bath)) (store1)
     do (vz=maxwell($bath)) (store1)

     dynamics cartesian
       nstep=$nstep
       timestep=$timestep
       tcoupling=true temperature=$bath
       nprint=$nstep
       cmremove=false
     end

     {- turn on all energy terms -}
     flags include dihe ? end

     {- set repel to ~vdw radii -}
     parameter
       nbonds
         repel=0.89
       end
     end

     minimize powell nstep=500 nprint=50 end

     flags exclude * include bond angl impr dihe vdw end

     {- return masses to something sensible -}
     do (mass=refy) (store1)

     do (vx=maxwell($bath)) (store1)
     do (vy=maxwell($bath)) (store1)
     do (vz=maxwell($bath)) (store1)

     dynamics cartesian
       nstep=$nstep
       timestep=$timestep
       tcoupling=true temperature=$bath
       nprint=$nstep
       cmremove=false
     end

     {- some final minimisation -}
     minimize powell
       nstep=500
       drop=40.0
       nprint=50
     end

     print thres=0.02 bonds
     print thres=5. angles

   end if
  
   fix selection=( none ) end

 end if

 set echo=false end
 show sum(1) (not(known))
 if ( $result < 100 ) then
   for $id in id (not(known)) loop print
     show (segid) (id $id)
     evaluate ($segid=$result)
     show (resname) (id $id)
     evaluate ($resname=$result)
     show (resid) (id $id)
     evaluate ($resid=$result)
     show (name) (id $id)
     evaluate ($name=$result)
     buffer message
       display unknown coordinates for atom: $segid[a4] $resname[a4] $resid[a4] $name[a4]
     end
   end loop print
 else
   buffer message
     display unknown coordinates for more than 100 atoms
   end
 end if
 set echo=true end

 if (&set_bfactor=true) then
   do (b=&bfactor) ( all )
 else
   show ave(b) (known and not(store1))
   do (b=$result) (store1 and (attr b < 0.01))
 end if

 if (&set_occupancy=true) then
   do (q=&occupancy) ( all )
 end if

 set echo=false end
 show sum(1) (store1)
 if ( $result < 100 ) then
   for $id in id (store1) loop print
     show (segid) (id $id)
     evaluate ($segid=$result)
     show (resname) (id $id)
     evaluate ($resname=$result)
     show (resid) (id $id)
     evaluate ($resid=$result)
     show (name) (id $id)
     evaluate ($name=$result)
     buffer message
       display coordinates built for atom: $segid[a4] $resname[a4] $resid[a4] $name[a4]
     end 
   end loop print
 else
   buffer message
     display coordinates built for more than 100 hundred atoms
   end
 end if
 set echo=true end

 set remarks=reset end

 buffer message
   to=remarks
   dump
 end
 
   write coordinates output=&coordinate_outfile end

stop