{+ file: anneal.inp +}
{+ directory: xtal_refine +}
{+ description: Crystallographic simulated annealing refinement +}
{+ authors: Axel T. Brunger, Luke M. Rice and Paul D. Adams +}
{+ copyright: Yale University +}

{+ reference: A.T. Brunger, J. Kuriyan and M. Karplus, Crystallographic
              R factor Refinement by Molecular Dynamics, Science
              235, 458-460 (1987) +}
{+ reference: A.T. Brunger, A. Krukowski and J. Erickson, Slow-Cooling 
              Protocols for Crystallographic Refinement by Simulated 
              Annealing, Acta Cryst. A46, 585-593 (1990) +}
{+ reference: A.T. Brunger, The Free R Value: a Novel Statistical 
              Quantity for Assessing the Accuracy of Crystal Structures, 
              Nature 355, 472-474 (1992) +}
{+ 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) +}
{+ reference: N.S. Pannu and R.J. Read, Improved structure refinement 
              through maximum likelihood, Acta Cryst. A52, 659-668 (1996) +}
{+ reference: P.D. Adams, N.S. Pannu, R.J. Read and A.T. Brunger, 
              Cross-validated Maximum Likelihood Enhances Crystallographic 
              Simulated Annealing Refinement, Proc. Natl. Acad. Sci. USA 
              94, 5018-5023 (1997) +}

{- 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 =========================}

{* molecular topology file *}
{===>} structure_infile="amy.mtf";

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

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

{====================== crystallographic data ========================}

{* 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;

{* anomalous f' f'' library file *}
{* If a file is not specified, no anomalous contribution will be included *}
{+ choice: "CNS_XRAYLIB:anom_cu.lib" "CNS_XRAYLIB:anom_mo.lib" "" user_file +}
{===>} anom_library="";

{* reflection files *}
{* specify non-anomalous reflection files before anomalous reflection files. *}
{* files must contain unique array names otherwise errors will occur *}
{===>} reflection_infile_1="amy.cv";
{===>} reflection_infile_2="";
{===>} reflection_infile_3="";

{* reciprocal space array containing observed amplitudes: required *}
{===>} obs_f="fobs";

{* reciprocal space array containing sigma values for amplitudes: required *}
{===>} obs_sigf="sigma";

{* reciprocal space array containing test set for cross-validation: required *}
{* cross-validation should always be used, with the possible exception
   of a final round of refinement including all data *}
{* cross-validation is always required for the maximum likelihood targets *}
{===>} test_set="test";

{* number for selection of test reflections: required for cross-validation *}
{* ie. reflections with the test set array equal to this number will be
       used for cross-validation, all other reflections form the working set *}
{===>} test_flag=1;

{* reciprocal space array containing weighting scheme for observed
   amplitudes: optional *}
{* only used for the "residual" and "vector" targets - this will
   default to a constant value of 1 if array is not present *}
{===>} obs_w="";

{* reciprocal space array containing observed intensities: optional *}
{* required for the "mli" target *}
{===>} obs_i="";

{* reciprocal space array containing sigma values for intensities: optional *}
{* required for the "mli" target *}
{===>} obs_sigi="";

{* reciprocal space arrays with experimental phase probability
   distribution: optional *}
{* Hendrickson-Lattman coefficients A,B,C,D *}
{* required for the "mlhl" target *}
{+ table: rows=1 "HL coefficients" cols=4 "A" "B" "C" "D" +}
{===>} obs_pa="";
{===>} obs_pb="";
{===>} obs_pc="";
{===>} obs_pd="";

{* complex reciprocal space array containing experimental phases: optional *}
{* required for the "mixed" and "vector" targets *}
{===>} obs_phase="";

{* reciprocal space array containing experimental figures of merit: optional *}
{* required for the "mixed" target *}
{===>} obs_fom="";

{* resolution limits to be used in refinement *}
{* the full resolution range of observed data should be used in refinement.
   A bulk solvent correction should be applied to allow the use of low
   resolution terms. If no bulk solvent correction is applied, data must
   be truncated at a lower resolution limit of between 8 and 6 Angstrom. *}
{+ table: rows=1 "resolution" cols=2 "lowest" "highest" +}
{===>} low_res=500.0;
{===>} high_res=2.0;

{* apply rejection criteria to amplitudes or intensities *}
{+ choice: "amplitude" "intensity" +}
{===>} obs_type="amplitude";

{* Observed data cutoff criteria: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs)/sigma < cutoff are rejected. *}
{===>} sigma_cut=0.0;

{* rms outlier cutoff: applied to amplitudes or intensities *}
{* reflections with magnitude(Obs) > cutoff*rms(Obs) will be rejected *}
{===>} obs_rms=10000;

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

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

{============ initial B-factor and bulk solvent corrections ==========}

{* overall B-factor correction *}
{+ choice: "no" "isotropic" "anisotropic" +}
{===>} bscale="anisotropic";

{* bulk solvent correction *}
{* a mask is required around the molecule(s). The region
   outside this mask is the solvent region *}
{+ choice: true false +}
{===>} bulk_sol=true;

{* bulk solvent mask file *}
{* mask will be read from O type mask file if a name is given
   otherwise calculated from coordinates of selected atoms *}
{===>} bulk_mask_infile="";

{* automatic bulk solvent parameter search *}
{+ choice: true false +}
{===>} sol_auto=true;

{* optional file with a listing of the results of the automatic bulk solvent grid 
   search *}
{===>} sol_output="";

{* fixed solvent mask parameters if the automatic option is not used *}
{+ table: rows=1 "bulk solvent" cols=2 "probe radius (A)" "shrink radius (A)" +}
{===>} sol_rad=1.0;
{===>} sol_shrink=1.0;

{* fixed solvent parameters if the automatic option is not used *}
{+ table: rows=1 "bulk solvent" cols=2 "e-density level (e/A^3)" "B-factor (A^2)" +}
{===>} sol_k=-1;
{===>} sol_b=-1;

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

{* select atoms to be included in refinement *}
{* this should include all conformations if multiple conformations are used *}
{===>} atom_select=(known and not hydrogen);

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

{* select atoms to be harmonically restrained during refinement *}
{===>} 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="";

{====================== refinement parameters ========================}

{* initial steps of minimization with x-ray term on *}
{===>} init_min=200;

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

{* number of minimization steps to regularize geometry before torsion md *}
{===>} geometry_min=100;

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

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

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

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

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

{* final steps of low temperature (300K) molecular dynamics *}
{===>} final_dyn=100;

{* final steps of minimization with x-ray term on *}
{===>} final_min=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;

{* refinement target *}
{+ list: mlf: maximum likelihood target using amplitudes
         mli: maximum likelihood target using intensities
        mlhl: maximum likelihood target using amplitudes
              and phase probability distribution
    residual: standard crystallographic residual
      vector: vector residual
       mixed: (1-fom)*residual + fom*vector
        e2e2: correlation coefficient using normalized E^2
        e1e1: correlation coefficient using normalized E
        f2f2: correlation coefficient using F^2
        f1f1: correlation coefficient using F +}
{+ choice: "mlf" "mli" "mlhl" "residual" "vector" "mixed" 
           "e2e2" "e1e1" "f2f2" "f1f1" +}
{===>} reftarget="mlf";

{* Wa weight for X-ray term *}
{* this will be determined automatically if a negative value is given.
   Note: wa can be very different depending on the target - if it is not
         determined automatically make sure an appropriate value is used *}
{===>} wa=-1;

{* number of bins for refinement target *}
{* this will be determined automatically if a negative value is given
   otherwise the specified number of bins will be used *}
{===>} target_bins=-1;

{* memory allocation for FFT calculation *}
{* this will be determined automatically if a negative value is given
   otherwise the specified number of words will be allocated *}
{===>} fft_memory=-1;

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

{* root name for output coordinate files *}
{* files  <output_root>_1.pdb ... <output_root>_n.pdb are written *}
{===>} output_root="anneal";

{* format output coordinates for use in o *}
{* if false then the default CNS output coordinate format will be used *}
{+ choice: true false +}
{===>} pdb_o_format=true;

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

 ) {- end block parameter definition -}

 checkversion 1.2

 evaluate ($log_level=quiet)

 structure @&structure_infile end

 coordinates @&coordinate_infile

 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

 xray

   @CNS_XTALLIB:spacegroup.lib (sg=&sg;
                                sgparam=$sgparam;)

   a=&a b=&b c=&c  alpha=&alpha beta=&beta gamma=&gamma

   @CNS_XRAYLIB:scatter.lib

   if ( &BLANK%reflection_infile_1 = false ) then
     reflection @@&reflection_infile_1 end
   end if
   if ( &BLANK%reflection_infile_2 = false ) then
     reflection @@&reflection_infile_2 end
   end if
   if ( &BLANK%reflection_infile_3 = false ) then
     reflection @@&reflection_infile_3 end
   end if

 end

 if ( &BLANK%anom_library = false ) then
   @@&anom_library
 else
   set echo=off end
   xray anomalous=? end
   if ( $result = true ) then
     display Warning: no anomalous library has been specified
     display          no anomalous contribution will used in refinement
   end if
   set echo=on end
 end if

 {- copy define parameters of optional arrays into symbols so 
    we can redefine them -}
    
 evaluate ($obs_i=&obs_i)
 evaluate ($obs_sigi=&obs_sigi)
 evaluate ($obs_w=&obs_w)
 xray
   @@CNS_XTALMODULE:checkrefinput (
                                  reftarget=&reftarget;
                                  obs_f=&obs_f;
                                  obs_sigf=&obs_sigf;
                                  test_set=&test_set;
                                  obs_pa=&obs_pa;
                                  obs_pb=&obs_pb;
                                  obs_pc=&obs_pc;
                                  obs_pd=&obs_pd;
                                  obs_phase=&obs_phase;
                                  obs_fom=&obs_fom;
                                  obs_w=$obs_w;
                                  obs_i=$obs_i;
                                  obs_sigi=$obs_sigi;
                                  )
 
   query name=fcalc domain=reciprocal end
   if ( $object_exist = false ) then
      declare name=fcalc domain=reciprocal type=complex end
   end if
   declare name=fbulk domain=reciprocal type=complex end
   do (fbulk=0) ( all )

   binresolution &low_res &high_res
   mapresolution &high_res

   if ( &obs_type = "intensity" ) then
     if ( &BLANK%obs_i = true ) then
       display  Error: observed intensity array is undefined
       display         aborting script
       abort
     end if
     evaluate ($reject_obs=&obs_i)
     evaluate ($reject_sig=&obs_sigi)
     show min (amplitude(&STRIP%obs_i)) (all)
     evaluate ($obs_lower_limit=$result-0.1)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
     evaluate ($obs_lower_limit=0)
   end if

   declare name=ref_active domain=reciprocal type=integer end
   declare name=tst_active domain=reciprocal type=integer end

   do (ref_active=0) ( all )
   do (ref_active=1) ( ( amplitude($STRIP%reject_obs) > $obs_lower_limit ) and
                       ( &low_res >= d >= &high_res ) )

   statistics overall
     completeness
     selection=( ref_active=1 )
   end
   evaluate ($total_compl=$expression1)

   show sum(1) ( ref_active=1 )
   evaluate ($total_read=$select)
   evaluate ($total_theor=int(1./$total_compl * $total_read))

   show rms (amplitude($STRIP%reject_obs)) ( ref_active=1 )
   evaluate ($obs_high=$result*&obs_rms)
   show min (amplitude($STRIP%reject_obs)) ( ref_active=1 )
   evaluate ($obs_low=$result)

   do (ref_active=0) ( all )
   do (ref_active=1)
                  ( ( amplitude($STRIP%reject_obs) >= &sigma_cut*$STRIP%reject_sig ) and
                    ( $STRIP%reject_sig # 0 ) and
                    ( $obs_low <= amplitude($STRIP%reject_obs) <= $obs_high ) and
                    ( &low_res >= d >= &high_res ) )

   do (tst_active=0) (all)
   if ( &BLANK%test_set = false ) then
     do (tst_active=1) (ref_active=1 and &STRIP%test_set=&test_flag)
   end if

   show sum(1) ( ref_active=1 and tst_active=0 )
   evaluate ($total_work=$select)
   show sum(1) ( ref_active=1 and tst_active=1 )
   evaluate ($total_test=$select)
   evaluate ($total_used=$total_work+$total_test)

   evaluate ($unobserved=$total_theor-$total_read)
   evaluate ($rejected=$total_read-$total_used)
   evaluate ($per_unobs=100*($unobserved/$total_theor))
   evaluate ($per_reject=100*($rejected/$total_theor))
   evaluate ($per_used=100*($total_used/$total_theor))
   evaluate ($per_work=100*($total_work/$total_theor))
   evaluate ($per_test=100*($total_test/$total_theor))

   associate fcalc ( &atom_select )

   tselection=( ref_active=1 )

   cvselection=( tst_active=1 )

   method=FFT          
   
 {- MODIFIED 2/15/06 -}
 end
 

 show min ( b ) ( &atom_select )
 evaluate ($b_min=$result)
 @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=auto;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
                            
 xray
 {- END MODIFICATION -}
      
   tolerance=0.0 lookup=false

   if ( &wa >= 0 ) then
      wa=&wa
   end if

 end                  

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

 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 -}
 
 @CNS_XTALMODULE:setupfixed (
                           mode="minimization";
                           atom_select=&atom_select;
                           atom_fixed=&atom_fixed;
                           atom_total_fixed=store6;
                           atom_multiplicity=rmsd;
                           mset=$mset;
                           )
                           
 fix selection=( store6 ) end

 fastnb grid end

 flags                                       
    exclude elec include pvdw xref                   
    ?                                        
 end                                         

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

 xray
    predict
      mode=reciprocal
      to=fcalc
      selection=(ref_active=1)
      atomselection=( &atom_select ) 
    end
 end

 {- BEGIN MODIFICATION -}
 @CNS_XTALMODULE:scale_and_solvent_grid_search (
                             bscale=&bscale;
                             sel=( ref_active=1 );
                             sel_test=( tst_active=1 );
                             atom_select=( &atom_select );
                             bulk_sol=&bulk_sol;
                             bulk_mask=&bulk_mask_infile;
                             bulk_atoms=( &atom_select );
                             
                             sol_auto=&sol_auto;
                             sol_k=&sol_k;
                             sol_b=&sol_b;
                             sol_rad=&sol_rad;
                             sol_shrink=&sol_shrink;
  
                             fcalc=fcalc;
                             obs_f=&STRIP%obs_f;
                             obs_sigf=&STRIP%obs_sigf;
                             obs_i=$STRIP%obs_i;
                             obs_sigi=$STRIP%obs_sigi;                             
                             fpart=fbulk;
                             
!
! Begin modification (6/28/06)                             
                             Baniso_11=$Baniso_11;
                             Baniso_22=$Baniso_22;
                             Baniso_33=$Baniso_33;
                             Baniso_12=$Baniso_12;
                             Baniso_13=$Baniso_13;
                             Baniso_23=$Baniso_23;
                             Biso=$Biso_model;
! End modification
! 
                             
                             sol_k_best=$sol_k_ref;
                             sol_b_best=$sol_b_ref;
			     solrad_best=$solrad_best;
			     shrink_best=$shrink_best;
                             
                             b=b;

                             low_b_flag=$low_b_flag;
                            
                             sol_output=&sol_output;
                             
                             )
 xray
   @@CNS_XTALMODULE:calculate_r (
                                 fobs=&STRIP%obs_f;
                                 fcalc=fcalc;
                                 fpart=fbulk;
                                 sel=( ref_active=1 );
                                 sel_test=( tst_active=1 );
                                 print=true;
                                 output=OUTPUT;
                                 r=$start_r;
                                 test_r=$start_test_r;)
 end

 {- check the gridding again since the minimum B-factor may have changed -}
 show min ( b ) ( &atom_select )
 evaluate ($b_min=$result)
 @@CNS_XTALMODULE:fft_parameter_check ( 
                             d_min=&high_res; 
                             b_min=$b_min;
                             grid=auto;
                             fft_memory=&fft_memory;
                             fft_grid=$fft_grid;   
                             fft_b_add=$fft_b_add; 
                             fft_elim=$fft_elim; 
                                      )
{- END MODIFICATION -}

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

 if ( &md_type = "cartesian" ) then
   evaluate ($start_temp=&temperature)
   evaluate ($time_step=0.0005)
   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

 xray
   @@CNS_XTALMODULE:refinementtarget (target=&reftarget;
                                     sig_sigacv=0.07;
                                     mbins=&target_bins;
                                     fobs=&STRIP%obs_f;
                                     sigma=&STRIP%obs_sigf;
                                     weight=$STRIP%obs_w;
                                     iobs=$STRIP%obs_i;
                                     sigi=$STRIP%obs_sigi;
                                     test=tst_active;
                                     fcalc=fcalc;
                                     fpart=fbulk;
                                     pa=&STRIP%obs_pa;
                                     pb=&STRIP%obs_pb;
                                     pc=&STRIP%obs_pc;
                                     pd=&STRIP%obs_pd;
                                     phase=&STRIP%obs_phase;
                                     fom=&STRIP%obs_fom;
                                     sel=(ref_active=1);
                                     sel_test=(tst_active=1);
                                     statistics=true;)
 end

 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
  
 if ( &wa < 0 ) then
   @@CNS_XTALMODULE:getweight (
                              selected=&atom_select;
                              fixed=(store6);
                              wa=$wa_init;
                             )
 end if

 if ( &init_min > 0 ) then
   minimize powell                             
     nstep=&init_min                                             
     nprint=5
     drop=10.0
   end                                            
 end if

 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)

   {- check isolated atoms and atoms at special positions and add to
      list of fixed atoms if needed - store6 will be used -}
 
   @CNS_XTALMODULE:setupfixed (
                             mode=&md_type;
                             atom_select=&atom_select;
                             atom_fixed=&atom_fixed;
                             atom_total_fixed=store6;
                             atom_multiplicity=rmsd;
                             mset=$mset;
                              )
                           
   fix selection=( store6 ) end

   xray
      predict
        mode=reciprocal
        to=fcalc
        selection=(ref_active=1)
        atomselection=( &atom_select ) 
      end
   end

   xray
     @@CNS_XTALMODULE:refinementtarget (target=&reftarget;
                                       sig_sigacv=0.07;
                                       mbins=&target_bins;
                                       fobs=&STRIP%obs_f;
                                       sigma=&STRIP%obs_sigf;
                                       weight=$STRIP%obs_w;
                                       iobs=$STRIP%obs_i;
                                       sigi=$STRIP%obs_sigi;
                                       test=tst_active;
                                       fcalc=fcalc;
                                       fpart=fbulk;
                                       pa=&STRIP%obs_pa;
                                       pb=&STRIP%obs_pb;
                                       pc=&STRIP%obs_pc;
                                       pd=&STRIP%obs_pd;
                                       phase=&STRIP%obs_phase;
                                       fom=&STRIP%obs_fom;
                                       sel=(ref_active=1);
                                       sel_test=(tst_active=1);
                                       statistics=true;)
   end

   if ( &wa < 0 ) then
     @@CNS_XTALMODULE:getweight (
                                selected=&atom_select;
                                fixed=(store6);
                                wa=$wa_dyna;
                               )
   end if

   if ( $md_steps > 0 ) then
     if ( &md_type = "torsion" ) then
       do (refx=x) (all)
       do (refy=y) (all)
       do (refz=z) (all)
       do (harm=10) (all)
       flags include harm exclude xref end
       if ( &geometry_min > 0 ) then
         minimize powell
           nstep=&geometry_min
           nprint=10
         end
       end if
       flags exclude harm include xref end
       do (harm=0) (all)
     end if
   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
  
   parameter
     nbonds
        repel ?  evaluate ($repel_old=$result)
        rcon ?   evaluate ($rcon_old=$result)
        if ( $repel_old = 1 ) then
           repel=1.  rcon=100.
        else
           repel=.75 rcon=50.
        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 )

   xray
     tolerance=0.2 lookup=true
   end

   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  

   do (store5=mass) ( all )
   do (mass=max(10,min(30,mass))) ( all )

   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
           vscaling=true  
           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
           vscaling=true
           temperature=$curr_temp
         end
       end if
       evaluate ( $curr_temp = $curr_temp - &cool_rate )
     end loop cool

   elseif ( &md_scheme = "constant" ) then
 
     if ( &md_type = "torsion" ) then
       dynamics torsion
         timestep=$time_step
         nstep=$md_steps
         nprint=5
         cmremove=false
         vscaling=true  
         temperature=&temperature
       end
     end if
     if ( &md_type = "cartesian" ) then
       dynamics cartesian
         timestep=$time_step
         nstep=$md_steps
         nprint=10
         vscaling=true
         temperature=&temperature
       end
     end if
   end if
  
   parameter
     nbonds
       repel=$repel_old rcon=$rcon_old
     end
   end
  
   if ( &md_type = "torsion" ) then
     dynamics torsion
       nstep=&final_dyn
       timestep=0.001
       nprint=10 
       cmremove=false
       vscaling=true 
       temperature=300
     end
   end if
   if ( &md_type = "cartesian" ) then
     dynamics cartesian
       nstep=&final_dyn
       timestep=0.001
       nprint=10
       vscaling=true
       temperature=300
     end
   end if

   if ( &md_type = "torsion" ) then
     dynamics torsion
       nstep = 0
       cmremove=false
       topology
         reset
       end
     end
   end if

   do (mass=store5) ( all )
   
   xray
     tolerance=0.0 lookup=false
   end
  
   {- check isolated atoms and atoms at special positions and add to
      list of fixed atoms if needed - store6 will be used -}
 
   @CNS_XTALMODULE:setupfixed (
                             mode="minimization";
                             atom_select=&atom_select;
                             atom_fixed=&atom_fixed;
                             atom_total_fixed=store6;
                             atom_multiplicity=rmsd;
                             mset=$mset;
                             )
                            
   fix selection=( store6 ) end

   if ( &final_min > 0 ) then
     minimize lbfgs
       nstep=&final_min
       nprint=10
     end
   end if
  
   xray
     @@CNS_XTALMODULE:refinementtarget (target=&reftarget;
                                       sig_sigacv=0.07;
                                       mbins=&target_bins;
                                       fobs=&STRIP%obs_f;
                                       sigma=&STRIP%obs_sigf;
                                       weight=$STRIP%obs_w;
                                       iobs=$STRIP%obs_i;
                                       sigi=$STRIP%obs_sigi;
                                       test=tst_active;
                                       fcalc=fcalc;
                                       fpart=fbulk;
                                       pa=&STRIP%obs_pa;
                                       pb=&STRIP%obs_pb;
                                       pc=&STRIP%obs_pc;
                                       pd=&STRIP%obs_pd;
                                       phase=&STRIP%obs_phase;
                                       fom=&STRIP%obs_fom;
                                       sel=(ref_active=1);
                                       sel_test=(tst_active=1);
                                       statistics=true;)
   end

   if ( &wa < 0 ) then
     @@CNS_XTALMODULE:getweight (
                                selected=&atom_select;
                                fixed=(store6);
                                wa=$wa_final;
                               )
   end if

   if ( &final_min > 0 ) then
     minimize lbfgs                             
       nstep=&final_min                                             
       nprint=10
     end                                            
   end if
  
   xray
      predict
        mode=reciprocal
        to=fcalc
        selection=(ref_active=1)
        atomselection=( &atom_select ) 
      end
      @@CNS_XTALMODULE:calculate_r (fobs=&STRIP%obs_f;
                                   fcalc=fcalc;
                                   fpart=fbulk;
                                   sel=(ref_active=1);
                                   sel_test=(tst_active=1);
                                   print=true;
                                   output=OUTPUT;
                                   r=$full_r;
                                   test_r=$full_test_r;)
   end
  
   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 simulated annealing refinement
   display REMARK refinement resolution: &low_res - &high_res A 
   if ( $total_test > 0 ) then
     display REMARK starting r= $start_r[f6.4] free_r= $start_test_r[f6.4]
     display REMARK final    r= $full_r[f6.4] free_r= $full_test_r[f6.4]
   else
     display REMARK starting r= $start_r[f6.4]
     display REMARK final    r= $full_r[f6.4]
   end if
   display REMARK rmsd bonds= $rmsd_bond[f8.6]  rmsd angles= $rmsd_angle[f8.5]
   xray wa=? end
   evaluate ($wa_print=$result)
   if ( &wa < 0 ) then
     display REMARK wa_initial= $wa_init  wa_dynamics= $wa_dyna  wa_final= $wa_print
   else
     display REMARK wa= $wa_print
   end if
   display REMARK target= &STRIP%reftarget  md-method= &STRIP%md_type  annealing schedule= &STRIP%md_scheme
   display REMARK starting temperature= &temperature  total md steps= $md_temp * $md_steps
   display REMARK sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
   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%anom_library = false ) then
     display REMARK anomalous f' f'' library: &STRIP%anom_library
   end if
   if ( &BLANK%reflection_infile_1 = false ) then
     display REMARK reflection file= &STRIP%reflection_infile_1
   end if
   if ( &BLANK%reflection_infile_2 = false ) then
     display REMARK reflection file= &STRIP%reflection_infile_2
   end if
   if ( &BLANK%reflection_infile_3 = false ) then
     display REMARK reflection file= &STRIP%reflection_infile_3
   end if
   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
   if ( &bscale # "no" ) then
     if ( $low_b_flag = true ) then
       display REMARK warning: B-correction gave atomic B-values less than zero
       display REMARK          they have been reset to zero
     end if
   end if
   ! 
   ! Begin modification (6/28/06)
   if ( &bscale = "anisotropic" ) then
    display REMARK Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
    display REMARK   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
    display REMARK   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
    display REMARK Isotropic component added to coordinate array B: $Biso_model[f8.3]
   elseif ( &bscale = "isotropic" ) then
    display REMARK B-factor applied to coordinate array B: $Biso_model[f8.3]
   else
    display REMARK initial B-factor correction: none
   end if
   ! End modification
   !
   
   {- MODIFIED 5/18/05 -}
   if ( &bulk_sol = true ) then 
     display REMARK bulk solvent: probe radius=$solrad_best, shrink value=$solrad_best
     display REMARK bulk solvent: density level= $sol_k_ref e/A^3, B-factor= $sol_b_ref A^2
   else
     display REMARK bulk solvent: false
   end if
   {- END MODIFICATION -}
 
   if ( &obs_type = "intensity" ) then
     display REMARK reflections with Iobs/sigma_I < &sigma_cut rejected
     display REMARK reflections with Iobs > &obs_rms * rms(Iobs) rejected
   else
     display REMARK reflections with |Fobs|/sigma_F < &sigma_cut rejected
     display REMARK reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
   end if
   xray anomalous=? end
   if ( $result = true ) then
     display REMARK anomalous diffraction data was input
   end if

   {- MODIFIED 2/15/06 -}
   display REMARK fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
   {- END MODIFICATION -}
 
   display REMARK theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
   display REMARK number of unobserved reflections (no entry or |F|=0): $unobserved[I6] ( $per_unobs[f5.1] % )
   display REMARK number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
   display REMARK total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
   display REMARK number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
   display REMARK number of reflections in test set:                    $total_test[I6] ( $per_test[f5.1] % )

   remark
   
   @CNS_XTALMODULE:write_pdb (pdb_o_format=&pdb_o_format;
                              coordinate_outfile=$coordinate_outfile;
                              sgparam=$sgparam;)

   set display=OUTPUT end

   evaluate ($trial=$trial+1)

 end loop main
 
 stop