{+ 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(

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

{* coordinate file *}
{===>} coordinate_infile="amy.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.mtf";

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

{* force field settings file *}
{===>} force_field_infile="";

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

{* 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 optimization for e-density level sol_k (e/A^3) and B-factor sol_b (A^2) *}
{+ choice: true false +}
{===>} sol_auto=true;

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

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

{* solvent mask parameters *}
{+ table: rows=1 "bulk solvent" cols=2 "probe radius (A) (usually set to 1)" "shrink radius (A) (usually set to 1)" +}
{===>} sol_rad=1.0;
{===>} sol_shrink=1.0;

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

{- begin modification, ATB, 11/29/08 -}
{* atom selections for non-default rigid groups during torsion angle dynamics *}
{* note: the selections must be non-overlapping *}
{===>} atom_rigid_1=(none);
{===>} atom_rigid_2=(none);
{===>} atom_rigid_3=(none);
{===>} atom_rigid_4=(none);
{===>} atom_rigid_5=(none);
{===>} atom_rigid_6=(none);
{===>} atom_rigid_7=(none);
{===>} atom_rigid_8=(none);
{===>} atom_rigid_9=(none);
{===>} atom_rigid_10=(none);
! to add more groups add more numbered entries:
!   {===>} atom_rigid_11=(none);
!   {===>} atom_rigid_12=(none);
!   {===>} atom_rigid_13=(none);
! etc

{* 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="cartesian";

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

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

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

{* torsion angle topology modification file *}
{===>} torsion_infile="CNS_TOPPAR:torsionmdmods";

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

 xray

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

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

   @CNS_XRAYLIB:scatter.lib

   evaluate ($counter=1)
   evaluate ($done=false)
   while ( $done = false ) loop read
    if ( &exist_reflection_infile_$counter = true ) then
      if ( &BLANK%reflection_infile_$counter = false ) then
       reflection
         @@&reflection_infile_$counter
       end
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
  end loop read

 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)
   else
     evaluate ($reject_obs=&obs_f)
     evaluate ($reject_sig=&obs_sigf)
   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) ( ( $STRIP%reject_sig # 0 ) 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                                       
    include xref                   
   ?                                        
 end      
 if ( &BLANK%force_field_infile = true ) then
    flags                                       
       exclude elec pele include vdw pvdw                 
      ?                                        
    end      
 else
     @&force_field_infile
 end if                                   

 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

 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 >>> 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 >>> 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 >>> 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 >>> 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 >>> custom patch file = &STRIP%patch_infile
   end if

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

 display >>> input coordinates: &STRIP%coordinate_infile
 if ( &BLANK%anom_library = false ) then
   display >>> anomalous f' f'' library: &STRIP%anom_library
 end if

 evaluate ($counter=1)
 evaluate ($done=false)
 while ( $done = false ) loop read
    if ( &exist_reflection_infile_$counter = true ) then
      if ( &BLANK%reflection_infile_$counter = false ) then
         display >>> reflection file $counter : &STRIP%reflection_infile_$counter       
      end if
   else
     evaluate ($done=true)
   end if
   evaluate ($counter=$counter+1)
 end loop read
 
 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):          $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 lbfgs                             
       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 lbfgs
           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

   {* modification: moved mass modification before velocity assignment, ATB 7/8/2008 *}
   do (store5=mass) ( all )
   do (mass=max(10,min(30,mass))) ( all )

   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=-1
         maxchain=-1
         maxtree=-1
         kdihmax = 95.
         
         {- begin modification, ATB, 11/29/08 -}
         @CNS_XTALMODULE:messagesilent(&message_old_silent=$message_old_silent; &echo_old_silent=$echo_old_silent; )     
         evaluate ($atr_count=1)
         evaluate ($atr_done=false)
         while ( $atr_done = false ) loop atrl
          if ( &exist_atom_rigid_$atr_count = true ) then
           fix group ( &atom_rigid_$atr_count ) 
           evaluate ($atr_count=$atr_count+1)
          else
           evaluate ($atr_done=true)
          end if
         end loop atrl
         @CNS_XTALMODULE:messagerevert(&message_old_silent=$message_old_silent; &echo_old_silent=$echo_old_silent; )     
         
         if ( &BLANK%torsion_infile = false ) then
            @&torsion_infile
         else         
            @CNS_TOPPAR:torsionmdmods
         end if
         
       end
       nstep=0
       cmremove=true
     end
   end if
  
   if ( &md_scheme = "slowcool" ) then
  
     evaluate ( $curr_temp = &temperature )
     while ( $curr_temp > 0.0 ) loop cool
       if ( &md_type = "torsion" ) then
         dynamics torsion
           timestep=$time_step
           nstep=$md_steps
           nprint=5
           cmremove=false
           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