{+ file: model_rigid.inp +}
{+ directory: general +}
{+ description: rigid body minimization with no
                experimental energy terms +}
{+ authors: Axel T. Brunger and Paul D. Adams +}
{+ copyright: Yale University +}

{- 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 store8 are available for general use -}

{- begin block parameter definition -} define(

{============================ coordinates ============================}

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

{==================== molecular information ==========================}

{* topology files *}
{===>} topology_infile_1="CNS_TOPPAR:protein.top";
{===>} topology_infile_2="CNS_TOPPAR:dna-rna.top";
{===>} topology_infile_3="CNS_TOPPAR:water.top";
{===>} topology_infile_4="CNS_TOPPAR:ion.top";
{===>} topology_infile_5="CNS_TOPPAR:carbohydrate.top";
{===>} topology_infile_6="";
{===>} topology_infile_7="";
{===>} topology_infile_8="";

{* linkage files for linear, continuous polymers (protein, DNA, RNA) *}
{===>} link_infile_1="CNS_TOPPAR:protein.link";
{===>} link_infile_2="CNS_TOPPAR:dna-rna-pho.link";
{===>} link_infile_3="";

{* parameter files *}
{===>} parameter_infile_1="CNS_TOPPAR:protein_rep.param";
{===>} parameter_infile_2="CNS_TOPPAR:dna-rna_rep.param";
{===>} parameter_infile_3="CNS_TOPPAR:water_rep.param";
{===>} parameter_infile_4="CNS_TOPPAR:ion.param";
{===>} parameter_infile_5="CNS_TOPPAR:carbohydrate.param";
{===>} parameter_infile_6="";
{===>} parameter_infile_7="";
{===>} parameter_infile_8="";

{* molecular topology file: optional (leave blank for auto generation) *}
{* 
   Auto generation of the molecular topology from the coordinates should only 
   be used if:
   (1) Each distinct protein, DNA, or RNA chain must have a separate segid 
       (or chainid if the chainid is non-blank). 
   (2) Each contiguous protein, RNA, or RNA chain must not be disrupted by 
       other types of residues or ligands.  Rather, these other residues 
       should be listed after protein, RNA/DNA chains. 
   (3) Disulphides are automatically detected based on distances between the sulfur atoms
      (must be less than 3 A apart).
   (4) Broken protein/RNA/DNA chains without terminii must be more than 2.5 A apart to be recognized as such.
   (5) N-linked glycan links are automatically recognized if the bonded atoms are less than 2.5 A apart.
   (6) Automatic generation cannot be used with alternate conformations. 
   For ligands, the user must make suitable topology and parameter files.
   For non-standard covalent linkages, the custom patch file should be used.
   Alternatively, the generate.inp or generate_easy.inp task files
   can be used to generated the mtf prior to running this task file.
    *}
{===>} structure_infile="amy_hydrogen.mtf";

{* for auto generation: extra linkages and modifications by custom patches *}
{===>} patch_infile="";

{===================== 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 in refinement *}
{* it is essential to include hydrogen atoms if anything other than
   regularization of geometry is being performed *}
{* this should include all conformations if multiple conformations are used *}
{===>} atom_select=(known);

{* select fixed atoms *}
{===>} atom_fixed=(none);

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

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

{* select atoms in group 1, use all known atoms to minimize
   everything as a rigid group *}
{===>} r_group_1=(all);

{* select atoms in group 2, use the selection none if not required *}
{===>} r_group_2=(none);

{* select atoms in group 3, use the selection none if not required *}
{===>} r_group_3=(none);

{* select atoms in group 4, use the selection none if not required *}
{===>} r_group_4=(none);

{* select atoms in group 5, use the selection none if not required *}
{===>} r_group_5=(none);

{* select atoms in group 6, use the selection none if not required *}
{===>} r_group_6=(none);

{* select atoms in group 7, use the selection none if not required *}
{===>} r_group_7=(none);

{* select atoms in group 8, use the selection none if not required *}
{===>} r_group_8=(none);

{* select atoms in group 9, use the selection none if not required *}
{===>} r_group_9=(none);

{* select atoms in group 10, use the selection none if not required *}
{===>} r_group_10=(none);

{* select atoms in group 11, use the selection none if not required *}
{===>} r_group_11=(none);

{* select atoms in group 12, use the selection none if not required *}
{===>} r_group_12=(none);

{* select atoms in group 13, use the selection none if not required *}
{===>} r_group_13=(none);

{* select atoms in group 14, use the selection none if not required *}
{===>} r_group_14=(none);

{* select atoms in group 15, use the selection none if not required *}
{===>} r_group_15=(none);

{* select atoms in group 16, use the selection none if not required *}
{===>} r_group_16=(none);

{* select atoms in group 17, use the selection none if not required *}
{===>} r_group_17=(none);

{* select atoms in group 18, use the selection none if not required *}
{===>} r_group_18=(none);

{* select atoms in group 19, use the selection none if not required *}
{===>} r_group_19=(none);

{* select atoms in group 20, use the selection none if not required *}
{===>} r_group_20=(none);

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

{=============== rigid body minimization parameters ==================}

{* number of minimization steps *}
{===>} rigid_nstep=30;

{* number of cycles *}
{===>} num_cycles=1;

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

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

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

{* output coordinate file *}
{===>} coordinate_outfile="model_rigid.pdb";

{===========================================================================}
{        things below this line do not normally need to be changed          }
{===========================================================================}

 ) {- end block parameter definition -}

 checkversion 1.3

 evaluate ($log_level=quiet)

  if ( $log_level = verbose ) then
   set message=normal echo=on end
 else
   set message=off echo=off end
 end if

 if ( &BLANK%structure_infile = true ) then
 
    {- read topology files -}
    topology
     evaluate ($counter=1)
     evaluate ($done=false)
     while ( $done = false ) loop read
      if ( &exist_topology_infile_$counter = true ) then
         if ( &BLANK%topology_infile_$counter = false ) then
            @@&topology_infile_$counter
         end if
      else
        evaluate ($done=true)
      end if
      evaluate ($counter=$counter+1)
     end loop read
    end
    
    @CNS_XTALMODULE:mtfautogenerate (
                                  coordinate_infile=&coordinate_infile;
                                  convert=true;
                                  separate=true;
                                  atom_delete=(not known);
                                  hydrogen_flag=true;
                                  break_cutoff=2.5;
                                  disulphide_dist=3.0;
                                  carbo_dist=2.5;
                                  patch_infile=&patch_infile;
                                  O5_becomes="O";
                                 )

 else

   structure @&structure_infile end
   coordinates @&coordinate_infile

 end if

 {- read parameter files -}
 parameter
  evaluate ($counter=1)
  evaluate ($done=false)
  while ( $done = false ) loop read
   if ( &exist_parameter_infile_$counter = true ) then
      if ( &BLANK%parameter_infile_$counter = false ) then
         @@&parameter_infile_$counter
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
  end loop read
 end

 set message=normal echo=on end

 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

 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

 igroup
   interaction ( &atom_select )
               ( &atom_select )
 end

 fix selection=( &atom_fixed ) end

 fastnb grid end

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

 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

 evaluate ($ngroup=1)
 evaluate ($group=1)
 evaluate ($done=false)
 while ( $done = false ) loop group
   if ( &exist_r_group_$group = true ) then
     show sum(1) ( &r_group_$group )
     evaluate ($size_$group=$result)
     evaluate ($ngroup=$ngroup+1)
   else
     evaluate ($done=true)
     evaluate ($ngroup=$ngroup-1)
   end if
   evaluate ($group=$group+1)
 end loop group

 evaluate ($cycle=1)

 while ($cycle <= &num_cycles) loop main

   if ( &rigid_nstep > 0 ) then

     minimize rigid

       evaluate ($counter=1)
       while ( $counter <= $ngroup ) loop sele
         if ( $size_$counter > 0 ) then
           group=(&atom_select and (&r_group_$counter))
         end if
         evaluate ($counter=$counter+1)
       end loop sele

       nstep=&rigid_nstep
       drop=10.0
     end
   end if

   evaluate ($cycle=$cycle+1)

 end loop main

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

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

 set display=&coordinate_outfile end

 display REMARK coordinates from rigid body minimization
 display REMARK rmsd bonds= $rmsd_bond[f8.6]  rmsd angles= $rmsd_angle[f8.5]
 display REMARK nonbonded cutoff= &nonb_cutoff Angstroms  dieletric= &dielectric
 display REMARK cycles= &num_cycles steps= &rigid_nstep
 if ( &use_cryst = true ) then
   display REMARK sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
 end if

 evaluate ($counter=1)
 evaluate ($done=false)
 while ( $done = false ) loop read
  if ( &exist_parameter_infile_$counter = true ) then
    if ( &BLANK%parameter_infile_$counter = false ) then
      display REMARK parameter file $counter  : &STRIP%parameter_infile_$counter
    end if
  else
   evaluate ($done=true)
  end if
  evaluate ($counter=$counter+1)
 end loop read
 
 if ( &BLANK%structure_infile = true ) then
   display REMARK molecular structure file: automatic

   evaluate ($counter=1)
   evaluate ($done=false)
   while ( $done = false ) loop read
    if ( &exist_topology_infile_$counter = true ) then
      if ( &BLANK%topology_infile_$counter = false ) then
        display REMARK topology file $counter  : &STRIP%topology_infile_$counter
      end if
    else
     evaluate ($done=true)
    end if
    evaluate ($counter=$counter+1)
   end loop read

   evaluate ($counter=1)
   evaluate ($done=false)
   while ( $done = false ) loop read
    if ( &exist_link_infile_$counter = true ) then
      if ( &BLANK%link_infile_$counter = false ) then
        display REMARK linkage file $counter  : &STRIP%link_infile_$counter
      end if
    else
     evaluate ($done=true)
    end if
    evaluate ($counter=$counter+1)
   end loop read

   if ( &BLANK%patch_infile = false ) then
      display REMARK custom patch file = &STRIP%patch_infile
   end if

 else
   display REMARK molecular structure file: &STRIP%structure_infile
 end if

 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 format=PDBO end

 stop