{+ file: model_anneal.inp +}
{+ directory: general +}
{+ description: simulated annealing/molecular dynamics +}
{+ authors: Axel T. Brunger, Luke M. Rice and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: L.M. Rice and A.T. Brunger, Torsion Angle Dynamics: 
              Reduced Variable Conformational Sampling Enhances 
              Crystallographic Structure Refinement, Proteins: Structure, 
              Function, and Genetics, 19, 277-290 (1994) +}

{- Guidelines for using this file:
   - all strings must be quoted by double-quotes
   - logical variables (true/false) are not quoted
   - do not remove any evaluate statements from the file 
   - the selections store1 through store3 are available for general use -}

{- begin block parameter definition -} define(

{======================= molecular structure =========================}

{* structure file *}
{===>} structure_infile="amy_hydrogen.mtf";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein.param";
{===>} parameter_infile_2="";
{===>} parameter_infile_3="";
{===>} parameter_infile_4="";
{===>} parameter_infile_5="";

{* coordinate file *}
{===>} coordinate_infile="amy_hydrogen.pdb";

{===================== crystallographic symmetry =====================}

{* use crystallographic symmetry *}
{+ choice: true false +}
{===>} use_cryst=true;

{* space group *}
{* use International Table conventions with subscripts substituted
   by parenthesis *}
{===>} sg="P2(1)2(1)2(1)";

{* unit cell parameters in Angstroms and degrees *}
{+ table: rows=1 "cell" cols=6 "a" "b" "c" "alpha" "beta" "gamma" +}
{===>} a=61.76;
{===>} b=40.73;
{===>} c=26.74;
{===>} alpha=90;
{===>} beta=90;
{===>} gamma=90;

{=================== non-crystallographic symmetry ===================}

{* NCS-restraints/constraints file *}
{* see auxiliary/ncs.def *}
{===>} ncs_infile="";

{========================== atom selection ===========================}

{* select atoms to be included *}
{* it is essential to include hydrogen atoms if a free MD simulation is
   being performed *}
{* this should include all conformations if multiple conformations are used *}
{===>} atom_select=(known);

{* select fixed atoms *}
{* note: isolated atoms and atoms are diatomic molecules are automatically 
   fixed during torsion angle dynamics. So, you don't have to 
   explicitly fix them here. *}
{===>} atom_fixed=(none);

{* select atoms to be harmonically restrained *}
{===>} atom_harm=(none);

{* harmonic restraint constant - for harmonically restrained atoms *}
{===>} k_harmonic=10;

{* select atoms to be treated as rigid groups during torsion angle dynamics *}
{===>} atom_rigid=(none);

{* select atoms in alternate conformation 1 *}
{===>} conf_1=(none);

{* select atoms in alternate conformation 2 *}
{===>} conf_2=(none);

{* select atoms in alternate conformation 3 *}
{===>} conf_3=(none);

{* select atoms in alternate conformation 4 *}
{===>} conf_4=(none);

{* additional restraints file *}
{* eg. auxiliary/dna-rna_restraints.def *}
{===>} restraints_infile="";

{====================== annealing parameters ========================}

{* type of molecular dynamics *}
{+ choice: "torsion" "cartesian" +}
{===>} md_type="cartesian";

{* annealing schedule *}
{+ choice: "slowcool" "constant" +}
{===>} md_scheme="constant";

{* starting temperature *}
{* used for both constant-temperature and slowcooling schemes *}
{===>} temperature=298;

{* temperature control method *}
{* either coupling to a temperature bath or velocity scaling *}
{+ choice: coupling scaling +}
{===>} tcontrol="scaling";

{* number of molecular dynamics steps *}
{* only used for constant-temperature annealing schedule *}
{===>} constant_steps=1000;

{* drop in temperature (K) per cycle of dynamics *}
{* only used for slowcooling annealing schedule *}
{===>} cool_rate=25;

{* molecular dynamics time step (ps) *}
{===>} time_step=0.0005;

{* nonbonded cutoff (Angstroms) *}
{===>} nonb_cutoff=13;

{* dielectric constant *}
{===>} dielectric=1;

{* number of trials to carry out with different initial velocities *}
{===>} num_trials=1;

{* frequency of writing trajectory (in steps) *}
{* this only applies to the constant temperature option *}
{* a trajectory will not be written if this value is 0 or less *}
{===>} traj_freq=100;

{* seed for random number generator *}
{* change to get different initial velocities *}
{===>} seed=82364;

{* torsion-angle MD parameters *}
{* increase these values if the program terminates with the message 
   that one of these parameters is exceeded *}

{* maximum unbranched chain length *}
{* increase for long stretches of polyalanine *}
{===>} torsion_maxlength=50;
{* maximum number of distinct bodies *}
{===>} torsion_maxtree=15;
{* maximum number of chains (increase for large molecules) *}
{===>} torsion_maxchain=1000;
{* maximum number of bonds to an atom *}
{===>} torsion_maxbond=6;

{=========================== output files ============================}

{* root name for output files *}
{+ list:
        coordinate files will be written: <output_root>_n.pdb 
        trajectory file will be written: <output_root>_n.crd +}
{===>} output_root="model_anneal";

{===========================================================================}
{         things below this line do not normally need to be changed         }
{         except for the torsion angle topology setup if you have           }
{         molecules other than protein or nucleic acid                      }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates @&coordinate_infile

 if ( &use_cryst = true ) then
   xray
     @@CNS_XTALLIB:spacegroup.lib (sg=&sg;
                                   sgparam=$sgparam;)
     a=&a b=&b c=&c  alpha=&alpha beta=&beta gamma=&gamma
   end
 end if

 parameter
   if ( &BLANK%parameter_infile_1 = false ) then
      @@&parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      @@&parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      @@&parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      @@&parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      @@&parameter_infile_5
   end if
 end

 if ( &use_cryst = true ) then
   flags exclude * include bond angle impr dihe vdw elec pvdw pele ? end
 else
   flags exclude * include bond angle impr dihe vdw elec ? end
 end if

 param nbonds tolerence=? end end
 evaluate ($toler=$result)

 evaluate ($ctofnb=&nonb_cutoff-(2*$toler))
 evaluate ($ctonnb=$ctofnb-1.0)

 param 
   nbonds 
     cutnb=&nonb_cutoff
     ctonnb=$ctonnb
     ctofnb=$ctofnb
     eps=&dielectric
     ?
   end
 end

 if ( &BLANK%ncs_infile = false ) then
    inline @&ncs_infile
 end if

 if ( &BLANK%restraints_infile = false ) then
     @&restraints_infile
 end if

 do (store6=0) (all)

 evaluate ($nalt=1)
 evaluate ($alt=1)
 evaluate ($done=false)
 while ( $done = false ) loop nalt
   if ( &exist_conf_$alt = true ) then
     show sum(1) ( &conf_$alt )
     if ( $result > 0 ) then
       evaluate ($nalt=$nalt+1)
     end if
   else
     evaluate ($done=true)
     evaluate ($nalt=$nalt-1)
   end if
   evaluate ($alt=$alt+1)
 end loop nalt

 evaluate ($alt=1)
 while ( $alt <= $nalt ) loop alt
   do (store6=$alt) ( &conf_$alt )
   evaluate ($alt=$alt+1)
 end loop alt

 igroup
   interaction ( &atom_select and not(attr store6 > 0))
               ( &atom_select and not(attr store6 > 0))
   evaluate ($alt=1)
   while ( $alt <= $nalt ) loop alcs
     interaction ( &atom_select and ( attr store6 = $alt or attr store6 = 0 ))
                 ( &atom_select and ( attr store6 = $alt ))
     evaluate ($alt=$alt+1)
   end loop alcs
 end

 {- check isolated atoms and atoms at special positions and add to
    list of fixed atoms if needed - store6 will be used -}
 
 do (store6=0) (all)

 connectivity
    selection=( &atom_select and not &atom_fixed ) 
    nsetto=store6
 end
 
 display
 display  list of isolated (non-covalently bonded) atoms:
 show element ( name ) ( attribute store6 = 1 )
 if ($select=0) then
       display  --none--
 end if

 display
 display  list of isolated (non-covalently bonded) di-atomic molecules:
 show element ( name ) ( attribute store6 = 2 )
 if ($select=0) then
       display  --none--
 end if
 
 if ( &md_type = "torsion" ) then
 
 {- for torsion angle dynamics we have to fix 
    isolated atoms and explicitly fixed atoms -}
    
    ident (store6) ((attribute store6 = 1) or
                    (attribute store6 = 2) or
                    ( not ( &atom_select ))  or
                            &atom_fixed ) 
                     
    display  $select isolated atoms, atoms in di-atomic molecules, 
    display  explicitly fixed atoms, and atoms not selected will be fixed. 

 else

    ident (store6) ((not ( &atom_select )) or
                           &atom_fixed )

 end if
    
 fix selection=( store6 ) end

 fastnb grid end

 show sum(1) (&atom_harm)
 if ( $result > 0 ) then
   evaluate ($harmonic=true)
 else
   evaluate ($harmonic=false)
 end if

 if ( &md_type = "torsion" ) then
   evaluate ($start_temp=&temperature)
   evaluate ($md_steps=6)
   evaluate ($fbeta=200)
 end if

 if ( &md_type = "cartesian" ) then
   evaluate ($start_temp=&temperature)
   evaluate ($md_steps=50)
   evaluate ($fbeta=100)
 end if

 if ( &md_scheme = "constant" ) then
   evaluate ($md_steps=&constant_steps)
 end if

 set seed=&seed end

 do (store7=x) (all)
 do (store8=y) (all)
 do (store9=z) (all)

 evaluate ($trial=1)

 while ( $trial <= &num_trials ) loop main

   do (x=store7) (all)
   do (y=store8) (all)
   do (z=store9) (all)

   if ( $harmonic = true ) then
     do (refx=x) (all)
     do (refy=y) (all)
     do (refz=z) (all)
     do (harm=0) (all)
     do (harm=&k_harmonic) (&atom_harm)
     flags include harm end
   end if
  
   parameter
     nbonds
        repel ?  evaluate ($repel_old=$result)
        rcon ?   evaluate ($rcon_old=$result)
        if ( $repel_old > 0 ) then
          if ( $repel_old = 1 ) then
             repel=1.  rcon=100.
          else
             repel=.75 rcon=50.
          end if
        end if
     end
   end

   do (fbeta=$fbeta) ( ( &atom_select ) and not store6 )

   do (vx=maxwell($start_temp)) ( ( &atom_select ) and not store6 )
   do (vy=maxwell($start_temp)) ( ( &atom_select ) and not store6 )
   do (vz=maxwell($start_temp)) ( ( &atom_select ) and not store6 )

   if ( &md_type = "torsion" ) then
     dynamics torsion
       topology
       
         maxlength=&torsion_maxlength
         maxchain=&torsion_maxchain
         maxtree=&torsion_maxtree
         maxbond=&torsion_maxbond
         kdihmax = 95.
         
         @CNS_TOPPAR:torsionmdmods
         
         fix group ( &atom_rigid )
         
       end
       nstep=0
       cmremove=true
     end
   end if  

   if ( &md_scheme = "slowcool" ) then
  
     evaluate ( $curr_temp = &temperature )
     while ( $curr_temp > 0.0 ) loop cool
       if ( &md_type = "torsion" ) then
         dynamics torsion
           timestep=&time_step
           nstep=$md_steps
           nprint=5
           cmremove=false
           if ( &tcontrol = "scaling" ) then
             vscaling=true  
           elseif ( &tcontrol = "coupling" ) then
             tcoupling=true
           end if
           temperature=$curr_temp
         end
       end if
       if ( &md_type = "cartesian" ) then
         dynamics cartesian
           if ($curr_temp=&temperature) then
              cmremove=true
           else
              cmremove=false
           end if
           timestep=&time_step
           nstep=$md_steps
           nprint=10
           if ( &tcontrol = "scaling" ) then
             vscaling=true  
           elseif ( &tcontrol = "coupling" ) then
             tcoupling=true
           end if
           temperature=$curr_temp
         end
       end if
       evaluate ( $curr_temp = $curr_temp - &cool_rate )
     end loop cool

   elseif ( &md_scheme = "constant" ) then
 
     evaluate ($traj_outfile=&output_root + "_" + encode($trial) + ".crd")

     if ( &md_type = "torsion" ) then
       dynamics torsion
         timestep=&time_step
         nstep=$md_steps
         nprint=5
         cmremove=false
         if ( &tcontrol = "scaling" ) then
           vscaling=true  
         elseif ( &tcontrol = "coupling" ) then
           tcoupling=true
         end if
         temperature=&temperature
         if ( &traj_freq > 0 ) then
           trajectory=$traj_outfile
           nsavc=&traj_freq
         end if
       end
     end if
     if ( &md_type = "cartesian" ) then
       dynamics cartesian
         timestep=&time_step
         nstep=$md_steps
         nprint=10
         if ( &tcontrol = "scaling" ) then
           vscaling=true  
         elseif ( &tcontrol = "coupling" ) then
           tcoupling=true
         end if
         temperature=&temperature
         if ( &traj_freq > 0 ) then
           trajectory=$traj_outfile
           nsavc=&traj_freq
         end if
       end
     end if
   end if
  
   parameter
     nbonds
       repel=$repel_old rcon=$rcon_old
     end
   end
  
   if ( &md_type = "torsion" ) then
     dynamics torsion
       nstep = 0
       cmremove=false
       topology
         reset
       end
     end
   end if

   if ( &md_scheme = "slowcool" ) then
     evaluate ($md_temp=(&temperature-0)/&cool_rate)
   else
     evaluate ($md_temp=1)
   end if

   print threshold=20.0 bond
   evaluate ($rmsd_bond=$result)

   print threshold=50.0 angle
   evaluate ($rmsd_angle=$result)

   evaluate ($coordinate_outfile=&output_root + "_" + encode($trial) + ".pdb")
   set display=$coordinate_outfile end

   display REMARK coordinates from molecular dynamics
   display REMARK rmsd bonds= $rmsd_bond[f8.6]  rmsd angles= $rmsd_angle[f8.5]
   display REMARK starting temperature= &temperature  temperature control= &STRIP%tcontrol 
   display REMARK nonbonded cutoff= &nonb_cutoff Angstroms  dieletric= &dielectric
   display REMARK total md steps= $md_temp * $md_steps  time step= &time_step ps
   if ( &use_cryst = true ) then
     display REMARK sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
   end if
   if ( &BLANK%parameter_infile_1 = false ) then
      display REMARK parameter file 1  : &STRIP%parameter_infile_1
   end if
   if ( &BLANK%parameter_infile_2 = false ) then
      display REMARK parameter file 2  : &STRIP%parameter_infile_2
   end if
   if ( &BLANK%parameter_infile_3 = false ) then
      display REMARK parameter file 3  : &STRIP%parameter_infile_3
   end if
   if ( &BLANK%parameter_infile_4 = false ) then
      display REMARK parameter file 4  : &STRIP%parameter_infile_4
   end if
   if ( &BLANK%parameter_infile_5 = false ) then
      display REMARK parameter file 5  : &STRIP%parameter_infile_5
   end if
   display REMARK molecular structure file: &STRIP%structure_infile
   display REMARK input coordinates: &STRIP%coordinate_infile
   if ( &BLANK%restraints_infile = false ) then
      display REMARK additional restraints file: &STRIP%restraints_infile
   end if
   if ( &BLANK%ncs_infile = false ) then
     display REMARK ncs= &STRIP%ncs_type  ncs file= &STRIP%ncs_infile
   else
     display REMARK ncs= none
   end if

   remark

   write coordinates output=$coordinate_outfile end

   set display=OUTPUT end

   evaluate ($trial=$trial+1)

 end loop main
 
 stop