REMARKS Grid search for Da, and rhombicity 
REMARKS MKF 2/18/03


{* read in the PSF file and initial structure *}

structure @protG.psf @axis_600.psf end

{* parameter files for protein,and axis for dipolars. *}
param @TOPPAR:protein.par @TOPPAR:axes.par end

coor @axis_600.pdb


{* set the weights for the experimental energy terms. *}
    


evaluate ($da = -30.0)
evaluate ($lowestrms = 1000.0)
evaluate ($lowestda = -30.5)
evaluate ($lowestr = -0.05)
evaluate ($dacnt = 0)
evaluate ($rhom = 0.00)
evaluate ($ksani = 0.01)
sani
     nres=1600
     class JNHp
     force $ksani
     potential harmonic
     coeff 0.0 $da $rhom
     @bicelles_new_nh.tbl   {* dipolar coupling restraints for protein amide NH. *}
end 
     
while ($dacnt < 60) loop da
   evaluate ($dacnt = $dacnt + 1)
   evaluate ($rcnt = 0)
   evaluate ($rhom = 0.00)
   while ($rcnt<20) loop rhombicity 
      evaluate ($rcnt = $rcnt +1)




 
{* Read Structures 1 -> 10 *}
   evaluate ($no_struct = 1)
   evaluate ($totalrms_this_da_r = 0.0)
   while ($no_struct < 4) loop structures
      evaluate ($file= "compare_" + encode($no_struct) + ".pdb")
      coor @$file
      sani class JNHp force $ksani 
            coeff  0.0 $da $rhom   end


      flags 
        exclude *
        include sani 
      end
    vector do (mass  = 100.0) (not (resid 600) )   
                                {* all atom masses are set the same so that *}
                                {* motion of each atom is equally likely in *}
                                {* the dynamics. Note that the dynamics is  *}
                                {* not a realistic model of the motion of   *}
                                {* the system; it is a method of minimizing *}
                                {* a highly complicated function.           *}

    vector do (mass  = 30.0) (resid 600)
    vector do (fbeta = 10.0) (all)         {* coupling to heat bath *}

  dynamics internal
    reset
    itype=powell
    stepsize = 0.1
    nstep=1000
    depred=1
    fix (resid 1:400) ! NOTE: change this to the correct range for your protein
    group = (resid 600)    {* Treat the protein and the axis as rigid bodies *}
    etol = 0.0000001
    gtol = 0.0000001
    nprint= 1
  end

 
  sani print threshold=0.0 class JNHp end
   evaluate ($totalrms_this_da_r=$totalrms_this_da_r+$result)
      evaluate ($no_struct=$no_struct+1)

end loop structures

if ($totalrms_this_da_r < $lowestrms) then
    evaluate ($lowestrms = $totalrms_this_da_r)
    evaluate ($lowestda = $da)
    evaluate ($lowestr =$rhom)
end if
evaluate ($rhom = $rhom +0.05)
end loop rhombicity
evaluate ($da = $da + 0.5)
end loop da

Remark Da = $lowestda  R = $lowestr

write coor output=test.pdb end

stop