!! ANNEAL.IN - SUITABLE ANNEALING FOR ANY HELIX
!! copied from IDK, msps 21.6.93

evaluate ($init = 1000) 	{ sets the initial temperature }
set message = on display = data.ann.log end
evaluate ($timetot = 0)
evaluate ($1="time     geom        impr         dihe        vdw")
display $1

! read the paramete file and structure file

parameter
	@/sansom/disk4/usr/xplor_3.0/toppar/param19x.proian
end


structure
	@m2i1.psf
end

! set up energy parameters for inclusion or exclusion
! note the exclusion of the electrostatic energy term. Do not want
! the growing helix to explode. Bang!

set seed = 1297 end

flag
	exclude *		{ turns all energy terms off }
	include bonds angle impr vdw dihe { except these }
end


! STEP ZERO ====================================================
! here is a loop that performs most of the functions 


for $struct in ("m2i1a") loop outer		{ end of loop is the last line }

	evaluate ($end_count = 5)
	evaluate ($count = 0)

	while ($count < $end_count) loop main	{ end loop is penultimate line }
! remember. Vector statements alter any character of the molecule
! here it is used to set the friction coefficient and uniform heavy
! masses for the molecular dynamics calculations. The use of a uniform
! heavy mass speeds up the calculation.

		vector do (fbeta=100) (all)	{ fbeta is in manual }
		vector do (mass=10) (all)

		evaluate ($count=$count+1)
		evaluate ($kba = 0.001)	{ bond scales }
		evaluate ($impsc = 0.001) {improper scale}
		evaluate ($vdwsc = 0.000)	{ vdw scale }
		evaluate ($dihsc = 0.001)	{ dihedral scale }

		evaluate ($filpdb = $struct  +  encode($count) + ".pdb")
		evaluate ($filpsf = $struct  +  encode($count) + ".psf")
		
		evaluate ($structin = "m2i1.pdb")
		coor @@$structin 

! now we need to set up the parameters for the molecular dynamics
! calculations. The following is annotated from section 8.2 of the
! xplor manual. Since the attractive term of the VDW energy is turned off,
! the atoms would tend to move further apart than they would normally be.
! Hence the repel parameter is set to 0.9 to scale down the repulsive force.
! The rexp value of 2 introduces a quartic repulsive term for the VDW 
! interactions allowing atoms to pass through one another.

parameter			{ invokes a parameter function }
	nbonds			{ the non-bonded parameters will be set }
		repel=1.0	{ see  above }
		rexp=2
		irexp=2	 	{ see above }
		rcon=1.0	{ energy constant for repel function }
		nbxmod=-2	{ exclusion list options. -2 = excludes 1-2 }
		wmin=0		{ close contact warning distance. }
		cutnb=4.0	{ non-bonded interaction cutoff distance } 
		ctonnb=2.99	{ r at which switching becomes active }
		ctofnb=3.0	{ r at which sw or sh function forces E = 0 }
		tolerance=0.5	{ r an atom may move before list is updated }
				{ cutnb should be > 2*tol + ctofnb }
	end
end

	
		cons 
			fix=(name ca) 
		end 		{ fix ca's to retain helix }


! PART ONE ==============================================================
!  AA..here is the initial minimisation step
! note line removed here that reduced furtther the energy term weights. Note
! also the inclusion of a dihedral scaling in the weights term.

!  			 constraints interaction
!				(all) (all) weights bond $kba angl $kba dihe $dihsc impr $impsc vdw $vdwsc end 
!  			 end
!
 		minimize powell 
 			nstep = 5
 			npri  = 5
 			drop  = 50.0 
		end

! This next section determines the "ideal" thermal energy of a system with
! 3*natom-6 degrees  of freedom at the initial temperature $init
! Ekin = (3*natom-6) * kboltz * temp / 2

		vector show ave ( x ) (all)
		evaluate ($ekin = (3*$select - 6) * $kboltz * $init / 2.0) 
		
! Now the bond and angle scales are set such that the bond and angle energy
! terms are equal to the kinetic energy at 1000K.

		energy end
		evaluate ($epot = $bond + $angl + $impr + $dihe)	{ pot energy }
		evaluate ($ba = $bond + $angl)
		evaluate ($escale = $ekin/$ba)
		evaluate ($kba = $escale)
		evaluate ($impsc = $kba / 1000)
		evaluate ($dihsc = $kba / 1000)

  		constraints interaction
  			(all) (all) weights bond $kba angl $kba dihe $dihsc impr $impsc vdw $vdwsc end 
  		end
! check now that $bond and $angle are equal to ekin

			energy end
! additional lines to set the improper and dihedral scales as follows and as
! detailed in Brungers paper


! PART TWO ==============================================================
! Now comes the interesting bit ie. the high temperature molecular dynamics
! vx,vy,vz are x,y, and z components of current molecular velocities
! maxwell is the Maxwellian random number distribution function. $init
! its argument. The two stage boiling is NOT identical to that of Nilges and 
! Brunger BUT contrasts that which Dave Doak supplied.
! in this section the terms for bond and angle are multiplied by 1.2 per cycle
! and the improper and dihedral terms by 1.1 squared

		vector do (vx=maxwell($init)) (all)
		vector do (vy=maxwell($init)) (all)
		vector do (vz=maxwell($init)) (all)

		while ($kba < 10.) loop stage1
	
				evaluate ($kba = $kba * 1.25)

! next three lines are new from dtox.annneal.in
			
			if ($impsc < $kba) then
				evaluate ($impsc = $impsc * 1.25 ** 2)
				evaluate ($dihsc = $dihsc * 1.25 ** 2)
			end if

! AA here is the set up section for the dynamics calculations. see xplor manual
! chapter 12 for more details. Iasvel can be set to maxwellian, uniform or current
 			evaluate ($a = $kba/1000)
 			evaluate ($i = $impsc/1000)
 			evaluate ($d = $dihsc/1000)

   			constraints interaction
 			    (all) (all) weights bond $a angl $a dihe $d impr $i end 
   			end

			dynamics verlet		
!				isvfrq=100 save=restart.file		
!				restart=restart.file
!				nstep=100-no. when crash occurred
				nstep=100	
				timestep=0.001	
				iasvel=current  
				tcoupling=true	
				tbath=$init
				nprint=100	
				iprfrq=100	
			end
		evaluate ($timetot = $timetot + 0.1)
		evaluate ($2 = " " + encode($impsc) + " " + encode($dihsc) + " " + encode($vdwsc))
		evaluate ($1 = " " + encode($timetot) + " " + encode($kba))
		display $1 $2 
		end loop stage1

! STEP TWO ==============================================================
! introducing a small vdw scaling factor
! in this section the terms for bond and angle are multiplied by 1.1 per cycle
! and the improper and dihedral terms by 1.1 squared
		evaluate ($vdwsc = 0.0001)
		while ($kba < 500.) loop stage1a
	
			evaluate ($kba = $kba * 1.1)
			evaluate ($impsc = $impsc * 1.1 ** 2)
			evaluate ($dihsc = $dihsc * 1.1 ** 2)
			if ($kba > 500) then
				evaluate ($kba = 500)
			end if

			if ($vdwsc < 0.003) then
				evaluate ($vdwsc = $vdwsc * 2.)
			else 
				evaluate ($vdwsc = $vdwsc * 1.25)
			end if

			if ($vdwsc > 0.1) then
				evaluate ($vdwsc = 0.1)
			end if

			if ($impsc > 200 ) then
				evaluate ($impsc = 200.0)
			end if
			if ($dihsc > 200 ) then
				evaluate ($dihsc = 200.0)
			end if

 			evaluate ($a = $kba/1000)
 			evaluate ($i = $impsc/1000)
 			evaluate ($d = $dihsc/1000)

   			constraints interaction
 			    (all) (all) weights bond $a angl $a dihe $d impr $i vdw $vdwsc end 
   			end


			dynamics verlet
!				isvfrq=800 save=restart.file		
!				restart=restart.file
!				nstep=800-no. when crash occurred
				nstep=800
				timestep=0.001	
				iasvel=current  
				tcoupling=true	
				tbath=$init
				nprint=800	
				iprfrq=800	
			end
		evaluate ($timetot = $timetot + 0.8)
		evaluate ($2 = " " + encode($impsc) + " " + encode($dihsc) + " " + encode($vdwsc))
		evaluate ($1 = " " + encode($timetot) + " " + encode($kba))
		display $1 $2 
		end loop stage1a

! STEP THREE ==============================================================
! cooling the system down again

		parameter
			nbonds			{ reset such that 1-2 and 1-3 interactions are }
				nbxmod=-3	{ not counted 1-3 determined }
				repel = 0.8	{s factor in Brunger}
			end			{ by bond angles }
		end

		evaluate ($tempstep=25)		
		evaluate ($ncycle=28)	
		evaluate ($bath = $init)
		evaluate ($icool = 0)
		evaluate ($vdwsc = 4)
		evaluate ($kba = 500)
		evaluate ($impsc = 200)
		evaluate ($dihsc = 200)

			while ($icool < $ncycle ) loop cooling
				evaluate ($icool = $icool + 1)

				evaluate ($bath = $bath - $tempstep)
 			evaluate ($a = $kba/1000)
 			evaluate ($i = $impsc/1000)
 			evaluate ($d = $dihsc/1000)

   			constraints interaction
 			    (all) (all) weights bond $a angl $a dihe $d impr $i vdw $vdwsc end 
   			end


				dynamics verlet
!				isvfrq=100 save=restart.file		
!				restart=restart.file
!				nstep=100-no. when crash occurred
					nstep = 300
					time = 0.001
					iasvel=current
					firstt=$bath
					tcoupling=true
					tbath=$bath
					nprint=150
					iprfrq=300
				end
			evaluate ($timetot = $timetot + 0.3)

			end loop cooling


! STEP FOUR ===================================================
! THE FINAL MINIMISATION STEP.
! the s factor in the vdw repulsive term is set to zero so switching on 
! a full lennard-jones potential

		parameter
			nbonds
				repel = 0
			end
		end

		constraints interaction
  			 	(all) (all) weights bond $a angl $a dihe $d impr $i vdw $vdwsc end 
  			 end

		minimize powell
			nstep = 2000
			drop = 10.0	{ expected initial drop in energy }
			nprint = 100
		end

! STEP FIVE ===================================================
! WRITE OUT THE FINAL STRUCTURE!!!!!!!!
! threshold sets the min difference of the property before it will be reported
! see xplor 10.1.1

		print threshold = 10. dihe	{ puts answers in $result and $violations }
		evaluate($rms_dihe = $result)
		evaluate($violations_dihe = $violations)
		print threshold = 0.05 bonds
		evaluate($rms_bonds = $result)
		print threshold = 10. angles
		evaluate ($rms_angles = $result)
		print threshold = 10. impropers
		evaluate ($rms_impropers = $result)


! what is the final filename 

		write coordinates output = $filpdb end
		write structure output = $filpsf end
	end loop main  { see several pages above }

end loop outer	{ see several pages above }
	
stop