!! FOR CONVEX.
!! Part 3 of the Nilges and Brunger method.
!! The electrostatic terms for side-chains is introduced slowly, while the
!! backbone electrostatics are constant throughout.
!! 5/11/92: scaled electrostatic on sidechains added
!! 6/11/92: outer loop of four as well, generating n*5 structures.
!! NOE for ca 18-24 replaced by 18-20 to again free c-termini
!! NOE FOR BACKBONE HELICITY AS WELL
!! 16/11/92: RESTART OPTION ADDED
!! CPU TIME ca. 4hrs per pentamer?
!! msps for a20 22.6.93
!! version with SIX noes per helix pair 28.6.93
!! for m2s1 19.8.93 

parameter           {* read the parameter file *}
   @/sansom/disk4/usr/xplor_3.0/toppar/param19x.proian
end

structure @m2i1.psf end

set echo = off message = all display = data.dyn.log precis = 4 end
display "temp" "tote" "totk" "ave_bond" "ave_angl" "ave_impr" "ave_vdw" "ave_elec" "ave_noe"

flags include noe end  

evaluate ($rand = 9731)            {* for generating seed *}
evaluate ($elecsc = 0.05)

! noe distance constraints ?????? hopefully. 1st effort 24/11/92
! if R < d then we want Enoe to be small, hence dminus = large
! if R > d then we want Enoe to be large, hence dplus = small
! pages 300-302 of the Xplor manual
noe
	reset
	average * center
	nres = 500
	potential * biharmonic
	scale * 2.5
	ceil = 50

end
	
noe @noe.intra.in end
parameter          
   nbonds
      repel=0                                    
      wmin=1.5
      rdie           {* constant/r-dependent dielectric, they say *}
      cutnb=9.5      {* cutnb >= ctofnb + 2(tole) *}
      switch         {* ??! *}
      ctonnb=5.0     {* switching *}
      ctofnb=9.0
      tolerance=0.25 {* quick nbon update default 0.5 *}
   end   
end

for $struct in ("m2i1a") loop outout

for $ann in (1 2 3 4 5) loop ann

for $dyn in (1) loop dyn

for $run in (1) loop run


evaluate ($structin = $struct + encode $ann + ".pdb"  )
evaluate ($structout = $struct + encode $ann + "d" + encode $dyn + ".pdb")
evaluate ($psfout = $struct + encode $ann + "d" + encode $dyn + ".psf")

   evaluate ($rand = $rand * 23)       {* set seed *}
   set seed $rand end 


evaluate ($init_t = 500)
evaluate ($cool = 50)
evaluate ($initr_ca = 20)
evaluate ($count = 0)

              

  coor init end
  coor disposition=reference init end

  coor @@$structin              {* read coors & ref for harmonic restraints *}
  coor disposition=reference @@$structin

! this vector statement moved here 5/11/92. Why reassign velocities every
! roun of dynamics?????

 vector do (vx = maxwell ($init_t)) (all)   {* assign velocties *}
 vector do (vy = maxwell ($init_t)) (all)
 vector do (vz = maxwell ($init_t)) (all)
 vector do (fbeta = 100) (all) {sets temperature frictional constant}

while ($count < 5) loop cooling

! vector statement moved inside dynamics loop 5/11/92, in order to gradually
! increase electrostatic effect of sidechain atoms
! setting these values to 0.3 had a disastrous effect 4/11/92

vector do (charge = ($elecsc * charge))
 ((resname ser) and (not (name ca or name c or name n or name o or name h )))
vector do (charge = ($elecsc * charge))
 ((resname his) and (not (name ca or name c or name n or name o or name h )))
vector do (charge = ($elecsc * charge))
 ((resname arg) and (not (name ca or name c or name n or name o or name h )))
! use the above for each polar sidechain when these are present

if ($count = 0) then
	evaluate ($elecsc = 1.7)
end if


 evaluate ($temp = ($init_t - ($count * $cool)))   {* temp *}
 evaluate ($rest_ca = ($initr_ca - ($count * 5.0)))  {* harmonic restraints *}

 vector do (harmonic=$rest_ca) (name ca and resid 1:29)  {* harmonic restraints on ca *}
 vector do (harmonic=0.0) (not (name ca and resid 1:29))
 constraints
 harmonic
 exponent=2
 end
!! NB "constraints" -> "restraints" in above for Xplor 21 -> 31
! flags include harmonic end - not actually needed, automatic with exp=2
! 9/11/92 now to strenthen the backbone-backbone hydrogen bonds we need to add a
! weak half-parabolic restraint with an energy constant of 2.5 kcals mol/
! Angstrom squared. A parameter-hbonds statement is not necessary

evaluate ($res = "res." + encode $count)
   dynamics  verlet
        isvfrq=800 		
!	restart=restart.file
!       nstep=800-no. when crash occurred
      nstep=800   timestep=0.001    iasvel=current   tcoupling=true 
      tbath=$temp  nprint=100  iprfrq=400 
   end

set message=on end
 display $temp $tote $totk $ave_bond $ave_angl $ave_impr $ave_vdw $ave_elec $ave_noe
set message = off end
 evaluate ($count = $count + 1)

end loop cooling

! release all harmonic restraints

flags exclude harmonic end
vector do (harmonic=0.0) (all)
 
! the Tconst dynamics follows. INstead of 5000 steps of 0.001 ps it is now
! 5 lots of 1000 steps of 0.001 ps.........

evaluate ($count=1)
	while ($count < 6) loop constant

evaluate ($res = "rest." + encode $count)
   dynamics  verlet
     isvfrq=1000 		
!	restart=restart.file
!       nstep=500-no. when crash occurred
      nstep=1000   timestep=0.001    iasvel=current   tcoupling=true 
      tbath=300    nprint=100  iprfrq=500 
   end
set message=on end
 display $temp $tote $totk $ave_bond $ave_angl $ave_impr $ave_vdw $ave_elec $ave_noe
set message = off end
evaluate ($count = $count + 1)

end loop constant

! here is the final   minimisation

 minim powell 
     nstep=1000 drop=40 nprint=100
 end

{====>}
  
write coor output=$structout end

write structure output = $psfout end

evaluate ($elecsc = 0.1197)

end loop run

end loop dyn

end loop ann

end loop outout
 stop