# Copyright (c) 2002-08 Peter Guntert. All rights reserved.
## 7MACROS: anneal - CYANA macro
##
## Parameters: thigh=<real>                   (default: 8.0)
##             steps=<integer>                (default: 10000)
##             highsteps=<integer>            (default: -1)
##             minsteps=<integer>             (default: 1000)
##             relax
##
## Performs simulated annealing on the current structure with a total of N MD
## steps, starting with Nhigh MD steps at temperature Thigh followed by slow
## cooling during N - Nhigh MD steps to a final temperature of Tend. Finally,
## n steps of conjugate gradient minimization are added. The temperature is
## measured in target function units per degree of freedom. Optionally, more
## minimization can be performed in order to relax strong overlaps and
## restraint violations prior to the start of the MD calculation. The relax
## option can be useful for larger (above 200 residues) proteins if otherwise
## the maximal length of the pair list would be exceeded.

var echo lev step i mr t0 t v ahigh aend it steps1 steps2 steps3

syntax $macro:thigh=0.0<@r=8.0 steps=@i=10000 highsteps=@i=-1 minsteps=@i=1000 relax

echo:=off
t0=cputime
if (highsteps.lt.0) highsteps=steps/5

# initialize rdc force constants
if (rdc_vts_on < 2) rdc_vts_on=-1
if (rdc_vas1_on < 2) rdc_vas1_on=-1
if (rdc_vas2_on < 2) rdc_vas2_on=-1
if (pcs_ws_on < 2) pcs_ws_on=-1

# Minimization to remove strong overlaps

atom set HYDROGEN radius=-1.0 info=none
atom set HEAVY radius=radius-0.2 info=none
atom set HBOUND radius=radius+0.15 info=none
weight_vdw=0.5
aco_type=2
mr=rnum(nr)-rnum(1)
if (relax) then
  lev:=1,3,10,20,50,100,150,200,300,600,1000,2000,100000
  do i 2 length('lev')
    if (lev(i).gt.mr) break
    minimize 100 level=lev(i-1)
  end do
else
  minimize 100 level=3
end if
level=mr
minimize 100

i=$val('structure','0')
if (i.gt.0) seed=17453+17*i
v=50

tmed=thigh*0.05
ahigh=0.005; aend=0.0001
#amed=(aend/ahigh)**(1.0/3.0)
#print "$ahigh , $amed , $aend"
steps3=200
steps1=nint((steps-steps3)*0.33)
steps2=steps-steps1-steps3
#print "$steps1 + $steps2 + $steps3 = $steps"

if (steps1.le.100 .or. steps2.le.100) error "Not enough timesteps."

# First simulated annealing stage

step=0.02
t=real(highsteps)/steps
md steps=steps1 \
   temperature=min(thigh,(thigh-tmed)*(1.0-(s-t)/(1.0-t))**4+tmed) \
   accuracy=ahigh \
   dt=step tau=20.0 angdev=36000.0 vdwupdate=v nprint=steps/20

# Second simulated annealing stage

atom set HEAVY radius=* info=none
minimize 50
step=0.5*timestep
t=0.67

if (rdc_vts_on < 2) rdc_vts_on=1
if (rdc_vas1_on < 2) rdc_vas1_on=1
#if (rdc_vas2_on < 2) rdc_vas2_on=1
if (pcs_ws_on < 2) pcs_ws_on=1
md steps=steps2/2 \
   temperature=tmed*(1.0-s*t)**4 \
   accuracy=ahigh*(aend/ahigh)**(s*t) \
   dt=step tau=20.0 angdev=36000.0 vdwupdate=v nprint=steps/20 continue
if (pcs_ws_on < 2) pcs_ws_on=0
if (rdc_vts_on < 2) rdc_vts_on=0
if (rdc_vas1_on < 2) rdc_vas1_on=0
#if (rdc_vas2_on < 2) rdc_vas2_on=0


atom set * radius=* info=none
minimize 50
step=0.5*timestep

if (rdc_vas2_on < 2) rdc_vas2_on=1
md steps=steps2/2 \
   temperature=tmed*(1.0-(t+s*(1-t)))**4 \
   accuracy=ahigh*(aend/ahigh)**(t+s*(1-t)) \
   dt=step tau=20.0 angdev=36000.0 vdwupdate=v nprint=steps/20 continue
if (rdc_vas2_on < 2) rdc_vas2_on=0

# Final MD and minimization with all atoms

aco_type=1
weight_vdw=2.0
minimize 50
step=0.5*timestep
md steps=steps3 dt=step temperature=0.0 accuracy=aend tau=20.0 \
   angdev=36000.0 vdwupdate=v nprint=100 continue
minimize minsteps

print "    Structure annealed in $nint(cputime-t0) s, f = $tf."