{+ file: optimize_wa.inp +}
{+ directory: xtal_refine +}
{+ description: Optimize Wa (xray/geometry weight) +}
{+ authors: Axel T. Brunger 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 store4 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 *}
{===>} test_set="test";

{* number for selection of test reflections: required *}
{* 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="";

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

{* wa scale factors to test *}
{===>} wa_scale=0.1 0.2 0.5 0.75 1 1.5 2 5 10 20 50;

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

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

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

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

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

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

{* final steps of minimization with x-ray term on *}
{===>} final_min=50;

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

{* final steps of low temperature molecular dynamics *}
{===>} final_dyn=50;

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

{* output listing file *}
{===>} list_outfile="optimize_wa.list";

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

   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)
   do (tst_active=1) (ref_active=1 and &STRIP%test_set=&test_flag)

   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                  

 xray
   @@CNS_XTALMODULE:checkcv (active=ref_active;
                             test=tst_active;)
 end

 evaluate ($best_test_r=1)

 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=&md_type;
                           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 -}


 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

 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

 if ( &md_scheme = "slowcool" ) then
   evaluate ($md_temp=&temperature/&cool_rate)
 else
   evaluate ($md_temp=1)
 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 ( &wa < 0 ) then
   @@CNS_XTALMODULE:getweight (
                              selected=&atom_select;
                              fixed=(store6);
                              wa=$wa_init;
                             )
 else
   evaluate ($wa_init=&wa)
 end if

 set display=&list_outfile end

 display >>> refinement resolution: &low_res - &high_res A  
 display >>> starting r= $start_r[f6.4] free_r= $start_test_r[f6.4]
 display >>> target= &STRIP%reftarget  md-method= &STRIP%md_type  annealing schedule= &STRIP%md_scheme
 display >>> starting temperature= &temperature  total md steps= $md_steps * $md_temp
 display >>> sg= &STRIP%sg a= &a b= &b c= &c alpha= &alpha beta= &beta gamma= &gamma
 if ( &BLANK%parameter_infile_1 = false ) then
    display >>> parameter file 1  : &STRIP%parameter_infile_1
 end if
 if ( &BLANK%parameter_infile_2 = false ) then
    display >>> parameter file 2  : &STRIP%parameter_infile_2
 end if
 if ( &BLANK%parameter_infile_3 = false ) then
    display >>> parameter file 3  : &STRIP%parameter_infile_3
 end if
 if ( &BLANK%parameter_infile_4 = false ) then
    display >>> parameter file 4  : &STRIP%parameter_infile_4
 end if
 if ( &BLANK%parameter_infile_5 = false ) then
    display >>> parameter file 5  : &STRIP%parameter_infile_5
 end if
 display >>> molecular structure file: &STRIP%structure_infile
 display >>> input coordinates: &STRIP%coordinate_infile
 if ( &BLANK%anom_library = false ) then
   display >>> anomalous f' f'' library: &STRIP%anom_library
 end if
 if ( &BLANK%reflection_infile_1 = false ) then
   display >>> reflection file= &STRIP%reflection_infile_1
 end if
 if ( &BLANK%reflection_infile_2 = false ) then
   display >>> reflection file= &STRIP%reflection_infile_2
 end if
 if ( &BLANK%reflection_infile_3 = false ) then
   display >>> reflection file= &STRIP%reflection_infile_3
 end if
 if ( &BLANK%restraints_infile = false ) then
   display >>> additional restraints file: &STRIP%restraints_infile
 end if
 if ( &BLANK%ncs_infile = false ) then
   display >>> ncs= &STRIP%ncs_type  ncs file= &STRIP%ncs_infile
 else
   display >>> ncs= none
 end if
 if ( &bscale # "no" ) then
   if ( $low_b_flag = true ) then
     display >>> warning: B-correction gave atomic B-values less than zero
     display >>>          they have been reset to zero
   end if
 end if
 ! 
 ! Begin modification (6/28/06)
 if ( &bscale = "anisotropic" ) then
   display >>> Anisotropic B-factor tensor Ucart of atomic model without isotropic component :
   display >>>   B11=$Baniso_11[f8.3] B22=$Baniso_22[f8.3] B33=$Baniso_33[f8.3]
   display >>>   B12=$Baniso_12[f8.3] B13=$Baniso_13[f8.3] B23=$Baniso_23[f8.3]
   display >>> Isotropic component added to coordinate array B: $Biso_model[f8.3]
 elseif ( &bscale = "isotropic" ) then
   display >>> B-factor applied to coordinate array B: $Biso_model[f8.3]
 else
   display >>> initial B-factor correction: none
 end if
 ! End modification
 !
 
 {- MODIFIED 5/18/05 -}
 if ( &bulk_sol = true ) then 
     display >>> bulk solvent: probe radius=$solrad_best, shrink value=$solrad_best
     display >>> bulk solvent: density level= $sol_k_ref e/A^3, B-factor= $sol_b_ref A^2
 else
     display >>> bulk solvent: false
 end if
 {- END MODIFICATION -}
 
 if ( &obs_type = "intensity" ) then
   display >>> reflections with Iobs/sigma_I < &sigma_cut rejected
   display >>> reflections with Iobs > &obs_rms * rms(Iobs) rejected
 else
   display >>> reflections with |Fobs|/sigma_F < &sigma_cut rejected
   display >>> reflections with |Fobs| > &obs_rms * rms(Fobs) rejected
 end if
 xray anomalous=? end
 if ( $result = true ) then
   display >>> anomalous diffraction data was input
 end if
 
 {- MODIFIED 2/15/06 -}
 display >>> fft gridding factor = $fft_grid, B factor offset = $fft_b_add A^2, Elimit = $fft_elim
 {- END MODIFICATION -}
 
 display >>> theoretical total number of refl. in resol. range:    $total_theor[I6] ( 100.0 % )
 display >>> number of unobserved reflections (no entry or |F|=0): $unobserved[I6] ( $per_unobs[f5.1] % )
 display >>> number of reflections rejected:                       $rejected[I6] ( $per_reject[f5.1] % )
 display >>> total number of reflections used:                     $total_used[I6] ( $per_used[f5.1] % )
 display >>> number of reflections in working set:                 $total_work[I6] ( $per_work[f5.1] % )
 display >>> number of reflections in test set:                    $total_test[I6] ( $per_test[f5.1] % )

 display
 display ------> initial wa estimate= $wa_init
 display

 set display=OUTPUT end

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

 for $scale in ( &wa_scale ) loop main

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

   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

   evaluate ($wa_current=$wa_init*$scale)
   xray wa=$wa_current end

   if ( &init_min > 0 ) then
     minimize powell                             
       nstep=&init_min                                             
       nprint=5
       drop=10.0
     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_temp;
                               )
     evaluate ($wa_current=$wa_temp*$scale)
   else
     evaluate ($wa_current=$wa_init*$scale)
   end if
   xray wa=$wa_current end

   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

   do (mass=store5) (all)

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

   xray
     tolerance=0.0 lookup=false
   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_temp;
                               )
     evaluate ($wa_current=$wa_temp*$scale)
   else
     evaluate ($wa_current=$wa_init*$scale)
   end if
   xray wa=$wa_current end

   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
  
   print threshold=20.0 bond
   evaluate ($rmsd_bond=$result)

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

   set display=&list_outfile end

   evaluate ($wa_print=$wa_init*$scale)

   display wa= $wa_print[f15.4] r= $full_r[f6.4] free_r= $full_test_r[f6.4] rmsd bonds= $rmsd_bond[f8.6]  rmsd angles= $rmsd_angle[f8.5]

   set display=OUTPUT end

   if ( $full_test_r <= $best_test_r ) then
     evaluate ($best_test_r=$full_test_r)
     evaluate ($best_wa=$wa_print)
   end if

 end loop main

 set display=&list_outfile end

 display
 display ------> lowest free_r for target &STRIP%reftarget with wa= $best_wa
 display
 display ------> The best value for wa is one that produces a low 
 display         free_r value and also good model geometry. As a general
 display         guide, good geometry is indicated by rms deviations from
 display         ideal values of less than:
 display              bonds -> 0.01 A
 display              angles -> 2.00 degrees
 display  
 stop