REMARK REFINE.INP
REMARK Author: Michael Nilges
REMARK Ref: Nilges, M. (1993), A calculation strategy for the structure
REMARK      determination of symmetric dimers by 1H NMR.
REMARK      Proteins 17, 297-309.
REMARK
REMARKS structure determination for dimeric proteins

evaluate ($seed  = 4367578) 
set seed $seed  end

!----------------------------------------------------------------------
! read in the PSF file and initial structure

@psf_bld_monomer.inp

duplicate
  SELEction=(name *)
  SEGId=B
end

coor disp=comp @tr2_min_avg.pdb

flags exclude ELEC DIHE include NOE CDIH NCS end


evaluate ($cool_steps = 10000)
evaluate ($init_t  = 2000.01)

vector do (mass  = 100.0) (all)         ! uniform mass for all atoms
vector do (fbeta = 10.0) (all)         ! coupling to heat bath


evaluate ($count = 0)
while ($count < 50) loop stru
  evaluate ($count = $count + 1)

evaluate ($file = "random_" + encode($count) + ".pdb")
coor @@$file

evaluate ($ini_rad  = 0.9)        evaluate ($fin_rad  = 0.78)
evaluate ($ini_con=  0.003)       evaluate ($fin_con=  4.0)
evaluate ($ini_ncs = 0.01)          evaluate ($fin_ncs = 10)
evaluate ($ini_ang = 0.4)         evaluate ($fin_ang = 1.0)
evaluate ($ini_imp = 0.1)         evaluate ($fin_imp = 1.0) 
evaluate ($ini_cdi = 10)          evaluate ($fin_cdi = 200) 
evaluate ($ini_noe = 2.0)           evaluate ($fin_noe = 25.0)

evaluate ($kcdi  = $ini_cdi)   ! torsion angles

evaluate ($knoe  = $ini_noe)   ! slope of NOE potential
evaluate ($kncs  = $ini_ncs)

ncs restraints
 initialize
 group 
  equi (segid A)
  equi (segid B)
  weight = $kncs
 end
 ?
end
noe
  reset
  nrestraints = 1500                   ! allocate space for NOEs
  ceiling 100
  
  class      all @@original_inter.tbl
                 @@original_intra.tbl
                 @@tr2_calcium.tbl

  averaging  all sum
  monomers   all 2
  potential  all soft
  scale      all $knoe                           
  sqconstant all 1.0
  sqexponent all 2
  soexponent all 1
  rswitch    all 1.0
  sqoffset   all 0.0
  asymptote  all 2.0
  

  class      symm @@symmetry.tbl

  potential  symm symmetry
  scale      symm 1.0                           
  sqconstant symm 1.0
  sqexponent symm 2
  soexponent symm 1
  rswitch    symm 0.5
  sqoffset   symm 0.0
  asymptote  symm 1.0

end

@@dihed.tbl                    

restraints dihed 
  scale $kcdi  
end

parameters
  nbonds
    atom
    nbxmod 3
    wmin  =   0.01  ! warning off
    cutnb =   4.5   ! nonbonded cutoff
    tolerance 0.5
    repel=    1.2   ! scale factor for vdW radii = 1 ( L-J radii)
    rexp   =  2     ! exponents in (r^irex - R0^irex)^rexp
    irex   =  2
    rcon=1      ! actually set the vdW weight
  end
end 


parameter nbonds nbxmod 4 cutnb 100 tolerance 45 end end

constraints 
        interaction (not name ca) (all) 
        weights * 1 angl 0.4 impr 0.1 vdw 0 elec 0 end 
        interaction (name ca) (name ca) 
        weights * 1 angl 0.4 impr 0.1 vdw 0.1 end 
end

set message off echo off end


dynamics  verlet
      nstep=10000           ! 
      timestep=0.005        ! 60 ps of dynamics
      iasvel=maxwell        firsttemp= $init_t
      tcoupling = true
      tbath = $init_t   
      nprint=50
      iprfrq=0
      ntrfr = 99999999
end






evaluate ($final_t = 100)     { K }
evaluate ($tempstep = 50)     { K }

evaluate ($ncycle = ($init_t-$final_t)/$tempstep)
evaluate ($nstep = int($cool_steps/$ncycle))



evaluate ($bath  = $init_t)
evaluate ($k_vdw = $ini_con)
evaluate ($k_vdwfact = ($fin_con/$ini_con)^(1/$ncycle))
evaluate ($radius=    $ini_rad)
evaluate ($radfact = ($fin_rad/$ini_rad)^(1/$ncycle))
evaluate ($ncs_fac = ($fin_ncs/$ini_ncs)^(1/$ncycle))
evaluate ($k_ncs = $ini_ncs)
evaluate ($k_ang = $ini_ang)
evaluate ($ang_fac = ($fin_ang/$ini_ang)^(1/$ncycle))
evaluate ($k_imp = $ini_imp)
evaluate ($imp_fac = ($fin_imp/$ini_imp)^(1/$ncycle))
evaluate ($noe_fac = ($fin_noe/$ini_noe)^(1/$ncycle))
evaluate ($knoe = $ini_noe)
evaluate ($cdi_fac = ($fin_cdi/$ini_cdi)^(1/$ncycle))
evaluate ($kcdi = $ini_cdi)


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

evaluate ($i_cool = 0)
while ($i_cool < $ncycle) loop cool
    evaluate ($i_cool=$i_cool+1)


    restraints dihed 
      scale $kcdi  
    end

    constraints interaction (all) (all) weights
        * 1 angles $k_ang improper $k_imp
    end end
    ncs restraints
        initialize
        group equi (segid A)  equi (segid B) weight = $k_ncs end
        ?
    end
    parameter nbonds nbxmod 3
        cutnb=4.5 tolerance 0.5 rcon=$k_vdw nbxmod=-3 repel=$radius
    end       end
    noe scale all $knoe end

    dynamics  verlet
        nstep=$nstep  timestep=0.005 iasvel=current  firsttemp= $bath
        tcoupling = true  tbath = $bath   nprint=$nstep  iprfrq=0
        ntrfr = 99999999
    end

    evaluate ($bath  = $bath  - $tempstep)
    evaluate ($k_vdw=min($fin_con,$k_vdw*$k_vdwfact))
    evaluate ($radius=max($fin_rad,$radius*$radfact))
    evaluate ($k_ncs = $k_ncs*$ncs_fac)
    evaluate ($k_ang = $k_ang*$ang_fac)
    evaluate ($k_imp = $k_imp*$imp_fac)
    evaluate ($knoe  = $knoe*$noe_fac)
    evaluate ($kcdi  = $kcdi*$cdi_fac)
end loop cool


mini powell nstep= 250 end
coor sele=(name ca and resid 128:160) fit end
coor sele=(name ca and resid 128:160) rms end
evaluate ($rmsdx = $result)
coor rms end
vector do (b=rmsd) (all) 

remarks weight on cdih $kcdi
remarks initial weight on NOE $ini_noe
remarks final rmsd to NMR $rmsdx
remarks final noe and cdi energies $noe $cdih
evaluate ($file = "refine_" + encode($count) + ".pdb")
write coor output= $file end

end loop stru
stop